Resolved-scale kinetic energy budget with implicit dissipation

[sdbrenner]

Apologies in advance for the long-winded question.

I am interested in energy budgets in a set of simulations I am running. In this question I am looking for some feedback about the resolved-scale kinetic energy budget. There are a number of threads about the Turbulent Kinetic Energy budget (e.g., #3159, #3710, #3875), which have been informative, but I think my question here is distinct, and I haven't been able to find relevant discussion threads.

Specifically, I am interested in how I might calculate the equivalent resolved-scale dissipation when using a high-order WENO advection scheme with no explicit closure scheme. Since the dissipation is implicit, I’m exploring whether it can be inferred as a residual from the kinetic energy budget. I've come up with an approach, but would love some input about whether my logic is sound and/or if there's some other better approach I should take instead.

For context: my simulations are in a 50 km x 50 km doubly-periodic domain, started from rest. I force the model with a heterogeneous convective surface buoyancy flux. The heterogeneity leads to lateral buoyancy gradients within the domain, which eventually develop into submesoscale and mesoscale features. The surface buoyancy flux varies in space based on, effectively, a binary-mask covering some sections of the domain, so regardless of resolution, there are sharp spatial gradients in the buoyancy flux at the scale of the grid. I can run the model using a Centered advection scheme with a hyper-viscosity closure (HorizontalScalarBiharmonicDiffusivity); for other closures, the sharp gradients lead to excess noise in the solutions. By using the explicit closure with the Centered advection scheme, I've been able to close the resolved-scale KE budget (using functions from Oceanostics). However, an implicit closure through the use of a high-order WENO advection scheme (together with closure=nothing) produces sharper, smaller-scale features that would be nice to capture. But without an explicit closure, there is no model diffusivity νₑ field, so I don't see a clear way to quantify the associated resolved-scale dissipation.

My thinking was that when using WENO, I could estimate the dissipation as the residual of other terms in the KE budget. For the domain-integrated KE budget:
∫ᵥ∂ₜKE = ∫ᵥwb - ∫ᵥε,
where ε here is the dissipation of resolved-scale KE, I use Oceanostics to calculate:

∫ᵥε = Integral( Field(BuoyancyProductionTerm(ocnModel) - KineticEnergyTendency(ocnModel)) )

This seems to work out (see figure, below). The results are (I think) sensible, and compare well enough to simulations performed with hyper-viscosity closure (values differ as I didn't try to optimize diffusivities to make them match, but temporal patterns make sense).

However, I am also interested in the structure of the dissipation within the domain. I thought I could take the same approach and reconstruct the dissipation as the residual of the other KE budget terms:

KE = KineticEnergy(ocnModel) 
KET = KineticEnergyTendency(ocnModel) 
ADV = AdvectionTerm(ocnModel) 
wb = BuoyancyProductionTerm(ocnModel) 
PR = PressureRedistributionTerm(ocnModel, pressure=ocnModel.pressures.pHY′) 
ε = Field( -PR + wb - KET - ADV ) 

But this didn't work out; the reconstructed ε values are near-zero (shown below as a snapshot of horizontally-averaged fields):

I've realized that this is because the implicit dissipation is somehow wound up in the Advection term (which now doesn't integrate to zero over the domain, and—as seen in the above figure—has a structure quite similar to the Biharmonic ε). In retrospect, this makes some sense due to the fact that the model dissipation is an implicit result of the advection scheme, and the AdvectionTerm calculation from Oceanostics uses info from the scheme.

With that in mind, the approach that I'm now taking is to assume that the output of Oceanostics AdvectionTerm is a combination of both a "true" advective flux and the implicit/numerical dissipation. I then reconstruct the dissipation as the difference between the AdvectionTerm and a directly-calculated version of KE advection:

ADV = Field( AdvectionTerm(ocnModel) )
uⱼ∂ⱼKE = Field( @at (Center,Center,Center) u*∂x(KE) + v*∂y(KE) + w*∂z(KE) )
ε = Field( ADV - uⱼ∂ⱼKE )

The results from this approach look promising, with the reconstructed advection and ε having similar vertical structures to those from the Biharmonic scheme. But does this make sense?

[glwagner]

You are onto something I think.

I think the reason you cannot close the budget pointwise is because you are not taking into account the temporal discretization. Consider a tracer equation discretized with a foward Euler method:

$$ c^{n+1} - c^n = \Delta t G_c^n $$

where $G_c$ is some tendency. To derive an equation for $c^2$, we need to multiply this equation not by $c^n$, but by $(c^{n+1} + c^n) / 2$. This is the key. This produces

$$ \frac{1}{2} \left [ (c^2)^{n+1} - (c^2)^n \right ] = \frac{\Delta t}{2} \left (c^{n+1} + c^{n} \right ) G_c^n $$

The right hand side is an approximation to $\partial_t c^2 / 2$.

It should be immediately clear that your method is making a discretization error. For example, the function BuoyancyProductionTerm depends only on the fields at time step n. In other words, you have implicitly derived a variance equation using $c^n$, ignoring the temporal discretization. This introduces an error and you can't close the budget as a result.

It is possible to exactly compute the change in tracer variance every time-step taking into account the temporal discretization correctly and thereby compute the numerical dissipation rate, by expanding upon this concept. You have to take into account the entire time discretization of your problem (note: time-stepping methods are described here: https://arxiv.org/abs/2502.14148 but it is actually best to check the source code).

@simone-silvestri has succeeded in computing numerical mixing of tracers using a technique of this kind. He can comment further. Extending his work to kinetic energy is a bit more effort. But it is clearly important. It'd be worth of a paper and would be a significant contribution to Oceananigans.

I also personally believe that it is important to correctly compute terms in the kinetic energy and turbulent kinetic energy budget, from a scientific accuracy perspective. Coming up with discretely-correct diagnostics does require work, because one has to design a diagnostic that makes computations on consecutive iterations. Still this is not really all that complicated (time averaging is way worse). I suspect that computing budget terms in a discretely-correct way will reduce the averaging required for high quality results, too.

Computing the difference between Centered(order=2) and WENO is an interesting approach for sure. This must be making a temporal approximation for the reasons outlined above, but it is possibly closer to the truth than the budget approach. @simone-silvestri has toyed with that technique, too.

You might be interested in this paper: https://www.sciencedirect.com/science/article/pii/S1463500314000754

Tagging @jm-c and @tomchor as well.

[tomchor]

Thanks @sdbrenner, this is a very interesting question. I had discussions about this with @zhihua-zheng and if I recall correctly we actually came to a similar quantification to what you did (i.e. WENO5 implicit advection minus explicit advection of KE) but we never looked too much into it. So it's nice to see your results. @zhihua-zheng might have more to add.

I also agree with @glwagner that ideally budgets would be taking the time discretization into account, although I don't know if that'll make a bit difference. An investigation of that matter along with algorithms to close budgets to machine-precision would indeed be a significant contribution!

I have a plan to soon start looking more seriously into budgets in Oceanostics. The plan is to add a couple of doc examples (or modify the existing ones) where we actually go over these questions and do a complete budget using WENO, and probably in the process develop some functionality that helps with that. If you're interested in helping out, please let me know!