Closing Kinetic Energy budget

[zhihua-zheng]

This is meant to be a discussion. I need some help to interpret my KE ($k = u_i u_i /2$) budget calculations, which follow this equation:

$$\frac{\partial k}{\partial t} = -\frac{\partial u_j k}{\partial x_j} - \frac{\partial u_i p}{\partial x_i} + u_i b_i - u_i\frac{\partial \tau_{ij}}{\partial x_j}$$

These 5 terms correspond to the KineticEnergyTendency, AdvectionTerm, PressureRedistributionTerm, BuoyancyProductionTerm, KineticEnergyStressTerm in Oceanostics' TKEBudgetTerms, respectively.

When integrated in volume for the entire domain, this becomes

$$\frac{\partial (\int k dV)}{\partial t} = -\int u_j k \cdot n_j dS - \int u_i p \cdot n_i dS + \int u_i b_i dV - \int u_i\frac{\partial \tau_{ij}}{\partial x_j} dV$$

For horizontally periodic, vertically bounded domain, the integrated AdvectionTerm and PressureRedistributionTerm should be zero.

This is the code I used:

using Oceananigans
using Oceananigans.Units: minutes, days
using Statistics
using Oceanostics
using Oceanostics.TKEBudgetTerms: KineticEnergyTendency, AdvectionTerm, KineticEnergyStressTerm,
                                  PressureRedistributionTerm, BuoyancyProductionTerm

grid = RectilinearGrid(GPU(), size=(32, 32, 32), extent=(128, 128, 64),
                       topology=(Periodic, Periodic, Bounded))

const hᵢ = 20
const Qᵘ = -3.72e-5
const N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient
b_boundary_conditions = FieldBoundaryConditions(bottom = GradientBoundaryCondition(N²))
u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

model = NonhydrostaticModel(; grid,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (b=b_boundary_conditions, u=u_boundary_conditions))

@inline Ξ(z) = randn() * exp(z / 4)
@inline stratification(z) = z < - hᵢ ? N² * z : N² * (-hᵢ)
@inline bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz
@inline uᵢ(x, y, z) = 1e-1 * Ξ(z)
@inline wᵢ(x, y, z) = 1e-1 * Ξ(z)
set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

simulation = Simulation(model, Δt=45, stop_time=2days)
wizard = TimeStepWizard(cfl=0.65, diffusive_cfl=0.65)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

## KE budget using resolved velocities
KE = Field(Integral(KineticEnergy(model)))
KE_tendency = Field(Integral(KineticEnergyTendency(model)))
KE_advection = Field(Integral(AdvectionTerm(model)))
KE_stress = Field(Integral(KineticEnergyStressTerm(model)))
KE_pressure = Field(Integral(PressureRedistributionTerm(model)))
KE_buoyancy = Field(Integral(BuoyancyProductionTerm(model)))
outputs = (; KE, KE_tendency, KE_advection, KE_stress, KE_pressure, KE_buoyancy)

output_interval = 5minutes 
simulation.output_writers[:fields] =
     NetCDFOutputWriter(model, outputs,
                     schedule = TimeInterval(output_interval),
                     filename = "mwe_KE_budget.nc",
                     overwrite_existing = true)
run!(simulation)

The result from the simulation is shown above. The residual is nearly 0. I mainly have two points of confusion:

1. The nonzero integral of AdvectionTerm conflicts with what I imagined from the equation.
2. The KineticEnergyTendency being mostly negative conflicts with the fact that the integral of KE grows with time.

I could have overlooked something here, please let me know what you think.

[zhihua-zheng]

On the second point, I think it is because KineticEnergyTendency does not include boundary tendency from prescribed surface fluxes?

[tomchor]

@zhihua-zheng thanks for opening this! Here are a couple of comments:

These 5 terms correspond to the KineticEnergyTendency, AdvectionTerm, PressureRedistributionTerm, BuoyancyProductionTerm, KineticEnergyStressTerm in Oceanostics' TKEBudgetTerms, respectively.

