Fixes ice area inconsistency under convergence in single column mode

[davidclemenssewall]

For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium, please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers

PR checklist

• [x] Short (1 sentence) summary of your PR: Changed Icepack such that when there is a closing forcing, ice with uniform properties is advected into the cell (so that closing does not create open water).
• [x] Developer(s): David Clemens-Sewall
• [x] Suggest PR reviewers from list in the column to the right.
• Tony Craig, Elizabeth Hunke, Dave Bailey
• [x] Please copy the PR test results link or provide a summary of testing completed below. Icepack base_suite passed 138/138 (including new checks toggling 'ice_advection_type') in ubuntu docker container. CICE quick_suite passed 16/16 on cheyenne.
• How much do the PR code changes differ from the unmodified code?
  • [ ] bit for bit
  • [ ] different at roundoff level
  • [x] more substantial
• Does this PR create or have dependencies on CICE or any other models?
  • [ ] Yes
  • [x] No
• Does this PR add any new test cases?
  • [ ] Yes
  • [x] 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/.)
  • [x] Yes
  • [ ] No, does the documentation need to be updated at a later time?
    • [ ] Yes
    • [ ] No
• [x] Please provide any additional information or relevant details below:

Problem: In the consortium main standalone configuration of Icepack, we read in a mechanical forcing file that contains and opening and a closing rate at each time step. When the forcing is net closing, the ridging scheme in Icepack produces open water (which is not physical) as a consequence of maintaining a constant cell area. In reality, and in CICE, this does not happen because it is ice advection from the surrounding area/grid cells that is driving the closing.

Proposed solution: I've resolved this erroneous open water creation by adding an ice advection step prior to the ridging step (in the same way that CICE calls step_dyn_horiz before calling step_dyn_ridge). I assume that our single column grid cell is surrounded by ice of identical properties, so this amounts to scaling up the ice area, volume, and tracers before ridging. I have set the area of ice that is advected in to be equal to the closing rate times the timestep. I'd be happy to hear thoughts if others think a different area of ice should be advected in.

I've added a namelist parameter 'ice_advc_type' to icepack_in where 'uniform' means apply this advection step and 'none' means do nothing (a.k.a. the behavior without this fix). I've set the default to be 'uniform' because I think that most people running a single column model would not expect closing to produce open water.

When running standalone cases with the SHEBA ice dynamics forcing, including the advection of ice into the cell on closing dramatically increases the ice thickness and aice. This is to be expected since we are adding ice to the cell on closing. Note, ice is always advected out of the cell on opening (both in the current version and this modified version).

[eclare108213]

Also, please, would you fix the formatting in the template above? Checkboxes need to be

- [x]

instead of

• [X ]

Thank you!

[apcraig]

Had a quick look, we need to be thoughtful about adding new public interfaces to Icepack's columnphysics. If icepack_step_advection_scm can be migrated out to the driver, then I support that effort. I assume that will require some additional queries for columnphysics "data", and is why it was put inside the columnphysics in the first place. Hopefully, all the data that is needed is accessible to the driver. Let me know if there are any questions.

[davidclemenssewall]

@eclare108213 @apcraig @dabail10, this PR is ready for review.

I have inlined the changes in step_advection_scm and removed all modifications to the columnphysics. Because tracers are stored a densities it no longer modify the tracer values (essentially the scaling I was doing to the tracer values was just undone by the compute_tracers call). I've also made the requested documentation, naming and testing changes. The current commit passes the icepack base_suite and CICE quick_suite.

[davidclemenssewall]

Sorry for the confusion. I had thought that if 'uniform_ice' is to be the default it would make sense for all tests that include a closing forcing to use 'uniform_ice'. If we don't want this PR to change all tests that use a closing forcing, then I think our two options are:

1) don't have 'uniform_ice' be the default (i.e., have the default be 'none' so a user has to specifically choose to advect in ice on closing). Replace set_nml.advcopenw with namelist setting option file to set ice_advection_type to 'uniform_ice' and add a few tests with that set.

2) leave 'uniform_ice' as the default and add the advcopenw setting option to each of the current tests that include a closing forcing (i.e., all of them except the bgc runs)

Of the two options, I think from a code testing perspective I'm more inclined to option 1. But I'm happy to set it to whatever you think is best.

[eclare108213]

