Discontinuous Cs when using dynamic Smagorinsky model

[zwei961120]

Hi all!

I am using the dynamic Smagorinsky model, and I am noticing some unexpected behavior in the computed dynamic coefficient Cs.

Below is an example of what I observe (horizontally-averaged Cs and update every 1 iteration):

What happens is:

1. At first, Cs looks physically reasonable.
2. Then, suddenly a discontinuity appears.
3. After that moment, the lower part of the water column collapses into an almost depth-constant value of Cs. This constant region does not seem to be physically meaningful and appears instantly rather than gradually.

What I currently suspect:
This behavior might be related to the variable time stepping. Since the Lagrangian averaging procedure uses the time step to compute the running average, changes in Δt can affect how the averaging is performed. As turbulence develops, the CFL-based time step typically decreases, which may introduce inconsistencies in the Lagrangian averaging and lead to the discontinuous Cs behavior I observed. So when I fix the time step, the problem seems much less severe or does not appear, and the evolution of Cs shows like:

So, my question is:
Is this behavior known or expected for the dynamic Smagorinsky model? I could not find references in previous literature reporting this, so I wonder if it may indicate an implementation issue or a model limitation.

Here is the code:

grid = RectilinearGrid(GPU();
                       size = (480, 480, 160), 
                       halo = (3, 3, 3),
                       extent = (240, 240, 80))

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4))

dTdz = 0.01 # K m⁻¹

T_bcs = FieldBoundaryConditions(bottom = GradientBoundaryCondition(dTdz))

u₁₀ = 10   
cᴰ = 1.25e-3 
ρₐ = 1.225
ρₒ = 1026.0
Qᵘ =  ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀)


u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

Qˢ = 0.0
S_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Qˢ))


S₀ = 35.0
T₀ = 11.85
T0(z) = T₀ + dTdz * (z)


gaussian_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)

u_sponge = v_sponge = w_sponge = Relaxation(rate=1/hour, mask=gaussian_mask)

T_sponge = Relaxation(rate = 1/hour,
                      target = LinearTarget{:z}(intercept=T₀, gradient=dTdz),
                      mask = gaussian_mask)
S_sponge = Relaxation(rate = 1/hour,
                      target = S₀,
                      mask = gaussian_mask)



const amplitude = 0.6 
const wavenumber = 0.1 
const frequency = sqrt(9.8 * wavenumber) 


const Uˢ = frequency * amplitude^2 * wavenumber 


uˢ(x, y, z) = frequency * amplitude^2 * wavenumber * exp(2 * wavenumber * z)
∂z_uˢ(z, t) = frequency * amplitude^2 * wavenumber * 2 * wavenumber * exp(2 * wavenumber * z)

# ## Model instantiation

model = NonhydrostaticModel(
                advection = WENO(order=5),
              timestepper = :RungeKutta3,
                     grid = grid,
                  tracers = (:T, :S),
                 coriolis = FPlane(f=0),
                 buoyancy = buoyancy,
                 closure = Smagorinsky(coefficient=DynamicCoefficient(averaging=LagrangianAveraging(), schedule = IterationInterval(1))),
      boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs),
                stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                  forcing = (u=u_sponge, v=v_sponge, w=w_sponge, T=T_sponge, S=S_sponge))


rng = MersenneTwister(1256) 

Ξ(z) = randn(rng) * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

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

uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)

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

simulation = Simulation(model, Δt=10.0, stop_iteration=30000000)

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

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

𝒥ᴸᴹ = model.diffusivity_fields.𝒥ᴸᴹ

𝒥ᴹᴹ = model.diffusivity_fields.𝒥ᴹᴹ

cs2 = Field(𝒥ᴸᴹ) / Field(𝒥ᴹᴹ)
cs = Field(cs2)^(1/2)
cs_mean = Average(cs, dims=(1, 2))

Any guidance or clarification would be greatly appreciated! Thanks!