An example with a non-trivial thermal-wind balance

[tomchor]

The issue of calculating a thermally-wind-balanced background has popped up a couple times before (more recently here) and I think it would be nice to have an example where b is computed from u for a general function. I myself don't know how to do that properly and efficiently (I think), so it would definitely be helpful for a larger audience I believe.

@glwagner @francispoulin

[glwagner]

Hmm. I think given general u, v we need to solve an elliptic equation to find psi, eg lap(psi) = dx v - dy u. Then the geostrophic buoyancy field is b = f * dz(psi). For some boundary conditions or assumptions of homogeneity, we can possibly simplify this procedure. One could in principle use the conjugate gradient solver to write a fully general utility, though. Hopefully this task will be taken up once the hydrostatic model is operational. Some care must be taken with the geostrophic mode / computation of the geostrophic free surface displacement for the hydrostatic model.

Or, given a geostrophic buoyancy field we can integrate in z (eg summation of b * dz). This would be the intent of IntegratedField, but nobody has implemented that yet.

Absent those techniques I think the best method for idealized problems is to use a geostrophic streamfunction. This is illustrated by the Eady turbulence example:

Oceananigans.jl/examples/eady_turbulence.jl

Lines 38 to 41 in 0153dec

	# The geostrophic basic state in the Eady problem is represented by the streamfunction,
	#
	# ```math
	# ψ(y, z) = - α y (z + L_z) \, ,

Oceananigans.jl/examples/eady_turbulence.jl

Lines 151 to 152 in 0153dec

	U(x, y, z, t, p) = + p.α * (z + p.Lz)
	B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z

It's also used in a number of validation cases, but those aren't always useful for general audiences.

[tomchor]

Or, given a geostrophic buoyancy field we can integrate in z (eg summation of b * dz). This would be the intent of IntegratedField, but nobody has implemented that yet.

That was more along the lines of what I was thinking. Basically showing users how to create a u, find out what b is, and then indicate if it's best to use that as BackgroundField, or an Array or something.

But it seems like IntegratedField needs to be in place for something as simple as what I'm envisioning?

[glwagner]

Or, given a geostrophic buoyancy field we can integrate in z (eg summation of b * dz). This would be the intent of IntegratedField, but nobody has implemented that yet.

That was more along the lines of what I was thinking. Basically showing users how to create a u, find out what b is, and then indicate if it's best to use that as BackgroundField, or an Array or something.

But it seems like IntegratedField needs to be in place for something as simple as what I'm envisioning?

Just to clarify, the quoted example suggests integrating b to find psi.

If u is a function of y only I think we can integrate the relation dy(psi) = u to find psi.

If the geostrophic streamfunction is a function of x, y then I think solving an elliptic problem is required, probably?

To summarize:

• To obtain psi(x, y, z) from b(x, y, z), we integrate the relation b = f * dz(psi). IntegratedField is useful.
• To obtain psi(x, y, z) from from u(x, y, z), v(x, y, z) I think (might be worth thinking about this more) we solve the horizontal elliptic problem hlap(psi) = dx(v) - dy(u) at every vertical level.
• To obtain a 2D psi(y, z) (for example), we might be able to integrate dy(psi) = u.

These techniques could be important when initializing from data, perhaps. I'm not sure. For example, if a 3D distribution of T, S were obtained from data without knowledge of the velocity field, it might make sense to initialize the problem with a geostrophic velocity field. It'd be great to have utilities for all of these cases eventually.

For idealized work I think its better to specify the streamfunction psi. This method is demonstrated by the Eady turbulence example.

[tomchor]

What is hlap?

I was thinking more about the 2D example. From my experience most thermal wind cases are rotated in a way that they are treated as a 2D problem. I agree that specifying Ψ (which I'm assuming you meant as the streamfunction) is the preferred way, but that's not always easy and it would be nice (IMO) to have a documented way to do a numerical approach.

I definitely think that coming up with something to solve the more general case of getting Ψ from a general u and v is overkill.

[glwagner]

What is hlap?

horizontal Laplacian; hlap(psi) = dx^2(psi) + dy^2(psi).

[glwagner]

I definitely think that coming up with something to solve the more general case of getting Ψ from a general u and v is overkill.

This is useful for GCM experiments; I think we'll want it eventually.

[glwagner]

To clarify: I think there is utility in a simple example. However, I also recommend using a streamfunction for idealized cases, since its easier to reproduce and understand. Idealized cases with simple geostrophic flows include experiments like the Eady turbulence problem, or LES cases that simulate turbulence in the presence of a simple geostrophic flow.

Another class of experiments (likely more prevalent and important than LES around geostrophic flows?) are realistic / quasi-realistic simulations on the sphere initialized from data or complicated buoyancy distributions. For these it's often important to find a flow field that's in discrete geostrophic balance with an arbitrary buoyancy field. We can show how to integrate the buoyancy field for this purpose; but having IntegratedField implemented will make it a bit cleaner.

The problem of finding a streamfunction from a velocity field is likewise be important for simulations on the cubed sphere; or cases in which the grid coordinate system does not align with, for example, a geographic coordinate system. This is the case that @francispoulin was attempting to solve in #1826, I think. @navidcy may illustrate how to do this in post processing using the conjugate gradient solver in his horizontal convection example, perhaps.

[francispoulin]

This all sounds very good and what @glwagner suggests will certainly work, and yes, that is exactly what I was trying to do in #1826.

As a question of style, introducing a streamfunction is convenient for those of us that like streamfunctions, but really we're computing the pressure. Since the governing equations that we state for IncompressibleModel and HydrostaticModel are in terms of pressure, then that might be good as an introductory example.

Oceananigans has a pressure solver already. Since the basic profiles that we consider have nonlinear terms, then I wonder if it's possible to use the built in pressure solver to determine p and then compute b using a derivative operator. Clearly the functionality is there and being able to solve for the pressure of your initial conditions, which is what we're really doing, might be helpful.

[glwagner]

I'll convert this to discussion on how to implement simulations in thermal wind balance...