default 'hs0' value different in code/doc (0.03) vs namelist (0)

[phil-blain]

While auditing our namelist for our project of transitioning to CICE6, I noticed that hs0 is set to 0 in the namelist, but its default value (in both ice_init and icepack_parameters ) is 0.03 m. It's also documented in both the CICE and Icepack documentation as defaulting to 0.03.

I'm not familiar with the Delta-Eddington code but I just thought I'd point that out because it seemed a little weird... From a quick search in https://github.com/CICE-Consortium/CICE-svn-trunk it seems this change in the namelist was done sometimes between 4.0 and 5.1.2...

[eclare108213]

Good catch, @phil-blain. We need to clarify the documentation.
Please, would you check something in the code? I suspect that hs0 is not used for every melt pond scheme, and in particular it's possible that it's not used for our default melt pond scheme (hence set to 0), but when it is used, we recommend the 0.03 value (i.e. 0.03 is the default for the non-default scheme(s)). I'm guessing based on vague memory, but if that's the case, we can document it as such, otherwise then we should decide which value to call 'default'.

[phil-blain]

It is used only in icepack_shortwave, at 2 places:

1. In subroutine shortwave_dEdd_set_snow
   https://github.com/CICE-Consortium/Icepack/blob/53ffce04d729a3e802b397340c4410126b5b0ac9/columnphysics/icepack_shortwave.F90#L3655-L3661

      ! set snow horizontal fraction
      hs = vsno / aice
      
      if (hs >= hs_min) then
         fs = c1
         if (hs0 > puny) fs = min(hs/hs0, c1)
      endif

