Residual ice

[eclare108213]

Very small volumes of ice appear and/or remain in grid cells with no mechanism for getting rid of them.

If there is already an issue in either CICE or Icepack repos for this, please point it out. Potentially related issues are #564, CICE-Consortium/Icepack#148 and CICE-Consortium/Icepack#149.

There are a few other potential contributing factors:

1. The dynamics has a minimum ice area and mass for being activated, below which the ice does not move. The assumption here is that the thermodynamics will either melt such small amounts of ice or they will continue to grow until large enough to move.

     a_min = p001         , & ! minimum ice area
     m_min = p01              ! minimum ice mass (kg/m^2)

2. The thermodynamics has a minimum thickness for being activated, due to CFL considerations (hi_min=0.01 or 0.1). The assumption here is that additional frazil ice will be sent from the ocean in future time steps, or that the dynamics will ridge the ice into thickness large enough for the thermodynamics to handle.

  hi_min = p01    ! minimum ice thickness allowed (m) for thermo
                  ! note hi_min is reset to 0.1 for kitd=0, below

These 2 factors are obviously a catch-22 and there could easily be a gap in those parameters where nothing will ever happen.

3. Under freezing conditions, frazil ice sent from the ocean fills the open ocean area of a grid cell to a certain thickness, after which it is put onto the bottom of existing ice. If there is not enough frazil ice volume to occupy the open ocean area with a thick enough (for what?) sheet of ice, then it is added with a smaller area coverage. Ice formed in open water is added to the first ice thickness category. This is in subroutine add_new_ice and is closely tied to the ITD.

! a note regarding hi_min and hin_max(0):
! both represent a minimum ice thickness. hin_max(0) is
! intended to be used for particular numerical implementations
! of category conversions in the ice thickness distribution.
! hi_min is a more general purpose parameter, but is specifically
! for maintaining stability in the thermodynamics.
! hin_max(0) = 0.1 m for the delta function itd
! hin_max(0) = 0.0 m for linear remapping

Factor 3 handles frazil ice appropriately upon freezing, provided the parameters are consistent with item 2 above. But the ice could partially melt in the next time step, dropping below the minimums in 1 and 2.

4. Under freezing conditions, a volume of ice could potentially be spread out over the whole grid cell at a thickness near or below the minimums in 1 and 2. The large concentration would limit the atmosphere-ocean exchanges that would allow the ocean to warm up enough to subsequently melt the ice (and could be affecting things like ocean circulation/AMOC). This scenario could be occurring in models with sub-daily coupling at times of the year when the ice forms during colder hours and should melt during warmer hours but doesn't due to inconsistencies between assumptions in the ice (e.g. factors 1 and 2) and the coupler (high ice concentration implies little solar flux to ocean). Some models have a parameter to control the division of ice volume between area and thickness; I do not think we do, and maybe we should.

5. Averaging in the history files also could contribute to this issue, if there's ice present only for a few timesteps and the rest of the month is ice-free. There is an "ice-present" history variable available to help understand this type of thing.

There may be other parameters or factors that I'm not thinking of at the moment.

[eclare108213]

@rgrumbine @TillRasmussen @proteanplanet
Please provide examples of this from your model configurations, showing ice concentration, thickness, volume and mass. Any other information such as time of year would be helpful.

I suspect it's a coupled thermo + dynamics problem that can not be solved just in Icepack, and it may be a more broad coupling issue. @apcraig @JFLemieux73 since you've both recently run QC cases, could you look for this in the output to see if it shows up in the stand-alone CICE model? I think it does.

[daveh150]

For discussion, I think we ran into this issue recently. Liz Douglass here at NRL is using our coupled HYCOM+CICE5.1.2 (yes, older CICE…) without data assimilation, and was investigating the different thermo schemes. She started from Jan 15, 2017. When she uses the ‘mushy’ thermo (ktherm=2), her run is successful. When she changes to ‘linear_salt’ (and ktherm=1), after 17 hours it crashed with a very small ice concentration, 2.85e-6. The log info is below (it only write aice, not vice). I’m looking through the mushy vs. linear_salt methods, but I don’t see if there is a difference in the minimum ice area or thickness. Is this what others are seeing? Thermo energy conservation error istep1, my_task, i, j: 16442727 72 100 100 Flux error (W/m^2) = 1.523171639483836E-003 Energy error (J) = 0.274170895107090 Initial energy = -454195.357306166 Final energy = -375892.611155715 efinal - einit = 78302.7461504512 fsurfn,flatn,fswint,fhocn, fsnow*Lfresh: 2.13535318803204 -18.5648469793598 0.000000000000000E+000 -414.313533052364 0.000000000000000E+000 Input energy = 78302.4719795561 fbot(i,j),fcondbot(ij): -414.313533052364 46.7704253736093 Intermediate energy = -462229.396128675 efinal - einter = 86336.7849729600 einter - einit = -8034.03882250877 Conduction Error = 5.489471579767269E-002 Melt/Growth Error = 0.329065610904763 Advection Error = 0.000000000000000E+000 dt*(fsurfn, flatn, fswint, fhocn, fsnow*Lfresh, fadvocn): 384.363573845766 -3341.67245628477 0.000000000000000E+000 -74576.4359494256 0.000000000000000E+000 0.000000000000000E+000 istep1, my_task, iblk = 16442727 72 1 category n = 1 Global block: 73 Global i and j: 899 603 Lat, Lon: 68.4469529589150 -18.7899398812531 aice: 2.856320136459712E-006 n: 1 aicen: 2.851444384406774E-006 ice: Vertical thermo error [edited to remove original email text]

