In valeriyafilipova/StochasticPQRUT: Functions for Stochastic PQRUT

---

title: "Stochastic PQRUT" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Stochastic PQRUT} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8}

---

This R package implements the Stochastic PQRUT model, which performs flood frequency analysis.The model simulates Precipitation, temperature and initial conditions and uses the deterministic rainfall-runoff model PQRUT to obtain the final discharge values. The return periods are estimated using plotting positions. The analysis is performed using seasonal values but currently the seasonal split is hard coded to season 1= "9-11",season 2="12-2",season 3="3-5",season 4="6-8".

library(StochasticPQRUT)

  data(Qd)
  data(P)
    data(SWE)
  data(sl)  
    data(p1)
  data(h1)
    data(Qh)
       data(Td)

1. Load the library,set seed and define number of simulations

library(StochasticPQRUT)
set.seed(567)
  Nsim=10000

1. Find critical duration - the duration of the storm event that results in maximum discharge. For this Precipitation and Discharge with daily timestep are used. Sample data for Hørte (stN 16.193) is included in this package.

  data(Qd);data(P)
   head(Qd); head(P)

POT events are extracted over the quantile, specified by qtT and seperated by number of days intEvent (using extRemes). This is performed for predifined seasons, as explained earlier. Then these events are correlated to the Precipition on 1,2,3 days before the peak.

a=criticalduration(Q=Qd,P=P,qtT=0.9,PDFplots=FALSE,intEvent=7,writeResults=FALSE)

The function returns a list, containing the critical duration and dataframes of the extracted POT.

1. Split the precipitation values into seasons and fit GPD or Exponential distribution to POT events using the extRemes package. The parameters qFtset,qWtset,qWtset,qSptset are the threshold quantiles that are used to extract the POT events. These events do not coincide with the events, extracted in step 1.

d=rainPOT(datarainfall=P,pathmain,durt=24,qFtset=0.98,qWtset=0.99,qSptset=0.99,qStset=0.98,distfunc="GP",
          writeResults=FALSE,PDFplots=FALSE)

4.Extract initial conditions (Snow water equivalent, soil moisture deficit and initial discharge) for the POT flood events, identified in step 1. In this case, the data has been generated using the hydrological model DDD.

  data(sl);data(SWE);data(Td)
   head(sl); head(SWE); head(Td)

The function returns a list of initial conditions for the given seasons.

b=initcond(Qobs=Qd,sl=sl,
           POTF=a$POTF,POTW=a$POTW,POTSp=a$POTSp,POTS=a$POTS,SWE=SWE,Td=Td,durt=a$d,PDFplots=FALSE,writeResults=FALSE)

5.Extract hyetographs for flood events. This function matches the dates for the POT events, extracted in step 3 to hourly precipitation dataset. It is assumed that the daily data is recorded at 6 am.

#load data data(p1)
head(p1)
h=stormp(prpFt=d$pF,prpWt=d$PW,prpSpt=d$PSp,prpSt=d$PS,p1 = p1,durt=24,writeResults=FALSE,PDFplots=FALSE)

1. Extract temperature sequence. Similarly to 4, the function matches the dates of the POT events to hourly Temperature data

#load data
data(h1)
head(h1)
#error due to NULL values in sequence
h1=temprsim(incondFt=b$ntp1F,incondWt=b$ntp1W,incondSpt=b$ntp1Sp,incondSt=b$ntp1S,durt=24,tm=h1,writeResults=FALSE,PDFplots=FALSE)

1. Simulate initial conditions. The values are simulated using truncated multivariate normal distribution. This ensures that the simulated values are within the limits of the observed series.

layout(matrix(1:4, nrow=2))
gincon=initconditions(incondFt=b$ntp1F,incondWt=b$ntp1W,incondSpt=b$ntp1Sp,incondSt=b$ntp1S, Nsim=Nsim,durt=a$d,writeResults=FALSE,PDFplots=FALSE)
#plot for season "09-11"
print(gincon$plot1)

8.Simulate precipitation and temperature values and disaggregate to 1 hour. This function uses the parameters of the GPD, obtained in step 2 and the results of step 5 and step 6.

 h=simulateP(Nsim=Nsim,durt=24,distfunc="GP",hyet=h,TempSeq=h1,Pexp=d$par,writeResults=FALSE,PDFplots=FALSE)

9.Calibrate the hydrological model PQRUT for selected events. In this case optim is used , however a better option is to use ppso package.