This is not important here, but just to clarify, the KETendency isn't exactly $\partial k / \partial t$. It is that, minus the contribution from the nonhydrostatic pressure. This is said in the docstring, but I just realized it's incomplete!:

Oceanostics.jl/src/TKEBudgetTerms.jl

Lines 129 to 138 in ffa54aa

	"""
	$(SIGNATURES)
	Return a `KernelFunctionOperation` that computes the tendency of the KE except for the nonhydrostatic
	pressure:
	```julia
	```
	"""
	function KineticEnergyTendency(model::NonhydrostaticModel; location = (Center, Center, Center))

In this case the integrated contributions from the nonhydrostatic pressure are zero, so once you integrate the tendency there's no difference.

1. The nonzero integral of AdvectionTerm conflicts with what I imagined from the equation.

This indeed is weird. My first guess is that this might come from the momentum boundary condition, although I thought that was implemented through the stress term. (I'm actually having a hard time finding where exactly that BC is applied, but I'll look more into it later.)

2. The KineticEnergyTendency being mostly negative conflicts with the fact that the integral of KE grows with time.

That is also weird. The first thing I'd compare is the KineticEnergyTendency with the actual calculated $\partial k / \partial t$ from KineticEnergy. I'm guessing there might be a sign wrong somewhere...?

The way I'd move forward is to remove the buoyancy forcing, keep only the momentum forcing for simplicity, and add Coriolis. That way the BL depth will equilibrate (at $\sim u_\star / f$ if I recall correctly) and the whole system reaches a steady state. Then you can plot the relevant terms as a function of $z$. This should give us insight into where the advection contribution is coming from. Then we can deal with the tendency later.

Also a note: once we understand this better, I think it would be great to convert this code into an example of KE budgeting for the docs!

[glwagner]

Can you link to the code where these things are implemented? AdvectionTerm must somehow be a numerical approximation, so perhaps there is a discrete discrepany. Ultimately to get discrete correctness we have to derive every term from the discretely exact kinetic energy equation. This can be a little onerous though...

As a sanity check and to shed some light on whether discreteness is playing a role --- you should try repeating the problem with different time-steps (perhaps 10x smaller) and also on different resolution grids.

[zhihua-zheng]

Can you link to the code where these things are implemented?

It is using the velocity multiplied by the advective flux divergence kernel:
https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L161C1-L166C4

[zhihua-zheng]

That is also weird. The first thing I'd compare is the KineticEnergyTendency with the actual calculated $\partial k/ \partial t$
from KineticEnergy. I'm guessing there might be a sign wrong somewhere...?

They don't match, the blue line shows the KE grows with time.

[zhihua-zheng]

Following your suggestions, I have done two more experiments, both with finer stretched vertical grid (~ 0.5 m at surface), Coriolis and wind forcing only. Different advection schemes are used, and results are sensitive to this choice (sign of numerical dissipation?).

Below I show two figures for each experiment: evolution of vertically integrated horizontal mean budget term; time averaged horizontal mean of budget term as a function of depth. Note that different than the figure posted above, I have included the sign of each RHS term here to indicate the contribution to local KE tendency. The top most grid point is neglected due to outstandingly large value, and this appears to bring the KineticEnergyTendency inline with actual change of KineticEnergy with time.

The first one uses CenteredSecondOrder().

With WENO(order=5):

[tomchor]

Okay, just to make I'm understanding the main takeaways here:

1. You observe qualitative differences when you change advection scheme. For example with the CenteredSecondOrder() the main balance is between the advection term creating KE and the stress term removing it. With WENO() the balance seems similar, but the sign of the contributions is flipped, no?
2. The first grid point of the AdvectionTerm is unreasonably large and ignoring it improved results (specifically dk/dt matches the tendency term). @glwagner can this be related to the halos not being filled? I don't really understand very well how halo filling works in the code, but I can't think of anything else.

