GOSI10: Investigate sea ice issues in the Obb Estuary

[emmafiedler]

This mass budget assessment plot shows erroneous frazil ice growth in June in NEMO5 (in both GOSI.beta.4 and GOSI.beta.5 configurations). This is occurring in the Obb Estuary, which is located along the north fold. This is not seen for GOSI runs using NEMO4.2.2.

These runs are for 1958-2022. The problem does not occur in the first three years of the run, but is seen in the 1992-2021 average used for the plot.

[cplltr]

Hi @emmafiedler

We (ECMWF) experienced similar issues in the same area (Ob river delta). In our case, the model would supercool the water (down to -20ºC) instead of forming sea ice in the winter, when losing surface heat. This is due to that river delta being very fresh (almost zero salinity), and SI3's thermodynamics being ill-defined in very fresh waters.

Sea-ice enthalpy (internal energy) is derived from eq. 5 of Bitz et al. 1999, and the last ratio (mu*S/T) gets close to zero divided by zero when S, T -> 0.
That enthalpy formulation is used in the open water sea-ice formation routine. If the model's cooled down at the surface, then the model will evaluate the sea-ice enthalpy, and compare that to the residual (here "residual" means, once the energy needed for bringing SSTs to freezing point) energy loss. If there's still energy to be lost, then the model should create sea ice. Our problem was that the sea-ice enthalpy was numerically bogus because of the 0/0 problem referred to above, so no sea ice got formed, and then the residual energy loss was consumed by further cooling down the ocean. Issue only occurred when using single precision NEMO though, which I think is consistent (the 0/0 issue is easier to hit with single precision).

I think @dcalve recently put (or reported) of a patch for this in the NEMO trunk. It looks like it would have addressed our issue, but in any case, that patch worked for addressing ours back then:

diff --git a/nemo/NEMOGCM_V40/src/ICE/icethd_do.F90 b/nemo/NEMOGCM_V40/src/ICE/icethd_do.F90
index fb7ff8d9a8..ca061e4d94 100644
--- a/nemo/NEMOGCM_V40/src/ICE/icethd_do.F90
+++ b/nemo/NEMOGCM_V40/src/ICE/icethd_do.F90
@@ -44,6 +44,7 @@ MODULE icethd_do
    REAL(wp) ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
    REAL(wp) ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
    REAL(wp) ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
+   REAL(wp) ::   rn_opw_fresh  ! minimum sss_1d under which opw sea ice formation is assumed to happen as if in freshwater
 
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -299,9 +300,14 @@ CONTAINS
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = ji1, ji2
             ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
-            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
-               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
-               &                      - rcp   *         ztmelts )
+
+            IF( sss_1d(ji) < rn_opw_fresh ) THEN
+               ze_newice(ji) = rhoi * rLfus
+            ELSE
+               ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
+                    &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+                    &                      - rcp   *         ztmelts )
+            ENDIF
          END DO
 
          ! --- Age of new ice --- !
@@ -477,7 +483,7 @@ CONTAINS
       !!-------------------------------------------------------------------
       INTEGER  ::   ios   ! Local integer 
       !!
-      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
+      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, rn_opw_fresh
       !!-------------------------------------------------------------------
       !
       REWIND( numnam_ice_ref )              ! Namelist namthd_do in reference namelist : Ice thermodynamics
@@ -498,6 +504,7 @@ CONTAINS
          WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom     rn_maxfraz = ', rn_maxfraz
          WRITE(numout,*) '      Threshold relative drift speed for collection of frazil   rn_vfraz   = ', rn_vfraz
          WRITE(numout,*) '      Squeezing coefficient for collection of frazil            rn_Cfraz   = ', rn_Cfraz
+         WRITE(numout,*) '      Minimum sss_1d under which opw ice formation is assumed to happen as if in fresh conditions rn_opw_fresh  = ', rn_opw_fresh
       ENDIF
       !
       IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )

which essentially says to assume "sea ice" is actually fresh under a certain SSS threshold (I think I used 0.1 PSU).

Not sure that's the same issue you're having, but it does look similar-ish. I think the key point is the 0/0 division is eq. 5 of Bitz 1999. It can create problems in different parts of the model, IMO.

Charles

[dcalve]

The patch @cplltr mentioned can be found here. It looks like the exact same expression as shown for icethd_do.F90, except that code already had the divide-by-zero fix.

icethd_do.F90

ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
&                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
&                      - rcp   *         ztmelts )

icevar.F90

e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
&                   + rLfus * ( 1._wp - ztmelts / MIN( t_i_1d(ji,jk) - rt0, zepsi ) )   &
&                   - rcp   * ztmelts )

In our case, this issue only manifested when we switched to our current HPC and compiler, and it led to infinite values everywhere.

I can see other occurrences of this same code pattern elsewhere in the code, but they have already been similarly fixed (with the exception of OCE/BDY/bdyice.F90, which would only affect regional configurations).