Zap residual ice

[eclare108213]

• Short (1 sentence) summary of your PR:
  Adds minimum area and mass parameters and zaps grid-cell volumes of ice less than those.
• Developer(s):
  @eclare108213, @DeniseWorthen
• Suggest PR reviewers from list in the column to the right.
• Please copy the PR test results link or provide a summary of testing completed below.
  Significant testing is documented in the comments below, including Icepack test suites, QC tests in CICE and ice-ocean coupled simulations. The results are not BFB but do change the results in regions where there was little ice to begin with. QC tests indicate these are not climate-changing differences.
• How much do the PR code changes differ from the unmodified code?
  • bit for bit
  • different at roundoff level
  • more substantial
• Does this PR create or have dependencies on CICE or any other models?
  • Yes
  • No
• Does this PR add any new test cases?
  • Yes
  • No
• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please document the changes in detail, including why the changes are made. This will become part of the PR commit log.
  This PR addresses the long-standing residual ice problem in CICE by zapping grid-cell volumes with area or mass less than the new namelist parameters min_area and min_mass. To work properly, these parameters should likely be set to the same values as the dynamics area and mass minimum values dyn_area_min and dyn_area_mass. Since those parameters in CICE are quite large (and everyone should decrease them), the default values in Icepack are currently set to 1.e-11 and 1.e-10, respectively, maintaining the area/mass=0.1 relationship of the dynamics parameters. The driver (e.g. CICE) should initialize them as needed via an icepack_init_parameters call.

[eclare108213]