This subroutine is called at the beginning of rundEdd, for all melt pond schemes. It is also called a second time in rundEdd, but only for the level melt pond scheme (https://github.com/CICE-Consortium/Icepack/blob/53ffce04d729a3e802b397340c4410126b5b0ac9/columnphysics/icepack_shortwave.F90#L1029-L1048)

2. It is also used directly in rundEdd , but only for the CESM meltpond scheme: https://github.com/CICE-Consortium/Icepack/blob/53ffce04d729a3e802b397340c4410126b5b0ac9/columnphysics/icepack_shortwave.F90#L1013-L1024

            ! set pond properties
            if (tr_pond_cesm) then
               ! fraction of ice area
               fpn = apndn(n)
               ! pond depth over fraction fpn 
               hpn = hpndn(n)
               ! snow infiltration
               if (hsn >= hs_min .and. hs0 > puny) then
                  asnow = min(hsn/hs0, c1) ! delta-Eddington formulation
                  fpn = (c1 - asnow) * fpn
                  hpn = pndaspect * fpn
               endif

[eclare108213]

Okay, that's more complicated -- I was only remembering item 2. Table 1 in the paper cited below shows that hs0=0 except for the CESM test, but that hs1=0.03 for almost all of the tests. hs0 is the snow patchiness factor on sea ice; hs1 is for snow on top of refrozen ponds. Need to be careful that hs0 might be used as a local variable in a subroutine but hs1 is being passed in... sorry, that's confusing I know.

E. C. Hunke, D. A. Hebert, O. Lecomte (2013). Level-ice melt ponds in the Los Alamos Sea Ice Model, CICE. Ocean Modelling, 71, 26–42. doi: 10.1016/j.ocemod.2012.11.008.

[eclare108213]

See also https://bb.cgd.ucar.edu/cesm/threads/snow-depth-transition-to-bare-sea-ice-icepack.7700/

[eclare108213]

I'm starting to remember what the deal is with hs0. I think the code is okay but admit that it's confusing and could be documented better. Here's what I think is true -- points (5) and (6) below should probably be checked by stepping through the radiation calculation in some runs in the various configurations:

1. hs0 sets the snow fraction fs on bare ice via a ramp based on snow thickness, which allows some bare ice to be available for a radiation calculation, i.e. radiation is computed for both bare ice and snow (and/or ponds) on each ice thickness category. The original reason for hs0 was numerical: To prevent a sudden change in the shortwave reaching the sea ice (which can prevent the thermodynamics solver from converging), thin layers of snow on pond ice are assumed to be patchy, thus allowing the shortwave flux to increase gradually as the layer thins.
2. hs1 does a similar thing for snow on refrozen ponds. See E.C. Hunke et al., Ocean Modelling 71 (2013) 26–42 (the melt pond paper with David Hebert and Olivier Lecomte) for a discussion of the impacts of varying hs1.
3. At each time step, fs=0 initially, then if the snow is deep enough, a snow fraction 0 < fs < 1 is computed, as long as hs0 > 0. If hs0 = 0, then fs = 1.
4. fs is used for the super-simple pond formulation in dEdd (which is just in there for testing) to compute the pond fraction fp. shortwave_dEdd_set_pond could be inlined, which might help reduce confusion.
5. fs is used --sometimes -- when topo ponds are on. If topo ponds are present, then their fraction fp is computed by the topo scheme and fp is then used to set fs. But if there aren't any topo ponds present (fp=0) then the value of fs from (3) is used instead.
6. fs=1 for the level-ice pond scheme initially, because the ponds are allowed to infiltrate snow over the level ice area until they get deep enough to show through, at which point fp < 1 has been computed. Like the topo scheme, fp is used to set fs, and in this case fs+fp=1. The ponds are only on the level ice area, and so there's still snow on the ridges even if the entire level ice area becomes filled with ponds. Maybe, when there aren't any ponds, the snow fraction could also be ramped down, but that would need to be tested. If hs0 /= 0 in namelist, the code currently resets hs0=0 with level-ice ponds. This is done because the infiltration of snow by pond water accomplishes the gradual radiative forcing transition for which the patchy-snow parameters were originally intended.
7. In the tests shown in the Hunke et al paper above, hs0=0 for all cases except when the cesm pond scheme is on. That indicates that hs0 was only intended to be used with that scheme, not with topo or level-ice ponds.
8. Since level-ice ponds ON is the default configuration for CICE and Icepack, we set hs0 = 0. But when level-ice ponds aren't used, then hs0 = 0.03 traditionally. This is what needs to be clarified in the documentation.

[eclare108213]

Actually, since the cesm pond scheme is being deprecated, we could set hs0=0 as the default in the code and the docs, and also explain the situation above in the docs, for better clarity.

[eclare108213]

To make things even more confusing, we also have a parameter snowpatch = 0.02 that functions the same way as hs0, but for the CCSM3 shortwave scheme.

[eclare108213]

I changed the default values of hs0 and added an explanation to the Icepack documentation (see the thermodynamics section 2.7 in the Icepack PR).
#787
CICE-Consortium/Icepack#411

[eclare108213]

I changed the default values of hs0 and added an explanation to the Icepack documentation (see the thermodynamics section 2.7 in the Icepack PR). #787 CICE-Consortium/Icepack#411

I'm wondering if it would be better to not change the hs0 defaults in these PRs, and clean that up separately. It's a slightly bigger can of worms than I originally thought. In particular, I'm not sure whether hs0=0 will make the topo scheme more prone to thermo crashes, since it "fills" the snow with melt water only one category at a time, unlike the level-ice scheme, which allows meltwater to infiltrate snow on level ice in all categories. For that matter, maybe we should allow nonzero hs0 for the level-ice scheme, to be able to check whether sudden radiation jumps is the issue if/when the thermo crashes with that one? The level-ice scheme has a different tapering between the snow and pond fractions, but if hs0=0 and there aren’t any ponds, then there’s no tapering between snow and bare ice. Could that be a problem? I’m not sure if the problem that hs0 (and likewise hs1 and snowpatch) fixes is only in the melting season or also other times of year. Snow grain metamorphosis in the new snow physics should help with this problem during the melting season — it’s a more realistic tapering mechanism.

[eclare108213]

CESM ponds are being deprecated, and hs0 was originally implemented for that scheme for numerical stability reasons. The physical rationalization is that it allows a representation of snow “patchiness” effects on the surface radiation budget, which smooths the transition between fully snow-covered ice and bare ice based on snow depth.

Part of the confusion is that depending on how we change the code and namelist settings, the BFB-ness of our tests and users’ configurations may be different. There are a number of ways to approach the problem (more than are listed here).

Request for Consortium team: Would someone please double-check the code and make sure that my logic below is correct? All, please weigh in on the choices below.

Question for Consortium team and the community: In your simulation configurations, does the thermo intermittently and mysteriously fail to converge when tr_pond_lvl=true? If so, maybe we should encourage hs0 > 0 all the time (D below).

Background on the current state of the code and scripts:

• hs0 = 0.03 at initialization (both CICE and Icepack)
• hs0 = 0 in ice_in and icepack_in
• nonzero hs0 has an effect when tr_pond_topo = T
• nonzero hs0 is reset to hs0 = 0 when tr_pond_lvl = T and so does not have an effect
• hs0 is not included in any options files and so hs0 = 0 for all of our tests (primary namelist value)

In all cases:

• Note hs0 values/behavior explicitly in the release notes.
• Explain hs0 uses/behavior in the Icepack documentation.
• Review the troubleshooting guidance and update appropriately. Ensure that the primary reason for thermo convergence errors is clear — inconsistent forcing (units, coastlines, etc) — and mention hs0 as a potential engineering fix for radiation-related failures.
• Move the initial call to shortwave_dEdd_set_snow into the tr_pond_lvl and ‘else’ portions of the ‘set pond properties’ block and change the comment to ‘set snow and pond properties’ (so that tr_pond_lvl is not calling it twice)

A. BFB:
Status quo

• Leave hs0=0.03 as the initial value in the code and hs0=0 in ice_in and icepack_in, but document the situation clearly.

B. Potentially nonBFB fix 1:
hs0=0 everywhere
Allows snow patchiness for level-ice ponds

• Set hs0=0 as the initial value in the code
• Add hs0=0 for namelist option files with tr_pond_topo=T
• Add hs0=0 for namelist option files with tr_pond_lvl=T
• Remove the initialization code that resets hs0 to 0 for tr_pond_lvl

tr_pond_topo=T: BFB for in our tests and BFB for users unless hs0 is not set in their namelist
tr_pond_lvl=T: BFB in our tests and BFB for users unless hs0 is nonzero or not set in their namelist

C. Potentially nonBFB fix 2:
As in B but change tr_pond_lvl behavior in icepack_shortwave.F90
Maintains physical consistency with level-ice ponds design, conceptually, by not allowing snow patchiness

• Set hs0=0 as the initial value in the code
• Add hs0=0 for namelist option files with tr_pond_topo=T
• Add hs0=0 for namelist option files with tr_pond_lvl=T
• Remove the initialization code that resets hs0 to 0 for tr_pond_lvl
• In subroutine shortwave_dEdd_set_snow, comment out
  if (hs0 > puny) fs = min(hs/hs0, c1)
  and add
  if (.not.(tr_pond_lvl) .and. hs0 > puny) fs = min(hs/hs0, c1).
• Note in the troubleshooting documentation that if the thermo aborts with a nonconvergence error and all other potential forcing problems have been ruled out, then uncommenting the original line that sets fs might help.

tr_pond_topo=T: BFB in our tests and BFB for users unless hs0 is not set in their namelist
tr_pond_lvl=T: BFB in our tests and BFB for users unless hs0 is not set in their namelist

D. NonBFB fix:
hs0=0.03 everywhere
Allows snow patchiness for level-ice ponds

• Keep hs0=0.03 as the initial value in the code
• Set hs0=0.03 in ice_in and icepack_in
• Add hs0=0.03 for namelist option files with tr_pond_topo=T
• Add hs0=0.03 for namelist option files with tr_pond_lvl=T
• Remove the initialization code that resets hs0 to 0 for tr_pond_lvl

tr_pond_topo=T: nonBFB in our tests and nonBFB for users unless hs0=0.03 or is not set in their namelist
tr_pond_lvl=T: nonBFB in our tests and nonBFB for users unless hs0=0.03 or is not set in their namelist

[DeniseWorthen]

@eclare108213

Question for Consortium team and the community: In your simulation configurations, does the thermo intermittently and mysteriously fail to converge when tr_pond_lvl=true? If so, maybe we should encourage hs0 > 0 all the time (D below).

I would say generally that yes, we see do see issues with seemingly random thermo convergence issues in UFS. We run only the tr_pond_lvl case. We use hs0 = 0. and hs1 = 0.03.

[eclare108213]

@proteanplanet found in RASM that with tr_pond_lvl, hs0 must equal 0. This was because of issues with exact restarts associated with the radiation being called twice each time step.

[dabail10]

I am wondering what the plan was here? Should I set up testing for the suggestions above?