Do we want to try well-balanced schemes?

[francispoulin]

A well-balanced numerical scheme is one that preserves the steady-state solution exactly (to numerical precision). Most schemes are not, and this holds true for Oceananigans.

Here are a couple of references that might be of interest.

This article shows one approach for low order finite volume. Looks pretty simple but I admit that I still need to sort out the details.

And this one here, presents an approach that establishes a high-order well-balanced finite volume WENO scheme for a shallow water model that allows for wetting and drying. This looks more complicated but since it uses WENO, I thought it might be of interest.

Is there any interest in trying a well-balanced scheme? If yes, I would suggest first trying this for ShallowWater and if it's not too difficult and has some benefits, maybe we could extend it to HydrostaticModel?

Thoughts?

[glwagner]

I have to admit that I don't totally understand what the problem is that this scheme solves. The idea is that our current numerical methods don't have the same geostrophic basic states as the numerics? Is there an example of a solution / setup which is clearly deficient in this respect and needs to be improved?

There's a related problem --- accurate calculation of the hydrostatic pressure gradient in a warped / non-geopotential coordinate system, which is well-known to be one of the most difficult aspects of implementing a general-vertical-coordinate hydrostatic model. The latest I've heard on this problem is https://www.sciencedirect.com/science/article/pii/S1463500317300665?casa_token=CTwbsGJFWusAAAAA:QTMbdHl-vWx59_kZwPXNutYbJ6zjVXqvdzvZaNMot-w6MD99zz8I55B2xcrwXaVXWY0tLDc

But that problem is about the hydrostatic pressure gradient / hydrostatic balance, rather than geostrophic balance. The "problem" for hydrostatic balance is well illustrated by simulating an ocean at rest in a warped coordinate system (which is the example used in the above paper).

I don't quite see how these papers are related though --- can you help?

I think the main downside is that this scheme looks like it requires a totally different time-stepping method. So, it might be better to investigate in a simpler setting than Oceananigans. Not sure.

[francispoulin]

I don't know of any examples in Oceananigans that is deficient because of this problem. I am not claiming this is a problem but it is a possible treak that might be an improvement? But I honestlyl don't know as I haven't tried it yet.

This is independent of computing the hydrostatic pressure, and can perhaps be best understood in the context of the ShallowWaterModel.

If one wanted to initialize the model with a geostrophically balanced jet the first thing one might think of doing is to initialize the velocity and height with analytical expressions, which is not wrong. However, if you do that, you will find that the solution will evolve slightly at the first time step, and maybe more, because the discretizastions of the analytical solutions are not well-balanced. Small changes changes are bound to happen. What they do given time depends on the stability of the system.

Another approach is to specify the height field and then compute the velocity numerically. This is what @glwagner, @fadaie91 and I have been doing recently. This, I beleive, gives a well-balanced solution since the time deperivatives will be zero if the solution is an exact solution. This approach is much simpler than what's talked about in these papers, and probably more attractive for that reason alone. I will continue to do this in the future and will put these papers on the backburner. If other problems arise then maybe I'll revist them.

Thanks for the discussion @glwagner and happy to put this to rest, unless others have other thoughts.

[glwagner]

If one wanted to initialize the model with a geostrophically balanced jet the first thing one might think of doing is to initialize the velocity and height with analytical expressions, which is not wrong. However, if you do that, you will find that the solution will evolve slightly at the first time step, and maybe more, because the discretizastions of the analytical solutions are not well-balanced. Small changes changes are bound to happen. What they do given time depends on the stability of the system.

Ah, I understand better now! Yes, we have the concept of "numerical geostrophic balance" (and also "numerically divergence-free"), both of which require defining a streamfunction rather than a velocity field.

Perhaps what is still unclear to me is whether that slight difference between numerical balance and analytical balance degrades the quality of solutions when integrated forward, or whether its purely a problem with defining balanced states.

Maybe one sample problem would be a non-trivial surface stress distribution that should correspond with a particular pattern of Ekman transport. The numerical solution might then differ from the analytical solution due to the difference between numerical geostrophic balance and analytical geostrophic balance; one could argue that this has an effect on real circulations (where either observed winds, or winds derived from an atmospheric model, forces near-surface oceanic momentum).

I think coming up with a clear test case that reveals the deficiency is the first thing to do (and perhaps worth a publication in its own right). Then with a strong motivation in hand we'd be well-positioned to invest some time to find a solution.

[francispoulin]

The scientific computing communittee has done a lot fo work in the area, so they think this is important, but I also share your cocnerns. We don't just want to start off with the steady solution and remain with the steady solution forever, as that would be pretty boring. Does this affect the time-evolving dynamics? I will glance more at the papers and see if I can gather some intel.

I like your suggestion about a non-trivial surface stress distribution but I don't know the problem well enough to say. I would be happy to discuss further if you like. I imagine an imbalance in the Ekman layer could cause some strange things to go on, maybe?

I would guess that in most cases the difference you get from a non-balanced scheme would be like noise in the forcing, if you have it, or maybe perhaps numerical error. I guess it depends on how big it is.

I do have an idea, which involves linear stabilty theory. When we are looking at the stability of a jet, and there is enough symmetry, often we can boil it down to a 1D eigenvalue problem, and that is something that can be solved easily. However, even 2D eigenvalue problems can be difficult. I know because I had a PhD student who did this a couple of years ago. As an alternative, what you can do is integrate the linearized equations forward in time and see what develops. If it converges to something than that is the most unstable mode. I believe Oceananigans does this with Kelvin-Helmholtz instability, but the perturbatoins are normalized instead of solving the equations. Both should yield the same thing.

I could imagine that in this case, the error you get could give rise to spurious modes, as it occurs at every time step. Or at the very least it could cause you to converge more slowly than if you didn't have this noise.

Anyhow, my two cents worth for now.