2015-06-18 Bias Checks

[lawinslow]

2015-06-18:Run future models, look at bias in past

We know there is bias in the driver input. We don't know how important it is. Here is an example of GENMOM annual average wind speed compared to NLDAS

Using GENMOM drivers, how does bias look in the water temp estimates? And can we remove some of the bias and improve the temp outputs?

Bias in water temp for model run between 1979 and 1999.

Model	Bias (deg C, mean(Observed - Modeled))
NLDAS	0.1246422
GENMOM	0.9924086
GENMOM Debiased	0.7389677

That doesn't get us far. GENMOM keeps the lakes a bit cold.

Simple Debias algorithm

    #Debias wind with a multiplier
    wnd_multip = 1/(mean(overlap$WindSpeed)/mean(overlap$nldas_WindSpeed))
    dbiased$WindSpeed = dbiased$WindSpeed*wnd_multip

    #debias airT with linear model
    air_lm = lm(nldas_AirTemp~AirTemp, overlap)

    dbiased$AirTemp = predict(air_lm, dbiased)

Only applied to air temp and wind speed. It was a lot less clear how to remove bias in the other drivers.

Shortwave NLDAS vs GENMOM

RelHum NLDAS vs GENMOM

LongWave NLDAS vs GENMOM

Rain and Snow were right out.

Can we improve the debias to the point we get near NLDAS bias?

This paper does a simple offset bias correction to some of their data. Maybe we try that. Hmm.

Let's add a offset debias to ShortWave.

  #Debias wind with a multiplier
  wnd_multip = 1/(mean(overlap$WindSpeed)/mean(overlap$nldas_WindSpeed))
  dbiased$WindSpeed = dbiased$WindSpeed*wnd_multip

  #debias airT with linear model
  air_lm = lm(nldas_AirTemp~AirTemp, overlap)

  dbiased$AirTemp = predict(air_lm, dbiased)

  #lets try shortwave if requested
  if(shortwave){
    dbiased$ShortWave = dbiased$ShortWave + (mean(overlap$nldas_ShortWave) - mean(overlap$ShortWave))
  }

Model	Bias (deg C, mean(Observed - Modeled))
NLDAS	0.1246422
GENMOM	0.9924086
GENMOM Debiased w/ shortwave	0.5396926

Improves things fairly substantially. Hmmm.

Also ran full GENMOM run, all past and all future

The results are pretty big. 3GB Rdata file.

[jiwalker-usgs]

Random question, have you tried using compress="gzip|bzip2" on the save? Might be able to get some good compression on that file.

[lawinslow]

Oh yeah, that is compressed already.

[jordansread]

@lawinslow is it worth starting at the point of http://onlinelibrary.wiley.com/doi/10.1111/gcb.12262/epdf with de-biasing applied relative to PRISM air temperatures? Using this to de-bias air temperatures, and then to apply that same simple offset correction to longwave (not doing a separate de-bias fit to LW, instead using the SB equation with a modified T^4 input) would be robust and pretty straightforward. To do this, you would need to assume that changes in RH and SW are negligible (and deal separately w/ the wind scaling).

Also, the same method could be applied to all three downscaled models, which are likely all important to bring along along for the simulations for comparison. In the AL-CHOKHACHY et al. 2013 paper, the 3 models had pretty different warming rates when applied to streams (also differences in seasonal timing of warming).

I can help w/ the PRISM de-biasing if you would like to give this a try.

[lawinslow]

Yeah, I was mainly trying to explore if de-biasing was worth the effort and had NLDAS at my fingertips. It seems like the answer is probably "yes".

What would work best is if you could get the PRISM airT data for the WiLMA lakes

Also, dove into Al-Chokhachy et al supplement. Take a look at this.

All the models show a lot of warming in the fall. GENMOM doesn't show much throughout the rest of the year (at least for their sites). Interestingly enough, it isn't as clear in our data. Here is weekly aggregated sens slope.

Hmmm, I'll have to dig into how they made their figure. I'm also very curious now what ECHAM and CM2.0 water temp model output will look like.

[lawinslow]

@jread-usgs Do you want to get together the PRISM data for the WiLMA lakes? It's easy to add a debias step in if I have the data available by ID.

[jordansread]

Yep. Is the SB WFS still up?

[lawinslow]

Yup, I think these are the droids you're looking for.

[lawinslow]

[jordansread]

Perfect, Luke

[jordansread]

@lawinslow I added #65, which includes data and code (code is in demo/airT_from_prism.R). the GDP output is there and I think you have a loader function, otherwise geoknife::parseTimeseries(file) will do it.

That Al-Chokhachy et al supplement looks like an R plot...wonder if they share code... GFDL CM2.0 looks quite hot in the summer.

[lawinslow]

@jread-usgs To calculate LW in from AirT/RH, are you ok with the Crawford and Duchon (1999) approach used in LakeMetabolizer?

[lawinslow]

To do this correction, use simplified SB correction. Crawford and Duchon are mainly concerned with a more nuanced bulk emissivity calc from cloud cover. We don't have that info, but we can back-calculate bulk emissivity.

Calculate bulk emissivity from current met data. builk_emis = biased$LongWave/(5.67E-8 * (biased$AirTemp + 273.13)^4)

Then bias-correct the temp using dbiased$AirTemp = biased$AirTemp + (mean(overlap$nldas_AirTemp) - mean(overlap$AirTemp))

The re-calc LW using bulk emissivity value and debiased air temp

dbiased$LongWave = builk_emis * 5.67E-8 * (dbiased$AirTemp + 273.13)^4

Sneak edit: Fixed blackbody radiation eqn.

[lawinslow]

EDIT: I originally used NLDAS instead of PRISM. Updated bias numbers. Different, but not really better.

@jread-usgs Using the below bias correction based on airT offset model and LW corrected to match, bias actually increased using GENMOM.

Sim	Bias	RMSE
GENMOM Uncorrected	0.9807225	3.38653
airT Offset and LW corrected	-0.9492572	4.598

I am double-checking the GENMOM uncorrected number now.

names(obs) = paste0('obs_', names(obs))
  names(obs)[1] = 'time'

  #prism is monthly avg, so avg dscale to monthly first

  mon_dscale = group_by(dscale, year = year(as.POSIXct(time)), month=month(as.POSIXct(time))) %>% 
    summarise(AirTemp=mean(AirTemp)) 
  mon_dscale$time = as.character(ISOdate(mon_dscale$year, mon_dscale$month, 1, hour = 0))

  overlap = merge(obs, mon_dscale , by='time')

  builk_emis = dscale$LongWave/(5.67E-8 * (dscale$AirTemp + 273.13)^4)

  #debias airT with offset model
  dbiased$AirTemp = dbiased$AirTemp + (mean(overlap$obs_AirTemp) - mean(overlap$AirTemp))

  dbiased$LongWave = builk_emis * 5.67E-8 * (dbiased$AirTemp + 273.13)^4

[lawinslow]

PRISM data does weird things to our model

[jordansread]

hmm. long live nldas