Also can you pelase share the new code along with the code used to plot the figures? I wanna do some investigating of my own.

[tomchor]

It boggles me the most that the contribution from the AdvectionTerm isn't zero even far from the surface. Inspecting the code, I don't see anything wrong with it. Like @zhihua-zheng mentioned, this is the relevant kernel:

Oceanostics.jl/src/TKEBudgetTerms.jl

Lines 161 to 166 in ffa54aa

	@inline function uᵢ∂ⱼuⱼuᵢᶜᶜᶜ(i, j, k, grid, velocities, advection)
	u∂ⱼuⱼu = ℑxᶜᵃᵃ(i, j, k, grid, ψf, velocities.u, div_𝐯u, advection, velocities, velocities.u)
	v∂ⱼuⱼv = ℑyᵃᶜᵃ(i, j, k, grid, ψf, velocities.v, div_𝐯v, advection, velocities, velocities.v)
	w∂ⱼuⱼw = ℑzᵃᵃᶜ(i, j, k, grid, ψf, velocities.w, div_𝐯w, advection, velocities, velocities.w)
	return u∂ⱼuⱼu + v∂ⱼuⱼv + w∂ⱼuⱼw
	end

And that's being called from here:

Oceanostics.jl/src/TKEBudgetTerms.jl

Lines 196 to 199 in ffa54aa

	function AdvectionTerm(model::NonhydrostaticModel; velocities = model.velocities, location = (Center, Center, Center))
	validate_location(location, "AdvectionTerm")
	return KernelFunctionOperation{Center, Center, Center}(uᵢ∂ⱼuⱼuᵢᶜᶜᶜ, model.grid, velocities, model.advection)
	end

What the code is doing, using index notation, is

$$\begin{align} u_i \partial_j(u_j u_i) &= \partial_j(u_j u_i^2) - u_j u_i \partial_j u_i\\\ &= \partial_j(u_j u_i^2) - u_j \partial_j u_i^2/2\\\ &= \partial_j(u_j u_i^2)/2 - u_i^2 \partial_j u_j/2 \end{align}$$

I think this is highly unlikely, but is it possible that $\partial_j u_j$ is not zero? i.e. can the flow be erroneously divergent here?

But also, once you remove the first grid point, it seems that the budget is approximately closing, indicating that the diagnostics are (apart from the first point) quantifying things correctly... or at least with compensating errors...

[zhihua-zheng]

Also can you pelase share the new code along with the code used to plot the figures? I wanna do some investigating of my own.

Sure. Code for simulation:

using Oceananigans
using Oceananigans.Units: minutes, days
using Oceananigans.Advection: div_Uc
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceanostics
using Oceanostics.TKEBudgetTerms: KineticEnergyTendency, AdvectionTerm, KineticEnergyStressTerm,
                                  PressureRedistributionTerm, BuoyancyProductionTerm

function  KineticEnergyAdvectionTerm(model::NonhydrostaticModel, ke; velocities = model.velocities, location = (Center, Center, Center))
    return KernelFunctionOperation{Center, Center, Center}(div_Uc, model.grid, model.advection, velocities, ke)
end

const Nz = 64
const Lz = 50
const refinement = 1.5 # controls spacing near surface (higher means finer spaced)
const stretching = 10  # controls rate of stretching at bottom
h(k) = (k - 1) / Nz
ζ₀(k) = 1 + (h(k) - 1) / refinement
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(GPU(), size=(32, 32, Nz), x=(0, 100), y=(0, 100), z=z_faces,
                       topology=(Periodic, Periodic, Bounded))

const hᵢ = 20
const Qᵘ = -2e-4
const Qᵇ = 0#1.3e-7
const N² = 6.4e-5 # s⁻², initial and bottom buoyancy gradient
b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
                                                bottom = GradientBoundaryCondition(N²))
u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