An accompanying CICE PR (CICE-Consortium/CICE#1055) initializes the parameters to the dynamics values. In my brief CICE tests so far, the Icepack default values produce BFB results but I expect they will not be BFB in longer tests. To recover perfectly BFB results, set them both to zero to prevent extra zapping. I have intentionally not used zeroes as the default values even though that might not be BFB -- I'm happy to hear arguments to the contrary and change it, if you think it would be better.

I have not completed full standalone test suites in either model. While I'm doing that, I hope that others of you will test this in your coupled systems to make sure that it really does work as advertised, and perhaps to test the sensitivity of the parameter values.

[DeniseWorthen]

I've got your residual zap branch up and running now and I'm going to do some longer comparison runs w/ our current develop branch and the ICs that start w/ residual ice.

For sensitivity, are there suggested 'non-large' values of dyn_area_min and dyn_area_mass that I should try?

[eclare108213]

For sensitivity, are there suggested 'non-large' values of dyn_area_min and dyn_area_mass that I should try?

I'd start with the two end points and one in the middle:

min_area = 1.e-11, min_mass = 1.e-10
min_area = 1.e-7, min_mass = 1.e-6
min_area = 1.e-3, min_mass = 1.e-2

[DeniseWorthen]

I've completed 4 runs with our 1/4deg configuration, which starts with large counts of aice_h<1.0e-6. In the figure below, the 'base' case is our current develop branch (c08a7a6) and 3 runs with the zap-residual modifications. I'm plotting using instantaneous history files every 6 hours.

I know it looks suspicious that the case w/ the smallest min parameters starts higher than the base case, but I verified that they have identical counts in the IC that is written at start up (iceh_ic.2021-05-01-00000.nc) but the first 6 hours results in ~7000 additional "residual ice" points which quickly get zapped away.

[eclare108213]

I've completed 4 runs with our 1/4deg configuration

Yay! This mostly makes sense to me. After it settles down, the green line indicates about how much ice you have with concentrations between 1.e-7 and 1.e-6. The larger criteria used for the orange line zaps everything less than 1.e-6 and then some.

By changing dyn_area_min and dyn_mass_min, in addition to zapping you're also changing where the dynamics calculation occurs for values less than those in your base case, so it's not surprising that the red line starts off differently. At least for these few months, it looks like your configuration doesn't crash with the puny ice area (red), but it's interesting that it grows after the initial drop, while the blue baseline stabilizes, like there is something weird going on in the red case. So maybe that's too small. It would be informative, I think, to see what the velocity values are around the ice edge. They might be quite big for the smallest criteria.

I think you should choose the parameters you use based on validation against whatever observational data you use. There might actually be 1.e-7 bits of ice in some grid cell at some time, but it shouldn't matter much for the simulation and I would guess that you can use larger cutoff values to better match your statistics.

Does this all make sense to you?

[NickSzapiro-NOAA]

There are some more places with "minimum" thresholds like hi_min=p01 for icepack_therm_itd, aice>p001 in neutral_drag_coeffs in icepack_atmo, dh_min=p001 in icepack_brine. Are these independent of the parameters here?

[DeniseWorthen]

I agree that it is satisfying that all three parameter sets, nothing crashes and the impacts seem rational and consistent.

For the red-line case, yes, the residual ice trends upward, but in way which is more physically satisfying (imho). This is what the residual ice looks like in the IC--large patches which have formed but have no way to be removed

and at the end of red-line run

Residual ice now tracks the actual ice edge; the parameters adjust how fast/hard you want to tighten that gradient.

[eclare108213]

@NickSzapiro-NOAA, good catch. We have tried to pull out all the hardwired parameters but clearly missed some of them.

hi_min is a namelist parameter that is intended to maintain stability in the vertical thermodynamics for thin ice (so mainly the thinnest category). Given an area, it would be related to the mass minimum (which is the total over all categories), but functionally they are pretty independent. Conceptually, if the ice thickness in a category is less then hi_min, then the vertical thermo calculation isn't done, but the ice might still be transported and ridged and end up greater than both hi_min and the zapping minima, allowing the thermodynamics to act on it in the next time step. This is the way the code was originally intended to work, and why we spent so much time trying to include hi_min in the zapping criteria.

The p001 criterion for aice in ice_atmo.F90 should not be hardwired. I think this could track dyn_area_min for the wind stress -- it's only used for the form drag parameterization so shouldn't affect configurations without form drag. This change will need to be tested by someone actually using form drag.

@njeffery, the brine height parameterization uses both dh_min=0.001 and thinS=0.05. Could/should those be related to hi_min or the dynamics minima, or are they truly independent?

[njeffery]

@eclare108213 : They should be larger than or equal to the thermodynamic minimum thickness, but not smaller.

[eclare108213]

@eclare108213 : They should be larger than or equal to the thermodynamic minimum thickness, but not smaller.

but hi_min=0.01 and dh_min=0.001?

[njeffery]

Then it shouldn't be making a difference. That's not a problem though. It just means we always calculate brine height changes rather than assuming its the same as ice thickness.

[apcraig]

Let me know if you want me to run a full test suite on Derecho with the updated implementation in Icepack and CICE, when ready, before we merge.

[eclare108213]

Let me know if you want me to run a full test suite on Derecho with the updated implementation in Icepack and CICE, when ready, before we merge.

That would be awesome, thank you. Since this is not BFB, a QC test will also be needed. But I think there also needs to be broader testing by the various modeling centers before this is merged.

AND I would really like to get opinions on the parameter values. Should we recommend reducing dyn_area_min and dyn_mass_min in the CICE PR? And if so, to what? I'd rather not use the large values we currently have for zapping, but that is how the code is currently set up in that PR.

Thinking out loud... maybe 2 QC tests could help decide, one with the parameters as currently set in the CICE PR, and a second one with parameters set at the other extreme, as they are currently in Icepack. Compare both test simulations with the present code and with each other. I wouldn't expect either to be climate changing in that test, but who knows? (Caveat: I'm not sure whether the QC configuration has residual ice.)

[apcraig]

OK, happy to help with testing. I can run the test suite when we have final parameters chosen. Do you want me to run the QC tests with the current PR. What values for the parameters should I use for the two tests. I will compare them to each other and to the baseline. Also understand there may be no zapping, so we'll see how it goes.

[eclare108213]

Running the QC tests with the current PR will capture the extremes of potential parameters, I think. The two tests would be

A
dyn_area_min = 1.e-3
dyn_mass_min = 1.e-2

and

B
dyn_area_min = 1.e-11
dyn_mass_min = 1.e-10

The baseline run will have the dynamics parameters as in A but no extra zapping.
Thanks!

[apcraig]

The QC tests all pass. The difference plots don't show much of anything. There are some differences in the Canadian archipelago and something in the southern ocean, but nothing much at the ice edge. I attach the three difference plots. zapbase is the current code. zapqc1 is A above with 0.001 and 0.01 parameter values. zapqc2 is B above with 1e-11 and 1e-10 values.

[eclare108213]

The QC tests all pass. The difference plots don't show much of anything.