[eclare108213]

@daveh150 I think that is a different issue - the thermodynamics is running (although it crashes). I suspect that it does not run at all for the residual ice amounts that are the subject of this issue, although I'm not sure about that.

[TillRasmussen]

@daveh150 This look a lot like what I saw in our system (also HYCOM CICE). Is she running with form drag?
The latent heat is numerically very high -3341.67245628477
@mhrib pushed a solution for that.
See #633.
@eclare108213 I will find a example. The low concentration appears as a band around the sea ice. It is places where ice should disappear and in reality has disappeared for a while.

[daveh150]

Thank you Till! In this case she is not using formdrag. I need to do some research on the effect of using formdrag, Would that be potential problem for the latent heat? Thanks! David

[TillRasmussen]

Thank you Till! In this case she is not using formdrag. I need to do some research on the effect of using formdrag, Would that be potential problem for the latent heat? Thanks! David

I dont know if this is limited to formdrag, but in our HYCOM-CICE model we experienced issues with the iteration in the function atmo_boundary_layer. Maybe try to use the constant boundary layer.
We often saw the issue when the 2 meter temperature was higher than the ice surface temperature. This may be the reason why the ice concentration is low. Not sure.

[eclare108213]

There are two separate issues associated with very small ice areas:
A. large regions of small ice mass and/or thickness, e.g. hi < hi_min (@TillRasmussen and others)
B. grid cells with small concentrations of 'actual' ice thickness vice/aice >>1 (@rgrumbine).

I'm tackling A first, initially via increased lateral melt. Unfortunately I am not finding these large areas of small ice concentration in the standard gx3 test case, so debugging is challenging. Please, would someone provide a standalone configuration where this problem is visible? Preferably gx3, to allow quicker runs and debugging.

As aice decreases, the relative importance of lateral melting should increase relative to bottom melting, given the same ice thickness. When the FSD is not in use, we assume a constant floe diameter of 300 m. However, when aice is very small, there is often not enough ice in a grid cell to form even one floe of that diameter, which is physically inconsistent. I've tried 2 things so far:

1. calculating the floe diameter of a single floe containing all of the ice, and using the minimum of that value and 300 m to compute the lateral melting flux. Note this approach introduces a grid dependency.
2. maximizing the lateral melt for floes smaller than 300 m by setting rside=1, which should provide enough heat to melt the entire floe. In practice rside will be less than 1 because fbot > 0 (they are scaled so that the total amount of heat available for melting is less than frzmlt).

(2) seems like a sledgehammer approach, but (1) doesn't seem to make much difference. This could possibly aggravate B, if thinner ice melts out first and small areas of the thickest ice remain.

@TillRasmussen @mhrib I would like to know whether either of these choices helps. My code is here:

https://github.com/eclare108213/cice -b residual
https://github.com/eclare108213/icepack -b iresidual

Currently, it is set up for item (1) above. The (2) configuration can be run by uncommenting the

rside=c1

line in icepack_therm_vertical.F90. Could you try it out in your problematic case, please?

If increasing lateral melting doesn't work, then we can introduce a namelist parameter to optionally zap these small areas for A. (I'm also curious as to whether turning on the FSD helps with this problem, but we can save that for later.)

[eclare108213]

Actually for case (2), uncomment these 3 lines:

if (floediameter < floediam) then
! print*,'floe diameter ',floediameter, aice
rside = c1
endif

[TillRasmussen]

I will look at it. We have not tested on gx3 only on our own grids. I will see what happens with your modifications and send a test example

[eclare108213]

