Problem with flow normalized values
khhydro opened this issue · 4 comments
Hello,
I am a graduate student looking to run EGRET and WRTDS as part of my thesis research to estimate sediment loading trends over a 26 year period based upon 1,965 observations of TSS. I have run into an issue where I am seeing very large flow normalized concentrations (~1x10^13 mg/l) and fluxes and I am having trouble determining what is going wrong. Below is my code and attached are my input files and results table. Any insight would be greatly appreciated and thanks in advance for your patience with my novice-ness!
Thank you!
#import synthetic daily discharge record
filePath <- ".../csv"
RedDaily <- readUserDaily(filePath, "RedDaily.csv")
#import TSS observations
RedSample <- readUserSample(filePath,"RedTSS.csv")
#import info
RedINFO <- readUserInfo(filePath, "RedInfo.csv")
#merge eList
eListRed <- mergeReport(RedINFO,RedDaily,RedSample)
#Discharge Record is 9496 days long, which is 26 years
#First day of the discharge record is 1985-10-01 and last day is 2011-09-30
#The water quality record has 1965 samples
#The first sample is from 1986-03-27 and the last sample is from 2010-07-15
#Discharge: Minimum, mean and maximum 0.102 1.78 18.1
#Concentration: Minimum, mean and maximum 0.38 23 2400
#Percentage of the sample values that are censored is 0 %
#Red WRTDS
eListRed <- modelEstimation(eListRed)
plotConcHist(eListRed)
I'll take a look at this afternoon!
I was looking at EGRET to see if anything was new and saw this post. I have an interest in models that go wrong, so I took a look at it.
First, I am an advocate for model calibration and validdation, so the first things that I did was to look at some of the plotting functions that address at least the linearity assumption that is critical for calibration:
plotResidPred(eListRed)
plotResidQ(eListRed)
plotResidTime(eListRed)
boxResidMonth(eListRed)
The first two show some interesting patterns, but nothing suspiscious. The residuals over time show typical variation above and below the 0 line over years, but what stood out was the grouping pattern in the samples. The boxplot confirmed the grouping pattern--samples were collected only in March through July with a few in August. That means that predictions are valid only for those months. There just is not enough seasonal information to make predictions outside of those months. That would be true no matter you used for load estimation.
That does solve the issue that you observed. The paStart value in the INFO component can be set to 3, March, and the paLong value can be set to 5 (months). The plotConcHiist plot looks much better with those settings. If you look at the monthly data, you will see that the months of November and December are those with the most crazy values.
But you should also look at the available validation data:
errorStats(eListRed)
fluxBiasStat(eListRed$Sample)
The first just gives a description of the expected variability, which is very large, but my experience is that sediment can have very large variability relative to other constituents. The second gives an estimate of bias. The reported value, about -0.27 indicates the the observed values are about 27 % larger than the estimated values, or that the estimated value are about 21 % less than the observed values. Either way, that is very large bias for load estimation and you may not trust your predicted values or the trend you observe. The plotConcQSmooth and plotConcTimeSmooth functions can be used to fine tune the windowQ and windowY arguments for WRTDS. I did not try that.
Yup! WRTDS doesn't use 0 concentrations, so there were 7531 rows of data that were deleted. So as Dave mentioned, I think the only thing you can do with this data set is to set the period of interest to be March through July.
library(EGRET)
filePath <- here::here()
RedDaily <- readUserDaily(filePath, "RedDaily.csv")
RedSample <- readUserSample(filePath,"RedTSS.csv")
RedINFO <- readUserInfo(filePath, "RedInfo.csv")
eListRed <- mergeReport(RedINFO,RedDaily,RedSample)
eListRed <- modelEstimation(eListRed)
plot(eListRed)
eListRed <- setPA(eListRed, paStart = 3, paLong = 5)
plotConcHist(eListRed)
Thank you both for your help!