That's excellent news! So at least in this configuration, the parameter cutoff values don't matter. Do the QC test diagnostics consider all the output, or do they cut it off at e.g. aice=0.15? Searching the QC directory for 0.15 didn't turn up anything for me, but a plotting or analysis cutoff could have a different value.

[apcraig]

I do not see any cutoff values in the QC test. The test reads only 'hi' as far as I can tell. There is a mask generated for gridcells that never have hi > 0.01 but that is not an aice criteria. I think the QC test considers all gridcells. It could be that zapping plays little role. What's interesting is that baseline - A has the most differences. This suggests that the baseline is operating a lot like B with, I assume, less zapping happening in both vs B. My conclusion is that the new implementation is zapping more than the current implementation. I think that's what we expected. Whether it's doing that mainly at the ice edge seems like an open question though.

[eclare108213]

There is a mask generated for gridcells that never have hi > 0.01 but that is not an aice criteria.

In cice.t-test.py:

    Calculate the difference and mask the points where the the difference at
    every timestep is 0, or the sea ice thickness for every timestep is < 0.01 meters
    for data_a or data_b

The hi > 0.01 criterion is probably masking out most small-area cells. Why this is done here is understandable, to prevent comparing cells that have a little bit of ice in one run with cells that never have any in the other run, which probably shouldn't matter for climate. But we are looking for exactly those cells! What happens if the criterion value is reduced? Could it be set to min(dyn_mass_min/aice) over the two runs assuming aice=1, or 1.e-10 here? Maybe QC isn't the best test... but it is good to know that there is a hardwired parameter in there.

EDIT: On the other hand, the criterion is for every timestep, and our small values may be moving around. So maybe some of the cells are being tested? I would guess that the "not climate-changing" tests will pass, regardless.

[apcraig]

Right. It's a cutoff for "every timestep". That means the test is probably assessing most of the ice edge for most of the year and would only exclude gridcells that only have very small amounts of ice probably just at the peak of the ice season at the far edge.

[anton-seaice]

Apologies if this turns into a hijack - I can run some tests, but I am not sure what to look for.

In a recent run with the default settings for minimum area, I can't identify a problem (based on the criteria I think is being used above).

This is cells with 0<aice<1e-6 in the last day of the run, the low concentration cells all look to be in reasonable locations ?

[eclare108213]

@anton-seaice thanks for running a test! Is this a one-year run or longer? Resolution? It looks like you have some larger areas of low-concentration ice developing. They are in reasonable locations; the problem is that they persist through the melting season and shouldn't. Do you see that behavior?

[anton-seaice]

@anton-seaice thanks for running a test! Is this a one-year run or longer? Resolution?

This is end of 52 years, Its 0.25deg global - data atmosphere, interactive ocean.

It looks like you have some larger areas of low-concentration ice developing. They are in reasonable locations; the problem is that they persist through the melting season and shouldn't. Do you see that behavior?

Ok thanks - I did a plot of the test above AND ice is present for the entire final year. (So ice is present in all daily diagnostics for one year, and is less than 1e-6 in the last day):

The same figure for ice present for several years is similar. So that looks like the same problem.

So I think I could do two tests with the build from CICE-Consortium/CICE#1055 with:

• default dyn_area_min and dyn_mass_min,
• very small dyn_area_min and dyn_mass_min

and re-run the last 2/4/10 years of the experiment?

[eclare108213]

@anton-seaice That would be fantastic! Is this a coupled configuration?

[anton-seaice]

Coupled with an ocean, and using a data-atmosphere forcing

[anton-seaice]

As far as I can tell this works ok.

I ran two experiments with a 25km global ocean + sea ice model forced with a repeat year forcing of a selected year from JRA55do, starting from a restart about ~40 years into an control experiment

• `control' - without this change
• small_min : ( dyn_area_min = 1.e-11, dyn_mass_min = 1.e-10 )
• large_min : (dyn_area_min = 1.e-3, dyn_mass_min = 1.e-2)

This change appears to zap ice sucessfully:

Although some results are counter-intuitive, e.g. the control (without the extra zapping) appears to have more ice in the 1.e-6 to 1e-3 concentration range than the change with the zapping using small_min. (Does this mean the extra zapping is preventing some ice growth ?)

The impact on Sea Ice Area is small (the change in maximum area is similar, not shown):

It appears to deal with residual ice ok:

Figure is cells where ice exists for all of the final year AND final concentration is <1e-6. There is a weird blip near Coronation Island off the end of the Antarctic Peninsula in the small_min run.

Global ocean temperature is marginally cooler in the large_min run, and marginally warmer in the small_min run:

[anton-seaice]

How about aborting if min_area or min_mass is less than `puny? Those setings mean the user has tried to do something which the model doesn't support.

