CE 569: Watershed Analysis
SUMMARY REPORT - LAB 5 - AWBM-MODEL ASSESSMENT
Professor: Dr. Tyler Smith
By Cássio Rampinelli
October, 25th, 2018
OBS: This script was done in R programming language
__________________________________________________________________________________________________________________________
INTRODUCTION
The goal of this script is summarizing and briefly discuss the outputs of the AWBM hydrologic model applied to a study site at the Gauging Station SEAUS_218007. The main idea is showing basic hydrologic model assessment procedures by presenting the programming code lines to do that with some plots. The guidelines of the assessment were provided as part of Lab 5 assignment of the course Watershed Analysis - Fall/2018.
Although the algorithms and basic code were provided in Matlab programming language, they were converted to R programming language. The conversion task not only helped a better understanding of the algorithms but also allow sharing a version of the algorithms in R.
The data (daily precipitation, evapotranspiration and observed flow) were provided for the study site.
LOADING THE DATA
The provided data was the file: “seaus_218007.mat” this file contains daily precipitation, evapotranspiration and observed flow data for an approximated period of 10 years. To load this type of file in R (.mat) we need to install a package called “R.matlab” and use a function of this file named readMat. The content of the file will be transformed in a type of object called “list” in R. The code lines below shows how to load the file and extract from them the series of precipitation, potential evapotranspiration and observed flows. These data will be saved in vectors named “p”, “pet” and “qobs” respectively. (Obs: the data should be in the same folder were the R script is. To guarantee that we set the working directory.)
###Installing package to load .mat data
#install.packages("R.matlab")
require("R.matlab")## Loading required package: R.matlab## Warning: package 'R.matlab' was built under R version 3.4.4## R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.##
## Attaching package: 'R.matlab'## The following objects are masked from 'package:base':
##
## getOption, isOpen#Setting working directory
setwd("E:/Users/Cassio/Docs/Clarkson University/Fall_2018/CE_569_Watershed Analysis/Labs/Lab5/Rcodes")
#Loading the data
dados=readMat("seaus_218007.mat")
#Accessing all the data and separating them in new variables.
dados=dados[[1]]
#Precipitation,potential evapotranspiration and observed flow
p=dados[[6]]
pet=dados[[7]]
qobs=dados[[5]]
Checking the plots of the time series
par(mfrow=c(1,1))
plot(p,ylab="Precipitation",xlab="days",type="l")
plot(pet,ylab="Potential Evapotranspiration",xlab="days",type="l")
plot(pet,ylab="Observed Flow",xlab="days",type="l")
BASIC FUNCTIONS AND ALGORITHMS
In this section, the R code of the basic functions and algorithms that will be used to perform the analysis are presented. The details about the metrics and equations should be verified in specific references and they are not the scope of this report.
AWBM Hydrologic Model
A simplified version of the Australian Water Balance Model (AWBM) was employed to perform the analysis.
More details about the model can be checked in Boughton(2004)(https://www.sciencedirect.com/science/article/pii/S1364815203002196) and are out of the scope of this summary report. Here we just present a picture of the basic framework and the function of the model coded in R.
#AWBM hydrologic Model Function
#input:
#P= Precipitation time serie
#PET=Potential Evapotranspiration time serie
#PARAMS=Vector with 3 parameters:S,K and BFI
#Output:
#Q=Simulated flow
awbm<-function(P,PET,PARAMS){
S = PARAMS[1]
K = PARAMS[2]
BFI = PARAMS[3]
n = length(P)
#Pre-allocate variable
AET = 0
SS = 0
BS = 0
ER = 0
SF = 0
BF = 0
Q = 0
for(i in 2:n){
#prevent AET > available water
AET[i] = min((SS[i-1]/S)*PET[i], SS[i-1]+P[i]);
SS[i] = SS[i-1] + P[i] - AET[i];
#prevent negative ER
ER[i] = max(0, SS[i]-S)
BS[i] = BS[i-1] + ER[i]*K #Eqn4
SF[i] = ER[i]*(1-K) #Eqn5
BF[i] = BS[i]*(1-BFI) #Eqn6
Q[i] = SF[i] + BF[i] #Eqn7
BS[i] = BS[i] - BF[i] #remove BF from BS
}
return(Q)
}
Objective Functions
The assessment was done based on three objective functions: Nash-Sutcliffe Efficiency (NSE), Deviation of Runoff Volume (DRV) and Kling-Gupta Efficiency (KGE). The equations for each objective functions are presented below followed by the codes in R. Specific questions about them can be checked in the specific references and are out of the scope of this report.
Figure 2: Objective functions expressions (CE_569_Smith_Lecture Note)
#Nash Sutcliffe Efficiency Function
of_nse<-function(obs,pred){
#CALCULATE NASH-SUTCLIFFE EFFICIENCY --------------------------------%%%
# obs = observed streamflow
# pred = predicted streamflow
O=obs
Obar=mean(obs)
P=pred
numer=(O-P)^2
denom=(O-Obar)^2
yfit= 1-(sum(numer)/sum(denom))
#DEFINE OUTPUTS
return(yfit)
}
#KLING-GUPTA Efficiency Function
of_kge<-function(obs,pred){
#CALCULATE KLING-GUPTA Efficiency --------------------------------%%%
# obs=observed streamflow
# pred=predicted streamflow
O=obs
P=pred
r= cor(O,P)
b= mean(P)/mean(O)
g= (sd(P)/mean(P))/(sd(O)/mean(O))
yfit= 1 - sqrt( (r-1)^2 + (b-1)^2 + (g-1)^2 )
output=list(KGE=yfit,Corr=r, Beta= b,Gamma=g)
return(output)
}
#DEVIATION OF RUNOFF VOLUME Function
of_drv<-function(obs,pred){
#CALCULATE DEVIATION OF RUNOFF VOLUME --------------------------------%%%
# obs=observed streamflow
# pred=predicted streamflow
O=obs
P=pred
yfit=sum(P)/sum(O);
return(yfit)
}
Hydrologic Signatures
In this report, the Runoff Ratio (RR) and the Flashiness Index (FI) were used as hydrologic signatures. The equations and codes are presented below.
Figure 3: Runoff Ratio and Flashiness Index equations (CE_569_Smith_Lecture Note)
#Runoff Ratio Function
hs_rr<-function(q,p){
#Runoff Ratio
#Inputs:
#q: time series (of streamflow)
#p: time series (of precipitation)
#Outputs:
#rr: runoff ratio
#rr = mean(q)/mean(p)
rr = mean(q)/mean(p)
return(rr)
}
#Flashiness Index Function
hs_fi<-function(q){
#Flashiness Index - refer to Baker et al. (2004), JAWRA.
#Inputs:
#q: time series of streamflow
#Outputs:
#fi: flashiness index of q
#fi = sum( abs( q(i)-q(i-1) ) ) / sum( q(i) )
n=length(q);
t1=0
for(i in 2:n) {
t1[i] = q[i]-q[i-1];
}
fi= sum( abs( t1 ) ) / sum ( q );
return(fi)
}Particle Swarm Optimization Algorithm
In addition to the requirements of this assignment, the Particle Swarm Optimization (PSO) algorithm will be also presented. The algorithm is a population-based stochastic optimization technique developed by Dr. Eberhart and Dr. Kennedy in 1995 (http://www.swarmintelligence.org/index.php). More information can be found in specific references and is out of the scope of this report. The algorithm was adapted to be applied in a generic hydrologic model as an additional tool to perform the calibration process and is presented below.
#PSO Function
#Input Objective Function (funcao)
#Parameters range (limites)
PSO<-function(funcao,limites)
{
n<-ncol(limites)#number of columns
npar<-100#total number of particles
maxger<-100#Number of iterations
POP<-matrix(0,nrow=npar,ncol=n)#Record the population
######################
###INITIAL POPULATION
######################
for (i in 1:npar)
{
for (j in 1:n)
{
POP[i,j]= (limites[2,j]-limites[1,j])*runif(1)+limites[1,j]
}
}
###Initializing constants###
criterio<-0#Stoped criteria
geracao<-0#Counter of the number of iterations
#Minimum value of the function for a particle
pbestAJ<-matrix(0,nrow=1,ncol=npar)+99999999999999999999
#Lowest particle value
pbestXX<-matrix(99999999999999999999,nrow=npar,ncol=n)
#Velocity
vant<-matrix(0,nrow=npar,ncol=n)
#Trust parameter
c1=2
c2=2
#Particle inertia.
#Higher values contributes to global otptimun.
w=1
passo=(w-0.4)/maxger
#maximum velocity
vmax=(limites[2,]-limites[1,])*0.2
v<-matrix(0,nrow=npar,ncol=n)
AJ<-0
x<-0
#######Iterative process#####################
#######
while (!criterio)
{
geracao=geracao+1 #update the iteration number.
w=w-passo #inertial parameter update
for (i in 1:npar)
{
for(j in 1:n)
{
x[j]<-POP[i,j]#particle position
}
AJ[i]=funcao(x)
if(AJ[i]<pbestAJ[i])
{
pbestAJ[i]=AJ[i] #best particle value
pbestXX[i,]=POP[i,] #best position
}
}
gbestAJindice<-which(pbestAJ==min(pbestAJ),arr.ind=TRUE)
gbestXX=pbestXX[gbestAJindice[1,2],]
for (i in 1:npar)
{
v[i,]=(w*vant[i])+(c1*runif(1)*(pbestXX[i,]-POP[i,]))+(c2*runif(1)*(gbestXX-POP[i,]))
for(j in 1:n)
{
if(v[i,j]>vmax[j])
{
v[i,j]=vmax[j]
}
}
POP[i,]=POP[i,]+v[i,]
vant[i]=v[i]
}
if (geracao==maxger)
{criterio=1}
}
result<-list("Par_Otimo:",gbestXX,"FO:",funcao(gbestXX))
return(result)
}
ASSESSMENT: PART 1
In this section the first part of the assignment requirements is presented.
a) setup the AWBM model (using seaus_218007.mat as forcing data) for a 6-yr calibration period followed by a 4-year validation period.
The data was split separating 6 years for calibration and 4 years for validation. To simplify, the number of days in a year is assumed 365 days ( it was assumed that we don’t have leap years are). The variable calibYears is used to set the number of years for calibration considering that the number of days will be the number of years times 365 days. The remaining days will be used for calibration. Those variables are used to refer the limits of the vector/matrix indexes as will be presented along the code.
#Setting the number of years for calibration
calibYears=6
#Days for calibration
calibDays=calibYears*365
b) draw 100,000 parameter sets.
#Sampling 100,000 parameters
n=100000 # of samples to draw
S=runif(n,0,500)
K=runif(n,0,1)
BFI=runif(n,0,1)
#Matrix to save the parameters
params=matrix(0,nrow=n,ncol=3)
params=cbind(S,K,BFI)
#Matrix to save the simulated flows for the calibration period
qcalb=matrix(0,nrow=calibDays,ncol=n)
#Matrix to save the simulated flows for the validation period
qvald=matrix(0,nrow=(length(qobs)-calibDays),ncol=n)
c) determine the NSE, KGE, and DRV for calibration and validation periods.
#Initializing NSE, KGE e DRV variables for calibration period
nse=0
kge=0
drv=0
#Running for Calibration Period
for (i in 1: n){
teste= awbm(p[1:calibDays],pet[1:calibDays],params[i,])
qcalb[,i]=awbm(p[1:calibDays],pet[1:calibDays],params[i,])
nse[i]=of_nse(qobs[1:calibDays],qcalb[,i])
temp=of_kge(qobs[1:calibDays],qcalb[,i])
kge[i]=temp[[1]]
temp=0
drv[i]=of_drv(qobs[1:calibDays],qcalb[,i])
}
#Initializing NSE, KGE e DRV variables for validation period
nseVal=0
kgeVal=0
drvVal=0
#Running for Validation Period
for (i in 1: n){
qvald[,i]=awbm(p[(calibDays+1):length(qobs)],pet[(calibDays+1):length(qobs)],params[i,])
nseVal[i]=of_nse(qobs[(calibDays+1):length(qobs)],qvald[,i])
temp=of_kge(qobs[(calibDays+1):length(qobs)],qvald[,i])
kgeVal[i]=temp[[1]]
temp=0
drvVal[i]=of_drv(qobs[(calibDays+1):length(qobs)],qvald[,i])
}
d)Create dotty plots of parameters for each of the above objective functions.
Dotty plots of the parameters considering each objective function for the calibration period.
Nash-Sutcliffe Efficiency (NSE)
plot(params[,1],nse,ylab="NSE",xlab="S")
plot(params[,2],nse,ylab="NSE",xlab="K")
plot(params[,3],nse,ylab="NSE",xlab="BFI")
Kling-Gupta Efficiency (KGE)
plot(params[,1],kge,ylab="KGE",xlab="S")
plot(params[,2],kge,ylab="KGE",xlab="K")
plot(params[,3],kge,ylab="KGE",xlab="BFI")
Deviation of Runoff Volume (DRV)
plot(params[,1],drv,ylab="DRV",xlab="S")
plot(params[,2],drv,ylab="DRV",xlab="K")
plot(params[,3],drv,ylab="DRV",xlab="BFI")
Dotty plots of the parameters considering each objective function for the validation period.
Nash-Sutcliffe Efficiency (NSE)
plot(params[,1],nseVal,ylab="NSE",xlab="S")
plot(params[,2],nseVal,ylab="NSE",xlab="K")
plot(params[,3],nseVal,ylab="NSE",xlab="BFI")
Kling-Gupta Efficiency (KGE)
plot(params[,1],kgeVal,ylab="KGE",xlab="S")
plot(params[,2],kgeVal,ylab="KGE",xlab="K")
plot(params[,3],kgeVal,ylab="KGE",xlab="BFI")
Deviation of Runoff Volume (DRV)
plot(params[,1],drvVal,ylab="DRV",xlab="S")
plot(params[,2],drvVal,ylab="DRV",xlab="K")
plot(params[,3],drvVal,ylab="DRV",xlab="BFI")
e) create a plot of the observed hydrograph, along with the optimal model prediction based on NSE, KGE, and DRV (for both calibration and validation periods).
First, the best set of parameters will be defined based on the mean of those set of parameters associated with the top 20% highest values for each objective function (NSE, KGE, and DRV ). A box plot of the top 20% parameters is also presented for each objective function. Then, a comparison between the observed and simulated flow for the best set of parameters for each objective function.
Nash-Sutcliffe Efficiency (NSE)
Calibrated period
Mean of the top 20% parameters
top20NseValues=tail(sort(nse),20)
top20NsePositions=which(nse %in% top20NseValues)
#top20NseValues
optimumParamsNSE=params[top20NsePositions,]
bestParameterNSE=colMeans(optimumParamsNSE)
bestParameterNSE## S K BFI
## 267.5779110 0.8485195 0.9513059Boxplot for the top 20% parameters
optimumParamsNSE=as.data.frame(optimumParamsNSE)
boxplot(optimumParamsNSE$S,xlab=c("S"))
boxplot(optimumParamsNSE[,2],optimumParamsNSE[,3],xlab=c("K BFI"))
Simulated X Observed flow
Qcalc= awbm(p[1:calibDays],pet[1:calibDays],bestParameterNSE)
plot(Qcalc,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for calibration period")
lines(qobs[1:calibDays],type ="l",col="black",lty=2)
legend(1500,22,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
Validation period
Mean of the top 20% parameters
top20NseValuesValid=tail(sort(nseVal),20)
top20NsePositionsValid=which(nseVal %in% top20NseValuesValid)
#top20NseValues
optimumParamsNSEValid=params[top20NsePositionsValid,]
bestParameterNSEValid=colMeans(optimumParamsNSEValid)
bestParameterNSEValid## S K BFI
## 380.0139234 0.5223576 0.6263451Boxplot for the top 20% parameters
optimumParamsNSEValid=as.data.frame(optimumParamsNSEValid)
boxplot(optimumParamsNSEValid$S,xlab=c("S"))
boxplot(optimumParamsNSEValid[,2],optimumParamsNSEValid[,3],xlab=c("K BFI"))
Simulated X Observed flow
QcalcVald= awbm(p[(calibDays+1):length(qobs)],pet[(calibDays+1):length(qobs)],bestParameterNSEValid)
plot(QcalcVald,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for Validation period")
lines(qobs[(calibDays+1):length(qobs)],type ="l",col="black",lty=2)
legend(1500,52,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
Kling-Gupta Efficiency (KGE)
Calibrated period
Mean of the top 20% parameters
top20KgeValues=tail(sort(kge),20)
top20KgePositions=which(kge %in% top20KgeValues)
#top20NseValues
optimumParamsKGE=params[top20KgePositions,]
bestParameterKGE=colMeans(optimumParamsKGE)
bestParameterKGE## S K BFI
## 491.2071483 0.3987526 0.9862826Boxplot for the top 20% parameters
optimumParamsKGE=as.data.frame(optimumParamsKGE)
boxplot(optimumParamsKGE$S,xlab=c("S"))
boxplot(optimumParamsKGE[,2],optimumParamsKGE[,3],xlab=c("K BFI"))
Simulated X Observed flow
Qcalc= awbm(p[1:calibDays],pet[1:calibDays],bestParameterKGE)
plot(Qcalc,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for calibration period")
lines(qobs[1:calibDays],type ="l",col="black",lty=2)
legend(1500,22,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
Validation period
Mean of the top 20% parameters
top20KgeValuesValid=tail(sort(kgeVal),20)
top20KgePositionsValid=which(kgeVal %in% top20KgeValuesValid)
#top20NseValues
optimumParamsKGEValid=params[top20KgePositionsValid,]
bestParameterKGEValid=colMeans(optimumParamsKGEValid)
bestParameterKGEValid## S K BFI
## 330.6106751 0.6490375 0.9914439Boxplot for the top 20% parameters
optimumParamsKGEValid=as.data.frame(optimumParamsKGEValid)
boxplot(optimumParamsKGEValid$S,xlab=c("S"))
boxplot(optimumParamsKGEValid[,2],optimumParamsKGEValid[,3],xlab=c("K BFI"))
Simulated X Observed flow
QcalcVald= awbm(p[(calibDays+1):length(qobs)],pet[(calibDays+1):length(qobs)],bestParameterKGEValid)
plot(QcalcVald,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for Validation period")
lines(qobs[(calibDays+1):length(qobs)],type ="l",col="black",lty=2)
legend(520,22,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
Deviation of Runoff Volume (DRV)
Calibrated period
Mean of the top 20% parameters
top20DrvValues=tail(sort(drv),20)
top20DrvPositions=which(drv %in% top20DrvValues)
#top20NseValues
optimumParamsDRV=params[top20DrvPositions,]
bestParameterDRV=colMeans(optimumParamsDRV)
bestParameterDRV## S K BFI
## 87.6737235 0.4973332 0.4774285Boxplot for the top 20% parameters
optimumParamsDRV=as.data.frame(optimumParamsDRV)
boxplot(optimumParamsDRV$S,xlab=c("S"))
boxplot(optimumParamsDRV[,2],optimumParamsDRV[,3],xlab=c("K BFI"))
Simulated X Observed flow
Qcalc= awbm(p[1:calibDays],pet[1:calibDays],bestParameterDRV)
plot(Qcalc,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for calibration period")
lines(qobs[1:calibDays],type ="l",col="black",lty=2)
legend(1500,50,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
Validation period
Mean of the top 20% parameters
top20DrvValuesValid=tail(sort(drvVal),20)
top20DrvPositionsValid=which(drvVal %in% top20DrvValuesValid)
#top20NseValues
optimumParamsDRVValid=params[top20DrvPositionsValid,]
bestParameterDRVValid=colMeans(optimumParamsDRVValid)
bestParameterDRVValid## S K BFI
## 79.8262720 0.4561368 0.4329308Boxplot for the top 20% parameters
optimumParamsDRVValid=as.data.frame(optimumParamsDRVValid)
boxplot(optimumParamsDRVValid$S,xlab=c("S"))
boxplot(optimumParamsDRVValid[,2],optimumParamsDRVValid[,3],xlab=c("K BFI"))
Simulated X Observed flow
QcalcVald= awbm(p[(calibDays+1):length(qobs)],pet[(calibDays+1):length(qobs)],bestParameterDRVValid)
plot(QcalcVald,type = "l",col="red",ylab="Flow",xlab="Days",main="Observed X Simulated Flow for Validation period")
lines(qobs[(calibDays+1):length(qobs)],type ="l",col="black",lty=2)
legend(1500,52,c("Simulated","Observed"),col=c("red","black"),lty=c(1,2),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
GLOBAL OPTIMUM
Complementary, the PSO algorithm was also applied to the calibration period to define the global optima. To do that the Nash Sutcliffe function was multiplied by (-1) and adjusted to the name FO1 so that finding optima should be converted in a minimization problem of the objective function.
#Nash Sutcliffe Efficiency Function FO1 for the calibration period
FO1<-function (para)
{
pen=0
if(para[1]>=limites[2,1]) pen=100000000000000000000
if(para[1]<=limites[1,1])pen=100000000000000000000
if(para[2]>=limites[2,2]) pen=100000000000000000000
if(para[2]<=limites[1,2]) pen=100000000000000000000
if(para[3]>=limites[2,3]) pen=100000000000000000000
if(para[3]<=limites[1,3]) pen=100000000000000000000
Qcalc<-awbm(p[1:calibDays],pet[1:calibDays],para)
#Calcula a função objetivo
NS1<-(-1)*(1-((sum((qobs[1:calibDays]-Qcalc)^2))/(sum((qobs[1:calibDays]-mean(qobs[1:calibDays]))^2))))+pen
return(ifelse(is.na(NS1), 100000000000000000000, NS1))
}
#Search limits for the parameters
limites<-rbind(c(0,0,0),c(500,1,1))
#Run PSO Function
resultadoCalibracaoPSO<-PSO(FO1,limites)
#Show the global optima paramaters
resultadoCalibracaoPSO[2]## [[1]]
## [1] 484.5697801 0.9550344 0.9953164
f) calculate RR, FI for both observed and model (based on KGE) scenarios for calibration and validation periods.
Runoff Ratio (RR) for the Calibration considering the observed flow
hs_rr(qobs[1:calibDays],p[1:calibDays])## [1] 0.1777375Runoff Ratio (RR) for the Calibration considering the simulated flow
Qcalc=0
Qcalc=awbm(p[1:calibDays],pet[1:calibDays],bestParameterKGE)
hs_rr(Qcalc,p[1:calibDays])## [1] 0.168435Flashiness Index (FI) for the Calibration considering the observed flow
hs_fi(qobs[1:calibDays])## [1] 0.4781125Flashiness Index (FI) for the Calibration considering simulated flow
hs_fi(Qcalc)## [1] 0.09463796
Runoff Ratio (RR) for the Validation considering the observed flow
hs_rr(qobs[(calibDays+1):length(qobs)],p[(calibDays+1):length(qobs)])## [1] 0.1867452Runoff Ratio (RR) for the Validation considering the simulated flow
QcalcVald=0
QcalcVald=awbm(p[(calibDays+1):length(qobs)],pet[(calibDays+1):length(qobs)],bestParameterKGEValid)
hs_rr(QcalcVald,p[(calibDays+1):length(qobs)])## [1] 0.283569Flashiness Index (FI) for the Validation considering the observed flow
hs_fi(qobs[(calibDays+1):length(qobs)])## [1] 0.4516941Flashiness Index (FI) for the Validation considering simulated flow
hs_fi(QcalcVald)## [1] 0.1033999