Option 1 is definitely simpler. These tests are mainly to exercise the code -- we don't necessarily expect them to have physically realistic outcomes, although we try to keep them from being grossly unphysical. But I see your point, that the SHEBA forcing (and therefore closing) is turned on by default in icepack_in. It's okay with me to make these consistent, i.e. set 'ice_uniform' in icepack_in (as you have it) and 'none' for some of the tests for code coverage. What I would ask for in that case is a full comparison of the base_suite runs, before and after the change, for sanity's sake and as documentation, since the results are not BFB. By "full comparison" I mean a set of timeseries plots (made with Icepack's timeseries script is fine) for each case, at least for the ice state (area, thickness) and ridging variables, and others if it's easy/convenient.
@apcraig @dabail10 Do you have a strong opinion either way?

[davidclemenssewall]

I wonder if we should stick with option 1 until we have a mutually consistent forcing dataset.

I've done a little bit of ad hoc analysis of the output and it definitely changes the ice state. However, since the default forcing doesn't correspond to a set of conditions that actually occurred anywhere, it's hard to evaluate whether the changes in the ice state are what we should expect. For example, with the defaults and no advection cells 2 and 3 (the two that start with ice in them in January) simulate more open water than SHEBA (peak around 30%, vs. 17% observed) and the timing is offset (SHEBA peak open water was in mid-August). Including ice advection the simulated maximum open water fraction is just 5% (for cells 2 and 3) but the timing is more consistent with the SHEBA observations. However, the ice thicknesses in the no advection case are more reasonable, and appear to be stable. In contrast, in the advection case the ice thickness is not stable but increases without limit if we run the simulation out several years. Which is not too surprising since we don't have any ice strength feedback (i.e., more ice is pushed into the grid cell and ridging happens even when the ice gets thick enough that it should resist further ridging). But I think this behavior might be unexpected for testing.

With the MOSAiC test case (and hopefully other future test cases with self-consistent atmosphere, dynamics, and ocean forcing) I think that it would make sense to turn on ice advection on closing as the default.

aice, vice, and vsno without advection (default settings):

aice, vice, and vsno with advection (default settings otherwise):

[eclare108213]

Thank you, this is very helpful. I agree that sticking with option 1 makes the most sense for now, considering the ice thickness behavior. In a simulation with horizontal extent, ice would be advected both in and out, while in this PR you're only adding ice to the grid cell when it's converging. Does it make sense to also remove some of it when diverging? That might be getting pretty complicated (like the real advection scheme), but it's worth thinking about -- seems like this will be a persistent issue for MOSAiC simulations in a single column.

[davidclemenssewall]

Hmm, I had thought that Icepack was advecting ice out when we had divergence. We certainly get a drop in aice when there is an opening forcing. Now that I go back and look at the code though, I'm not totally certain how that is happening. Line 1289 in icepack_mechred (subroutine ridge_shift) increases aice0 by opning*dt. If there were no closing forcing, I think that this would be the only change in the ice thickness distribution at this step. However, then when asum_ridging is called (line 410) asum should be greater than 1. Which would cause iterate_ridging to be true and closing_net should then be equal to the opening rate from the first iteration. On the next iteration of ridge_shift I would think this would just close the open water that was opened on the first step. So I'm not sure how it is that Icepack is producing open water at all.

Maybe in step_advection_scm we should set expansion_ratio = c1 + (closing(i) - opening(i)) * dt? That would at least be consistent with the assumptions that ridge_ice currently makes in line 299 (in icepack_mechred).

[eclare108213]

The code distinguishes between advection, which moves ice between grid cells, and open water closing/opening within a grid cell via ridging. Advection is necessarily a grid-dependent process and so is not included in Icepack's column physics, but the rates of divergence/convergence that result from the advection are used in the column physics to ridge ice in a single grid cell. And as ice ridges, it can cause the area of open water to increase -- that's what it should do, in the absence of other ice being advected in. The expert on this part of the sea ice model is there at NCAR: @whlipscomb. Maybe the two of you can put your heads together and propose a solution for Icepack with the MOSAiC data. This is definitely one of the challenges of working with Icepack as a Lagrangian model! Thanks for digging in.

[davidclemenssewall]

@whlipscomb and I had a good discussion about this yesterday.

In short, our conclusion was that for the mechanical redistribution in single-column Icepack to be consistent with the mechanical redistribution of Icepack inside CICE, there needs to be an "advection" step that changes the total grid cell area by an amount consistent with the net opening and closing.