e.g. could implementing as max(min_area, puny) leave people confused about why there is no ice with concentration less than puny, greater than min_area

[eclare108213]

I guess someone might try to set min_area and min_mass both equal to zero in an attempt to turn off zap_residual, which wouldn't work. We can certainly implement a flag to turn it off, as @NickSzapiro-NOAA suggests. Then under the covers we could do min_area = min(min_area, puny) at initialization to enforce the puny minimum, which I do think we should keep, although I'd rather not obscure what's happening. What do you suggest, @anton-seaice ?

[anton-seaice]

Yes the namelist flag to enable/disable zap_residual is ok - (I think it should be enabled by default though, but we could revisit the default in a year or two)

Does setting min_area < puny means no ice would be zapped ? (as thermodynamics wouldn't be active?)

I think abort if the user tries to set it less than puny, we could write in the error message to use the "zap_residual" flag to disable zapping if needed ?

[eclare108213]

Setting min_area < puny means that ice areas less than puny would still be zapped, which (if I remember correctly) is the current behavior. In that case, we could write a warning that areas less than puny will still be zapped and offer the namelist flag as an alternative to turn off zapping completely, with no need to abort. I agree that having zap_residual on as the default is what we want.

[NickSzapiro-NOAA]

To be able to keep current answers and test change, can this come in via a "zap_residual" namelist option?

[eclare108213]

To be able to keep current answers and test change, can this come in via a "zap_residual" namelist option?

Yes, good thought.

[eclare108213]

As far as I can tell this works ok.

@anton-seaice thank you for running this test. It looks like it's working as expected and doesn't do weird things to the ocean, which is what I was concerned about. But you're right, the results are a bit counterintuitive.

Does this mean the extra zapping is preventing some ice growth ?

I think ice in the runs with zapping have a higher bar to reach before they can survive zapping, whether the low-concentration ice arrives in a cell through growth/melt or advection. In that sense, the extra zapping could be preventing some ice growth. But the difference between the control and the low-min runs is both zapping on/off and the cutoff values for ice velocity (and therefore advection) to be calculated, so it's complicated. In your plot with areas between 1.e-6 and 1.e-3, what time of year are the main peaks occurring? (Are those plots SH?) There is a notable adjustment in the first month or two of the low-min test, which I'm not sure how to interpret.

What's also counterintuitive to me is that the annual minimum area in the control and small-min runs are so similar. I guess that could be driven more by the ocean than internal sea ice physics, and the large-min run simply zaps way too much ice.

[anton-seaice]

Thanks @eclare108213

• the extra zapping does cool the ocean slightly. I guess it accelerates the transfer of heat from ocean -> ice -> atmosphere by extracting the heat from the ocean to do the zapping

In your plot with areas between 1.e-6 and 1.e-3, what time of year are the main peaks occurring? (Are those plots SH?)

This is global, peaks are march / september - one for each hemisphere

There is a notable adjustment in the first month or two of the low-min test, which I'm not sure how to interpret.

Maybe this (extra low conc cells) is related to the extra points / lower limits for advection ?

[apcraig]

Testing in Icepack and CICE indicates this is ready to merge.

[dabail10]

Just an update. I finally did a sensitivity run in CESM with these changes. I think it might have more of an impact.

https://webext.cgd.ucar.edu/BLT1850/b.e30_alpha08b.B1850C_LTso.ne30_t232_wgx3.299_cice/ice/html/infrastructure/index.html

[anton-seaice]

Interesting - it might be neatest to make a new issue ?

So the plots show an increase in sea ice area in both hemispheres in a coupled (with interactive atmosphere) model with zap residual ice enabled?

[dabail10]

Interesting - it might be neatest to make a new issue ?

So the plots show an increase in sea ice area in both hemispheres in a coupled (with interactive atmosphere) model with zap residual ice enabled?

I think this might just be interannual variability. Also, if you look at the spatial plots, the differences are mainly around the ice edge.

[NickSzapiro-NOAA]

should this be vice*rhoi (instead of aice*rhoi) @eclare108213 ?