What happens when ice is at minimum thickness in CICE/UM?

[ofa001]

Spencer reported that when Ice was dropping below the minimum thickness there was some odd behaviour as ice concentration in the mediator was not being mirrored in the UM,

[ofa001]

hi @bimlim @anton-seaice In Icepack_therm_itd.F90 there is a remapping occuring in relation to ice volume when the ice drops below the ice volume, which also adjusts the pond variables.
As I said on the other thread the memory of switching betwen multi-layer and zero layer schemes goes further back to the "CSIRO" models which happened in both ice and snow schemes.

Here is the code.

! Make sure hice(1) >= minimum ice thickness hi_min.
!-----------------------------------------------------------------

     if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then

        da0 = aicen(1) * (c1 - hicen(1)/hi_min)
        aicen(1) = aicen(1) - da0
        hicen(1) = hi_min

        if (tr_pond_topo) &
           fpond = fpond - (da0 * trcrn(nt_apnd,1) &
                                * trcrn(nt_hpnd,1))
     endif

I will next check the flux issue to the ocean which will be covered under icepack therm vertical or icepack_therm_bl99.

Then its a case of what the mediator to UM was seeing in Spencer's plots..

[ofa001]

hi @blimlim @kieranricardo
This routine is in the original GSI8.1 set up Kieran reworked it a bit and turned it into a function, so it may be what you
are talking about in relation to the conductive fluxes and extra being used to alter the thickness by melt or growth this is from ice_therm_vertical.F90 in cice_gsi8.1 branch on cice5 repo.

! This routine is only called if UM-style coupling is being used, with top conductive flux as forcing.
! Check top conductive flux, ice thickness and top layer temperature.
! If the ratio of flux to thickness is too high, remove some of the flux and put it into fcondtopn_extra.
! If the top layer temperature is getting too low, and the flux is negative, also put some into fcondtopn_extra.
! The remainder, fcondtopn_solve, goes to the thermodynamic solver. fcondtopn_extra is added to the energy balance
! at the bottom of the ice in thickness_changes, and is thus used to grow / melt ice at the bottom.
!
! author Alex West, MOHC

  subroutine cap_conductive_flux(nx_block,ny_block,my_task,icells,indxi,indxj,fcondtopn,fcondtopn_solve,fcondtopn_extra,hin,zTsn,zTin,hslyr)  

etc...

[ofa001]

HI @blimlim, The excess flux passing to the ocean comes from ice/icepack therm bl99 using an array enum. This happens when the temprature solver gets into a loop where it cant solve and it needs assistance to help it to solve.

          if ((l_brine) .and. zTsn(m,k)>c0) then

! Alex West: return this energy to the ocean

              dqmat_sn(m,k) = (zTsn(m,k)*cp_ice - Lfresh)*rhos - zqsn(m,k)

              ! Alex West: If this is the second time in succession that Tsn(1) has been
              ! reset, tell the solver to reduce the forcing at the top, and
              ! pass the difference to the array enum where it will eventually
              ! go into the ocean
              ! This is done to avoid an 'infinite loop' whereby temp continually evolves
              ! to the same point above zero, is reset, ad infinitum  
       its not depending on the hi_min  setting at 0.2m  and hs_min at 0.1 min  though it may be sensitive to that  and its on both the snow and ice  temperature calculations.

For example here is some of the snow loop.

     if (l_snow(m) .AND. k == 1) then
                 if (Top_T_was_reset_last_time(m)) then
                    fcondtopn_reduction(i,j) = fcondtopn_reduction(i,j) + dqmat_sn(m,k)*hslyr(m) / dt
                    Top_T_was_reset_last_time(m) = .false.
                    enum(m) = enum(m) + hslyr(m) * dqmat_sn(m,k)            
                 else
                    Top_T_was_reset_last_time(m) = .true.
                 endif
              endif

              zTsn(m,k) = min(zTsn(m,k), c0)

You can look at code to follow through on enum and this code is form Gsi8.1 hence the (i,j)'s in the code.

Still puzzled what you were seeing in the mediator v the UM, this is simply some of the flux being reduced to help the temperature solve solve, its the same as if it went into the open water (ledas), or if it had gone in as a penetrative shortwave flux, through the ice which we still need to fix and reached the ocean that way, its small.

[blimlim]

Hi @ofa001, thanks for writing this up. I think it will take some time for me to properly understand the CICE code! I can add some details on what I think is happening on the UM side though.

I've attached the graph I shared yesterday below. It shows the thinnest category ice concentration exported by the mediator to the UM, the concentration output by the UM, and their difference, all at the first timestep:

A large part of the thinner ice is missing in the UM output. This impacts the total ice area of this category:

Meanwhile, the thicker categories have much better agreement in their areas at different stages:

It looks like the ice disappears when we process the imports to the UM. We're first clipping the ice fractions based on their thickness:

 ! clip ice fractions
      if (ice_fract_cat(i,j,k) >  0.0) then
        if (ice_thick_cat(i,j,k) / ice_fract_cat(i,j,k) < hi_min) THEN
          ice_fract_cat(i,j,k) = 0.0  ! other fields will be zeroed next
        end if
      end if

where hi_min=0.2, ice_fract_cat is the imported aicen, and ice_thick_cat is the imported vicen. This occurs in um_cap_import_export.F90 here.

If I take the mediator output and manually apply the above clipping to it, we get something much closer to the UM's output:

Applying the second processing step (clipping fractions < 0.01) then recovers the UM's output almost exactly:

[blimlim]

I'm not too sure what the different effects of this are UM side. Presumably it would reduce the albedo in grid cells with a large fraction of thin ice, and it might also impact the surface temperature calculated in um_cap_import_export.F90 using the aggregate ice fraction.

As @kieranricardo mentioned, this comes from the oasis coupling. There's a little bit of explanation there, which I don't properly understand – apparently when the zero layer scheme is used, the ice concentrations will be adjusted but ice volume won't be deleted, but when multilayer ice is used, thin ice is just deleted.

This also matches what was done in CM2 which has hi_min=0.2 hardcoded in oasis_geto2a.F90. It's interesting because 0.2m ice would still be given a much higher albedo in the jules_ssi_albedo routine if it were to exist!

The 0.2m hi_min is in the recommendations in the rose-meta.conf file for UM13.0:

[namelist:coupling_control=hi_min]
compulsory=true
description=Minimum sea ice thickness to couple
help=The minimum sea ice thickness to couple. Sea ice thicknesses below this have either:
    =  1) if multi layer sea ice = sea ice areas set to zero
    =  2) if single layer sea ice = sea ice area reduced to thickness/hi_min
    =Recommended settings are:
    =  * CICE with multi layers (l_sice_multilayers = True)  hi_min = 0.2
    =  * CICE with single layer (l_sice_multilayers = False) hi_min = 0.01
    =  * SI3 with multi layers  (l_sice_multilayers = True)  hi_min = 0.05
    =  * SI3 with single layer  (l_sice_multilayers = False) hi_min = 0.01
ns=namelist/Coupling/Main coupling controls
range=0.0:
sort-key=df12
type=real

[ofa001]

