Unable to Replicate MITGCM Simulation

[kenflat2]

I am trying to replicate the channel model simulation from the paper "Submesoscale Vertical Velocities Enhance Tracer Subduction in an Idealized Antarctic Circumpolar Current" (Balwada 2016) in Oceananigans. Their setup is run in the MITGCM and captures baroclinic equilibration in a channel model of the Southern Ocean with a ridge at the bottom.

The largest discrepancy between their results and mine is much stronger eddies and somewhat noisy vorticity. These eddies can be seen in the surface vorticity field. Their version looks like this

while ours looks like this

and is greatly intensified in the South and somewhat intensified in the North. This is problematic because our simulation seems to have much stronger baroclinic eddy heat transport than the MITGCM simulation reports.

For reference, I compiled a table that compares the parameters of the MITGCM simulation to my version.

We use a 20km resolution here.

The main difference is in parameterizations - they use KPP and a Leith closure. We parameterize the mixed layer using an enhanced near surface vertical viscosity and diffusivity and use convective adjustment. I don't impose extra horizontal diffusivity outside of what WENO does on its own. I have tried alternate versions of these parameterizations but all version exhibit the strange turbulence near the Southern boundary.

The script is at this link: https://gist.github.com/kenflat2/a3ab79fe0588937987305e3e6b61b1e1

Any help would be greatly appreciated.

[glwagner]

Converting this to a discussion because at first glance this is a discussion about a set up, rather than a bug in the source code that needs to be fixed.

I think the main difference is probably using WENO rather than Leith and biharmonic. With WENO you have a higher effective resolution and therefore resolve more features, including stronger (you described it as "noisy") vorticity. Also the no-slip boundary condition has no effect if you don't have an explicit viscosity. So you need to add bottom drag, and possibly side drag (eg on north/south) as well if you want to replicate the effect of drag on side walls. To test this, you might try adding a HorizontalScalarBiharmonicDiffusivity. @simone-silvestri might know how to replicate the "A4GridMax" MITgcm parameter.

I think also that using a static vertical diffusivity rather than dynamic KPP could have the effect of producing weaker mixing. You'll miss stronger / more intense mixing events that will dissipate strong features by using "climatological mixing". You might want to try using CATKEVerticalDiffusivity.

[kenflat2]

I added a Biharmonic horizontal viscosity as follows

biharmonic_viscosity = (Δx ^ 4) / 7.0days
horizontal_closure = ScalarBiharmonicDiffusivity(TurbulenceClosures.HorizontalFormulation(), ν=biharmonic_viscosity)

and so far the simulation looks better. Here is a snapshot from the surface 180 days in

The strong vortices near the Southern boundary are greatly damped, as desired.

Does the biharmonic viscosity value above look reasonable? I attempted to set the parameter so that grid-scale biharmonic diffusion happens on a timescale of 7 days. This timescale was chosen to match the timescale for surface baroclinic eddies to form, which is about a week. The eddies look quite smooth in the above simulation so this may be overkill, please let me know what you think.

I didn't base my choice of parameter on the viscA4GridMax=0.8 parameter because, from my understanding of the MITGCM documentation, the viscA4GridMax parameter sets a weak CFL bound on viscosity and only constrains the viscosity to be less than (20kilometers^4)/(32*1200s) = 4.16e12m^4s^{-1}, which is high compared to what is in the literature.

[glwagner]

You've discovered that these motions should exist, and are artificially / incorrectly damped by the biharmonic viscosity, revealing a possible flaw in the cited study. I don't think any value about 0 is justified when you have sufficient dissipation for stability from WENO. Perhaps the motions are "undesirable", but that is an issue with the setup, or the way the simulation is forced, rather than with physics (eg physics is physics, and is neither desirable nor undesirable...)

[glwagner]

I suggest removing the biharmonic viscosity and experimenting with reducing the WENO order for WENOVectorInvariant, including simulations with order 5 and order 9. In that case you may have a smooth solution that at least partially resolves the motions. Another option is to change the physics of the setup to remove the motions, or to increase your resolution so they are more properly resolved.

[glwagner]

Simone's paper has additional helpful info about the unphysically large dissipation associated with explicit schemes like biharmonic and Leith, which can be mitigated using WENO:

https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2023MS004130

[kenflat2]

Our study isn't interested in these specific boundary conditions but in the Southern Ocean thermocline, which was impacted in a way not matching reality by the strong Southern eddies.

I modified the model physics to damp out these eddies by adding a quadratic drag with coefficient $\mu=1.2e-3$ to the immersed bottom boundary and the North and South walls and then disabled the biharmonic viscosity. Advection is 9th order WENO with no horizontal parameterizations.

This resulted in the above vorticity and temperature field at the surface (after 1 year of simulation time), which is the result we were hoping for.

Thanks for the help!

[glwagner]

Our study isn't interested in these specific boundary conditions but in the Southern Ocean thermocline, which was impacted in a way not matching reality by the strong Southern eddies.

I would challenge the conclusion that biharmonic viscosity makes the solution more realistic --- it's obviously the opposite. You should look for the flaw in the setup in the forcing; the flaw is merely masked by overdamping eddies. Plus, I still see evidence that energy is accumulating at your Southern boundary. It's important to understand why you see activity at the southern boundary, and address the underlying cause. The biharmonic viscosity just produces fluid dynamics that is not physical. The important progress is that you have identified why the MITgcm solution behaves the way it does. Moreover, you have showed that the MITgcm simulation has a flaw.