Add conductive flux limiting for vertical thermodynamics when calc_Tsfc=.false.

[kieranricardo]

PR checklist

• Short (1 sentence) summary of your PR:
  This limits the conductive flux for the vertical thermodynamic solve when calc_Tsfc=.false. See: Add conductive flux limiting for vertical thermodynamics when calc_Tsfc=.false. #556
• Developer(s):
  @kieranricardo @blimlim @anton-seaice
• Suggest PR reviewers from list in the column to the right.
• Please copy the PR test results link or provide a summary of testing completed below.
• How much do the PR code changes differ from the unmodified code?
  • bit for bit
  • different at roundoff level
  • more substantial
• Does this PR create or have dependencies on CICE or any other models?
  • Yes
  • No
• Does this PR add any new test cases?
  • Yes
  • No
• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please document the changes in detail, including why the changes are made. This will become part of the PR commit log.

When the vertical thermodynamic solver is forced by prescribed surface fluxes (calc_Tsfc = .false.), it can become numerically unstable. This instability appears to arise in the vertical thermodynamic solve when the conductive flux is too large (which depends on ice thickness and time step). Flux limiting is needed to stabilise the solution in this configuration.

An existing implementation of this flux limiting is available in an ACCESS fork of CICE5. This implements a similar approach, enabling stable coupling with the UM.

This conductive flux limiting occurs in 2 places:

1. Before the thermodynamic solve. If the conductive flux to ice thickness ratio exceeds a certain value (default 1000 w/m^2) the conductive flux is reduced such the ratio matches this value. The excess energy is used for bottom melt.

2. During the thermodynamic solve. If any of the layer temperatures exceed 0C, the layer temperature is reduced to 0C. The excess energy is deducted from the conductive flux and passed to the ocean.

[anton-seaice]

This limiting is noted in the icepack docs here:

Icepack/doc/source/science_guide/sg_boundary_forcing.rst

Lines 83 to 85 in 222ab41

	returned to the atmosphere model via the coupler. Since the ice surface
	temperature is treated explicitly, the effective conductivity may need
	to be limited to ensure stability. As a result, accuracy may be

I can run icepack standalone tests when this is ready

[apcraig]

Is this change bit-for-bit with calc_Tsfc=.false.? Is that because flux limiting is rarely triggered? When it is triggered, it changes answers, correct? We probably should run this update in CICE as well as Icepack.

[anton-seaice]

should these be within a if ( .not. calc_Tsfc) ?

[kieranricardo]

It looks like there's no other if checks in icepack_write_parameters, so I think all parameters are output regardless of whether they're used. It might be inconsistent to add an if statement here.

[eclare108213]

this is okay as long as they are not used elsewhere except when calc_Tsfc = .false.

[anton-seaice]

Is this change bit-for-bit with calc_Tsfc=.false.? Is that because flux limiting is rarely triggered? When it is triggered, it changes answers, correct? We probably should run this update in CICE as well as Icepack.

Looks like alt03 sets calc_Tsfc=.false., so quite possibly could be answer changing in that test.

https://github.com/CICE-Consortium/CICE/blob/783a867d24520403cb35cec5f07cda8b275448b1/configuration/scripts/options/set_nml.alt03#L12

[kieranricardo]

@apcraig flux limiting only happens if calc_Tsfc=.false., so this shouldn't be answer changing for calc_Tsfc=.true.

[anton-seaice]

@kieranricardo has added these as configurable - but we don't plan to add them to the cice namelist (as we dont foresee changing them).

Are we happy to change answers when calc_Tsfc = .false. - @apcraig ? If so, should we not make these two new parameters configurable through icepack_init_parameters ?

[eclare108213]

Are there values of these parameters that would preserve BFBness when calc_Tsfc = .false.? If so, I recommend setting the code defaults to those values and allow them to be changed via namelist.

[kieranricardo]

We could pick extreme value to ensure the cap_conductive_fluxroutine doesn't do anything, but because the flux is also limited in bl99 with e_num this will still change answers when calc_Tsfc = .false..

[eclare108213]

Overall, this change seems fine to me. I would prefer that it be provided as an option whose default behavior does not change answers. In terms of the coding, it's usually better to use pre-existing variables rather than introducing new ones to do more-or-less the same thing. I'm not sure how easy it is in this case, but if possible, I recommend trying to use einex instead of e_num and doing it in a way that minimizes code invasiveness, which might be a bit of a balancing act.

[eclare108213]

this is okay as long as they are not used elsewhere except when calc_Tsfc = .false.

[eclare108213]

How are einex and e_num related? Could you modify einex instead, just above, so that this block of code does not change?

[kieranricardo]

einex and e_num are both excess energy that's been removed from the ice to bring the temperature below melting. They're calculated slightly differently, e_num is accumulated over solver iterations, whereas einex isn't.