If the mods to lateral melting are insufficient, then I will add a method to remove all ice less than some input concentration, conservatively. I already modified the code to do that using the zap_small_areas routine, but haven't added namelist etc. Either method will require a new parameter in Icepack, which has to be defined and passed back and forth through the driver/Icepack interface. I'd prefer to let the physics do the job instead of brute force removal, so will wait for your tests before fully implementing either option. @proteanplanet if you have other ideas for how to do this, which I'm missing, please pass them along.

[eclare108213]

I have updated the code to make testing a little easier, although namelist options still are not implemented, pending @TillRasmussen's feedback on what's helpful. It's in the same place: my fork, residual branch for cice and iresidual branch for icepack.

• in ice_step_mod.F90, search for 1==0 and change the conditional to allow floe diameters less than 300 m for calculating lateral melt, when there's not enough ice in a cell to make a floe that big
• in icepack_therm_vertical.F90, there is also a small code block commented out that would apply the maximum available heat for lateral melting, as noted above (rside = c1). This provides an upper bound for how much lateral melting could occur with this basic parameterization, and is not intended to be kept permanently.
• in icepack_mechred.F90, search for 1==0 and change the conditional to allow ice concentrations larger than puny to be zapped. This version will conserve. It's currently hard-coded to act for aice < 1.e-8 but that parameter can be put in namelist if we want to keep this option.

Note: I also changed a conditional in icepack_therm_vertical.F90 from aice > puny to aice > 0, which could potentially change the answers via lateral melting, but it's BFB in my 1-year gx3 smoke test. I'm not sure whether this is because there are never ice concentrations < puny at that point in the timestep (probably, having been cleaned up at the end of the prior timestep), or if calculations are not performed on those grid cells later, or both.

[TillRasmussen]

I put the test example on ftp.dmi.dk user anonymous folder tarftp, file testNAAf.tgz.
Extract the tar ball, create a executable independently of this. Modify the runme.sh with the path to the executable and potentially the rundir name. The submission job in cfg/cice.job probably needs to be altered according to the machine as well. Otherwize it should run the test.
The run started on the 1/9 1990 and has run within the coupled system until the 1/6. The example runs 1 day from the 1/6 as standalone. Note that if you view the iceh_ic file then you can see the low concentration with ncview, when you change linear to low.

Running @eclare108213 's branch with the modifcation as seen in
diff_lowcon.odt

The low conentration values are removed. It makes sense that there is a scaling issue between the 300 meter floe size and the concentration of a grid cell when the grid size changes. This setup is a 10km pan arctic run.

[eclare108213]

@TillRasmussen I'm having trouble connecting to ftp.dmi.dk. Error says "the share doesn't exist on the server". If you can send a complete link to the file via email, maybe it would let me in.

Glad to know the low concentration values disappear -- the last option is a pretty big hammer, so not surprising. The next thing is to understand which of the three modifications are needed, and which ones aren't helpful. If I'm not able to get your test problem, could you test these cases individually? Or maybe @apcraig can?

[eclare108213]

@TillRasmussen I've downloaded the tarball and now understand your test problem a bit better. I think that to truly test the options that I've implemented, the runs will need to be done in the coupled system. The standalone ice test serves to check that the changes can get rid of existing residual ice, but I think that's probably because we're bludgeoning it (item 3 below). I'd like to know whether changing the lateral melting implementation provides a more physically based solution, and ice-ocean interactions probably will be critical for whether or not this works. I suggest doing 3 separate tests:

1. in ice_step_mod.F90, test the new floediameter option by itself. I doubt it will get rid of existing residual ice in 1 day, so you might need to run longer, or start from a point without existing residual ice and see if it appears. If that's not enough, then
2. in icepack_therm_vertical.F90, uncomment the code block that uses floediameter to see if lateral melting can ever be a solution to this problem.
3. test the change in icepack_mechred.F90 separately to make sure it works as expected.

Thanks for helping with this!

[TillRasmussen]

After running all 8 combinations of possible test the low ice concentrations are removed when using bullet 1 and 2 or all 3 bullets.

My first conclusion is therefore that using bullet 1 and 2 should do it.
!!!! NOTE CONCLUSION IS EDITED
I will rerun the coupled test for a year to see the difference and if the conclusion holds.

[eclare108213]

Great! So bullet 3 alone doesn't work? I was thinking that we could keep that one (fix it if necessary) to provide more flexibility for operational uses. It should be equivalent to @mhrib 's approach but conservative.

[TillRasmussen]

