WENO reconstruction of energy density produces diabatic numerical fluxes

[glwagner]

We have formulated an anelastic equation set that prognoses moist static energy density, which we call $\rho e$. To compute the advection of $\rho e$, we reconstruct it with WENO and advect with model.velocities:

Breeze.jl/src/AtmosphereModels/dynamics_kernel_functions.jl

Line 189 in 71a4ae2

return ( - div_Uc(i, j, k, grid, advection, velocities, energy_density)

As pointed out on #214 this leads to discrete dynamics that dissipate the variance of energy density --- leading to the homogenization of energy density. This behavior, I believe, is the reason why the bomex case develops a spurious cold layer at the bottom. To illustrate this I doubled the resolution, turned off the forcing and boundary fluxes, and ran for one hour. Below I am plotting the horizontal averages of energy density, specific energy $e$, and potential temperature:

The energy density has a strong gradient from the initial condition, decrease from its maximum value at the bottom of the domain. After one hour, due to noise in the initial condition, the energy density has become slightly homogenized (left plot). This is the correct response to WENO advection / dissipation. However, this homogenization in turn produces a spurious cold layer as seen in potential temperature (rightmost panel).

This situation should be nearly adiabatic, which means that it is potential temperature (or equivalently, entropy) which should remain constant --- not energy. That's why the cold layer is spurious.

The best fix for this as far as I can tell is to prognose entropy or potential temperature instead of energy. Evolving entropy for sure is an equivalent application of the first law of thermodynamics (which our moist static energy equation is also designed to respect a la Pauluis). Potential temperature probably can do the same if formulated correctly.

Advecting specific energy rather than energy density could be a hot fix because specific energy has weaker gradients than energy density. I will explore this and post results here. But ultimately I think we probably want to have at least the possibility of evolving potential temperature.

[glwagner]

Reflecting on this more I am convinced that 1) we need an entropy or potential temperature formulation (or both) and 2) we need to do tendency computations in terms of specific quantities. For example with a potential temperature formulation we would evolve $\rho \theta$, but we want no mixing for constant $\theta$. This means we need to reconstruct $\theta$, not $\rho \theta$, during advection. A similar issue holds for diffusive fluxes.

Because we need specific quantities to compute tendencies (rather than densities, for the most part), we need to either 1) store all specific quantities paired up with their density counterparts or 2) convert densities to specific quantities immediately prior to computing tendencies.

I am inclined to do the second (except for velocities/momentum -- we should probably keep velocities around in addition to momentum)? I actually think it is simpler just to temporarily convert densities to specific quantities and then convert them back, rather than having to keep things around...

I think we should also build a user interface wherein users can universally deal with specific rather than density quantities. Basically a density quantity is one prefixed by :ρ.

As for microphysical velocities/momentum --- there is actually no issue for tracer advection, because we do not reconstruct momentum / velocities on the c-grid. We can add density to an advection term in place no problem. It will look similar to what our diffusive fluxes look like right now.

[simone-silvestri]

Since this behavior is restricted to the boundary, it might be the fault of the very low-order scheme that we default to near the boundary. Something very similar is happening in the internal_tide example in oceananigans, where the stratification is destroyed at the top and bottom boundary. That example is fixed by setting

tracer_buffer_scheme = WENO(order=5, buffer_scheme=Centered())
tracer_advection     = WENO(order=7, buffer_scheme=tracer_buffer_scheme)

(I am yet to open an issue / PR on this)
Even if the problem might be more about the formulation of what quantity is advected, maybe using this setting also helps for this case?

[glwagner]

That probably will help indeed.

[glwagner]

Why don’t we just make the default buffer scheme for order 3 or 5 into a Centered?

[glwagner]

Also posted on #214

Reconstructing specific quantities and using a Centered scheme near the boundary does seem to help:

However it looks like we still get the best results with an explicit turbulence closure for the anelastic bomex case:

[glwagner]

Worth stating, if we implement @simone-silvestri's suggestion --- in addition to making the change that advects specific energy rather than energy density --- with the unforced, no-boundary-condition case, we get

specifically I wrote