The main difference is what's done with this excess energy. einex is added to fcondbot so it affects bottom melt/growth. And e_num is added to fhocn, so it's just passed straight to ocean.

I'm actually not sure what the best approach is for calc_Tsfc=.false., the einex approach seems more physical, and would melt additional ice over summer which would be helpful. Thoughts @anton-seaice @ofa001? ACCESS-CM2 did the e_num approach for context.

[anton-seaice]

Are you suggesting not using e_num at all and seeing if it's stable using the already available einex ? Give it a try I guess, it looks like einex should be calculated fine when calc_Tsfc is false?

[kieranricardo]

I meant use the current approach for calculating e_num (I think we'll need this for stability), but use this for bottom melt/growth instead of passing this to the ocean.

In the code this would add extra checks around einex, calculating it normally if calc_Tfsc = .true., and calculating einex using the e_num approach if calc_Tfsc = .false..

[kieranricardo]

After some discussion it doesn't seem like this approach will work. With calc_Tfsc=.false. the extra numerical flux can be much larger than the calc_Tfsc=.true., so adding this to fcondbot could cause unrealistic bottom melt/growth. @eclare108213 so it looks like we won't be able to use einex approach here.

[eclare108213]

Thanks for looking into this. Then perhaps a compromise will work. In the PR comments at the top:

This conductive flux limiting occurs in 2 places:

1. Before the thermodynamic solve. If the conductive flux to ice thickness ratio exceeds a certain value (default 1000 w/m^2) the conductive flux is reduced such the ratio matches this value. The excess energy is used for bottom melt.
2. During the thermodynamic solve. If any of the layer temperatures exceed 0C, the layer temperature is reduced to 0C. The excess energy is deducted from the conductive flux and passed to the ocean.

Since the excess in 1. is used for bottom melt, that's essentially supplementing einex, isn't it? And then 2. handles it differently. So perhaps einex could be adjusted for 1. instead of using a new variable, and the new variable still used for 2. If that works, then I would rename the two variables so they are more intuitive, e.g. 1. einex_melt and 2. einex_ocnflux. Otherwise, we can keep the current functionality but rename the variables accordingly, e.g. 1. einex_sfc_calc and 2. einex_sfc_flux or something like that (implicit / explicit? Will need to look at the GEOS names to try to be consistent). einex translates to 'excess ein' for both cases.

[kieranricardo]

Good point! But that first approach could be a little involved so I think the second approach sounds good. The first approach would involve adding an extra inout argument to temperature_changes, being careful not to update einex_melt when calc_tsfc = .false., and would complicate the energy accounting.

[kieranricardo]

I've updated the variable name to einex_sfc_calc and einex_sfc_flux, but not sure which GEOS variable names this should match.

[eclare108213]

Are there values of these parameters that would preserve BFBness when calc_Tsfc = .false.? If so, I recommend setting the code defaults to those values and allow them to be changed via namelist.

[kieranricardo]

Thanks for the review @eclare108213!

Yep ratio_Wm2_m and cold_temp_flag cold temp are only used when calc_tsfc=.false.. Happy to enable these to be changed via a namelist, but what do you mean by BFBness?

[eclare108213]

Thanks for the review @eclare108213!

Yep ratio_Wm2_m and cold_temp_flag cold temp are only used when calc_tsfc=.false.. Happy to enable these to be changed via a namelist, but what do you mean by BFBness?

BFB = bit-for-bit. We try to maintain exact reproducibility to the extent possible, except when an answer-changing bug is being fixed. This often requires adding a flag to turn the new feature on and off or to select it.

It would be great if the einex functionality could work for your case!

[kieranricardo]

@eclare108213 ah gotcha, this will be BFB when calc_Tfsc=.true., but will change answers when it's false.

[eclare108213]

Initial code review. In addition to the formatting changes noted below, this functionality needs to be documented.

The GEOS approach is described in sections 2.1 and 2.7 of the Icepack documentation. It seems to be quite a bit different from what you are doing (especially 2.7). The last paragraph of section 2.1 generally covers this situation, and I think your approach is now the third one. This needs to be checked: only the GEOS approach applies to both BL99 and mushy thermo, and the other two only affect BL99. This needs to be clarified in the docs, along with how the functionality is accessed (e.g. it's on by default in your case).

[eclare108213]

@kieranricardo @anton-seaice: just checking on the status of this PR. Let me know if you'd like to discuss it.

[kieranricardo]

@eclare108213 thanks for your review, I've now added description and units for the all the new variables (and some related existing variables as well).

Shall I edit the documentation in this PR? I can add a description of the limiting in 2.7.4.1.

[anton-seaice]

It is covered somewhat in the second paragraph in 2.1, so probably update that ?