I have now run two simulations from the 1/9 1990 to the 1/7 1991. Conclusion is the same. Use the first two bullet points then the bug should be fixed. Output on the 1/ 1991 can be found on
ftp://ftp.dmi.dk/tarftp/ehunke/base_iceh_inst.1991-07-01-00000.nc
ftp://ftp.dmi.dk/tarftp/ehunke/cor_iceh_inst.1991-07-01-00000.nc
ftp://ftp.dmi.dk/tarftp/ehunke/diff.nc

base is the original code, cor is the run using the corrected code and diff.nc is the difference between the two.
Using ncview the issue can best be seen by changing the mode from linear to low

[eclare108213]

@TillRasmussen could you run another test? I found a bug in my alternative floediameter calculation, which will increase the amount of heat used for lateral melting by a factor of 4. I'd like to see if this is sufficient using just bullet 1, not 2 or 3.

[TillRasmussen]

I did run the updated code with only bullet 1. This do not remove the low concentration ice.

[eclare108213]

Okay, thank you. I need to think about this some more.

[eclare108213]

@TillRasmussen which thermodynamics option are you using? mushy?

[TillRasmussen]

Yes it is with the mushy thermodynamics.

[rgrumbine]

At long last, I'm working routinely on this. The issue is still present and discoverable in gx3 runs with the current ice_in and then option file with:

npt_unit = 'd'
npt = 365
dumpfreq = 'd'
dumpfreq_n = 365
histfreq = 'd','x','x','x','x'
f_aice = 'd'
f_hi = 'd'
f_uvel = 'd'
f_vvel = 'd'
f_Tsfc = 'd'

At the end of 2005, there are over 5 million km^2 of ice extent with concentration between 1e-2 and 1e-4

The initial conditions have much more ice extent than it finishes the year with, by about 4 million km^2 (and this is pretty independent of whether I use 1e-2, 1e-4, or 0.15 as the concentration criterion).

So I'm going to look at removing low concentration or thin ice from the initial conditions and rerun. Some perhaps large fraction of the issue is simply initial conditions. Can't be all of it as the model retains large areas of low concentration even after running for a year.

[rgrumbine]

Following up: I was using the 2017 gx3 initial conditions. Will rerun with the 2021 conditions.

[rgrumbine]

2017 was the gx3 v5 conditions.
2021 (starting from v6 for 2005-01-01) is identical

so back to seeing about editing the initial conditions. Alas, simply zeroing aice and hi doesn't work due to the other places affected by nonzero values for them.

[rgrumbine]

Finally a clean experimental narrowing of the issue. tl;dr -- much of the issue derives from the transport scheme, the rest seems to be failure to melt away the last little bit of ice.

• For initial conditions:
• • zero out all velocities
• • where either aice < 0.01 or ice thicknesses < 0.01, zero aice, vice, qiceN

This occurs for both 'upwind' and 'remap' advection. 'upwind' suffers more from the problem than 'remap', but only ca. 25%

advection.ts:

Test Grid PEs Sets BFB-compare

smoke gx3 1x1 med3,yr_out,nodyn
smoke gx3 1x1 med3,yr_out,nodyn,upwind
smoke gx3 1x1 med3,yr_out,nodyn,notrans
smoke gx3 1x1 med3,yr_out,nodyn,notrans,upwind
smoke gx3 1x1 med3,yr_out
smoke gx3 1x1 med3,yr_out,upwind
smoke gx3 1x1 med3,yr_out,notrans
smoke gx3 1x1 med3,yr_out,notrans,upwind

set_nml.notrans:
ktransport = -1

set_nml.nodyn
kdyn = -1

set_nml.upwind:
advection = 'upwind'

set_env.med3: (1 processor on gaea takes about 2.5 hours for a gx3 year)
setenv ICE_RUNLENGTH 3

set_nml.yr_out:
diagfreq = 1
npt_unit = 'd'
npt = 365
dumpfreq = 'd'
dumpfreq_n = 365
histfreq = 'd','x','x','x','x'
f_aice = 'd'
f_hi = 'd'
f_uvel = 'd'
f_vvel = 'd'
f_Tsfc = 'd'

When no dynamics are allowed (velocities remain zero throughout) results are not identical, which surprised me:
Date crit Area Extent Area Extent Area Extent
remap, transport: 2005-12-31 0.000 NH 15.760 15.921 SH 14.196 17.866 Global 29.956 33.787
upwind, transport: 2005-12-31 0.000 NH 15.760 15.921 SH 14.199 17.783 Global 29.959 33.704
remap, notrans: 2005-12-31 0.000 NH 15.760 15.921 SH 14.169 17.783 Global 29.929 33.704
upwind, notrans: 2005-12-31 0.000 NH 15.760 15.921 SH 14.169 17.783 Global 29.929 33.704

