Numerical instability in density_h20

[davidorme]

Sorry about these bugreps! I'm chasing down some numerical warnings in my python implementation that cropped up using real world inputs with more extreme ranges in test inputs.

There is particular behaviour in the equations for density_h20 that create aberrant results. Specifically, temperatures between about -44 and -46 show extreme numerical instability. The code below shows how the mean and sd of the predicted density go haywire for small changes in temperature in this region:

t <- seq(-88, 58, length=1001)
r <- runif(100, -0.5,0.5)
sd_t <- sapply(t, function(xx) sd(density_h2o(xx + r, p=101325)))
plot(sd_t ~ t, type='l')
rho_t <- sapply(t, function(xx) mean(density_h2o(xx + r, p=101325)))
plot(rho_t ~ t, type='l')

This isn't really a flaw in the equation - it is just being extrapolated way outside the intended range of use and happens to have behaviour that causes all sorts of downstream numerical issues.

I think the best fix here is probably a more general constraint on temperature inputs to the pmodel - so that no predictions are made for sub-zero temperatures, given how little the freezing point changes under environmental temperature and pressures.

[stineb]

Hi David,
Your inputs and raised issues are highly appreciated. Keep coming! We had a minimum temperature if in there earlier and should obviously revert to it.
In this case, the function density_h2o() applies to water, that is H2O in liquid phase. So we may just limit the function to treating only temperatures above 0 degrees C to begin with. Hence also rpmodel() the same way.
thanks
b

[davidorme]

You're welcome! I think it is true that nearly all of the functions with temperature as an input should only make predictions for temperatures above zero. It would be irritating (and slowing) to have a filter in all of them - I can't immediately think of an efficient mechanism to enforce that limit universally.