Errors in calculation of net solar radiation about crhmcode

[jhs507]

Are you gathering data using the same project as for #10?

[KevinShook]

Pretty much, although I used a SWE threshold of 1 mm.
This one calcualted the radiation components
radiation_model.zip
This one had calculaed the net radiation
Evap_model.zip

[jhs507]

What HRU's are you seeing this difference in?

[KevinShook]

I was looking at HRU 11

[jhs507]

I am not seeing differences in net, netD, or RnD that are outside of our tolerance of 5% or 0.0001.

I checked hru's 2, 3, 11 for the first two years of the simulation.

What years are you examining?

[KevinShook]

Sorry, the examples I plotted above are from 1995.
You can see that in many years, we get very good agreement between both models.
However in some years, chrmcode is underestimating runoff, probably because the evaporation
is too large

[jhs507]

And you believe that the evaporation is too large because our radiation values are too large as well?

[KevinShook]

Probably - haven't checked it yet, but it's the most likely scenario.
I'll do a more complete water balance to be sure.

[jhs507]

Alright I will investigate the overflow, evap, and rad, variables around 1960-1968 and then see if I can figure out what is causing it.

[KevinShook]

I plotted the daily evap for HRU 11 over the whole run period.
Most of the values are very similar, but there are quite a few outliers

The crhmcode mean daily evap is slightly smaller than the Borland evap for HRU 11

      date               Borland          crhmcode     
 Min.   :1960-01-01   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:1971-08-06   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :1983-03-11   Median :0.7606   Median :0.7227  
 Mean   :1983-03-11   Mean   :2.1168   Mean   :2.1071  
 3rd Qu.:1994-10-14   3rd Qu.:4.4962   3rd Qu.:4.4766  
 Max.   :2006-05-20   Max.   :8.0753   Max.   :8.0753

[jhs507]

I have found that the albedo is not matching for hru 11 I assume that this will affect the evap?

[KevinShook]

yes it will - this is important.

[KevinShook]

I got Logan to output the values of some of the internal variables from the Borland netall for HRU 11
for May 21-28, 1995.
netall_HRU11_May_21-28_1995_run2.txt

[KevinShook]

I output the daily evaporation, radiation and albedo for days where the Broaldn evaporation was non zero and > 1.1 x the crhmcode evaporation and vice versa.
The biggest thing that stands out is the number of days where one of the models was setting the albedo to 0.11, i.e. to soil rather than snow.
HRU11_daily_evap_netrad_albedo.zip

I wonder if this might be due to a differing precision when evaluating

    if(SWE[hh] > 5.0) {
      Albedo[hh]    = Albedo_snow[hh];
      winter[hh]    = 1;
    }
    else{
      Albedo[hh]    = Albedo_bare[hh];
      winter[hh]    = 0;
    }

in classalbedo

[jhs507]

That certainly is possible. Anytime we have comparisons with floating points we are going to have issues as comparisons on floating points are never exact.

Especially with moving from single to double precision.

[jhs507]

I implemented the instrumentation that Logan did and analyzed the same time period.
instrumented_results.zip

Results reveal that the instrumented variables are identical save for precision error on that time period but the radiation values do not match.

[jhs507]

I ran a test for the year of 1995 to see if that swe comparison was not being evaluated the same. I used the winter varaiable as a proxy to see what branch was evaluated.

For that year I did not see any times where borland and crhmcode differed I will try with a longer timeline on Monday.

[KevinShook]

I think I might have it.
I output meltflag for both models and crhmcode is seting meltflag to 1 in the interval 1995-06-21 00:00 through
1995-06-26 23:00, while the Borland version doesn't. It's being set by this condition which also sets the
albedo:

   if(SWE[hh] <= 0.0) {
        Albedo[hh] = Albedo_bare[hh];
        winter[hh] = 0;
        meltflag[hh] = 0;
      }

Looking at the SWE values, the crhmcode values are, again, very slightly greater than zero,
while the Borland versions are not.

We need to use a threshold for the SWE test. Why don't we try using the existing SWE threshold?

[KevinShook]

I think we should also check out the ratio of rainfall to precipitation calculated in classobs

  hru_rain[hh] = hru_p[hh]*ratio;
  hru_snow[hh] = hru_p[hh]*(1.0-ratio);

We are getting some very small values of SWE, which really should be taken care of

If we set the ratio to 1.0 and 0.0 for values >= 0.999 and <= 0.001 (or similar values) then we should avoid the roundoff error SWE values

[jhs507]

We need to use a threshold for the SWE test. Why don't we try using the existing SWE threshold?

We can certainly do this. My only concern would be that the parameter as we had specified it was a specific threshold for preventing evaporation. If you feel using the same number here makes sense in terms of the physical process we can certainly apply it.

We will have to get Logan to change the newly added parameter to a shared parameter.

However I think your second point is more important. If we can fix the root of the very small amount of snow problem this would possibly be the best in general. I will see if I can come up with reasonable numbers that have the desired effect.

[KevinShook]

Those are very good points. I agree that it's best to fix the root of the problem and to keep the code as similar as possible between the two versions

[jhs507]

