Seeking for help in making the TKE budget equation closed

[zwei961120]

Hello everyone,

Hope you all are doing well! I am writing to seek your help regarding a problem I encountered while calculating the terms of TKE budget equation and making the equation closed. For a simple pure convection case, the TKE budget equation always fails to close (although I don't expect it to close perfectly, there is always a relatively large residual term). I am not sure if it is due to the model configuration (maybe the time stepper or advection scheme issues?) or some calculation errors in the TKE budget terms. Therefore, I would appreciate it if you could take a look at my code and point out any possible issues. I also hope that my code can clearly show the potential problems.

The TKE budget equation that my calculation follows is shown as below (note that p is the total kinematic pressure and there's no need for ρ₀ in code, the calculation of some terms uses Oceanostics.jl):

using Print
using JLD2
using Oceananigans
using Oceananigans.Units: minute, minutes, hour, hours
using Oceanostics: IsotropicKineticEnergyDissipationRate, TurbulentKineticEnergy, ZShearProductionRate
using Oceananigans.BuoyancyModels: g_Earth

# ## The grid
const H = 120
grid = RectilinearGrid(GPU();
                       size = (128, 128, 288),
                       halo = (3, 3, 3),
                          x = (0, 100),
                          y = (0, 100),
                          z = (-H, 0))
# Buoyancy
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4))


# ## Boundary conditions

# surface boundary conditions: surface cooling
Qʰ = 200.0  # W m⁻², surface _heat_ flux
cᴾ = 3991.0 # J K⁻¹ kg⁻¹, typical heat capacity for seawater
Qᵀ = Qʰ / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux

T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ))

# ## Coriolis
f₀ = 1e-4
coriolis = FPlane(f=f₀)

# ## Background current                                                                                                                                                                             
U∞ = 0.25 # m s⁻¹


# ## Set up sponge layer
damping_rate = 2/hour
bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)
uvw_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target = 0)

model = NonhydrostaticModel(; grid,
                advection = WENO(grid; order=9),
              timestepper = :RungeKutta3,
                  tracers = (:T, :S),
                  coriolis = coriolis,
                 buoyancy = buoyancy,
                  closure = SmagorinskyLilly(),
              boundary_conditions = (;T = T_bcs),
                  forcing = (;u = uvw_sponge, v = uvw_sponge, w = uvw_sponge))
                  
# ## Set up initial conditions
S₀ = 35.0
T₀ = 20.0
dTdz = 0.1 # K m⁻¹
T0(z) = T₀ + dTdz * z
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

## Horizontal velocity initial condition: geostrophic flow with random noise superposed
uᵢ(x, y, z) = U∞ + 1e-5 * Ξ(z)

## Vertical velocity initial condition: random noise scaled by the friction velocity.
wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)

Tᵢ(x, y, z) = T0(z) + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=S₀)

# ## Setting up a simulation
simulation = Simulation(model, Δt=10.0, stop_time=240hours)

wizard = TimeStepWizard(cfl=0.5, max_change=1.1, max_Δt=1minute)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

## Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s, advective CFL: %.2e, diffusive CFL: %.2e\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w),
                                prettytime(sim.run_wall_time),
                                AdvectiveCFL(sim.Δt)(sim.model),
                                DiffusiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(10))

### Calculate TKE budget equation terms
νₑ = model.diffusivity_fields.νₑ
u = @at (Center, Center, Center) model.velocities.u
v = @at (Center, Center, Center) model.velocities.v
w = @at (Center, Center, Center) model.velocities.w

U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))
W = Field(Average(w, dims=(1, 2)))
up =  @at (Center, Center, Center) Field(u - U)
vp =  @at (Center, Center, Center) Field(v - V)
wp =  @at (Center, Center, Center) Field(w - W)

# Calculate TKE using Oceanostics
k = TurbulentKineticEnergy(model, u, v, w; U, V, W)
k_mean = Average(k, dims = (1, 2))

# Calculate TKE transport
ue = Field(Average(wp*k, dims=(1, 2)))
ue_mean = @at (Center, Center, Center) ∂z(ue)
ue_mean = Average(ue_mean, dims = (1, 2))

# Calculate buoyancy term 
b = BuoyancyField(model)
B = Field(Average(b, dims=(1, 2)))

wp = @at (Center, Center, Center) Field(w-W)
bp = @at (Center, Center, Center) (b-B)
b_mean = Average(wp*bp, dims = (1, 2))

