Skip to content

Fix aEVP rheology: restore compressive resistance for stationary ice (delta_min floor in stress_tensor_a)#946

Draft
patrickscholz wants to merge 2 commits into
mainfrom
workbench_fix_aevp_delta_min
Draft

Fix aEVP rheology: restore compressive resistance for stationary ice (delta_min floor in stress_tensor_a)#946
patrickscholz wants to merge 2 commits into
mainfrom
workbench_fix_aevp_delta_min

Conversation

@patrickscholz

Copy link
Copy Markdown
Contributor

Summary

A one-word fix to stress_tensor_a in the adaptive EVP (aEVP) ice rheology: the isotropic stress term r1 was missing the delta_min regularisation floor that the mEVP variant (stress_tensor_m) already applies correctly. As a result, stationary sea ice in the aEVP scheme has zero compressive resistance — it goes mechanically limp the moment it stops deforming.


Background: the role of delta and delta_min

In EVP rheology, delta is the total deformation rate invariant:

delta = sqrt(eps1² + vale*(eps2² + 4*eps12²))

It measures how fast the ice is deforming (divergence + shear). When ice is stationary, all strain rates are zero and delta = 0.

delta appears in the denominator of the ice pressure term:

pressure = P* * exp(-c*(1-A)) / (delta + delta_min)

delta_min (typically ~2×10⁻⁹ s⁻¹) is a small regularisation floor with two purposes:

  1. Prevent division by zero in the pressure denominator when delta = 0
  2. Maintain a residual compressive stress at rest (see below)

The bug: isotropic stress collapses to zero at rest

The stress tensor update for the isotropic (compressive/divergent) component is:

r1 = pressure * (eps1 - delta)

where eps1 = div(u) is the divergence. This term drives the isotropic part of the internal ice stress. For converging ice (eps1 < 0), r1 < 0 — the ice resists compression. For diverging ice, r1 > 0.

When ice is stationary: eps1 = 0 and delta = 0, so:

r1 = pressure * (0 - 0) = 0

The isotropic stress contribution vanishes entirely. Over the aEVP sub-cycling iterations, the stress relaxation term det1 * sigma then damps the existing stress state toward zero with no restoring force — the ice loses all memory of its compressive stress state without any physical forcing causing the relaxation.

In practice this means:

  • Land-fast ice that is held stationary by coastline geometry offers no internal resistance to wind or ocean-current forcing on the first sub-cycle. Stress has to be rebuilt from scratch each time the ice momentarily stops.
  • Converging pack ice that decelerates to rest momentarily loses its compressive rigidity, allowing unphysical yielding.
  • The severity scales with how often ice transitions through zero velocity, which is frequent in winter pack dynamics and wherever ice abuts land boundaries.

Why mEVP does not have this problem

stress_tensor_m (line 171) already applies the delta_min floor:

r1 = pressure * (eps1 - max(delta, ice%delta_min))

When ice is stationary: max(0, delta_min) = delta_min, so:

r1 = pressure * (0 - delta_min) = -pressure * delta_min

This is a small but nonzero compressive stress that is maintained even at rest. It ensures the ice "remembers" it should resist compression, consistent with the plastic yield curve. The aEVP variant was missing this guard — likely an oversight when stress_tensor_a was implemented separately from stress_tensor_m.


The fix

One-word change in stress_tensor_a (src/ice_maEVP.F90, line 1095):

! Before (wrong)
r1 = pressure * (eps1 - delta)

! After (correct)
r1 = pressure * (eps1 - max(delta, ice%delta_min))

This makes stress_tensor_a consistent with stress_tensor_m and restores physically correct compressive resistance for stationary ice in the aEVP scheme.


Scope

  • Only affects runs using the aEVP rheology variant (stress_tensor_a)
  • No effect on mEVP runs (stress_tensor_m was already correct)
  • No namelist changes required
  • Effect is most visible in land-fast ice regions and during winter convergence events where ice periodically decelerates to rest

In stress_tensor_a, r1 = pressure*(eps1 - delta) collapses to zero when
ice is stationary (eps1=0, delta=0), removing all compressive resistance.
stress_tensor_m already guards against this with max(delta, delta_min).
Apply the same fix to keep the aEVP variant consistent.
@patrickscholz patrickscholz added this to the FESOM 2.8 milestone Jun 26, 2026

@suvarchal suvarchal left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Bug hat on!

@JanStreffing

JanStreffing commented Jun 26, 2026

Copy link
Copy Markdown
Collaborator

This sounds like it needs a test as well. Would be good to quantify. I have no feeling for how of vel goes through zero. I convert to draft to make the shiny green merge button go away for now. Esp. cause aevp is widely used, e.g. in our CMIP7 and all sorts of high resolution projects.

@JanStreffing
JanStreffing marked this pull request as draft June 26, 2026 12:08
@patrickscholz

patrickscholz commented Jun 26, 2026

Copy link
Copy Markdown
Contributor Author

@JanStreffing: Did you not use mEVP in CMPI7 ? Than i will also do some bias runs for that one!

@JanStreffing

Copy link
Copy Markdown
Collaborator

We used it for much of the spinup, but then switched to aEVP. mEVP had proved to be the reason for max ice thickness unphysically high in some Greenland fjords.

@patrickscholz

Copy link
Copy Markdown
Contributor Author

done for aEVP

  • difference in Summer/Winter NH/SH sea ice concentration
a_ice_NH_nps_y2000-2019_mDJF a_ice_NH_nps_y2000-2019_mJJA a_ice_SH_sps_y2000-2019_mDJF a_ice_SH_sps_y2000-2019_mJJA
  • maximum NH/SH MLD2
MLD2_NH_nps_y2000-2019_mMar MLD2_SH_sps_y2000-2019_mSep
  • temperature/salinity bias
salt_rob_y2000-2019_dep10 0m temp_rob_y2000-2019_dep10 0m

Fazit: bugfix has very small effect !

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

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants