diff --git a/Analysis/Code/ggplot_time_series_with_trends.R b/Analysis/Code/ggplot_time_series_with_trends.R
index 986a51e..4982636 100644
--- a/Analysis/Code/ggplot_time_series_with_trends.R
+++ b/Analysis/Code/ggplot_time_series_with_trends.R
@@ -19,18 +19,40 @@ setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/P
 #Set inputs
 ######################
 Case.Study <- "CASESTUDY2"
-Estuary <- "HUNTER"
-Climvar <- 'tasmean'
-ggplotGAM.k <- 14
+Estuary <- "RICHMOND"
+Riv.Gauge.loc <- "OAKLANDROAD"     #GRETA 
+Est.Gauge.loc <-  "CORAKI"          #"RAYMONDTERRACE" # "HEXHAMBRIDGE"
+logtransformFlow <- TRUE
+ggplotGAM.k <- 15
+rivTempGAM.k <- 20
+Version <- 'V2'
 ######################
 
+######################
+Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/Riv_', Riv.Gauge.loc,'_Est_',Est.Gauge.loc,'_GAMk', ggplotGAM.k, '/', sep="")
+if (file.exists(Output.Directory)){
+  print('output folder already existed and was not created again')
+  } else {
+  dir.create(file.path(Output.Directory))
+  print('output folder did not exist and was created')
+  }
+######################
 
 ######################
 #Set input file paths
 ######################
-AirT_CSV_Path <- "./Data/SILO/CASESTUDY2/HUNTER/SILO_climdata_HUNTER_Catchment_-32.162479_150.5335812.csv"
-RivT_CSV_Path <- "./Data/River_Gauge_Data/HUNTER_210064_20180615/210064_Temp.csv"
-SST_CSV_Path <- './Data/NARCLIM_Site_CSVs/CASESTUDY2/HUNTER_32.751_151.690/sstmean_NNRP_HUNTER_33.034_152.156_NARCliM_summary.csv'
+
+pattern = paste('SILO_climdata_', Estuary,'*', sep="")
+AirT_CSV_Path <- list.files(paste("./Data/SILO/",Case.Study, '/',sep=""), pattern, full.names=T, recursive=T)
+
+pattern = paste(Estuary,'@', Riv.Gauge.loc, '.*.ALL.csv', sep="")
+RivT_CSV_Path <- list.files("./Data/River_Gauge_Data/", pattern, full.names=T)
+
+pattern = paste(Estuary,'@', Est.Gauge.loc, '.*.ALL.csv', sep="")
+EstT_CSV_Path <- list.files("./Data/River_Gauge_Data/", pattern, full.names=T)
+
+pattern = paste('sstmean_NNRP_', Estuary,'*', sep="")
+SST_CSV_Path <- list.files(paste("./Data/NARCLIM_Site_CSVs/",Case.Study, '/', sep=""), pattern, full.names=T, recursive=T)
 ######################
 
 
@@ -62,9 +84,21 @@ AirT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 35
 
 ############River temp
 #Load a daily (no gaps) time series as a time serie baseline for other time series used here
+#Here we use the raw DPI CSV format that comes with a bunch of metadata
 RivT.df <- data.frame(read.csv(RivT_CSV_Path))
-RivT.full.TS <- zoo(RivT.df$Temp, order.by= as.Date(RivT.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
-RivT.TS <- window(RivT.full.TS, start=as.Date("1995-01-01"), end=as.Date("2018-01-01"))
+char.df <- data.frame(lapply(RivT.df[2,], as.character), stringsAsFactors=FALSE)
+#dat <- data.frame(lapply(RivT.df[(4:nrow(RivT.df)),(2:ncol(RivT.df))], as.numeric), stringsAsFactors=FALSE)
+#dat <- RivT.df.num[!is.na(as.numeric(as.character(RivT.df.num))),]
+dat <- RivT.df[(4:nrow(RivT.df)),]
+colnames(dat) <- lapply(RivT.df[2,], as.character)
+dat$Date <- gsub(x=dat$Date,pattern="00:00:00",replacement="",fixed=T)
+colnames(dat) <- gsub(x=colnames(dat),pattern="Water Temp(C)",replacement="Temp",fixed=T)
+RivT.df <- dat
+rm(dat,char.df)
+RivT.full.TS <- zoo(as.numeric(as.character(RivT.df$Temp)), order.by= as.Date(RivT.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+RivT.full.TS <- window(RivT.full.TS, start=as.Date("1995-07-01"), end=as.Date("2018-01-01"))
+RivT.full.TS <- na.approx(RivT.full.TS)
+RivT.TS <- RivT.full.TS
 RivT.full.df <- data.frame(RivT.TS)   ### This is only done because 
 RivT.df <- data.frame(RivT.TS)
 colnames(RivT.df) <- 'rivTmean'
@@ -80,6 +114,37 @@ RivT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
 RivT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
 ############River temp
 
+############River flow
+#Load a daily (no gaps) time series as a time serie baseline for other time series used here
+#Here we use the raw DPI CSV format that comes with a bunch of metadata
+RivQ.df <- data.frame(read.csv(RivT_CSV_Path))
+char.df <- data.frame(lapply(RivQ.df[2,], as.character), stringsAsFactors=FALSE)
+dat <- RivQ.df[(4:nrow(RivQ.df)),]
+colnames(dat) <- lapply(RivQ.df[2,], as.character)
+dat$Date <- gsub(x=dat$Date,pattern="00:00:00",replacement="",fixed=T)
+colnames(dat) <- gsub(x=colnames(dat), pattern="Discharge (ML/d)",replacement="Flow",fixed=T)
+RivQ.df <- dat
+rm(dat,char.df)
+#RivQ.full.TS <- zoo(log10((as.numeric(as.character(RivQ.df$Flow))) +1) , order.by= as.Date(RivQ.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+RivQ.full.TS <- zoo(as.numeric(as.character(RivQ.df$Flow)), order.by= as.Date(RivQ.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+
+RivQ.TS <- window(RivQ.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
+RivQ.full.df <- data.frame(RivQ.full.TS)   ### This is only done because 
+RivQ.df <- data.frame(RivQ.TS)
+colnames(RivQ.df) <- 'RivQmean'
+colnames(RivQ.full.df) <- 'RivQmean'
+############trends
+RivQ.df$Julday1 <- seq(1,length(RivQ.df[,1]),1)
+linear.trend.model_EC_all <- lm(RivQmean ~ Julday1, RivQ.df)
+RivQ.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+RivQ.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+RivQ.full.df$Julday1 <- seq(1,length(RivQ.full.df[,1]),1)
+linear.trend.model_EC_all <- lm(RivQmean ~ Julday1, RivQ.full.df)
+RivQ.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+RivQ.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+############River Flow
+
+
 ############ SST
 #Load a daily (no gaps) time series as a time serie baseline for other time series used here
 SST.df <- data.frame(read.csv(SST_CSV_Path))
@@ -102,6 +167,42 @@ SST.lintrend <- summary(linear.trend.model_EC_all2)$coefficients[2,1] * 356
 ############ SST
 
 
+############Estuary temp
+#Load a daily (no gaps) time series as a time serie baseline for other time series used here
+#Here we use the raw DPI CSV format that comes with a bunch of metadata
+EstT.df <- data.frame(read.csv(EstT_CSV_Path))
+char.df <- data.frame(lapply(EstT.df[2,], as.character), stringsAsFactors=FALSE)
+#dat <- data.frame(lapply(EstT.df[(4:nrow(EstT.df)),(2:ncol(EstT.df))], as.numeric), stringsAsFactors=FALSE)
+
+#dat <- EstT.df.num[!is.na(as.numeric(as.character(EstT.df.num))),]
+dat <- EstT.df[(4:nrow(EstT.df)),]
+
+colnames(dat) <- lapply(EstT.df[2,], as.character)
+dat$Date <- gsub(x=dat$Date,pattern="00:00:00",replacement="",fixed=T)
+colnames(dat) <- gsub(x=colnames(dat),pattern="Water Temp(C)",replacement="Temp",fixed=T)
+EstT.df <- dat
+rm(dat,char.df)
+#replace negative values with NA
+EstT.df$Temp <- replace(EstT.df$Temp, which(as.numeric(as.character(EstT.df$Temp)) < 11), NA)
+EstT.full.TS <- zoo(as.numeric(as.character(EstT.df$Temp)), order.by= as.Date(EstT.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+plot(EstT.full.TS)
+EstT.TS <- window(EstT.full.TS, start=as.Date("2013-06-01"), end=as.Date("2018-01-01"))
+EstT.full.TS <- EstT.TS
+EstT.full.df <- data.frame(EstT.TS)   ### This is only done because of poor data at beginning
+EstT.df <- data.frame(EstT.TS)
+colnames(EstT.df) <- 'EstTmean'
+colnames(EstT.full.df) <- 'EstTmean'
+############
+EstT.df$Julday1 <- seq(1,length(EstT.df[,1]),1)
+linear.trend.model_EC_all <- lm(EstTmean ~ Julday1, EstT.df)
+EstT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+EstT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+EstT.full.df$Julday1 <- seq(1,length(EstT.full.df[,1]),1)
+linear.trend.model_EC_all <- lm(EstTmean ~ Julday1, EstT.full.df)
+EstT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+EstT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+############Est temp
+
 
 ######################
 #Plot
@@ -110,57 +211,91 @@ SST.lintrend <- summary(linear.trend.model_EC_all2)$coefficients[2,1] * 356
 ##################################### Full Time Period
 #Air temp Full period 
 p1air <- ggplot(AirT.full.df, aes(y=tasmean, x=index(AirT.full.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in catchment airT (SILO) lin trend was ",
+  ggtitle(paste(Estuary," Catchment Centroid - Linear and smooth trend in catchment airT (SILO) lin trend was ",
                 round(AirT.full.lintrend,3), ' C°/year with p=', round(AirT.full.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
   ylab("Air Temperature [C°]") + xlab("Time")
 
-#Riv temp Full period 
 p1riv <- ggplot(RivT.full.df, aes(y=rivTmean, x=index(RivT.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
+  ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
                 round(RivT.full.lintrend,3), ' C°/year with p=', round(RivT.full.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
   ylab("River Temperature [C°]") + xlab("Time")
 
+
+if(logtransformFlow ==TRUE){
+p1rivQ <- ggplot(RivQ.full.df, aes(y=log10(RivQmean+2), x=index(RivQ.full.TS))) + geom_line(alpha=0.5) + 
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
+                  round(RivQ.full.lintrend,3), ' ML/day /year with p=', round(RivQ.full.pvalNCV_ECall,10), sep=" ")) +
+    theme(plot.title=element_text(face="bold", size=9)) +
+    geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+    stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+    ylab(expression(paste("log10(River Flow [ML/day] + 2)", sep=""))) + xlab("Time")
+  } else {
+  p1rivQ <- ggplot(RivQ.full.df, aes(y=RivQmean, x=index(RivQ.full.TS))) + geom_line(alpha=0.5) + 
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
+                  round(RivQ.full.lintrend,3), ' ML/day with p=', round(RivQ.full.pvalNCV_ECall,10), sep=" ")) +
+    theme(plot.title=element_text(face="bold", size=9)) +
+    geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+    stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+    ylab(expression(paste("River Flow [ML/day]", sep=""))) + xlab("Time") 
+    }
+
+
 #Sea temp Full period 
 p1sst <- ggplot(SST.full.df, aes(y=SSTmean, x=index(SST.full.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
+  ggtitle(paste(Estuary, "NNRP (NARCLIM reanalysis) - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
                 round(SST.full.lintrend,3), ' C°/year with p=', round(SST.full.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
   ylab("Sea Surface Temperature [C°]") + xlab("Time")
 
+p1Est <- ggplot(EstT.full.df, aes(y=EstTmean, x=index(EstT.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary,'@', Est.Gauge.loc, " - Linear and smooth trend in Estuary temperature (GAUGE) lin trend was ",
+                round(EstT.full.lintrend,3), ' C°/year with p=', round(EstT.full.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+  stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+  ylab("Estuary Temperature [C°]") + xlab("Time")
+
+
 #export to png
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_tasmean_full_period_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_tasmean_full_period_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(p1air,ncol=1)
 dev.off()
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_rivTmean_full_period_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_rivTmean_full_period_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(p1riv,ncol=1)
 dev.off()
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_SSTmean_full_period_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_RivQmean_full_period_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1rivQ,ncol=1)
+dev.off()
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_SSTmean_full_period_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(p1sst,ncol=1)
 dev.off()
+png.name <- paste(Output.Directory, Estuary, '@', Est.Gauge.loc, '_Trends_estTmean_full_period_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1Est,ncol=1)
+dev.off()
 ##################################### Full Time Period
 
 
-
-
 ######################### 1990-present
-combined.TS <- window(merge(AirT.full.TS, window(RivT.full.TS, start=as.Date("1995-01-01"), end=end(RivT.full.TS)), SST.full.TS, all=T), start=as.Date("1990-01-01"), end=end(AirT.full.TS))
+combined.TS <- window(merge(AirT.full.TS, window(RivT.full.TS, start=as.Date("1995-01-01"), end=end(RivT.full.TS)), SST.full.TS,EstT.full.TS,RivQ.full.TS,  all=T), start=as.Date("1990-01-01"), end=end(AirT.full.TS))
 combined.df <- data.frame(combined.TS)
-colnames(combined.df) <- c('tasmean','rivTmean', 'SSTmean')
+colnames(combined.df) <- c('tasmean','rivTmean', 'SSTmean', 'EstTmean', 'rivQmean')
 
 #Air temp 
-p1air <- ggplot(combined.df, aes(y=tasmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in catchment airT (SILO) lin trend was ",
+p2air <- ggplot(combined.df, aes(y=tasmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary," Catchment centroid - Linear and smooth trend in catchment airT (SILO) lin trend was ",
                 round(AirT.lintrend,3), ' C°/year with p=', round(AirT.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
@@ -168,49 +303,109 @@ p1air <- ggplot(combined.df, aes(y=tasmean, x=index(combined.TS))) + geom_line(a
   ylab("Air Temperature [C°]") + xlab("Time")
 
 #Riv temp 
-p1riv <- ggplot(combined.df, aes(y=rivTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
+p2riv <- ggplot(combined.df, aes(y=rivTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
                 round(RivT.lintrend,3), ' C°/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
   ylab("River Temperature [C°]") + xlab("Time")
 
+#Riv flow
+if(logtransformFlow ==TRUE){
+  p2rivQ <- ggplot(combined.df, aes(y=log10(rivQmean+2), x=index(combined.TS))) + geom_line(alpha=0.5) + 
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
+                  round(RivQ.full.lintrend,3), ' ML/day /year with p=', round(RivQ.full.pvalNCV_ECall,10), sep=" ")) +
+    theme(plot.title=element_text(face="bold", size=9)) +
+    geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+    stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+    ylab(expression(paste("log10(River Flow [ML/day] + 2)", sep=""))) + xlab("Time")
+  } else {
+  p2rivQ <- ggplot(combined.df, aes(y=rivQmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
+                  round(RivT.lintrend,3), ' C°/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
+    theme(plot.title=element_text(face="bold", size=9)) +
+    geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+    stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+    ylab(expression(paste("River Flow [ML/day]", sep=""))) + xlab("Time")
+    }
+
 #Sea temp 
-p1sst <- ggplot(combined.df, aes(y=SSTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary, " - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
+p2sst <- ggplot(combined.df, aes(y=SSTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary,"NNRP(NARCLIM reanalysis) - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
                 round(SST.lintrend,3), ' C°/year with p=', round(SST.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
   ylab("Sea Surface Temperature [C°]") + xlab("Time")
 
+#sst temp 
+p2Est <- ggplot(combined.df, aes(y=EstTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary,'@', Est.Gauge.loc, " - Linear and smooth trend in Estuary temperature (GAUGE) lin trend was ",
+                round(EstT.lintrend,3), ' C°/year with p=', round(EstT.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+  #stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
+  ylab("Estuary Temperature [C°]") + xlab("Time")
+
 #export to png
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_tasmean_1990-present_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_tasmean_1990-present_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
-grid.arrange(p1air,ncol=1)
+grid.arrange(p2air,ncol=1)
 dev.off()
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_rivTmean_1990-present_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_rivTmean_1990-present_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
-grid.arrange(p1riv,ncol=1)
+grid.arrange(p2riv,ncol=1)
 dev.off()
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_SSTmean_1990-present_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_rivQmean_1990-present_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
-grid.arrange(p1sst,ncol=1)
+grid.arrange(p2rivQ,ncol=1)
+dev.off()
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_SSTmean_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p2sst,ncol=1)
 dev.off()
 
 #export an overview plot as png
-gA <- ggplotGrob(p1air)
-gB <- ggplotGrob(p1riv)
-gC <- ggplotGrob(p1sst)
-gA$widths <- gB$widths
-gC$widths <- gB$widths
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_Summary_1990-present_', Sys.Date(),".png", sep="")
+gA <- ggplotGrob(p2air)
+gB <- ggplotGrob(p2riv)
+gC <- ggplotGrob(p2sst)
+gD <- ggplotGrob(p2Est)
+gE <- ggplotGrob(p2rivQ)
+gA$widths <- gE$widths
+gC$widths <- gE$widths
+gD$widths <- gE$widths
+gB$widths <- gE$widths
+
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_SST_RivT_AirT_1990-present_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(gA,gB ,gC,ncol=1)
 dev.off()
 
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_SST_RivT_EstT_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gB,gD,gC,ncol=1)
+dev.off()
+
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_RivQ_RivT_EstT_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gE,gB,gD,ncol=1)
+dev.off()
+
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_RivQ_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gA,gB,gE,ncol=1)
+dev.off()
+
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gA,gB,ncol=1)
+dev.off()
 
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_RivT_RivQ_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gB,gE,ncol=1)
+dev.off()