Just pulled this from oasis3_geto2a which is just before the code that Kieran has duplicated across.

          ! Here we ensure that regardless of what has happened to the ice
          ! fraction and GBM ice depth during the coupling,
          ! the minimum local ice depth is 1cm (as in CICE). This is done by
          ! squashing the ice to a smaller area whilst conserving volume (so
          ! the GBM ice depth is unchanged).
          !
          ! Note that freezen is unchanged unless the local ice depth (i.e.
          ! hicen/freezen is less than 1.0e-02, in which case freezen is
          ! decreased so that hicen/freezen=1.0e-02

          ! If using multilayers....
          ! **CHANGE the way min local ice depth is enforced and increase
          ! min ice depth to hi_min, set at top.  Get rid of ice below that
          ! value - do not squash
          ! **INCREASE min ice concentration
          !
          ! Do not apply these limits to first-order ice concentration,
          ! which is used to produce an energy-conserving field of local
          ! fluxes going to the ocean

The thing is we did "squash" in the remap step in ice_therm_itd when we conserved volume on the CICE side to retain the thickness but reduce the area.

One of the other issues isMet office used the two step "time travelling" approach which mean they did have two time levels of ice concentration which as they say above they used to do produce the energy conserving fluxes to the ocean, Kieran has decide he doesn't need to do this, but is now following this step, which means he is not energy conserving after all.

What I also noticed Spencer in your 30 day graphs is that the result was noticeable in the first few days but adjusted away, by the end of the month, now that may be because of the initial condition which came from a fixed initial set up in the ACCESS-OM3, using ocean SST's and an ice thickness profile, and arbitrary AMIP run in the UM so well out of sync, and once the system has adjusted, this issue is going to be a minor, difference, possibly above round off error. I was surprised to see concentration errors of up to 0.5 "as the remapping step=squashing" had taken place in CICE and if this ice was less than hi_min it would reduce, it looked like these errors were outside the region set by the CICE/OM-3 initial condidtion, that means ice was advecting there, but probably thin ice, and some of it would be numerical diffusion at the very limit. After a month or so the two models would be more in sync.

Alos as you can see the hi-min limits are fairly arbitrarily chosen (trial and error at a guess). Its likely that CICE at 0.1 or 0.15 will mostly work its just they may have had problems early on

[ofa001]

I guess one question is to run a test on how far out it is if you do a restart in year2 or year 10, when the early mismatch between models is no longer an issue.

But the other issue in relation to the fluxes, yes thats where the 2 time levels of the ice concentration came in to get better fluxes, its slightly different because, Kieran is doing the ocean time step ahead of the atmosphere one, the Met office had the UM leading the ocean, I will send Spencer the document.

[blimlim]

Hi @ofa001, that's a good point about the large amount of ice being deleted in the UM despite us "squashing" it in CICE:

I was surprised to see concentration errors of up to 0.5 "as the remapping step=squashing" had taken place in CICE and if this ice was less than hi_min it would reduce, it looked like these errors were outside the region set by the CICE/OM-3 initial condidtion, that means ice was advecting there, but probably thin ice, and some of it would be numerical diffusion at the very limit. After a month or so the two models would be more in sync.

I couldn't work out what value of hi_min we were using to squash the ice in CICE, as there were a few values with the same name appearing in different places. I tried printing out what we actually use just before the squashing in icepack_therm_itd.F90:

      !-----------------------------------------------------------------
      ! Make sure hice(1) >= minimum ice thickness hi_min.
      !-----------------------------------------------------------------
         write(6,*) "SPENCER hi_min: ", hi_min
         if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then

            da0 = aicen(1) * (c1 - hicen(1)/hi_min)
            ...

Which printed out:

job.out: SPENCER hi_min:   1.000000000000000E-002

And so I think CICE is enforcing thickness > 1cm and the UM is enforcing thickness > 20cm – hence the large areas being deleted. I'm running a two year simulation now, and will upload some results when their done.

[ofa001]

Thats worrying, perhaps the way Kieran has implemented hi_min in icepack hasn't worked, in CICE5 it went into ice_itd.F90 perhaps in CICE6 it needs to be in icepack_parameters? I Will check what he has done. It also the aicenmin (met office) changes need to be passed to some of the subroutines cicecore, not sure about hi_min. I should have checked more carefully how he put these met office fixes in.

[blimlim]

@MartinDix pointed out we are currently trying to set hi_min=0.2 in cice here however this is appearing under an if block that's only active when kcatbound==0 and kitd!=1.

In our current config, we set kcatbound=0, and leave kidt to take the default value which I think is 1, which might be why the hi_min setting isn't taking effect.

[ofa001]

Hi @bimlim, I think I told Kieran that I didnt think that wasnt a good idea and he should put it in more general code with an access flag or in icepack parameters, as those settings could change.

[ofa001]

HI again I checked kitd the other day in the namelist (ice_in) but I just did it again its set to 1 so this code is not being reached. I guess Kieran put it there because it was logical next to the cesm coupled setting but I suggested he follow the Met office approach and have it which is right outside the kcatbound loops. and was not linked to any kitd and kcatbound settings.

Their code looks like this

! AEW: (based on Alison McLaren's vn4 modifications) Set a higher value
! of aicenmin if we're using multilayers with UM-style coupling.
! Also allow higher values of hi_min, hs_min to be set (this is a
! bit ad-hoc).
! Set aicenmin - this may want to be done via a namelist in future
!-----------------------------------------------------------------

  if (heat_capacity) then
    ! Set higher values to help with stability
    aicenmin = aicenmin_ml    ! 1.e-5.  Changed from 1.e-2
    hi_min   = p2     ! 0.2m
    hs_min   = p1     ! 0.1m
  else
    aicenmin = puny   ! Standard CICE setting
  endif

I am looking to see if Kieran has set aicenmin anywhere, as that is useful across as I said in cicecore (ice history) area and Met office also had a setting on hs_min again for stability.

'heat-capacity' was removed in CIC6.5 if it was the flag that if false switched to zero layer so no, I was really wrong when I was put on the spot the other day, the model would not jump back form multi-layer to zero layer (a single layer) in any 2024 CICE version. OOPs :(

[ofa001]

Hi @blimlim I found the change when Kieran put these on the branch no he did the absolute minimum he changed the hs_min in icepack_parameters and didn't introduce aicenmin at all though I think he did introduce it on the UM side as he copied their part of the code more closely there.

So he didnt do what was advised, and luckily you have found the consequences of the bug, and its shown up right at the beginning of the run when the 2 models are out of sync. The real issue now though, does not using the semi-implicit method, and the lack of the intermediate ice concentration which the Met office where using for the fluxes, mean the system is non conserving, as you are worrying about conservation through the mediator, this may be a good time to test it.

[ofa001]

Hi @bimlim this probably belong under a separate topic but by Kieran not including the changes to aicenmin which was used by the Met office in the zap small areas subroutine, and in some of the ice_history criteria, but copying the same name numbers from the Um into his "cap" again conservation has been broken, as this routine accounted for some of the conservation on the CICE side which was is not handled in the UM.

I should have been checking his code branchs more carefully not just that he did implement as following the cm2 approach and their approach when we spotted he had an error there not just the partial changes he has put in.