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.

[blimlim]

Hi everyone, I just have a small update on the treatment of thin ice. I roughly set the hi_min=0.2 in this branch to match the UM.

Running for 3 months gives the following total areas for the thinnest ice category in the UM and CICE:

Compared to our previous run which had hi_min=0.01:

The update reduces the discrepancy between CICE and the UM's ice areas, but doesn't remove it completely. The UM sees less ice especially at the start of the run:

Plotting the concentrations, the UM still ignores a lot of the ice around 65N on at the start of the simulation:

Quite a bit of the ice in this area ends up having a thickness slightly below 0.2, despite the resizing done in icepack_therm_itd.F90

I'm wondering if any of the CICE steps between after resetting the thickness to hi_min=0.2 in linear_itd but before the coupling could reduce the thickness below this 0.2 limit?

[ofa001]

Hi @bimlim In the icepack_step_therm2 the calls after linear_itd are in order "add new ice", "lateral melt" (which can reduce area) then "cleanup itd".

Look its highly likely some of the difference several weeks into the run once you are past the "initialization shock" could be coming from the "add new ice" which handles frazil ice which is way less than 20cm in the field it is set to hi0new which will normally lie between these settings in the lowest ice category, but is also set above hfrazilmin =0.05m
! hin_max(0) < new ice thickness < hin_max(1), but is very occasionally distributed across all the categories.
from icepack_itd for kcatbound=0,kitd=1
hin_max(0) = 0 and
do n = 1, ncat
x1 = real(n-1,kind=dbl_kind) / rncat
hin_max(n) = hin_max(n-1) &
+ cc1 + cc2*(c1 + tanh(cc3*(x1-c1)))
enddo

In "lateral melt" the lateral melt parametrization does change the volume and area, but does not specifially change the thickness.

The "cleanup itd" routine is the one that has the zap small area routines where aicenmin=1e-2 or 1e-5 matching the value in the UM should be set as in CM2 set up, I would suggest start with 1e-5 its greater than puny 1e-11 but not down at 1% ice cover levels, and test and see how sensitive the system is, switch value on the UM side to match, and will clean up numerical noise levels of ice left behind, it will mostly effect the lowest ice category as thicker categories melt down, normally.

There are other calls to "cleanup itd", however, from the dynamic parts of the code in particular the icepack mechred, the ice thickness can change through ridging/rafting so its highly possible that that may need to be checked, the call to the dynamics routines follows thermodynamics in the code.

I am beginning to wonder if we should be duplicating the hi_min=0.2 and aicenmin=1_e-2/1e-5 settings on the UM side of the coupler and its regridding with so many additional adjustments in CICE after they have been imposed taking into account those values are only approximations to the correct settings for stability. What happens if we temporarily switch them off and just have the hi_min and aicenmin (zap small area higher settings on) what do your curves look like after a 3 month run.

[ofa001]

hi @bimlim Was just checking the code in the UM in ESM1.5/6 on this issue the had hi_min at 0.01 as we knew but aicenmin on the atmosphere side was 2e-4 whilst it was 'puny' on the CICE side, and they accept there is squashing (maintaining volume) going on in the comments, unlike what they say in the more upto date models.

I am planning implementing at least a value of 2e-4 on the CICE side for aicenmin ESM1.6.

[blimlim]

Hi @ofa001, apologies on the lack of updates on my end! I've been a bit preoccupied with ESM1.5 work the last couple of weeks and so haven't had a chance to look at this in much detail yet.

I just have a couple of plots from some runs I'd done earlier which I'll upload for discussion in the next meeting. The first is a continuation of the run with hi_min=0.2 in both CICE and the UM, where I've run it for an additional year:

The difference in the thinnest category ice area remains fairly stable at around 40km^2 for most of the run, with a jump around late October. However the increases are still much smaller than at the very start of the run.

I then tried commenting out the sections of the UM which cull the ice based on thickness or ice fraction. In CICE, I'd left hi_min=0.2 and aicemin=puny, though will try again with a larger aicemin.

We see much smaller differences in the thinnest category ice fractions than before.

@kieranricardo mentioned that the UM sees ice that's less than 0.2m most likely as a result of the regridding, which doesn't conserve the mininum thickness, rather than additional steps in CICE. I'm currently running an extra simulation with timestep output to understand this better.

[ofa001]

Hi @blimlim, @kieranricardo is right the re-gridding is also going to reduce the ice below this minimum thickness, but I think the large change you got in the initial start up is due to frazil and dynamics, which came about in the mismatch between the initial condition in the ice-ocean set up and the atmosphere set up, the atmosphere was cold over a region of ocean that didn't have any ice on it in the CICE set up, and it need to form or advect there.

[blimlim]

