Spurious velocities in MOM6 ice shelf cavities

[cyung]

Hi MOM6 folk, I've got some idealized MOM6 ice shelf cavities running, based on the MOM6-examples ISOMIP examples. With thermodynamics turned on it works quite nicely and gives nice melt patterns and circulations. However, when thermodynamics is turned off I still see quite large spurious velocities emerge in multiple types of vertical coordinates.

We've tried turning KV_ML_INVZ2 = 0, which was an issue that got resolved last year, but unfortunately, we still have quite large (order cm/s) velocities emerging in both the ISOMIP+ case and other configurations I have set up. See this repo for configs and this notebook for some plots. The following plot shows the maximum horizontal speed in the unforced ice shelf cavity, where we expect velocities created from the imperfect ice shelf initialization to decay, and instead, they grow, particularly for the models run in sigmazstar and sigma (terrain or ice-shelf draft following) coordinates.



Does anyone have any ideas about how to reduce these spurious velocities or what could be going wrong (either in this config or in the model)?

 

[MattHarrison]

Hi Claire,

Suggest the following parameters:

FIXED_DEPTH_LOTW_ML = True
LOTW_VISCOUS_ML_FLOOR = True
HMIXED_FIXED = 20.0
KV_ML_INVZ2 = 1.e-4

When the model is in adiabatic mode with the bulk mixed layer disabled, unrealistic viscous coupling across massless layers in contact with the ice shelf was the source of problems reconciling large barotropic fluxes after initialization with layers which were being locked in place by the erroneous viscous coupling.

 

[cyung]

Hi Matt, thank you very much for your advice and explanation!

Unfortunately, the changes you suggested (turning LOTW_VISCOUS_ML_FLOOR on and making KV_ML_INVZ2 nonzero) only made very small changes to the velocities. See below a plot of the maximum horizontal speed over time in the zstar configuration with different KV_ML_INVZ2 values, which do not change much even with a very large KV_ML_INVZ2. The other coordinates (sigma_shelf_zstar and sigma) have similarly small changes to their order cm/s velocities.




What vertical coordinate and vertical resolution did you expect this change to work best? Is there anything else that might be wrong in my MOM_input (taken from the MOM6-examples ISOMIP layer example) and MOM_override?

Thanks!

 

[gmarques]

Do you get reasonable spurious velocities when you remove the ice shelf while keeping everything else the same? This is equivalent to the seamount test but using the ISOMIP topography. Before adding the ice shelf, even with the thermodynamics turned off, it would be good to make sure there isn't anything else in the configuration that might result in spurious velocities.

 

[cyung]

Thanks Gustavo for your suggestion! Indeed, your suspicion was correct and the velocities still grow to be quite large when I remove the ice shelf (by setting ICE_SHELF=False and TRIM_IC_FOR_P_SURF = False). However, the velocities are zero when I remove the topography and it is a rectangular box.

I have tested a few different parameters and find that the parameter that affects the magnitude of velocities the most, with and without an ice shelf, is the vertical stratification. Attached are some plots showing the different maximum horizontal velocities after 30 days as a function of the initial stratification (set by either a temperature T or salinity S gradient), in the three different coordinates and in linear and log plots to show the scaling. Here, I have set KV_ML_INVZ2 = 1.0e-4 as Matt suggested. Dots show the results without the shelf (but still with the ISOMIP+ topography) and crosses are with the shelf included.

As you can see, the velocities increase with density difference. Interestingly the relationship is approximately linear with the shelf present, especially for the zstar case. If I use a small enough stratification then the velocities are probably acceptably small (10^(-4) or 10^(-3)m/s), but they are not small for the stratification of the Asay-Davis (2016) ISOMIP+ configurations (the default, dark blue crosses are the ISOMIP+ Ocean1 experiment initial conditions), and especially not for the sigma_shelf_zstar and sigma coordinates with an ice shelf.

With the knowledge of this (approximately) linear scaling of velocity with stratification, does anyone have any ideas of what could be going wrong?

Are there any parameter settings I could use to improve the situation?

 

[gmarques]

Hi Claire,

The seamount tests should serve as reference for your "noshelf" tests. You might want to run the seamount tests to verify what's the "acceptable" maximum speed for each vertical coordinate in this kind of test.

You can also try to use the initial T&S and thicknesses from the seamount test case in your "noshelf" test. Does that help?
THICKNESS_CONFIG = "seamount"
TS_CONFIG = "seamount"

Another test is to use the topography from the seamount test and keep the T&S and thicknesses from the ISOMIP.

 

[cyung]

Thank you Gustavo for your suggestions. The seamount tests have very low velocities (10^-12 m/s), even when I modify the topography to be the same as the ISOMIP topography, and add the ISOMIP layer and TS configurations.

By comparing with the seamount test, I have modified the ISOMIP no shelf zstar configuration to have minimal velocities. The minimal parameters required to “fix” the configuration (starting from the MOM6-examples/layer repository) are different for zstar and sigma coordinates, suggesting there are multiple things going on, though I am yet to find the key problem parameters in the sigma case. The MOM_override for zstar coordinates is:

#override TS_CONFIG = "ISOMIP"
#override USE_REGRIDDING = True
#override REGRIDDING_COORDINATE_MODE = "ZSTAR" ! default = "LAYER", alternative "SIGMA_SHELF_ZSTAR"

#override ICE_SHELF = False ! Turn shelf off
#override TRIM_IC_FOR_P_SURF = False ! Also required once shelf off

#override DRAG_BG_VEL = 0.05 ! Default is 0.00
#override KH = 6.0 ! ISOMIP Asay-Davis protocol, original is 0
#override MIN_THICKNESS = 0.001 ! Original is 1E-12
#override MASS_WEIGHT_IN_PRESSURE_GRADIENT = False ! Original = True

This gives me the red dashed line the below maxmimum speed plot, which has a similar order of magnitude of maximum speed to the seamount test case. However, when I turn the shelf back on with these parameters (returning to ICE_SHELF = True and TRIM_IC_FOR_P_SURF = True), the zstar and sigma_shelf_zstar cases have very large velocities (purple and brown dashed lines). If I return to MIN_THICKNESS =1E-12, the ice shelf configurations improve a lot (orange and green lines for zstar and sigma_shelf_zstar), but then the case with no shelf (blue line) has exponentially growing velocities (note I also set MASS_WEIGHT_IN_PRESSURE_GRADIENT = True, it is worse if both this parameter is False and there is a very small MIN_THICKNESS) .

So there doesn’t seem to be a winning scenario here - either we can choose to have very small velocities with no shelf and very large velocities with the shelf (dashed lines), or exponentially growing velocities with no shelf and similar order of magnitude velocities with a shelf. (Solid lines).



Does anyone have any insight into why a smaller MIN_THICKNESS would cause spurious velocities to grow without a shelf but the converse is true when a shelf is added? And any other thoughts on this asymmetry between ice shelf and topography?

 

[gmarques]

Hi Clare,

Great progress! You are getting close...

My next suggestion is to compute the 2D momentum budget. Here is an example of how to do that. The time series of each term in a location where spurious velocities occur might help us understand what's happening and how to contain the pressure gradient errors. What's the balance when spurions velocities occur near the bottom? what's the balance when spurions velocities occur near the ice shelf?

Vertical friction and bottom stress are proportional to h^-2 and h^-1, respectively, where h is the layer thickness. I think we want to get a balance between pressure acceleration and vertical friction/boundary stress acceleration.
Changing MIN_THICKNESS mostly affects thicknesses near the bottom and that's why you see spurious velocities grow without a shelf when using a smaller MIN_THICKNESS -- my guess is that the combination of DRAG_BG_VEL = 0.05 + small thicknesses near the bottom increases the bottom stress and, therefore, accelerates the flow. In other words, too much bottom stress can destabilize the flow and lead to exponentially growing velocities (via feedback). We need to find the "right" stress near the bottom and under the ice shelf to stabilize the pressure errors in these locations.

 

[cyung]

Thanks so much for your advice, Gustavo! The momentum budget helped to discover a problem with the ice shelf initialization in ALE coordinates, where the initial conditions are not quite what was expected - I have submitted a MOM6 issue here describing the strange nonlinear initialization. I was able to fix this initialization for a 2D case and get very low spurious velocities in sigma_shelf_zstar ALE coordinates (though have not worked out a clean way to fix the initialization bug). However, unfortunately, when topography is added with a grounding line, or a more complex 3D ice shelf geometry such as ISOMIP+ is used, we still have reasonably large (mm/s) currents, and I haven't been able to work out why.

I received advice from GFDL that the layered isopycnal coordinate (rather than ALE) works better, and have had success with that -- in layer mode, I get low spurious velocities in a resting case, and reasonable circulation with melt on (though this requires some modification of input parameters, e.g. turning on bulk mixed layer). So I'd recommend to anyone wanting to get a simple MOM6 ice shelf working to use layer coordinates. However, it will be necessary to use ALE coordinates in the future for more realistic models and more flexibility to test vertical coordinates, so I'd be keen to hear if anyone works out how to get low spurious velocities in ALE ice shelves in the future.