# Calculate shear production term using Oceanostics
zshear = ZShearProductionRate(model, u, v, w, U, V, W)
s_mean = Average(zshear, dims = (1, 2))

# Calculate Z-direction pressure redistribution
p = model.pressures.pHY′ + model.pressures.pNHS
P = Field(Average(p, dims=(1, 2)))
pp = @at (Center, Center, Center) (p-P)
pp = Field(Average(pp*wp, dims=(1, 2)))
pp = @at (Center, Center, Center) ∂z(pp)
p_mean = Average(pp, dims = (1, 2))

# Calculate dissipation term using Oceanostics
diss = Field(Average(IsotropicViscousDissipationRate(model; U, V, W), dims = (1, 2))) # Follow the Yan et al., 2022

# Calculate SGS transport term d(uᵢ'τᵢ₃')/dz

xu =  @at (Center, Center, Center) Field(∂z(u))
yv =  @at (Center, Center, Center) Field(∂y(v))
zw =  @at (Center, Center, Center) Field(∂z(w))
yu =  @at (Center, Center, Center) Field(∂y(u))
xv =  @at (Center, Center, Center) Field(∂x(v))
xw =  @at (Center, Center, Center) Field(∂x(w))
zu =  @at (Center, Center, Center) Field(∂z(u))
zv =  @at (Center, Center, Center) Field(∂z(v))
yw =  @at (Center, Center, Center) Field(∂y(w))

ν = @at (Center, Center, Center) Field(model.diffusivity_fields.νₑ)

tau13 = Field(2*ν*(zu+xw)/2) 
tau23 = Field(2*ν*(yw+zv)/2)
tau33 = Field(2*ν*(zw+zw)/2)

t_new = Field(tau13*up + tau23*vp + tau33*wp)

sgs_trans_operataion_new = Field(Average(t_new, dims = (1, 2)))

sgs_trans = @at (Center, Center, Center) ∂z(sgs_trans_operataion_new)

sgs_trans =   Field(Average(sgs_trans, dims = (1,2)))

output_interval = 1minutes
simulation.output_writers[:averages] =
    NetCDFOutputWriter(model, (; k_mean, b_mean, s_mean, p_mean, ue_mean, sgs_trans, diss),
                     schedule = AveragedTimeInterval(5minute, window=5minute),
                     filename = "tke_buoyancy_weno_9_new.nc",
                     overwrite_existing = true)

run!(simulation)

And the result shows like (averaged over an inertial oscillation period in quasi steady state):

As you can see, there is a relatively large residual term (the sum of the right hand side terms of TKE budget equation) which can not be balanced by the relatively small tendency term (the left hand side term of TKE budget equation). So, do you have any ideas or suggestions about this issue? Thanks in advance!

[glwagner]

Can you remake the plot showing fewer terms? I can't figure out which line is the residual.

I'd also test dependence on resolution, and choice of closure (AnisotropicMinimumDissipation vs SmagorinskyLilly) and advection scheme (WENO(order=9) vs Centered(order=2) for example).

Some people that might be able to help: @tomchor @qingli411 @simone-silvestri

[simone-silvestri]

Also, you have to ensure that

1. You are in a statistically steady state
2. You are averaging for long enough to remove the time derivative

I would suggest a couple of changes

• you can try increasing the average time (maybe 5 minutes is not long enough).
• calculate the temporal average of the material derivative of the TKE, and make sure that the "missing budget" is not included in that physical term.
• (more of a good practice point) you define a fluctuation as the difference between your instantaneous value and the spatial average only. Then you perform a temporal average of the correlation terms. If you are in a statistically steady state, it is more natural to define a fluctuation as the difference between the instantaneous value and your time and space average, which corresponds to the temporal averaging window that you are choosing. In that way, you are sure that you are not removing temporal variations of your spatial mean in the TKE definition

[glwagner]

But here the tendency is included in the budget. Doesn't that mean we can define the average purely as a spatial average? I think it'd simpler to avoid time averaging.

[tomchor]

Tying up to @simone-silvestri's last point. It might be useful to start tackling this by making sure the kinetic energy budget (not the turbulent kinetic energy budget; i.e. with no Reynolds decomposition) closes first. It's much simpler and doesn't rely on any definition of average/background. That way you can exclude possible sources of error. For example, if the KE budget closes, then you know the issue is linked to your Reynolds decomposition.