library(hydroGOF)
#Prepare data
b1=b[[1]]
p1$date=as.Date(p1$date,format="%d/%m/%Y")
b2=match(as.character(b1$date+1),as.character(p1$date))
data(h1)
#extract Precipitation and temperature values
p_ev=lapply(1:length(b2), function(x) p1[(b2[x]-42):(b2[x]+5),] )

b3=match(as.character(b1$date+1),as.character(as.Date(Qh$Date)))
q_ev=lapply(1:length(b2), function(x) Qh[(b3[x]-42):(b3[x]+5),] )
h1$date=as.Date(h1$date,format="%d/%m/%Y")
b4=match(as.character(b1$date+1),as.character(as.Date(h1$date)))
t_ev=lapply(1:length(b2), function(x) h1[(b4[x]-42):(b4[x]+5),] )
PT_ev=lapply(1:length(b2), function(x) cbind(p_ev[[x]][3],t_ev[[x]][3]))
# 
bn=data.frame(b1[17,2:4])

#then take the medium of all values
  #in this case the model performance is poor
fn1=function(par,...){
qsim=PQRUT(int1=bn,tm1 = PT_ev[[17]],kd=3,durt =48,param.station =c(0.1,0.04,25),Area1 = 157, snpSpt=0.3,ttsnow=-3,Tmax=0.5,Tmin=0.5 )
KGE1=-KGE(qsim[[1]],q_ev[[17]][,2])

return(KGE1)
}
par1=optim(par = c(0.1,0.04,25),fn = fn1,method = "L-BFGS-B",lower=c(0.01,0.01,1),upper = c(0.5,0.5,100))
# 
# 
qsim=PQRUT(int1=bn,tm1 = PT_ev[[17]],kd=3,durt =48,param.station =par1$par,Area1 = 157, snpSpt=0.3,ttsnow=-3,Tmax=0.5,Tmin=0.5)
plot(q_ev[[1]][,2],type="l")
points(qsim[[1]],type="l",col="green")

1. Simulate discharge values. The model PQRUT is used to simulate discharge values. The model was calibrated using method described in Filipova 2016.

gFT= hydrolsim(seasn="Ft",ncl=2,param.station=c(0.1138291,0.028429,19.85680),Nsim=Nsim,int1=gincon$SWEFt,Pt=as.matrix(h[[1]]),E=h[[5]],durt=24,Area1=157,kd=2.8,snpSpt=0.3,ttsnow=-1,Tmax=0.5,Tmin=0.5,writeResults=FALSE,PDFplots=FALSE)
gWT= hydrolsim(seasn="Ft",ncl=2,param.station=c(0.1138291,0.028429,19.85680),Nsim=Nsim,int1=gincon$SWEFt,Pt=as.matrix(h[[2]]),E=h[[6]],durt=24,Area1=157,kd=2.8,snpSpt=0.3,ttsnow=-1,Tmax=0.5,Tmin=0.5,writeResults=FALSE,PDFplots=FALSE)
gSpT= hydrolsim(seasn="Ft",ncl=2,param.station=c(0.1138291,0.028429,19.85680),Nsim=Nsim,int1=gincon$SWEFt,Pt=as.matrix(h[[3]]),E=h[[7]],durt=24,Area1=157,kd=2.8,snpSpt=0.3,ttsnow=-1,Tmax=0.5,Tmin=0.5,writeResults=FALSE,PDFplots=FALSE)
gST= hydrolsim(seasn="Ft",ncl=2,param.station=c(0.1138291,0.028429,19.85680),Nsim=Nsim,int1=gincon$SWEFt,Pt=as.matrix(h[[4]]),E=h[[8]],durt=24,Area1=157,kd=2.8,snpSpt=0.3,ttsnow=-1,Tmax=0.5,Tmin=0.5,writeResults=FALSE,PDFplots=FALSE)

1. Calculate the return levels. This function uses plotting positions to estimate the return periods and returns a list of return levels.

gsim=list(Ft=gFT,Wt=gWT,Spt=gSpT,St=gST)
  d1=d[2:5]
  m1=annualfreqplot(qtT=0.9,Nsim=Nsim,durt=24,qobs=Qh,Pint=h[10:13],incond=gincon,nsy=d1,Qsim=gsim,writeResults=FALSE,PDFplots=FALSE)

valeriyafilipova/StochasticPQRUT documentation built on May 26, 2019, 5:34 a.m.