A few ideas for new diagnostics for wave-mean interaction problems

[glwagner]

Here are some ideas for additional diagnostics inspired by discussions with @jonathanlilly and @JeffreyEarly.

Eigensolver for vertical modes

I think it would be nice to have an eigensolver to compute vertical modes from a horizontally-averaged / reduced buoyancy profile (or taking in the vertical derivative of buoyancy directly). This means solving the eigenproblem

$$ \partial_z \left ( \frac{f^2}{N^2} \partial_z \phi_n \right ) + \frac{1}{R_n^2} \phi_n = 0 \qquad \text{with} \qquad \partial_z \phi = 0 \qquad \text{at} \qquad z = -H, 0 $$

where $1/R_n^2$ are the eigenvalues, $R_n$ is the Rossby radii of deformation of mode $n$, $f$ is the Coriolis parameter (a constant) and $N^2(z)$ is the vertical derviative of buoyancy (a function of $z$). The eigenfunctions $\phi_n(z)$ are the "vertical modes".

Finite volume details: $\phi$ is located at Center in the vertical, while N^2 is located at Face. This works out naturally so no shenanigans seem to be required here.

I think this can be implemented by meshing Oceananigans.AbstractOperations with Krylov.jl. In other words I am imagining a user interface wherein we write something like

B = Field(Average(b, dims=(1, 2)))
N² = ∂z(B)
ϕ = Field{Nothing, Nothing, Center}(grid)
L = ∂z(f^2 / N² * ∂z(ϕ))
sol = eigensolve(L)

where b would be the buoyancy field and f the Coriolis parameter. Many permutations on the above syntax might make sense too.

Diagnostics for wave-mean decomposition

Related to vertical modes would be diagnostics that allow one to decompose solutions into an internal wave and mean (or quasi-geostrophic) component. I think @JeffreyEarly knows how to do this well.

We can also support a direct solver for the Omega equation; eg solving

using a solver similar to the FourierTridiagonalPoissonSolver (combining horizontal-FFT with vertical tridiagonal). This would require modest (not heroic) efforts. The screenshot is from Danioux et al 2016.

[glwagner]

@tomchor we need an issue label for "new diagnostics idea" or something. Presumably ideas for new diagnostics are a prime thing to discsuss here.

[glwagner]

A related initiative is using Krylov to solve the Poisson equation in Oceananigans: CliMA/Oceananigans.jl#3812

@amontoison knows how to build matrix-free eigensolvers, I think.

[JeffreyEarly]

I think the omega-equation approach that you are suggesting is higher-order than what we are doing, but I don't understand the implications.

What we are doing is simply inverting QGPV using the projection operator in equation 78 here. To do the PV inversion with our lower order expression but using your approach you sketched out, I think you could compute,
$$\nabla^2 \psi +\frac{\partial }{\partial z} \left( \frac{f^2}{N^2} \frac{\partial \psi}{\partial z} \right) = Q$$
where Q is just QG PV computed from,
$$Q=\partial_x v - \partial_y u - f \partial_z \eta$$
which is equation 42 in the above paper.

[glwagner]

I think the objective is different, the omega equation obtains the balanced vertical velocity, while the QG inversion obtains the balanced streamfunction. So they are complementary / consistent approaches to obtain different variables perhaps?

It seems like the lower-order inversion could be a warm up to doing the higher-order inversion? Or should we ignore the low-order inversion.

It seems like it could make sense to solve a cleverly-forced tracer equation for the isopycnal displacement, rather than attempting to diagnose it. Does that make sense? I think that could be done fairly easily. The only fudging (associated with the staggered grid we use) is that the interface heights $\eta_e$ are at $w$-locations whereas a tracer would be shifted by half a cell. But another option could be to write a totally new dynamics that evolves the interfaces properly, at $w$ locations.

[JeffreyEarly]

My take on this is that you probably want to do QGPV inversion as a first pass, because it is really the only order at which everything is diagnostically crystal clear. The geostrophic streamfunction gives you $(u_g,v_g,\eta_g)$ from which you can deduce the wave portion of the flow without ambiguity---a separation of which is both complete and energetically orthogonal.

At the next order things get very murky: when you invert the $\Omega$ equation, you now have a geostrophic vertical velocity $w_g$ as well as corrections to $(u_g,v_g,\eta_g)$. But now there are many choices of equations of motion that result in these solutions which conserve (or don't) a variety of different expressions for PV and energy... so as a diagnostic it is more challenging to interpret. Remmel and Smith 2009 have a nice introduction summarizing various results, but maybe there's a newer take on this where things are more clear.

Seems like once you have the QGPV inversion the $\Omega$ equation would be fairly straightforward---so yeah! You should totally do both if it is not too hard.

For the QGPV inversion you can compute $\eta_e$ directly from density perturbation, you'll just need to scale it by $\bar{\rho_z}$. This is how we do it in the wave-vortex model. To compute the true APV, which requires the true isopycnal deviation, we start with a $\bar{\rho}$, and then use the bisection method using some code I borrowed from chebfun. It seems pretty fast, but I haven't fully validated its accuracy.