Without this step, Icepack does not handle closing or opening consistently with how CICE handles it. We've discussed the issues with closing at some length, but the opening issues are a new realization. In effect, what is happening in the current default configuration of Icepack is that when we provide it an opening forcing, it accomplishes the required amount of opening by ridging ice within the grid cell, not by advecting ice out of the grid cell. In other words, opening and closing have exactly the same effect in the current default version of Icepack. You can see in the following figure where we have opening followed by a similar magnitude of closing. Both the opening and closing decrease aice by a similar amount and neither has any impact on the trajectory of vice.

By adding including an advection step before this (where we have expansion_ratio = c1 + (closing(i) - opening(i)) * dt) we get more intuitive results where opening decreases the ice area and volume and closing increases both:

I think that the gradual decline in ice area with opening followed by an abrupt rise is due to the mixed layer response. The mixed layer starts a little above the freezing temperature and it takes a few hours of opening before enough heat is extracted to drop sst to the freezing temperature. Once it cools enough, frazil ice is rapidly added reducing the open water area.

With applying advection for opening and closing (the current change, assuming the grid cell is surrounded by uniform ice), the system appears stable and produces reasonable snow and ice thicknesses (see 5 year run below). I don't think that advection necessarily produces a generally stable system but perhaps this case is okay because the total closing and opening at SHEBA were very similar in magnitude. Additionally, the amount and timing of open water is now much more consistent with SHEBA (where the dynamic forcing comes) than the default case. In the simulation the maximum open water fraction 15% on August 5. At SHEBA the max was 18% on August 7. SHEBA flights only measured open water every 2 weeks. The atmospheric and oceanic forcing is different anyways so I think this is a helpful gut check but nothing more.

Because the current version does not handle opening consistently with CICE, I think we need to change the default and do a full comparison of all tests. I am happy to have the default either be assuming the simulated grid cell is surrounded by uniform ice of the same properties or by open water (but now actually advecting in the open water or uniform ice and advecting out ice on opening to be consistent with CICE). I have a slight preference for the default to be uniform ice but I'm happy either way. @eclare108213 @apcraig @dabail10 let me know what you think and I will implement it and do the full test comparison.

[eclare108213]

Nice, I like what you've done here. My vote would be for assuming that the grid cell is surrounded by ice that is the same as in the central cell, so that you aren't imposing some sort of spatial dependence on your results, and then have a test or two for the 'none' option (as you do now), for code coverage purposes.

[davidclemenssewall]

Apologies for my delay here, we had some family stuff come up.

I have finished the changes and compared the output from all of the tests before and after this fix. None of the state variables in the test output changes dramatically or unreasonably after the fix (e.g., ice thickness changes by 10s of cm). So I don't think it is likely to disrupt other testing.

In general, the largest impact is on aice, as discussed above. There are also downstream impacts. Specifically, because the change reduces the open water fraction in the summer, there is less solar heating of the mixed layer and less basal melt. Also, having higher ice area in September increases the snow depth as less snowfall is lost to the ocean.

Here are plots of output for all tests without the fix: plots_main_20230421.zip

And here with the fix: plots_convergence_open_water_6.zip

[davidclemenssewall]

I'm unsure why Read the Docs is now failing. I don't think that I have changed anything related to that. Otherwise I think this should be all set to review. It passes 138/138 on the base_suite.

[apcraig]

Seems to be an issue with readthedocs. Can you make a trivial change to your branch and push the commit. Lets trigger the checks again and see if the problem persists.

[apcraig]

OK, failed again. I then forced a rebuild and it worked. I think it's readthedocs at the moment. At any rate, there is an updated readthedocs for this PR, so I think we're OK for the moment. Let me know if you notice anything funny again. Thanks.

[davidclemenssewall]

Just checking in here, is there anything more that you'd like me to do on this PR? Otherwise, I think it is ready.

[apcraig]

Just to be clear, are all results in Icepack and CICE expected to change? Has anyone done a QC test in CICE? I can run a full CICE test suite on cheyenne and do a QC test if it's needed/helpful. Just let me know.

Quick review on my part shows no problems. There is an empty line with just spaces added in icepack_mechred, but I can fix that later. I have an automated tool that I should run on the code again. It's also likely readthedocs is broken, but not due to this PR. That's another thing we'll need to fix asap. For now, I think once we feel testing is complete and @eclare108213 gives approval, we can probably merge.

[apcraig]

Looked at the code again, this does not change anything in the Icepack columnphysics, so cannot affect CICE. Is that what we want/expect?

[eclare108213]

Looked at the code again, this does not change anything in the Icepack columnphysics, so cannot affect CICE. Is that what we want/expect?

Yes, this only changes how Icepack is forced, so shouldn't affect CICE at all.