Hi @ofa001, the following are from a simulation with hi_min=0.2 and aicemin=1.0e-5 on the CICE side, while the UM is set to use all the ice it receives. This is the same as for the last two plots my last comment, but with CICE's aicemin increased from puny.

The following shows the total ice area for the thinnest category, in CICE and in the UM:

The following shows the difference between the two:

It looks similar to the previous run with aicemin=puny on the CICE side, and I still get a small difference in area.

CICE's aicemin was previously increased from puny to 1e-2 and then reverted back to puny due to concerns about water conservation. I'll try to check with these latest runs if there are any impacts on water conservation.

[ofa001]

Hi @bimlim you say this above "CICE's aicemin was previously [increased from puny to 1e-2 due to concerns about water conservation" where did you find that, I saw them testing between 1-e-2 and 1-e-5 in the gsi8.1 version, that we picked up but we stayed with 1e-5 for cm2. In relation to water conservation, in the 'zapping small areas section" we save all the the fresh water, snow melt terms involved.

I spotted it was Kieran who suggested that it was a water conservation issue, but there were other water conservation issues in the code as well at the time, this wasn't it. I first thought you might have been talking about the "amin" term in the ice dynamics" which I think @anton-seaice said they now think they can run at a much smaller value in CICE6, I haven't considered the conservation issues in relation to this. I will think about it get back to you.

[blimlim]

Thanks @ofa001. Just to clarify, I was referring to the CM3 tests that @kieranricardo had done.

[ofa001]

Thanks @blimlim Didnt @kieranricardo have another water conservation issue in his model runs. Could he pin it down to coming from this ice change, the 'fresh' and 'salt' budgets in to 'zap small areas' are linked to aicenmin, and then fed through to the ocean terms, should work, I should check where he did it.

[ofa001]

@bimlim yes the water conservation is handled as done as the increments dfresh dfsalt which then get accumulated through the different model subroutines, the value of aicenmin wont impact it in that subroutine.

[ofa001]

Hi @bimlim, it looks like the CIC6 a_min (p001) in the ice_dyn_shared is the same as in CICE5 so hasnt been changed but from Anton they may be discussing the value.

[blimlim]

Thanks @ofa001,

I thought I'd write up some tables of the parameters used in the different models based on your earlier email, in case it is useful to have them in a central place for discussion in the next meeting.

CM3

	CICE	UM
hi_min	0.01	0.2
aicenmin	puny	0.01

The CICE hi_min=0.01 setting is due to a bug, and will be updated to 0.2.

CM2

	CICE	UM
hi_min	0.2	0.2
aicenmin	1.0e-5	0.01

ESM1.5

	CICE	UM
hi_min	0.01	0.01
aicenmin	puny	2.0e-4

[blimlim]

@kieranricardo mentioned that there was an additional conservation issue due to some rivers not connecting to the coastline, however increasing the minimum ice fraction for CICE's zapping did cause separate conservation errors. I'm currently running some simulations to reproduce this.

[ofa001]

Thanks @bimlim, I havent run a value for ESM1.6 yet but I was going to recommend aicenmin at 2e-4 on the CICE side, so

puny, it could still be 1 e-5.

[blimlim]

Hi @ofa001,
I've just put some of my simulations through a notebook for calculating water conservation from @kieranricardo. The following shows the change in total water mass across the models over time, for a few different values of aicenmin set in the zap_small_areas subroutine.

As @kieranricardo found earlier, there is a big impact on conservation when using a larger (i.e. 1e-2) value.

[ofa001]

Hi @bimlim that maybe why the comment statement in gsi8.1 said between 1e-2 and 1e-5 and the actual value used was 1e-5, and 1e-2 on the UM side, though I think that is too high.
We used 1e-5 in CICE in CM2, and I guess 1e-2 on the UM side.

I am suprised, it was so large though, it must represent additional, terms in the new ice growth and dynamics part of the code that are not being accounted for correctly, by just the zap small area code. 1e-2 is 3 orders of magnitude larger, and 1% of the ice cover though.

[blimlim]

Hi @ofa001, just adding a summary of the last meeting to keep track of tasks for me to do:

• Add a PR to set hi_min = 0.2 and aicenmin = 1.0e-5 in icepack, matching the values used in CM2.
• Run experiments varying hi_min simultaneously in icepack and the UM. Try e.g. hi_min=0.1 in both models.
• Try reducing the aicenmin in the UM. Try UM values 1e-3, 1e-4, while leaving icepack's aicenmin=1e-5. Investigate the sensitivity of the fluxes from the UM to these changes (I might need to get back to you to ask about the best things to look at).

Let me know if there was anything I've missed here, or if there's any other tests you would be interested in seeing.

[ofa001]

Hi @blimlim No that sounds reasonable summary of what we suggested, I will have a think on the test on the fluxes and how we can check it.