model = NonhydrostaticModel(; grid, coriolis = FPlane(f=1e-4),
                            advection = WENO(),#CenteredSecondOrder(),#WENO(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (b=b_boundary_conditions, u=u_boundary_conditions))

@inline Ξ(z) = randn() * exp(z / 4)
@inline stratification(z) = z < - hᵢ ? N² * z : N² * (-hᵢ)
@inline bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz
@inline uᵢ(x, y, z) = 1e-1 * Ξ(z)
@inline wᵢ(x, y, z) = 1e-1 * Ξ(z)
set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

simulation = Simulation(model, Δt=5, stop_time=4days)
wizard = TimeStepWizard(cfl=0.65, diffusive_cfl=0.65)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

# KE budget using resolved velocities
ke = Field(KineticEnergy(model))
KE = Field(Average(ke, dims=(1,2)))
KE_tendency = Field(Average(KineticEnergyTendency(model), dims=(1,2)))
KE_advection = Field(Average(AdvectionTerm(model), dims=(1,2)))
#KE_adv = Field(Average(KineticEnergyAdvectionTerm(model, ke), dims=(1,2)))
KE_stress = Field(Average(KineticEnergyStressTerm(model), dims=(1,2)))
KE_pressure = Field(Average(PressureRedistributionTerm(model), dims=(1,2)))
KE_buoyancy = Field(Average(BuoyancyProductionTerm(model), dims=(1,2)))

iKE = Field(Integral(KE))
iKE_tendency = Field(Integral(KE_tendency))
iKE_advection = Field(Integral(KE_advection))
#iKE_adv = Field(Integral(KE_adv))
iKE_stress = Field(Integral(KE_stress))
iKE_pressure = Field(Integral(KE_pressure))
iKE_buoyancy = Field(Integral(KE_buoyancy))

outputs = (; KE, KE_tendency, KE_advection, KE_stress, KE_pressure, KE_buoyancy,
           iKE, iKE_tendency, iKE_advection, KE_stress, iKE_pressure, iKE_buoyancy)

output_interval = 5minutes
simulation.output_writers[:profile] =
     NetCDFOutputWriter(model, outputs,
                     schedule = TimeInterval(output_interval),
                     filename = "mwe_KE_budget.nc",
                     overwrite_existing = true)
run!(simulation)

Code for plotting in python:

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

ds = xr.open_dataset('mwe_KE_budget.nc')
ds.close()
ds['timeTf'] = ds.time/np.timedelta64(int(np.around(2*np.pi/1e-4)), 's')

### Time series
plt.close()
plt.figure(figsize=(8,4), constrained_layout=True)

deltaz = np.diff(ds.zF)
ds['IKE_tendency'] = (ds.KE_tendency*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_advection'] = (ds.KE_advection*deltaz).isel(zC=slice(None,-1)).sum('zC')
#ds['IKE_adv'] = (ds.KE_adv*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_pressure'] = (ds.KE_pressure*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_buoyancy'] = (ds.KE_buoyancy*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_stress'] = (ds.KE_stress*deltaz).isel(zC=slice(None,-1)).sum('zC')

plt.plot(ds.timeTf, ds.IKE_tendency, 'k');
plt.plot(ds.timeTf, -ds.IKE_advection, 'C1');
plt.plot(ds.timeTf, -ds.IKE_pressure, 'C2');
plt.plot(ds.timeTf, ds.IKE_buoyancy, 'C3');
plt.plot(ds.timeTf, -ds.IKE_stress, 'C4');
plt.plot(ds.timeTf, (ds.IKE_tendency + ds.IKE_advection + ds.IKE_pressure - ds.IKE_buoyancy + ds.IKE_stress), ':C5', lw=1)
plt.legend(['tendency term', 'advection term', 'pressure term', 'buoyancy term', 'stress term', 'residual'],
           frameon=False, loc='lower right', ncol=3)
#plt.plot(ds.timeTf, -ds.IKE_adv, ':b', lw=1);
plt.xlabel(r'$T_{inertial}$')
plt.ylabel('Integrated KE budgets')
plt.ylim(-3.6e-5, 4.6e-5)
plt.ticklabel_format(axis='y', scilimits=(0,0))
plt.xlim(0,5.5)

ax1 = plt.gca().twinx()
color = 'C0'
ax1.set_ylabel(r'$\int k dV$', color=color)
ax1.plot(ds.timeTf, ds.iKE, color=color)
ax1.set_yticks([0, 1, 2, 3, 4])
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_ylim(-3.6, 4.6);
# plt.savefig('mwe_iKE_budgets', dpi=250)

### Profiles
plt.close()
plt.figure(figsize=(7.3,5), constrained_layout=True)

plt.subplot(121)
time_interval = (ds.timeTf >= 4) & (ds.timeTf < 4.5)
ke = ds.KE.where(time_interval, drop=True).isel(time=[0,-1]).T
dkedt = ke.diff('time')/ke.time.diff('time').dt.seconds
plt.plot(ds.KE_tendency.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'k')
plt.plot(dkedt, ds.zC, '--k')
plt.plot(-ds.KE_advection.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C1')
plt.plot(-ds.KE_pressure.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C2')
plt.plot(ds.KE_buoyancy.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C3')
plt.plot(-ds.KE_stress.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C4')
plt.plot((ds.KE_tendency + ds.KE_advection + ds.KE_pressure - ds.KE_buoyancy +
         ds.KE_stress).where(time_interval).mean('time')[:-1], ds.zC[:-1], ':C5', lw=1)
#plt.plot(-ds.KE_adv.where(time_interval).mean('time')[:-1], ds.zC[:-1], ':b', lw=1)
plt.text(0.02,0.42, r'$T_{inertial} \in [4,4.5)$', fontsize=12, transform=plt.gca().transAxes)
plt.grid('on', ls='--', lw=0.4)
plt.xlim(-3e-6, 3e-6)
plt.ylim(-40,0)
plt.xlabel('Time averaged KE budget term')

ax1 = plt.gca().twiny()
color = 'C0'
ax1.set_xlabel(r'$k$', color=color)
ax1.plot(ds.KE.where(time_interval).mean('time'), ds.zC, color=color)
ax1.set_xticks(np.array([0, 0.5, 1, 1.5])*1e-2)
ax1.tick_params(axis='x', labelcolor=color)
plt.ticklabel_format(axis='x', scilimits=(0,0))
ax1.set_xlim(-0.019, 0.019)

plt.subplot(122)
time_interval = (ds.timeTf >= 4.5) & (ds.timeTf < 5)
ke = ds.KE.where(time_interval, drop=True).isel(time=[0,-1]).T
dkedt = ke.diff('time')/ke.time.diff('time').dt.seconds
plt.plot(ds.KE_tendency.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'k')
plt.plot(dkedt, ds.zC, '--k')
plt.plot(-ds.KE_advection.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C1')
plt.plot(-ds.KE_pressure.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C2')
plt.plot(ds.KE_buoyancy.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C3')
plt.plot(-ds.KE_stress.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C4')
plt.plot((ds.KE_tendency + ds.KE_advection + ds.KE_pressure - ds.KE_buoyancy +
         ds.KE_stress).where(time_interval).mean('time')[:-1], ds.zC[:-1], ':C5', lw=1)
plt.legend(['tendency term', r'$\Delta k/\Delta t$', 'advection term', 'pressure term', 'buoyancy term',
            'stress term', 'residual'], frameon=False, loc='lower left', borderpad=0.1)
#plt.plot(-ds.KE_adv.where(time_interval).mean('time')[:-1], ds.zC[:-1], ':b', lw=1)
plt.text(0.02,0.42, r'$T_{inertial} \in [4.5,5)$', fontsize=12, transform=plt.gca().transAxes)
plt.gca().set_yticklabels([])
plt.grid('on', ls='--', lw=0.4)
plt.xlim(-3e-6, 3e-6)
plt.ylim(-40,0)
plt.xlabel('Time averaged KE budget term')

ax2 = plt.gca().twiny()
color = 'C0'
ax2.set_xlabel(r'$k$', color=color)
ax2.plot(ds.KE.where(time_interval).mean('time'), ds.zC, color=color)
ax2.set_xticks(np.array([0, 0.5, 1, 1.5])*1e-2)
ax2.tick_params(axis='x', labelcolor=color)
plt.ticklabel_format(axis='x', scilimits=(0,0))
ax2.set_xlim(-0.019, 0.019);
# plt.savefig('mwe_KE_budgets', dpi=250)

[glwagner]

https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L161C1-L166C4

This is not discretely correct. The discrete kinetic energy equation is obtained by multiplying the momentum equation by

$$ \frac{u^{n} + u^{n+1}}{2} $$

To understand this consider the time-discretized "momentum" equation (considering only $u$ for simplicity)

$$ u^{n+1} - u^{n} = \Delta t \int_{t_n}^{t_{n+1}} G_u \mathrm{d} t $$

If you multiply the above equation by $u^n$ you're not going to get a kinetic energy equation. You have to multiply it by $(u^{n} + u^{n+1}) / 2$, which yields

$$ \frac{1}{2} \left ( u^{n+1} \right )^2 - \frac{1}{2} \left ( u^n \right )^2 = \frac{1}{2} (u^{n} + u^{n+1}) \Delta t \int_{t_n}^{t_{n+1}} G_u \mathrm{d} t $$

The formulas in the code don't incorporate the velocity field from two time-steps so they just can't be correct. It also doesn't look like it's consistent with the time-stepping scheme, which involves tendencies from both $n$ and time-step $n-1$... ?

What you can try to show to get around this is that the budget doesn't change when you change the time-step. This would indicate that the error you are making wrt temporal discretization is negligible. But it's definitely not going to be negligible in general. On the contrary, we usually use time-steps close to the stability boundary.

[glwagner]

@glwagner can this be related to the halos not being filled? I don't really understand very well how halo filling works in the code, but I can't think of anything else.

The halos are filled within update_state!:

https://github.com/CliMA/Oceananigans.jl/blob/ed777806eea34ef02220e80676b3282f1f5a7f3f/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl#L30

which is called just before completing a time-step or substep.

Therefore whenever diagnostics are performed the halos for all prognostic variables should be filled.

[zhihua-zheng]

• You observe qualitative differences when you change advection scheme. For example with the CenteredSecondOrder() the main balance is between the advection term creating KE and the stress term removing it. With WENO() the balance seems similar, but the sign of the contributions is flipped, no?

With WENO() it is more like the time-averaged contribution from AdvectionTerm is ~ 0, while KineticEnergyStressTerm is positive, reflecting the wind work wins over dissipation in the mean? In this case, dissipation is indeed much smaller than the one with CenteredSecondOrder().

Ideally just based on the KE equation, I would only expect KineticEnergyStressTerm as the sole driver for KE tendency. In other words, the partition between advection and stress term does not seem right to me in either case.

[jhdong2016]

Hi Zhihua,

I am also doing some energy budget analysis and also have trouble in closing the budget. Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.

[zhihua-zheng]

Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.

Haven't tried a purely buoyancy forced case for KE budget yet, but it's on my list. I have tried that for a TKE budget though, the dissipation is similarly smaller with WENO().

[jhdong2016]

Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.

Haven't tried a purely buoyancy forced case for KE budget yet, but it's on my list. I have tried that for a TKE budget though, the dissipation is similarly smaller with WENO().

I am Jihai Dong. Jacob and Baylor told me about you. Do you have wechat? I think it will be convenient for us to discuss. You can send it to my email, jihai_dong@nuist.edu.cn.