I added some code so I could see what ratio was set to at various points and looked at the first two years of the model.

Perciptiation_Ratio.obs.txt

As you observed there are many points where this is leading to extremely small values for rain or snow. However, simply setting the ratio to 0 or 1 when it approaches that value will not solve our problem. As the data I posted shows there are points where the ratio value is set as 1 but the nature of floating points and the fact that snow is calculated as ppt* (1.0 - ratio) means that we get extremely small amounts of snow.

I will see if I can come up with a more comprehensive fix.

[KevinShook]

I was afraid of that. What would happen if we just set snow = 0.0 or snow = 1.0, instead of setting the ratio?
Would the roundoff error still occur?

[jhs507]

I think that is a more correct approach and what I was going to explore first.

  hru_snow[hh] = 0.0;

  hru_rain[hh] = 0.0;

  if(hru_p[hh] > 0.0) //rain or snow determined by ice bulb ratio

    hru_rain[hh] = hru_p[hh]*ratio;

  hru_snow[hh] = hru_p[hh]*(1.0-ratio);

In a related area the above code is the calucualtion we are discussing and I think it might be having other unintended behavior but I would like to run it by you.

As it is written now when hru_p is less than zero we are guaranteed to have a value of 0 for rain however their can still be a non zero value for snow is this intended?

[KevinShook]

No, both rain and snow should never be less than zero. More bad coding (although IIRC the values are checked somewhere else).

[jhs507]

But more specifically if hru_p is 0.0 then both hru_rain and hru_snow should be 0.0 as well?

Is it also safe to assume that hru_p = hru_rain + hru_snow?

[KevinShook]

Exactly. Conservation of mass requires that rain + snow = p.
Neither can ever be negative.

[jhs507]

Alright if that is the desired outcome that if statement should have both the setting of rain and snow within. As it stands now only rain is within the block of that if. I will fix that. It is unlikely that this causes issues but you may want to have it corrected in the borland code as well.

[KevinShook]

Thanks. I'll see if we get different values between the versions and will let Logan know about the changes.

[jhs507]

I added some logic to consider all precipitation as rain when ratio is close to 1 and all as snow when it is close to zero.

This did have the intended effect of zeroing out snow values where previously they were very small values on the order of 10^-10 ect...

Unfortunately this did not bring evaporation, or radiation into alignment. Albedo, melftlag, and winter are not being set correctly on some steps still. I don't know yet if this is because our new calculation is not catching all near zero snow cases or if it is because my new solution alters things too much from the borland code. I am going to see if I can figure that out.

Here is the data if you want to look yourself.

evap_with_ratio_fix.zip

Here is the altered code for reference.

double epsilon = 0.0001;

 hru_snow[hh] = 0.0;

 hru_rain[hh] = 0.0;

 if (hru_p[hh] > 0.0)  //rain or snow determined by ice bulb ratio
 {
     
     if ( abs(1.0 - ratio) <= epsilon) {
         //ratio is close enough to 1 that it should 100% rain
         hru_rain[hh] = hru_p[hh];
         hru_snow[hh] = 0.0;
     }
     else if (ratio <= epsilon) {
         //ratio is close enough to 0 that it should be 100% snow
         hru_rain[hh] = 0.0;
         hru_snow[hh] = hru_p[hh];
     }
     else {
         //calculate rain and snow normaly
         hru_rain[hh] = hru_p[hh] * ratio;
         hru_snow[hh] = hru_p[hh] * (1.0 - ratio);
     }
 }

[KevinShook]

Thanks. I suspect that we'll still have to use a zero SWE threshold as well.
I you can push the code to your branch, I'll check it out.

[jhs507]

I have pushed my code.

[KevinShook]

Thanks!

[KevinShook]

The albedo seems to be incorrect due to differences in the depth of SWE. This is occurring for depths >> 0.
It could be an error in the albedo decay calculations, or it could be an error in the melt rate.
It also seems that the program is not setting the Albedo to the soil value immediately when the snow depth is 0.
I'll start taking a look at the albedo and melt calculations.
OTOH, there has been some improvement in some of the previously incorrect values for net radiation as you can see for
May 1995.

When I looked at the annual discharges, some years are improved but surprisingly some are worse.

[KevinShook]

I just looked through Classabedo. It's filled with floating point tests.
This is going to take a while.

[jhs507]

Yes, usually in cases like this it is standard practice to implement some sort of tolerance when making comparisons with floating point values.

We could do so relatively easily however the difference between single and double precision would still lead to slight differences of outcome even if we did back fill these comparisons into the borland code. Although I would guess they would be less.

However the question becomes what is the relevant tolerance for each comparison. In some cases it might make sense to have it parameterized when the tolerance represents a physical process that could be varied but if we are just trying to account for floating point issues we can just assume a small tolerance.

[jhs507]

I have created a separate branch to explore implementing tolerances for all of the floating point comparisons. I will work on that there as I am not sure it is the best way forward but would like to explore it.

[KevinShook]

Sounds like a good idea.
Looking at the Classalbedo code, many of the comparisons are more empirically-based than physically-based.

[jhs507]

