Skip to content

fix: correct physics bugs in 3-equation ice cavity melt model#875

Open
JanStreffing wants to merge 2 commits intomainfrom
fix/cavity-3eq-issues
Open

fix: correct physics bugs in 3-equation ice cavity melt model#875
JanStreffing wants to merge 2 commits intomainfrom
fix/cavity-3eq-issues

Conversation

@JanStreffing
Copy link
Copy Markdown
Collaborator

Five physics/correctness fixes in cavity_heat_water_fluxes_3eq:

  1. Conduction term (ep31): use rhor^2 instead of rhor, because the ice draft (zice) must be converted to full ice thickness via D_ice = D_draft * (rhow/rhoi), and the heat balance is normalised by rhow. Verified against Holland & Jenkins (1999). ~11% correction.

  2. Reynolds number length scale: replace hardcoded 10 m with actual top-layer thickness hnode(nzmin,node), restoring the velocity- and geometry-dependent formulation from the original BRIOS code.

  3. Depth-to-pressure conversion: pass abs(zice)*1.025 [dbar] instead of abs(zice) [m] to potit() for the in-situ temperature calculation.

  4. Latent heat lhf: correct value from 3.33e5 to 334000 J/kg to match the L constant declared in the same subroutine.

  5. Discriminant guard: add max(..., 0.) before sqrt() to prevent NaN from negative discriminant under numerical extremes.

Also fix comment: tdif is thermal diffusivity [m2/s], not conductivity.

Five physics/correctness fixes in cavity_heat_water_fluxes_3eq:

1. Conduction term (ep31): use rhor^2 instead of rhor, because the ice
   draft (zice) must be converted to full ice thickness via
   D_ice = D_draft * (rhow/rhoi), and the heat balance is normalised by
   rhow. Verified against Holland & Jenkins (1999). ~11% correction.

2. Reynolds number length scale: replace hardcoded 10 m with actual
   top-layer thickness hnode(nzmin,node), restoring the velocity- and
   geometry-dependent formulation from the original BRIOS code.

3. Depth-to-pressure conversion: pass abs(zice)*1.025 [dbar] instead of
   abs(zice) [m] to potit() for the in-situ temperature calculation.

4. Latent heat lhf: correct value from 3.33e5 to 334000 J/kg to match
   the L constant declared in the same subroutine.

5. Discriminant guard: add max(..., 0.) before sqrt() to prevent NaN
   from negative discriminant under numerical extremes.

Also fix comment: tdif is thermal diffusivity [m2/s], not conductivity.
@JanStreffing JanStreffing self-assigned this Apr 2, 2026
@JanStreffing JanStreffing added the bug Something isn't working label Apr 2, 2026
@JanStreffing JanStreffing requested a review from xihao001 April 2, 2026 09:53
@JanStreffing
Copy link
Copy Markdown
Collaborator Author

Not sure if any of these will be major changes, perhaps not. But can't hurt to be accurate?

@JanStreffing JanStreffing requested a review from GMHuettner April 2, 2026 09:54

real(kind=WP),parameter :: pr = 13.8 !Prandtl number [dimensionless]
real(kind=WP),parameter :: sc = 2432. !Schmidt number [dimensionless]
real(kind=WP),parameter :: ak = 2.50e-3 !dimensionless drag coeff.
Copy link
Copy Markdown
Collaborator

@FinnHeu FinnHeu Apr 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The drag coefficient is the most effective tool to reduce the melt rates which are always too large. Numbers down to 0.5e-3 are reasonable. We could consider reducing the default value in this merge.

!rt s_surf_flux(i,j)=gas*(sf-(s(i,j,N,lrhs)+35.0))

heat_flux(node) = rhow*cpw*gat*(tin-tf) ! [W/m2] ! positive for upward
water_flux(node) = ep5*(sf-sal)/sf ! [m/s] !
Copy link
Copy Markdown
Collaborator

@FinnHeu FinnHeu Apr 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to Ralph Timmermann the fix in #855 was wrong and line 365 must be replaced with water_flux(node) = rhow / rhofw * gas * (sal-sf)/sf, with rhofw=1000kg/m3 #866 . Replacing ep5 with rhow / rhofw * gas seems correct to me based on Ralph's explanation. Replacing sf-sal with sal-sf would however switch the sign of the term and would require according modification of the code.

Copy link
Copy Markdown
Collaborator

@FinnHeu FinnHeu left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are good fixes, @JanStreffing. In my opinion, we should consider to additionally implement Losch 2008 (https://doi.org/10.1029/2007JC004368), which solves the problem of missing boundary layers right underneath the ice due to dilution of the heat and freshwater fluxes in coarse (or unevenly thick) z-coordinates. The general idea is to introduce an imaginary boundary layer of a certain thickness (e.g. 10 m) under the ice tongue no matter the thickness of the first active layer. The heat and freshwater fluxes are computed according to the temperature, salinity, and velocity of this boundary layer (which is relaxed to the full cells T and S to follow its evolution) and ultimately applied to both, the full cell and the boundary layer itself. This way, a cooler and fresher boundary layer develops and is used to compute the fluxes, which should further reduce the fluxes. I wanted to properly try this anyway and would rank it higher on my agenda.

@GMHuettner
Copy link
Copy Markdown
Collaborator

Trying these fixes in AWI-ESM2 with fesom-2.6.11, these fixes seem to have very little effect on the overall freshwater fluxes. Over 10 year test runs there is little difference to before, both in the total flux and specific cavity behaviour

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants