[MITgcm-support] Anomalous 'rebound' circulations at poles in 2D setup leading to numerical instability

[Flynn Ames]

Dear MITgcm Community,

Hello - I've been developing a 2D (latitude-depth) model setup for the MITgcm that extends pole to pole, for the subsurface ocean of Enceladus. I've been running into numerical stability issues at the poles. I have some ideas of how to overcome the problem but I'm hoping to get some thoughts on if any of them are particularly bad/good ideas, or anyone else has an alternative idea?

The numerical instability appears to arise from a 'rebound effect' of sorts where, owing to the 2D domain, the flow upon reaching the poles has nowhere to go and rebounds off of it. These circulations are always present but when the vertical viscosity is low (and for some reason, when along isopycnal mixing is considered), these anomalous circulations intensify and eventually blow up.
I've found that increasing the vertical viscosity can prevent the solution blowing up, but this is not ideal because it removes interesting features of the ocean circulation that I want to maintain and leads to an unreasonably high Prandtl ratio in the vertical. Similarly, I would like to consider the role of along-isopycnal mixing. Reducing the timestep dramatically is also not ideal as these simulations are computationally expensive owing to long equilibration times.

It's worth noting these anomalous circulations don't appear to alter the mean solution when running stably. I've attached snapshots of meridional velocity (hopefully it sends!), showing what these anomalous circulations look like when the model runs stably (with a reasonable vertical viscosity BUT no consideration of along-isopycnal mixing). My thinking is that these circulations are likely unphysical, and so some of the below approaches (NOTE: not ideal solutions) for overcoming the issue could potentially be justified (?). They are as follows:

1.) Try masking the grid cell directly over each pole with land and setting NO_SLIP_SIDES=.TRUE. This may damp the anomalous circulations, but comes at the cost of losing part of the model domain
2.) Use the RBCS package to create a sponge layer of sorts to relax meridional velocity to zero at the poles (with a few grid points next to the poles having weaker restoring). This would let me keep the whole domain. I don't think it should affect the mean solution significantly, but I wonder if this is an acceptable thing to do?
3.) Try implement a feature in the MITgcm that would let me specify a meridional profile for vertical viscosity, so I can keep it lower over most the domain, but increase it where numerically necessary (i.e., at the poles). This could let me hold onto features I want to maintain over most the domain. But once again is this an acceptable thing to do?
4.) An idea my project supervisor suggested is applying higher viscosity on divergent flows vs convergent flows. This could damp these anomalous circulations but would apply over the whole domain so could affect the mean solution.

I'd be very grateful for any thoughts or suggestions. If anyone has any questions about this, I can definitely answer them!

Cheers,
Flynn

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20230307/3aa64b32/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Vvel_snapshots_anom_circ.gif
Type: image/gif
Size: 1792751 bytes
Desc: Vvel_snapshots_anom_circ.gif
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20230307/3aa64b32/attachment-0001.gif>