Possible bug in SoilHydrologyMod due to bad variable names

[ekluzek]

Leo van Kampenhout < l.vankampenhout > - 2015-07-23 09:58:36 -0600
Bugzilla Id: 2193
Bugzilla CC: dlawren, oleson, rfisher, swensosc,
Bugzilla Attachment: snippet.F90 - variable is mistakingly used, latent heat, soil hydrology

Created attachment 57
variable is mistakingly used, latent heat, soil hydrology
(Bill S: Edit: See #118 (comment) for the contents of the old bugzilla attachment.)

revision clm4_5_1_r110

First of all, I do not know too much about CLM hydrology to be totally confident to say this is a bug, so please consider this report with a grain of salt. I think the impact if it WERE a bug would only be small but relevant.

The issue involves latent heat fluxes for snow covered areas. These are computed in SoilFluxesMod with the following variable names:

qflx_evap_grnd (liquid evaporation)
qflx_sub_snow (solid evaporation = sublimation)
qflx_dew_snow (solid deposition = frost)
qflx_dew_grnd (liquid deposition = dew)

See attached snippet of source code (Bill S: Edit: See #118 (comment) for the contents of the old bugzilla attachment) to see how this is done exactly (this involves a dependency on flux sign and ground temperature). Observe that the total of the four variables above always equals qflx_ev_snow, the so-called "evaporation flux from snow".

My first concern is the naming and commenting of these variables. The names are ambiguous and could easily be mistaken for something broader than what they represent (that is, latent heat to snow). I think this is in fact what happened to the suspected bug.

In SoilHydrologyMod.F90 the variable qflx_evap_grnd is used in a suspected wrong way. The variable is used in the urban column loop to determine surface water runoff, as if it represented evaporation in absence of snow.

 207             ! If there are snow layers then all qflx_top_soil goes to surface runoff
 208             if (snl(c) < 0) then
 209                qflx_surf(c) = max(0._r8,qflx_top_soil(c))
 210             else
 211                xs(c) = max(0._r8, &
 212                     h2osoi_liq(c,1)/dtime + qflx_top_soil(c) - qflx_evap_grnd(c) - &
 213                     pondmx_urban/dtime)

Possibly the coder meant to use qflx_evap_soi or qflx_evap_tot here, but again I'm by no means an expert in this code (especially urban).

If this turns out to be a bug, I request renaming of the variables in question to something more descriptive (e.g. qflx_evap_grnd and qflx_dew_grnd should have 'snow' in their names).

[ekluzek]

Erik Kluzek < erik > - 2017-11-03 14:57:33 -0600

I took an initial look into this, and thought it might be an issue. So I talked to Keith about it, and he agrees it needs more looking into. We think possibly the definition of these variables changed from applying everywhere to only applying over snow. So Keith will look into it more, and I'll look into the history a bit.

[ekluzek]

Erik Kluzek < erik > - 2017-11-03 16:10:14 -0600

OK, I traced it back. The definition of qflx_evap_snow changed when CLM4.5 came in. Specifically it came in on the h2osfc branch from Sean Swenson with this change to Biogeophysics2Mod.F90:

svn diff -r25822:r26395 https://svn-ccsm-models.cgd.ucar.edu/clm2/branches/h2osfc/models/lnd/clm/src/biogeophys/Biogeophysics2Mod.F90

@@ -409,22 +461,29 @@
        qflx_dew_snow(p) = 0._r8
        qflx_dew_grnd(p) = 0._r8
 
-       if (qflx_evap_soi(p) >= 0._r8) then
+!scs: instead of partitioning evap_soi, partition ev_snow for use in SnowHydrology
+!       if (qflx_evap_soi(p) >= 0._r8) then
+       if (qflx_ev_snow(p) >= 0._r8) then
           ! for evaporation partitioning between liquid evap and ice sublimation, 
          ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
          if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0.) then
-             qflx_evap_grnd(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
+!             qflx_evap_grnd(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
+             qflx_evap_grnd(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
          else
             qflx_evap_grnd(p) = 0.
          end if

[ekluzek]

Erik Kluzek < erik > - 2017-11-03 16:12:26 -0600

The change was presumably fine for other column types, but urban got merged in at this same time, and it was assuming the previous definition. In terms of trunk tags this comes in at clm4_0_60.

[ekluzek]

Keith Oleson < oleson > - 2017-11-03 16:34:33 -0600

It looks to me as if the original definition of the variable qflx_evap_grnd changed at some point since it's original use by the urban code. The urban code itself has not changed.

The original code was:

   if (qflx_evap_soi(p) >= 0._r8) then
      ! for evaporation partitioning between liquid evap and ice sublimation,
      ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
      if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0.) then
         qflx_evap_grnd(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
      else
         qflx_evap_grnd(p) = 0.
      end if
      qflx_sub_snow(p) = qflx_evap_soi(p) - qflx_evap_grnd(p)

So, it was defined as the liquid part of soil OR snow evaporation (qflx_evap_soi).

At some point that variable seems to have been co-opted as representing the liquid water part of evaporation from snow only:

     if (qflx_ev_snow(p) >= 0._r8) then
        ! for evaporation partitioning between liquid evap and ice sublimation,
        ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
        if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0.) then
           qflx_evap_grnd(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
        else
           qflx_evap_grnd(p) = 0.
        end if
        qflx_sub_snow(p) = qflx_ev_snow(p) - qflx_evap_grnd(p)

However, in the new code there is also this statement:

     ! set ev_snow, ev_soil for urban landunits here
     l = patch%landunit(p)
     if (lun%urbpoi(l)) then
        qflx_ev_snow(p) = qflx_evap_soi(p)
        qflx_ev_soil(p) = 0._r8
        qflx_ev_h2osfc(p) = 0._r8

Here, qflx_ev_sno is being assigned to qflx_evap_soi for urban. So, I think the new code is actually using qflx_evap_soi, which is correct.

I have no idea how this code got in there, but it appears to be ok I think in terms of how xs is derived.

However, it seems like the qflx_ev_snow variable on the history file (QSNOEVAP) will have evaporation from urban classified as snow evaporation which will be wrong when there is no snow. And possibly other variables derived from qflx_ev_snow.

I also completely agree that these variables names are confusing and that they should have snow included in their variable names and in their long_names in the history file output.

[billsacks]

Here are the contents of the attachment from bugzilla: Snippet from SoilFluxesMod.F90

Lines 128-132 (in that old version of the code):

qflx_evap_grnd          => waterflux_inst%qflx_evap_grnd_patch     , & ! Output: [real(r8) (:)   ]  ground surface evaporation rate (mm H2O/s) [+]
qflx_sub_snow           => waterflux_inst%qflx_sub_snow_patch      , & ! Output: [real(r8) (:)   ]  sublimation rate from snow pack (mm H2O /s) [+]
qflx_dew_snow           => waterflux_inst%qflx_dew_snow_patch      , & ! Output: [real(r8) (:)   ]  surface dew added to snow pack (mm H2O /s) [+]
qflx_dew_grnd           => waterflux_inst%qflx_dew_grnd_patch      , & ! Output: [real(r8) (:)   ]  ground surface dew formation (mm H2O /s) [+]
qflx_ev_snow            => waterflux_inst%qflx_ev_snow_patch       , & ! In/Out: [real(r8) (:)   ]  evaporation flux from snow (W/m**2) [+ to atm]

Lines 317-340 (in that old version of the code):

! Assign ground evaporation to sublimation from soil ice or to dew
! on snow or ground

qflx_evap_grnd(p) = 0._r8
qflx_sub_snow(p) = 0._r8
qflx_dew_snow(p) = 0._r8
qflx_dew_grnd(p) = 0._r8

if (qflx_ev_snow(p) >= 0._r8) then
   ! for evaporation partitioning between liquid evap and ice sublimation,
   ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
   if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0.) then
      qflx_evap_grnd(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
   else
      qflx_evap_grnd(p) = 0.
   end if
   qflx_sub_snow(p) = qflx_ev_snow(p) - qflx_evap_grnd(p)
else
   if (t_grnd(c) < tfrz) then
      qflx_dew_snow(p) = abs(qflx_ev_snow(p))
   else
      qflx_dew_grnd(p) = abs(qflx_ev_snow(p))
   end if

[billsacks]

From reading back through the above comments, it sounds like @olyson found that the initial suspected bug is not actually a bug, due to a compensating (confusing) correction in the code. But there are possibly two outstanding issues here:

(1) Over urban points, QSNOEVAP and possibly other diagnostic fields are wrong: see #118 (comment) for details. It sounds to me like the code should be revised to remove the confusing, compensating correction that leads to incorrect diagnostic fields (code snippet called out in #118 (comment) ) , and then the original issue identified by @lvankampenhout should similarly be revised accordingly. But I'm not clear enough on the details at this point to understand the exact changes needed.

(2) The four variables qflx_evap_grnd, qflx_sub_snow, qflx_dew_snow and qflx_dew_grnd may have misleading names (or at least, some of them might). I haven't come to understand this piece well enough, but if true, then a separate code cleanup issue should be opened to rename those variables.

[billsacks]

@olyson Based on input from @dlawrenncar I'm assigning this to you, but it's low priority.

[olyson]

FYI, I started working on this Friday.

[olyson]

I've added code to have urban columns explicitly use qflx_evap_soi in these equations instead of assigning qflx_ev_snow to qflx_evap_soi. This makes the code less confusing.
However, I'm finding that even in cases when there is no snow, qflx_ev_snow is non-zero, in fact it is equal to qflx_ev_soil. The code appears to be designed this way. And the sum of the variables qflx_evap_grnd, qflx_sub_snow, qflx_dew_snow, and qflx_dew_grnd will always be equal to qflx_ev_snow (and qflx_ev_soil) even in the absence of snow. So for example, qflx_sub_snow will either represent sublimation from snow ice in the presence of snow, or sublimation of soil ice in the absence of snow. On the other hand, maybe everybody knew this except me.
All of these variables have misleading names in some context then.
I'm not sure what to suggest for new variables names.
Maybe something like:

qflx_liqevap_from_soil_or_snow
qflx_solidevap_from_soil_or_snow
qflx_liqdew_on_soil_or_snow
qflx_soliddew_on_soil_or_snow

I know, I don't like them either.
Happy New Year!

[billsacks]

Thanks for looking into this @olyson . I agree that this is fairly confusing no matter what. I do find your variables better than what is currently there. Would an alternative be to replace from_soil_or_snow with from_top_layer, or maybe just top for short? Is that correct that these terms are all additions to / removals from the top layer, which might be soil or might be snow? (I'm not sure which I like better between my suggestions and yours, but just wanted to put this alternative out there for consideration).

[olyson]

Yes, they are all additions/removals to/from the top layer.
I think I like top_layer more than just top.
So I propose the following variable renaming and descriptions:

qflx_evap_grnd -> qflx_liqevap_from_top_layer (rate of liquid water evaporated from top soil or snow layer (mm H2O /s) [+])
qflx_sub_snow -> qflx_solidevap_from_top_layer (rate of ice evaporated from top soil or snow layer (sublimation) (mm H2O /s) [+])
qflx_dew_snow -> qflx_liqdew_on_top_layer (rate of liquid water deposited on top soil or snow layer (dew) (mm H2O /s) [+])
qflx_dew_grnd -> qflx_soliddew_on_top_layer (rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+])

Another minor suggestion would be to use "from/to" instead of "from/on".
I will give this a couple of days for any more comments and then proceed on with changes to my branch.

[billsacks]

I like your suggestions. Good point, too, about from/to being more symmetrical than from/on.

[olyson]

I've renamed the variables as discussed throughout the code.
I've also made a couple of minor changes as described in the commit logs.
These changes are all bfb (with the exception of QSNOEVAP) for a 1 month global 2deg simulation, but may not be bfb under all conditions. I tested for bfb before changing the history field variables themselves. QSNOEVAP is different because I've added in qflx_ev_snow for lakes.
So I think I can issue a pull request unless there are any objections. Sean S. has taken a look at this also.
My branch is https://github.com/olyson/ctsm/tree/issue118.

[billsacks]

Thanks @olyson ! Yes, please go ahead and open a PR.