Plot & compare the time series of ET fluxes
PlotFluxCompare(obsFlux1d.usnr1, "LE", modLsm.allrt.niw, "LE",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="ET: Niwot Ridge")
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")Set a data path to the Fourmile Creek test case.
dataPath <- '~/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0/'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).
lsmFiles <- list.files(path=paste0(dataPath, '/run.FullRouting'),
pattern=glob2rx('*LDASOUT*'), full.names=TRUE)
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 )We need to figure out where in the domain the Niwot tower sits. We have the Ameriflux site metadata (amfMeta) built into rwrfhydro, so we can look for the site by querying that metadata.
head(amfMeta)## name id_txt state cntry lon
## 1 Aiken US-Akn South Carolina USA -81.5656
## 2 Anaktuvuk River Moderate Burn US-An2 Alaska USA -150.2100
## 3 Anaktuvuk River Severe Burn US-An1 Alaska USA -150.2800
## 4 Anaktuvuk River Unburned US-An3 Alaska USA -150.2700
## 5 ARM SGP Burn US-ArB Oklahoma USA -98.0402
## 6 ARM SGP Control US-ArC Oklahoma USA -98.0401
## lat tz_us tz_glob UTC_offset UTC_DST_offset
## 1 33.3833 America/New_York America/New_York -05:00 -04:00
## 2 68.9500 America/Anchorage America/Anchorage -09:00 -08:00
## 3 68.9900 America/Anchorage America/Anchorage -09:00 -08:00
## 4 68.9300 America/Anchorage America/Anchorage -09:00 -08:00
## 5 35.5497 America/Chicago America/Chicago -06:00 -05:00
## 6 35.5465 America/Chicago America/Chicago -06:00 -05:00
## UTC_offset_int
## 1 -5
## 2 -9
## 3 -9
## 4 -9
## 5 -6
## 6 -6nrMeta <- subset(amfMeta, amfMeta$id_txt=="US-NR1")
nrMeta## name id_txt state cntry lon lat tz_us
## 96 Niwot Ridge US-NR1 Colorado USA -105.5464 40.0329 America/Denver
## tz_glob UTC_offset UTC_DST_offset UTC_offset_int
## 96 America/Denver -07:00 -06:00 -7Now we know the site lat/lon coordinates, so we can use the GetGeogridIndex function to calculate which coarse grid (LSM) cell that point sits in.
geoFile <- paste0(dataPath, "/DOMAIN/geo_em_d01.Fourmile1km.nlcd11.nc")
GetGeogridIndex(xy=nrMeta[,c("lat","lon")], ncfile=geoFile)## we sn
## 1 2 3Based on those indices, we can now setup the cell we want to extract and the stat we 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 3 cores using the doMC package.
library(doMC)
registerDoMC(3)
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 (US-NR1). We have already downloaded the BASE format data for this site from the Ameriflux website (note that data download now requires login/form access, so Ameriflux data is no longer supported as an automatic download within rwrfhydro). To properly format the timestamp within R, we need to know the local timezone (most Ameriflux data uses local timestamps). We can get this from the Ameriflux metadata record that we queried above.
nrMeta## name id_txt state cntry lon lat tz_us
## 96 Niwot Ridge US-NR1 Colorado USA -105.5464 40.0329 America/Denver
## tz_glob UTC_offset UTC_DST_offset UTC_offset_int
## 96 America/Denver -07:00 -06:00 -7We plug this timezone into the ReadAmeriflux function along with the downloaded file (BASE format).
obsFlux30min.usnr1 <- ReadAmeriflux('~/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0/OBS/AMF_US-NR1_BASE_HH_9-1.csv',
tz='America/Denver')Our LDASOUT time step was 1 day, so we need to aggregate the observations to a UTC daily time step for comparison. We apply the aggregations only to numeric variables to avoid a bunch of R warnings.
obsFlux1d.usnr1 <- aggregate(obsFlux30min.usnr1[sapply(obsFlux30min.usnr1, is.numeric)],
by = list(as.POSIXct(round(as.POSIXct(format(obsFlux30min.usnr1$POSIXct, tz="UTC"), tz="UTC"),
"days"))),
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")
CalcModPerf(modLsm.allrt.niw, obsFlux1d.usnr1, flxCol.mod="LE", flxCol.obs="LE")## Warning in log(m + 1e-04): NaNs produced## Warning in log(o + 1e-04): NaNs produced## Warning in log(m + 1e-04): NaNs produced## Warning in log(o + 1e-04): NaNs produced## Warning in log(m + 1e-04): NaNs produced## Warning in log(o + 1e-04): NaNs produced| n | nse | nselog | cor | rmse | rmsenorm | bias | msd | mae | errcom | errmaxt | errfdc | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| units | count | unitless | unitless | unitless | flux units | unitless | percent | flux units | flux units | hours|days | hours|days | flux units |
| t | 303 | 0.12 | 0.04 | 0.58 | 31.7292 | 20.99 | -9.9 | -4.2456 | 22.2789 | NA | NA | -3.98 |
| daily | 303 | 0.12 | 0.04 | 0.58 | 31.7292 | 20.99 | -9.9 | -4.2456 | 22.2789 | 0 | 0 | -3.98 |
| monthly | 11 | 0.62 | 0.11 | 0.9 | 15.1267 | 20.89 | -12 | -5.6063 | 13.148 | -0.55 | 2.55 | NA |
| yearly | 1 | NA | NA | NA | NA | NA | -9.9 | -4.2456 | 4.2456 | 26 | 81 | NA |
| max10 | 31 | -11.98 | -57.67 | -0.21 | 52.8217 | 98.61 | -27.9 | -30.0933 | 35.0188 | NA | NA | NA |
| min10 | 31 | -23.78 | -0.84 | -0.26 | 17.1013 | 132.62 | 322.7 | 10.7022 | 11.7519 | 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.
lsmFiles <- list.files(path=paste0(dataPath,'/run.FluxEval'),
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')
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)")
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)")
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)")
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$SWnet <- with(obsFlux30min.usnr1, SW_IN - SW_OUT)
obsFlux30min.usnr1$LWnet <- with(obsFlux30min.usnr1, LW_IN - LW_OUT)
obsFlux30min.usnr1$EnResid <- with(obsFlux30min.usnr1,
SWnet + LWnet - LE - H - G)Plot & compare the time series of incoming shortwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "SW_IN", modLsm.allrt.may.niw, "SWFORC",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Incoming Shortwave Radiation: Niwot Ridge (+ to surface)")
Plot & compare the time series of net shortwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "SWnet", modLsm.allrt.may.niw, "FSA",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Net Shortwave Radiation: Niwot Ridge (+ to surface)")
Plot & compare the time series of net longwave radiation
PlotFluxCompare(obsFlux30min.usnr1, "LWnet", modLsm.allrt.may.niw, "LWnet",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Net Longwave Radiation: Niwot Ridge (+ to surface)")
Plot & compare the time series of total net radiation
PlotFluxCompare(obsFlux30min.usnr1, "NETRAD", modLsm.allrt.may.niw, "Rnet",
labelObs="Observed Niwot Ameriflux",
labelMod1="All Routing", title="Total Net Radiation: Niwot Ridge (+ to surface)")
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)")
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))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")
Now we convert specific humidity to relative humidity to compare to the flux 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")