The deviations are largest in the summer hemisphere, as if ice melts to low concentration/thickness, but the last bits never go away. NH minimum extent about 8.6 million km^2, SH min about 12.5 million km^2 (in late April!)

With dynamics turned on, the disparities are larger:
remap, transport: 2005-12-31 0.000 NH 15.908 21.742 SH 8.445 28.326 Global 24.353 50.068
upwind, transport: 2005-12-31 0.000 NH 15.964 22.573 SH 6.236 35.243 Global 22.200 57.817
remap, notrans: 2005-12-31 0.000 NH 15.736 15.921 SH 9.281 17.403 Global 25.017 33.324
upwind, notrans: 2005-12-31 0.000 NH 15.736 15.921 SH 9.281 17.403 Global 25.017 33.324

Since dynamics, upwind, transport is the worst-behaving, I'll continue with that.

[eclare108213]

@rgrumbine did you test the approaches described in the comment above (#645 (comment))? I believe @TillRasmussen found that doing bullets 1 and 2 together took care of the residual ice issue. This must not have made it into the code, formally. Another change that has helped in E3SM (but not entirely removed the problem) is to reduce the parameter values for a_min and m_min to 1.e-11 and 1.e-10, respectively, something we should also test/fix in the Consortium codebase.

[rgrumbine]

Thanks for the pointer. While I've been scrupulous about keeping my codes up to the consortium code base, I haven't always kept them up to the discussions. I have a couple of things running, and will add these to my trials.

[proteanplanet]

I'm not sure if this is useful to the discussion, but I'm attaching this to share our experience in E3SM with regard to residual ice:

2024_Andrew_Roberts_CAMAS_Poster.pdf

[TillRasmussen]

Hi Andrew, Bob and Elizabeth I can add that at least in CICE the criterion for a_min is based on two different ice concentration. Aice for tmask and aice_init for umask. This can allow umask to be active without tmask being aktive. I think that is a bug but I am not sure. Till ons. 29. maj 2024 17.33 skrev Andrew Roberts ***@***.***>:

…

I'm not sure if this is useful to the discussion, but I'm attaching this to share our experience in E3SM with regard to residual ice: 2024_Andrew_Roberts_CAMAS_Poster.pdf <https://github.com/CICE-Consortium/CICE/files/15486715/2024_Andrew_Roberts_CAMAS_Poster.pdf> — Reply to this email directly, view it on GitHub <#645 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH7N7DVGU7Y2BN3N7VB5LUTZEXYOFAVCNFSM5GCEABQKU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJTG43TANRQGU4A> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[rgrumbine]

All:
I've checked, and my current (which should have all the floe diameter updates) runs still have the issues. With dynamics and upwind, the end of the year global extent in gx3 is over 50 million km^2

Following Andrew, I tried dropping a_min and m_min to 1.e-11 and 1.e-10, respectively, in ice_dyn_shared.F90. That made year end extent over 75 million km^2

Note: I'm using aice > 0 in determining extents. I don't really expect high quality match to observations, but do think even gx3 should be closer to year end (2005) observed extents than this.

[eclare108213]

Further thoughts:

• The bug that Till mentions above could be a factor -- we haven't investigated its effects in detail yet.
• Till's and and my tests relied on thermodynamics to remove residual ice via lateral melting, by altering the floe size diameter in the standard (not FSD) configuration. Bob, which thermo are you using? Is FSD turned off? Are you using the code mods from my branch? (note this branch is quite out of date by now with the main code base)

https://github.com/eclare108213/cice -b residual
https://github.com/eclare108213/icepack -b iresidual

• I wonder if there is some sort of coupling issue in your runs, which is reducing ice melt. E.g. are the ocean-ice heat fluxes being sent to/from the ice acting the way you'd expect?

[rgrumbine]

• I'm running the bog standard stand-alone CICE. All parameters as in the standard ice_in.
• Thermo is 'bubbly'
• Definitely not using your branches (I checked yours out and compared)
• Stand-alone CICE, so no coupled fluxes

Does the model rely on ocean fluxes to remove the trace amounts of ice?

One test I've completed is to confirm that state_to_work is cleanly inverted by work_to_state. Not a surprise, but good to know that the transport-related issue is strictly with the upwind (and remap) routines.

[eclare108213]

The model should rely on ocean fluxes to remove the trace amounts of ice. In stand-alone mode these should come from a forcing file. Sounds like you aren't testing with the changes from my 'residual' and 'iresidual' branches, so you wouldn't see the same results that Till did. We'll need to pull the commits onto that branch into a new one. When melting, the ice takes the heat that it needs from the ocean, and the changes in my branches allowed it to take more.