buffer_scheme = Centered(order=2)
buffer_scheme = WENO(order=3; buffer_scheme)
buffer_scheme = WENO(order=5; buffer_scheme)
buffer_scheme = WENO(order=7; buffer_scheme)
advection = WENO(order=9; buffer_scheme)

[bischtob]

@glwagner, question about the advection of specific variables vs. densities here: since we are using this thermodynamic formulation for energy conservation properties (among other things), I am wondering if this change plays nicely with those properties.

[glwagner]

So, the difference here is between these two options:

$$ \rho u e \approx \rho u [ e ] $$

versus

$$ \rho u e \approx u [ \rho e ] $$

where $[\phi]$ indicates a reconstruction. For second-order reconstruction $[\phi]$ is a simple average, something like $[\phi] = \left ( \phi_i + \phi_{i+1} \right ) / 2$. For WENO we do a biased, weighted reconstruction which depends on the smoothness of $\phi$ and also the sign of $u$ / $\rho u$.

The first option reconstructs the specific energy and produces a numerical flux of the form

$$ \rho u e \approx \rho u [e] \approx \rho u \langle e \rangle + C \Delta_x^n |\rho u| \partial_x^n e $$

for some unknown constant $C$ and an $n$ th order scheme (this is why $n$ must be odd-ordered to produce diffusion rather than dispersion), where $\langle e \rangle$ is centered reconstruction.

Global conservation properties are not affected because we still use a flux form, ie the momentum flux divergence is still

$$ \nabla \cdot \left ( \rho \boldsymbol{u} e \right ) $$

either way. The difference is just the properties of the flux $\rho \boldsymbol{u} e$.

There is more to say here though. Something that I can't find an appreciation for in the literature regarding ILES is that the numerical fluxes are not the same as the analytical / continuous fluxes of quantities. Take for example Romps 2008 energy equation:

the last three terms represent energy transport due to differential fluxes (denoted $d_v$ for vapor, etc) of the vapor, liquid, ice components relative to the bulk.

$\boldsymbol{d}_k$ are confusingly named "diffusion fluxes" but also include sedimentation velocities, so they aren't all diffusive. But they should also include things like SGS turbulent diffusion, since the governing equations for the water species are

where $e$ and $m$ represent phase change but not diffusion.

Here's the tricky part. For the budget to close numerically, I don't think we can conceive of $\boldsymbol{d}_k$ as encompassing just prescribed sedimentation velocities and differential SGS fluxes. In an ILES there is also a differential numerical flux if we use independent, diffusive (eg WENO) reconstructions of water species and energy, separately.

Basically we need to construct the energy equation from the discrete equations, not the analytical equations, for this to work correctly. If we try to build the energy equation from the continuous equations but then apply inconsistent numerical approximations, we don't end up with discrete global energy conservation.

You can see this is going to be pretty tricky. Or rather, most pressingly, we can't use schemes that strictly diffuse $e$ while also conserving it, because we need to use an advection for $e$ that is consistent with advection of $q$ instead. Would that be stable? Maybe. It's also a bit of effort (but definitely doable, that said). Romps (2008) later describes a code with diffusive advection scheme but doesn't mention how discrete total energy conservation might be achieved.

Maybe it depends how careful you want to be? Most models don't even use an analytical equation set that conserves energy, so issues of inconsistent numerical approximations are totally in the weeds. This is discussed by Bryan and Fritsch 2002. Also you can see in their figure 7, even with "wrong equations" we get energy conservation to within 0.001%. But they also argue that micropercentages do matter...

The issue is not escaped by using an entropy scheme. For entropy formulations we have to do some detailed bookkeeping to include entropy sources for every irreversible process. For example, Pressel et al 2015 include the entropy source due to SGS diffusion:

Well... if we have numerical diffusion of momentum, then don't we also need to include the numerical entropy source due to momentum diffusion in an ILES? (There is a corresponding term associated with SGS diffusion of moisture species.)

On the basis of all of this I think that energy or entropy conservation with diffusive numerics might be a wild goose chase. With Herculean effort we might get there eventually...

[glwagner]

#228 shows the advantage of a potential temperature (and entropy as well I think) formulation for ILES.