We are using WRF-Hydro to predict streamflow for Fourmile Creek at the Orodell USGS gage for the 2013 snowmelt period. We ran WRF-Hydro with NoahMP as the LSM for a 3-year spinup period and then did a daily output run for 5 months starting March 1, 2013. We want to evaluate the ET fluxes at a local Ameriflux tower site.
Load the rwrfhydro package.
library("rwrfhydro")
## To check rwrfhydro updates run: CheckForUpdates()
Set a data path to the Fourmile Creek test case.
dataPath <- '~/wrfHydroTestCases/Fourmile_Creek'
Import LSM water fluxes at the cell where the Ameriflux tower is located. Since we do NOT want a basin-wide average, but instead a time series of values for a specified cell, we will use the GetMultiNcdf tool directly. First, we do a bit of setup.
Generate the list of filenames for the GetMultiNcdf tool by providing the pathname to the model OUTPUT directory and the filename pattern to match (in this case we want the LDASOUT files).
<<<<<<< HEADlsmFiles <- list.files(path=paste0(dataPath, '/RUN.RTTESTS/OUTPUT_ALLRT_DAILY'),
pattern=glob2rx('*LDASOUT*'), full.names=TRUE)
=======
lsmFiles <- list.files(path=paste0(dataPath, '/RUN.RTTESTS/OUTPUT_ALLRT_DAILY'), pattern=glob2rx('*LDASOUT*'), full.names=TRUE)
>>>>>>> upstream/devBranch
flList <- list(lsm=lsmFiles)
Now setup which LDASOUT variables we want to bring in. We will be working with the accumulated ET fluxes at our daily time step.
lsmVars <- list('ACCECAN','ACCEDIR','ACCETRAN','ACCPRCP')
names(lsmVars) <- lsmVars
varList <- list( lsm = lsmVars )
Setup the cell you want to extract and the stat you want to calculate. We are only running a single cell (i=2, j=3, time=1) so we specify the same start and end cell and use “CalcMeanMinrm” as the stat. “CalcMeanMinrm” is a utility function that calculates a mean with an enforced minimum valid value (-1e+30). This is useful to remove intended “no data” values that were not caught by NetCDF checks.
niw <- list(start=c(2,3,1), end=c(2,3,1), stat='CalcMeanMinrm')
Match the variables to the desired cell and stat (here we are using the same for all variables).
lsmInds <- list( niw, niw, niw, niw )
names(lsmInds) <- names(lsmVars)
indList <- list( lsm=lsmInds )
Run the tool. To speed things up, we will run this over 8 cores using the doMC package.
library(doMC)
<<<<<<< HEAD
## Warning: package 'doMC' was built under R version 3.1.2
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 3.1.2
##
## Attaching package: 'foreach'
##
## The following object is masked from 'package:chron':
##
## times
##
## Loading required package: iterators
## Warning: package 'iterators' was built under R version 3.1.2
## Loading required package: parallel
=======
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
>>>>>>> upstream/devBranch
registerDoMC(8)
modDf <- GetMultiNcdf(file=flList,var=varList, ind=indList, parallel=TRUE)
Now we can “flatten” the output to make it work better with certain tools.
modLsm.allrt.niw <- ReshapeMultiNcdf(modDf)
rm(modDf)
Some of the LDASOUT fields we exported are accumulated totals, so we need to break them out into individual (per time step) fluxes.
modLsm.allrt.niw <- CalcNoahmpFluxes(modLsm.allrt.niw)
We want to examine total ET, so we need to calculate a sum column for all 3 evaporation and transpiration fluxes (in units of “mm”). We then use the density of water (1mm H2O * m2 = 1 kg), the latent heat of vaporization (assuming evaporation and not sublimation during this summer period, LHVAP = 2.5106e+6 J/kg), and the time step (1 day = 86400 sec) to calculate total daily latent heat flux in W/m2.
modLsm.allrt.niw$ET_mm <- with(modLsm.allrt.niw, DEL_ACCECAN + DEL_ACCETRAN + DEL_ACCEDIR)
modLsm.allrt.niw$LE <- with(modLsm.allrt.niw, ET_mm * 2.5106e+6 / 86400)
Import Ameriflux station data for Niwot Ridge (Level 2 Standardized). We specify “America/Denver” as the time zone (Ameriflux timestamps are in local time).
<<<<<<< HEADobsFlux30min.usnr1 <- ReadAmerifluxCSV(paste0(dataPath, '/OBS/AMF_USNR1_2013_L2_GF_V008.csv'),
"America/Denver")
Our LDASOUT time step was 1 day, so we need to aggregate the observations to a UTC daily time step for comparison.
obsFlux1d.usnr1 <- aggregate(obsFlux30min.usnr1[names(obsFlux30min.usnr1) != c("POSIXct")],
by = list(as.POSIXct(round(as.POSIXct(format(obsFlux30min.usnr1$POSIXct, tz="UTC"), tz="UTC"),
"days"))),
=======
obsFlux30min.usnr1 <- ReadAmerifluxCSV(paste0(dataPath, '/OBS/AMF_USNR1_2013_L2_GF_V008.csv'), "America/Denver")
Our LDASOUT time step was 1 day, so we need to aggregate the observations to a UTC daily time step for comparison.
obsFlux1d.usnr1 <- aggregate(obsFlux30min.usnr1[names(obsFlux30min.usnr1) != c("POSIXct")],
by = list(as.POSIXct(round(as.POSIXct(format(obsFlux30min.usnr1$POSIXct, tz="UTC"), tz="UTC"), "days"))),
>>>>>>> upstream/devBranch
CalcMeanNarm)
names(obsFlux1d.usnr1)[1] <- "POSIXct"
PlotFluxCompare(obsFlux1d.usnr1, "LE", modLsm.allrt.niw, "LE",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="ET: Niwot Ridge")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
CalcModPerf(modLsm.allrt.niw, obsFlux1d.usnr1, flxCol.mod="LE", flxCol.obs="LE")
<<<<<<< HEAD
## Warning: package 'pander' was built under R version 3.1.2
=======
>>>>>>> upstream/devBranch
n | nse | nselog | cor | rmse | rmsenorm | bias | mae | errcom | errmaxt | errfdc | |
---|---|---|---|---|---|---|---|---|---|---|---|
units | count | unitless | unitless | unitless | flux units | unitless | percent | flux units | hours|days | hours|days | flux units |
t | 152 | -0.44 | -1.14 | 0.43 | 36.56 | 26.06 | -0.7 | 27.01 | NA | NA | 0.07 |
daily | 152 | -0.44 | -1.14 | 0.43 | 36.56 | 26.06 | -0.7 | 27.01 | 0 | 0 | 0.07 |
monthly | 5 | -0.06 | -3.94 | 0.85 | 19.53 | 39.32 | -1.2 | 16.3 | -0.4 | 5 | NA |
yearly | 1 | NA | NA | NA | NA | NA | -0.7 | 0.43 | 4 | 74 | NA |
max10 | 16 | -11.34 | -62.95 | -0.28 | 42.47 | 90.8 | -20.6 | 25.8 | NA | NA | NA |
min10 | 16 | -54.84 | -6.96 | 0.07 | 32.72 | 208.25 | 94.4 | 21.75 | NA | NA | NA |
Let’s take a closer look at the month of May where our ET fluxes are significantly higher than the observed. We ran an hourly output run across the full month.
We will follow a similar setup as the daily fluxes, but we will grab the instantaneous water and energy fluxes from the LSM output and the wind inputs from the forcing data.
Provide the new pathname to the hourly output. We will limit the LDASIN data to just May 2013 since the FORCING directory contains a full year’s worth of data.
<<<<<<< HEADlsmFiles <- list.files(path=paste0(dataPath,'/RUN.FLUXCOMP/OUTPUT_ALLRT_MAY13_HOURLY'),
pattern=glob2rx('*LDASOUT*'), full.names=TRUE)
forcFiles <- list.files(path=paste0(dataPath, '/FORCING'), pattern=glob2rx('201305*LDASIN*'),
full.names=TRUE)
flList <- list(lsm=lsmFiles, forc=forcFiles)
Specify the variables for output
lsmVars <- list('ECAN', 'EDIR', 'ETRAN', 'RAINRATE', 'SWFORC', 'LWFORC',
'FIRA', 'FSA', 'LH', 'HFX', 'GRDFLX')
=======
lsmFiles <- list.files(path=paste0(dataPath,'/RUN.FLUXCOMP/OUTPUT_ALLRT_MAY13_HOURLY'), pattern=glob2rx('*LDASOUT*'), full.names=TRUE)
forcFiles <- list.files(path=paste0(dataPath, '/FORCING'), pattern=glob2rx('201305*LDASIN*'), full.names=TRUE)
flList <- list(lsm=lsmFiles, forc=forcFiles)
Specify the variables for output
lsmVars <- list('ECAN', 'EDIR', 'ETRAN', 'RAINRATE', 'SWFORC', 'LWFORC', 'FIRA', 'FSA', 'LH', 'HFX', 'GRDFLX')
>>>>>>> upstream/devBranch
names(lsmVars) <- lsmVars
forcVars <- list('T2D', 'Q2D', 'U2D', 'V2D', 'PSFC')
names(forcVars) <- forcVars
varList <- list(lsm=lsmVars, forc=forcVars)
Specify the indices and stats. We’ll use the same “niw” location and stat as above.
lsmInds <- list( niw, niw, niw, niw, niw, niw, niw, niw, niw, niw, niw )
names(lsmInds) <- names(lsmVars)
forcInds <- list( niw, niw, niw, niw, niw )
names(forcInds) <- names(forcVars)
indList <- list(lsm=lsmInds, forc=forcInds)
Run the tool and flatten results.
modDf <- GetMultiNcdf(file=flList,var=varList, ind=indList, parallel=TRUE)
modLsm.allrt.may.niw <- ReshapeMultiNcdf(subset(modDf, modDf$fileGroup=="lsm"))
modForc.allrt.may.niw <- ReshapeMultiNcdf(subset(modDf, modDf$fileGroup=="forc"))
rm(modDf)
Latent heat:
PlotFluxCompare(obsFlux30min.usnr1, "LE", modLsm.allrt.may.niw, "LH",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Latent Heat: Niwot Ridge (+ to atmosphere)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Sensible heat:
PlotFluxCompare(obsFlux30min.usnr1, "H", modLsm.allrt.may.niw, "HFX",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Sensible Heat: Niwot Ridge (+ to atmosphere)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Let’s look at total net latent and sensible heat fluxes.
obsFlux30min.usnr1$LE_H <- with(obsFlux30min.usnr1, LE + H)
modLsm.allrt.may.niw$LE_H <- with(modLsm.allrt.may.niw, LH + HFX)
PlotFluxCompare(obsFlux30min.usnr1, "LE_H", modLsm.allrt.may.niw, "LE_H",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Net Latent & Sensible Heat: Niwot Ridge (+ to atmosphere)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
We see that the model is overestimating peak LE fluxes and underestimating SH fluxes, but the net heat exchange is fairly well represented. So we now examine our radiation fluxes to see if we are over-estimating available energy.
We do a few calculations to get net radiation and the overall energy balance.
modLsm.allrt.may.niw$LWnet <- with(modLsm.allrt.may.niw, -FIRA)
modLsm.allrt.may.niw$Rnet <- with(modLsm.allrt.may.niw, FSA - FIRA)
modLsm.allrt.may.niw$EnResid <- with(modLsm.allrt.may.niw, Rnet - LE_H - GRDFLX)
obsFlux30min.usnr1$RgNet <- with(obsFlux30min.usnr1, Rg - RgOut)
obsFlux30min.usnr1$RglNet <- with(obsFlux30min.usnr1, Rgl - RglOut)
obsFlux30min.usnr1$EnResid <- with(obsFlux30min.usnr1,
(Rg - RgOut) + (Rgl - RglOut) - LE - H - FG)
Plot & compare the time series of incoming shortwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "Rg", modLsm.allrt.may.niw, "SWFORC",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Incoming Shortwave Radiation: Niwot Ridge (+ to surface)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Plot & compare the time series of net shortwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "RgNet", modLsm.allrt.may.niw, "FSA",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Net Shortwave Radiation: Niwot Ridge (+ to surface)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Plot & compare the time series of net longwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "RglNet", modLsm.allrt.may.niw, "LWnet",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Net Longwave Radiation: Niwot Ridge (+ to surface)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Plot & compare the time series of total net radiation
PlotFluxCompare(obsFlux30min.usnr1, "Rn", modLsm.allrt.may.niw, "Rnet",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Total Net Radiation: Niwot Ridge (+ to surface)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
Plot & compare the time series of energy residuals. The Ameriflux observations for ground heat flux do not cover this period, so we compare the observation residuals against the modelled ground heat flux.
PlotFluxCompare(obsFlux30min.usnr1, "EnResid", modLsm.allrt.may.niw, "GRDFLX",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Energy Residual: Niwot Ridge (+ to surface)")
<<<<<<< HEAD
=======
>>>>>>> upstream/devBranch
We do a few calculations to get the net wind speed.
modForc.allrt.may.niw$Wind <- with(modForc.allrt.may.niw, sqrt(U2D^2 + V2D^2))
<<<<<<< HEAD
Plot & compare the time series of wind speeds.
PlotFluxCompare(subset(obsFlux30min.usnr1,
POSIXct >= as.POSIXct("2013-05-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC") &
POSIXct < as.POSIXct("2013-06-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")),
"WS", modForc.allrt.may.niw, "Wind",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Wind Speed: Niwot Ridge")
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
Now we convert specific humidity to relative humidity to compare to the flow tower observations.
modForc.allrt.may.niw$RH <- with(modForc.allrt.may.niw,
0.263*PSFC*Q2D*(exp((17.67*(T2D-273.16))/(T2D-29.65)))^(-1))
Plot & compare the time series of relative humidity.
PlotFluxCompare(subset(obsFlux30min.usnr1,
POSIXct >= as.POSIXct("2013-05-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC") &
POSIXct<as.POSIXct("2013-06-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")),
"RH", modForc.allrt.may.niw, "RH",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Relative Humidity: Niwot Ridge")
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
=======
Plot & compare the time series of wind speeds. (This throws a warning, however the POSIXct timezones are handled correctly.)
PlotFluxCompare(subset(obsFlux30min.usnr1, POSIXct>=as.POSIXct("2013-05-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")
& POSIXct<as.POSIXct("2013-06-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")), "WS", modForc.allrt.may.niw, "Wind",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Wind Speed: Niwot Ridge")
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
Now we convert specific humidity to relative humidity to compare to the flow tower observations.
modForc.allrt.may.niw$RH <- with(modForc.allrt.may.niw, 0.263*PSFC*Q2D*(exp((17.67*(T2D-273.16))/(T2D-29.65)))^(-1))
Plot & compare the time series of relative humidity. (Again, timzones are correctly handled.)
PlotFluxCompare(subset(obsFlux30min.usnr1, POSIXct>=as.POSIXct("2013-05-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")
& POSIXct<as.POSIXct("2013-06-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC")), "RH", modForc.allrt.may.niw, "RH",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Relative Humidity: Niwot Ridge")
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
## Warning in check_tzones(e1, e2): 'tzone' attributes are inconsistent
>>>>>>> upstream/devBranch