I am noticing that this bit of code in Classalbedo.cpp line 117 is possibly causing problems:

if(net_snowD[hh] > 0.25) { // SF = SWE*2 if density 0.005 (0.5cm)
          Albedo[hh] = Albedo[hh] + net_snowD[hh]*0.1*2.0; // daily value
          if(Albedo[hh] > Albedo_snow[hh])
            Albedo[hh] = Albedo_snow[hh];
          continue;
        }

There are many times where the value of net_snowD[hh] appears to be very close to 0.25 ie (0.2499999999999999) could you explain a bit about what is being tested here?

[KevinShook]

This appears to be one of the empirical relationships/tests that I was referring to.
As far as I can tell (the comment makes zero sense), it's checking the daily snowfall,
and adding a very small value to the albedo.

[jhs507]

In the case where the daily snow is 0.2499999999999999 would it make sense to enter this branch or skip it?

[KevinShook]

I think we should be entering it

[jhs507]

Alright, That is currently what my "fix" code is doing. However we will need to move this new style of comparison back to the borland code if we want to continue to compare the two versions.

[KevinShook]

Agreed. However, we can test the results against the old version to see if it improves things.

[jhs507]

I did do such a test it seems that by adding a tolerance to all the comparisons in the albedo module that we moved towards closer to the borland in terms of results but I have indications that there are steps where we are moving in the opposite direction. The net effect seems to be getting closer to the borland results but there are still some places I am worried that we are "over correcting"

albedo with comparison tolerance.zip

[KevinShook]

Oof. Ok, thanks.
This module looks to be a real pain.

[jhs507]

I don't think there is a way to get a perfect match to the borland code without accounting for these floating point differences in the borland code as well. I am going to try a few different tolerances though to see if I can get closer without that.

[KevinShook]

Thanks. It looks like there are a lot of these tests in this module

[KevinShook]

Since you are creating a new branch, it would be great to use epsilon for all of the floating-point tests in Classalbedo.
This code is critical because of the positive feedback it sets up. If the albedo is too low, then snowmelt is increased, which then further decreases the albedo.
Could we try doing this?
Thanks

[jhs507]

That is what I did above, the data I posted is what happens when we place a comparison with tolerance in every floating point comparison in the albedo module.

[KevinShook]

Excellent. Can I try your code?

[jhs507]

Yes just give me a moment to push to that branch.

[KevinShook]

Thanks

[jhs507]

My code is committed on the comparison-tolerance branch. There is a logging line on common.cpp 343, and 375 that you may want to comment out before testing. Or send your output to file.

[KevinShook]

Thanks very much

[jhs507]

I just had a meeting with Kevin Schneider and we discussed our approach to these issues dealing with floating point precision. We think it would be productive to have a conversation with you about them.

What would be a good time for you? Kevin Schneider suggested that anytime after 2:00 today or after 1:00 tomorrow would work for him.

[KevinShook]

That would be great. I have a conference call from 1-2 today, so after 2 would be good for me.

[jhs507]

How about 2:15 then?

[KevinShook]

Works for me, thanks

[KevinShook]

Looking at Classebsm.cpp, line 205 is weird in that it is terminated with line break
double eamean = Common::estar(tmean[hh]) * rhmean[hh] / 100.0; \
This was also in the original Borland code. Any idea why it's there?

[jhs507]

That is weird but the compiler should just ignore it as the statement is properly ended with a ';'

[jhs507]

I did a bit more looking at the switches in Classalbedo. As we discussed yesterday at a certain point what path execution takes is a bit arbitrary but I focused in on the comparison on line 131.

Because it is a zero threshold I thought it would be particularity sensitive and correction would be easier. Also it is looking at SWE which we know often contains very small values.

After a bit of digging I found that in a 2 year period that by implementing a threshold of 1e-5 that there are 8967 instances where the SWE is above zero but under 1e-5. I wanted to see if they were evenly distributed on the HRU's and came up with this chart.

It seems HRU 5, and 10 are outliers. I don't know if that is relevant.

As you suggested above I think this really should use the same or similar threshold parameter as the pbsm module. It seems to me a place where we should change borland crhm as well.

[KevinShook]

I've been talking to John, and I think we should leave all of the thesholding alone for now.
There are too many of these things. Also, there are some issues with how we have been forcing the model
with daily precip, but hourly air temps and humiidities.
I'll send an update - an in the middle of soemthing right now.

[KevinShook]

I'm going to try re-running the model, using the old version of the code - the master branch.
This time I'm not going to use the Harder method for partitioning the precipitation into rain and snow.
I'm going to use air temperatures, which should be less susceptible to problems.
As I said above, I think that the problem is that the daily precipitation data is causing precipitation during
intervals when the RH is << 100%, which is not possible.
This is very likely the cause of snowfalls when they cannot occur.

I have run the Borland code with this configuration, so I can compare the crhmcode output against that run.

[jhs507]

Sounds good, let me know how it goes.

[KevinShook]

I changed the precipitation partitioning to use the air temperature rather than the Harder method.
As you can see, both models now the same evaporation values:

There are some issues with other values, but I think that we can consider this issue to be closed.