Inclusion of sea water content in sea ice basal mush for congelation growth

[proteanplanet]

Following from https://github.com/E3SM-Seaice-Discussion/E3SM/pull/3, the new branch https://github.com/E3SM-Project/E3SM/tree/proteanplanet/seaice/congelation has been created to include the contribution of liquid sea water in the energy of congelation ice growth for mushy thermodynamics.

Please consider the changes from

if (ktherm == 2) then

 qbotm = enthalpy_mush(Tbot, sss)
 qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
 qbot0 = qbotm - qbotp

 dhi = ebot_gro / qbotp     ! dhi > 0

 hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
 hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
 emlt_ocn = emlt_ocn - qbot0 * dhi

else

change to:

if (ktherm == 2) then

 qbotm = enthalpy_mush_liquid_fraction(Tbot, phi_i_mushy)
 qbot0 = enthalpy_brine(sst)
 qbotp = qbotm - qbot0

 dhi = ebot_gro / qbotp     ! dhi > 0

 hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
 hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
 emlt_ocn = emlt_ocn - qbot0 * dhi

else

The issue of whether or not congelation ice formation is correctly treated has previously been raised by @njeffery and @akturner. Now, others in the sea ice community have found that mushy as it's currently formulated can't replicate landface ice growth. One solution offered has been to decrease the liquid component of freshly-formed congelation ice. The approach intended in the the code provided here, proposed by @akturner is to include the full energy of sea water in the mush.

At present the branch runs with DEBUG switched on using Chrysalis and Anvil, but not without for a B-case. It runs successfully for D- and G-cases.

Please comment on this and alternate approaches. This has come to the fore because congelation ice growth is one of the main contributors to the Labrador Sea bias in winter, as shown here:

Congelation Ice Growth:

Snow Ice Formation:

Frazil Growth:

[njeffery]

@proteanplanet : If we fix a porosity for the new ice, phi_i_mushy, and a temperature, Tbot, then we are over prescribing by choosing bulk salinity to be sss. One solution would be to set the new salinity = phi_i_mushy * liquidus_brine_salinity_mush(Tbot). I think that would keep everything consistent.

The broader question is whether the changes in porosity from phi = 1 to phi_i_mushy come from 1) a change in temperature, 2) a change in salinity, or 3) a change in both. The first implementation looks like option 1). You're implementation looks like option 3). Is option 2) a possibility? Does it even matter?

[njeffery]

Option 2 would have qbotm = enthalpy_mush_liquid_fraction(sst,phi_i_mushy) and the bulk salinity value for the hstot would be phi_i_mushy * liquidus_brine_salinity_mush(sst). Also qbot0 = enthalpy_brine(sst)

[proteanplanet]

Here are results from a 10-year test using the above changes, which are superimposed on the hfdefault branch that was recently accepted into master:

[proteanplanet]

Also see https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.afroberts/E3SMv3/B12.junebaseline.congelation.E3SM-Project.anvil/mpas_analysis/ts_0001-0010_climo_0001-0010

[eclare108213]

I agree with @njeffery that hstot should not increase by dhi*sss if the mush enthalpy includes an ice fraction via the porosity phi_i_mushy. Her solution

set the new salinity = phi_i_mushy * liquidus_brine_salinity_mush(Tbot)

looks right to me. Also, emlt_ocn * dt is added to fhocn, which is the available ocean heat left over after use by ice and is sent back to the ocean. I do not think the increment to emlt_ocn should be the entire enthalpy of a dhi thickness of seawater.

This is my suggestion:

! qbotm = enthalpy of the new congelation/mush at temperature Tbot and porosity phi_i_mushy
! qbot0 = enthalpy of the brine fraction within the new mush at temperature Tbot
! qbotp = enthalpy of the solid ice fraction within the new mush
! qbotm = qbotp + qbot0

 qbotm = enthalpy_mush_liquid_fraction(Tbot, phi_i_mushy)
 qbot0 = phi_i_mushy * enthalpy_brine(Tbot)
 qbotp = qbotm - qbot0

 dhi = ebot_gro / qbotm     ! dhi > 0

 hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm            ! d(hqtot) = (fcondbot - fbot)*dt
 hstot = dzi(nilyr)*zSin(nilyr) + dhi*phi_i_mushy*liquidus_brine_salinity_mush(Tbot)
 emlt_ocn = emlt_ocn - qbot0 * dhi                     ! add heat associated with only the brine fraction

Note that this computes dhi using qbotm instead of qbotp, which means that qbotp is not used. Using this definition of qbotp to compute dhi would imply that the new ice contains no brine. In the original, qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy), which does include brine, but at which salinity and temperature?

I don't completely understand why emlt_ocn should be changed here -- I left it in because comments in the code indicate that the brine fraction is otherwise assumed to have zero enthalpy. If the ocean temperature is at the freezing point, then removing extra heat could cause additional frazil ice formation, but presumably the exact amount of water carrying this extra heat is also removed, so that the ocean temperature wouldn't change. Wouldn't this be captured in the ice enthalpy change? Conservation would need to be checked carefully for models using virtual water/salt fluxes as well as those that don't.

Not sure whether this all makes sense, just throwing it out there for further discussion.