Remove infiltration excess runoff for VIC

[billsacks]

NOTE: This PR was originally created in the temporary CTSM repository, Sep 13, 2017, where it was NCAR/clm-ctsm#17. I am moving over most of the comments from there: both my own and those of @martynpclark .

Here was the original comment:

Martyn Clark reviewed the VIC implementation, and felt that the current
implementation of infiltration excess runoff is inconsistent with the
standard VIC implementation. It appears that what was being called VIC's
infiltration excess runoff was actually just an attempt to give a better
numerical approximation to the solution for saturated surface excess
runoff. So deleting this leaves only a first-order approximation to
VIC's saturated surface excess runoff.

Eventually we may want to put in place a more accurate solution for
VIC's saturated surface excess runoff. But Martyn's feeling is that this
can come in with other changes we want to make regarding numerical
solutions in CTSM.

[billsacks]

Comment from @billsacks Sep 13, 2017

When we talked about this a few weeks ago, there were suggestions of using a global constant here.

Reasons I went with a locally-defined qinmax_huge rather than a globally-defined constant:

• Globally-defined is useful to eliminate the need for, "If I change this here, then I also need to change this over there". That doesn't seem to be the case here:
  • Nobody needs to check this
  • I couldn't see any benefits to having this shared between different parameterizations
• Problems with globally-defined
  • A different value might make more sense in other contexts
  • If another parameterization needs a different value and so changes the global constant, it would change this and other previous parameterizations that use it, which could inadvertently change behavior

[billsacks]

Comment from @martynpclark Sep 13, 2017

Agreed. This seems like a good approach.

[billsacks]

Comment from @martynpclark Sep 13, 2017

@billsacks -- good description -- the VIC implementation was actually the analytical solution to saturation-excess runoff (which depends on saturated area), not infiltration-excess runoff (which depends on the flux at the upper boundary).

[billsacks]

Comment from @martynpclark Sep 13, 2017

@billsacks -- we have already got to the point where we are making changes to the numerical solutions in CTSM. The VIC implementation provides an analytical solution, in the same way that we use analytical solutions for the ponded water store.

[billsacks]

Comment from @martynpclark Sep 13, 2017

The code for the VIC analytical solution is by no means obvious. The derivation is in this pdf file:

VIC_surfaceRunoff.pdf

[billsacks]

Comment from @martynpclark Sep 13, 2017

It is correct to delete this code here, and now that we have options for analytical numerical solutions, the code block below can be moved to the saturation excess runoff routine:

       max_infil = (1._r8+b_infil(c)) * top_max_moist(c)
       i_0 = max_infil * (1._r8 - (1._r8 - fsat(c))**(1._r8/b_infil(c)))
       if(qflx_in_soil(c) <= 0._r8) then
          rsurf_vic = 0._r8
       else if(max_infil <= 0._r8) then
          rsurf_vic = qflx_in_soil(c)
       else if((i_0 + qflx_in_soil(c)*dtime) > max_infil) then             !(Eq.(3a) Wood et al. 1992)
          rsurf_vic = (qflx_in_soil(c)*dtime - top_max_moist(c) + top_moist(c))/dtime
       else                                                                !(Eq.(3b) Wood et al. 1992)
          basis = 1._r8 - (i_0 + qflx_in_soil(c)*dtime)/max_infil
          rsurf_vic = (qflx_in_soil(c)*dtime - top_max_moist(c) + top_moist(c)    &
               + top_max_moist(c) * basis**(1._r8 + b_infil(c)))/dtime
       end if

We can use the same logic as in the ponded water (options for explicit Euler [current code] and the analytical solution [this code]).

Again, the pdf document above provides derivations.

[billsacks]

Comment from @martynpclark Sep 13, 2017

We can add the analytical solution here.

[billsacks]

Comment from @billsacks Sep 14, 2017

@martynpclark I feel like there are still two somewhat unresolved issues, both really boiling down to how much complexity we want to support:

1. I recall that moving this code into the saturated surface excess runoff calculation would mean that different parameterizations for saturated surface excess runoff now potentially differ in more than just fsat. That's fine in principle, but my sense from the Aug 24 ctsm meeting was that nobody was very happy with what that would mean in terms of code complexity, and that we only wanted to go there if we really had to. It's unclear to me whether the vic option is important enough to justify this extra complexity itself, or if this is only worth doing if some other parameterizations also require this.

2. I feel like we haven't totally worked out how we want to select between different solution methods. In this case, I'd personally be more inclined to just always include the more exact analytical solution, rather than providing an option to include it or exclude it, in order to avoid extra levels of complexity. Of course, I'm open to having this be an option, but I feel we should be careful to only include options like this where they're really warranted.

We can continue this discussion here if you want. Or we can defer this discussion to a later date if (as is my impression) this vic stuff isn't the highest priority.

[billsacks]

Comment from @martynpclark Sep 14, 2017

Hi @billsacks -- I agree. Let's defer the decision until once we converge on the overall structure that we are happy with. Yes, I agree for the operator splitting implementation that we should (ideally) always use the more accurate method if it is available -- one question is if there is a desire to have what was present before, since different solutions can lead to different model behavior.

[billsacks]

Comment from @billsacks Sep 14, 2017

@martynpclark

I agree. Let's defer the decision until once we converge on the overall structure that we are happy with.

So any feelings about what I should do in the meantime?:

1. Merge this PR as is (once I have tested it)

2. Hold off on merging this PR until we decide on the final solution

one question is if there is a desire to have what was present before, since different solutions can lead to different model behavior.

But here, isn't the more accurate method what was there before? (Just now we'd be applying it to saturated surface excess runoff rather than infiltration excess runoff?)

[billsacks]

Comment from @martynpclark Sep 14, 2017

Hi @billsacks

So any feelings about what I should do in the meantime?

My preference is to merge this PR as is

But here, isn't the more accurate method what was there before? (Just now we'd be applying it to saturated surface excess runoff rather than infiltration excess runoff?)

Exactly. My comments were intended to be more general. There may be cases where we add a more accurate solution (as in the case for ponded water), and we still may have the desire to retain what was present before (since different solutions can lead to different model behavior). In these cases we may not always want to use the more accurate solution, and we retain capabilities in order to compare against previous model versions.

[billsacks]

Comment from @billsacks Sep 14, 2017

@martynpclark thanks for the clarifications - this all makes sense now. I'll run testing on this as soon as I'm done with #13

[billsacks]

@martynpclark in case it wasn't clear: no action / replies are needed; I just wanted to move over our old conversation to save it in the new repository.