Large-Eddy Simulation Closures

When the grid resolves the energy-containing turbulent eddies, the subgrid-scale flux can be carried by a local eddy-viscosity closure rather than by PROPHET. ClimaAtmos provides three such closures: Smagorinsky–Lilly [10, 68, 69], anisotropic minimum dissipation [67, 70], and a constant horizontal diffusivity. They are used in box and plane configurations — large-eddy simulation, cloud-resolving radiative-convective equilibrium, density-current tests — and they act on the grid-mean variables only.

The horizontal tendencies are explicit. The vertical Smagorinsky–Lilly tendency follows implicit_diffusion, like the prescribed vertical diffusion; the vertical AMD tendency currently is explicit.

The eddy-viscosity form

Each closure supplies an eddy viscosity $\nu_t$ and an eddy diffusivity $D$. The momentum flux is proportional to the resolved strain rate,

\[\boldsymbol{\mathcal{T}} = -2 \nu_t \boldsymbol{S} , \qquad \boldsymbol{S} = \tfrac{1}{2} \left[ \nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^{T} \right] ,\]

and the scalar flux is down-gradient,

\[\boldsymbol{\mathcal{F}}_\psi = -\rho D \nabla \psi .\]

The strain rate is assembled on both centers and faces from the full three-dimensional velocity gradient, and its norm is $|\boldsymbol{S}| = \sqrt{2 \boldsymbol{S} : \boldsymbol{S}}$, evaluated after projection onto whichever axes the closure uses.

The energy flux is split, as in the other diffusive terms of the model (the prescribed vertical diffusion, the PROPHET fluxes, hyperdiffusion, and the viscous sponge), into a dry-static-energy part and the enthalpy carried by the diffusing water,

\[\boldsymbol{\mathcal{F}}_h = -\rho D \left( \nabla s_d + (h_{\mathrm{eff}} + \Phi) \nabla q_t^{\mathrm{eff}} \right) ,\]

with $h_{\mathrm{eff}}$ the mass-weighted enthalpy of the suspended water and $q_t^{\mathrm{eff}}$ the total water less rain and snow. The split makes the flux invariant to the reference temperature $T_0$ and ties the energy the water flux carries to the water flux itself (see Thermodynamics and the Working Fluid). Anisotropic minimum dissipation forms a diffusivity per scalar, so the two parts of the split use the diffusivity of the gradient each acts on. Smagorinsky–Lilly and the constant horizontal diffusion use a single scalar diffusivity for both parts of the split and for all grid-scale tracers: Smagorinsky–Lilly divides its eddy viscosity by the turbulent Prandtl number (below), and the constant horizontal diffusion prescribes the diffusivity outright. Because total water is part of the working fluid, its diffusive tendency is also added to the density tendency. All water species, precipitation included, diffuse with the same diffusivity and keep their own tendencies; the enthalpy that rain and snow carry is left out of the energy flux, since $q_t^{\mathrm{eff}}$ excludes them.

Smagorinsky–Lilly

The eddy viscosity is the product of a squared mixing length and the strain-rate norm, and the diffusivity follows from a constant turbulent Prandtl number,

\[\nu_t = L^2 |\boldsymbol{S}| , \qquad D = \frac{\nu_t}{Pr_t} .\]

The mixing length is the filter scale times the Smagorinsky coefficient $c_s$, and it is reduced in stable stratification by the Lilly correction

\[f_b = \left( 1 - \frac{Ri}{Pr_t} \right)^{1/4} , \qquad Ri = \frac{N^2}{|\boldsymbol{S}_v|^2} ,\]

with $f_b = 1$ where $Ri \le 0$, and with $N^2$ computed from the vertical gradient of virtual potential temperature. The Richardson number always uses the vertical strain-rate norm, whatever axes the closure acts on. The correction vanishes at $Ri = Pr_t$, which switches the closure off in strongly stable layers.

The smagorinsky_lilly configuration key selects which axes the closure acts on, and that choice also sets the filter scale:

Here, $\Delta x$ is the horizontal node spacing and $\Delta z$ the local layer thickness. In the coupled "UVW" case, one cube-root filter width and one strain-rate norm are used for both directions; in the other three cases the horizontal and vertical viscosities are formed from separate lengths and from strain-rate norms projected onto the corresponding axes. The stratification correction enters the vertical length only, except under "UVW" where there is a single length to correct.

Anisotropic minimum dissipation

The minimum-dissipation closure chooses the smallest eddy viscosity that keeps the resolved kinetic energy from growing at the grid scale, which makes it sensitive to grid anisotropy without an explicit filter-width prescription. Writing $\hat{\partial}_i = \Delta_i \partial_i$ for the grid-scaled derivative, the viscosity is

\[\nu_t = \max \left[ 0, \, -c_{amd} \frac{(\hat{\partial}_k u_i)(\hat{\partial}_k u_j) S_{ij}}{(\partial_l u_m)(\partial_l u_m)} \right] ,\]

and each scalar gets its own diffusivity,

\[D_\psi = \max \left[ 0, \, -c_{amd} \frac{(\hat{\partial}_k u_i)(\hat{\partial}_k \psi)(\partial_i \psi)}{(\partial_l \psi)(\partial_l \psi)} \right] ,\]

in which the scalar derivatives $\partial_i \psi$ are the horizontal ones in the horizontal tendency and the vertical one in the vertical tendency, so each scalar has a separate horizontal and vertical diffusivity. The velocity gradients in the viscosity are the full three-dimensional ones. The numerators use the grid-scaled derivatives and the denominators do not, so the ratio has the dimensions of a diffusivity. Clipping at zero removes backscatter. There is no Prandtl number: momentum and each scalar are treated separately.

The closure has a single switch: amd_les: true applies both the horizontal and the vertical tendency.

The coefficient $c_{amd}$ is set by the c_amd configuration key, default 0.29. It is the one coefficient among these closures that is read from the YAML configuration rather than from ClimaParams; the default configuration has a TODO to move it to parameters.toml.

Constant horizontal diffusion

The simplest option applies a fixed horizontal diffusivity to the scalars — the split energy flux and the grid-scale tracers — and leaves momentum alone. It is switched on with constant_horizontal_diffusion: true, and the diffusivity comes from the D_horizontal_diffusion parameter, default $5 \times 10^5$ m² s⁻¹. No configuration in the repository currently enables it.

Boundary conditions and surface fluxes

The vertical scalar fluxes vanish at the surface and at the model top: the divergence operator applied to energy and the tracers imposes a zero-flux boundary condition at both ends. Momentum is handled differently. The face strain rate is built with a zero vertical velocity gradient at the boundaries, but the horizontal gradients still contribute to the boundary stress, and the horizontal-momentum divergence applies that stress as computed; the vertical-momentum operator imposes zero divergence at the boundaries rather than zero flux. Either way the closure supplies no wall-layer stress of its own. Surface fluxes do not enter through these boundary conditions. They are added as a separate explicit tendency confined to the surface-adjacent cell, so surface drag and the surface enthalpy and moisture fluxes are applied whether or not an LES closure is active; see Surface Conditions.

Relation to the other diffusion schemes

Prescribed vertical diffusion (vert_diff) is a boundary-layer closure with a prescribed diffusivity profile, and it differs in two ways beyond the diffusivity itself: it interpolates the diffusivity to faces as a harmonic rather than an arithmetic mean, and it diffuses the effective total water $q_t - q_{rai} - q_{sno}$ and redistributes the resulting mass among the cloud species. See Diffusion for both closures and the water-species bookkeeping.

PROPHET. The LES closures and the PROPHET horizontal diffusive flux both apply horizontal subgrid-scale diffusion to the same fields, so combining edmfx_sgs_horizontal_diffusive_flux with smagorinsky_lilly or amd_les is rejected at configuration time; see PROPHET Horizontal Diffusion.

Hyperdiffusion is a numerical filter rather than a turbulence closure, and nothing in the code couples the two: whether it runs alongside an LES closure is left to the configuration. The Smagorinsky and AMD configurations set hyperdiff: ~, since the closure already damps the grid scale. See Hyperdiffusion.

Implicit vertical diffusion. With implicit_diffusion: true, the vertical Smagorinsky–Lilly tendency moves into the implicit tendency, and the eddy viscosity is refreshed on every Newton iterate so that it follows the current state. The Jacobian linearizes the enthalpy, tracer, and horizontal-momentum flux divergences with the eddy viscosity frozen; the $u_3$ term, the horizontal-gradient part of $\boldsymbol{\mathcal{T}}$, and the per-species condensate diagonals have no Jacobian contribution, which slows Newton convergence but leaves the tendency exact. A vertically acting Smagorinsky closure satisfies the configuration check for implicit diffusion on its own; the AMD closure has no Jacobian block and stays explicit. See Implicit Solver.

Configuration and diagnostics

Worked examples: les_isdac_box.yml and box_density_current_test.yml use "UVW", rcemipii_box_CRM_1M.yml uses "UV_W", and plane_density_current_test_amd.yml and baroclinic_wave_equil_amd.yml use AMD.

The Smagorinsky closure exposes its diffusivities and strain-rate norms as the Dh_smag, Dv_smag, strainh_smag, and strainv_smag diagnostics. AMD exposes none.

Symbols

Where this is implemented

See Configuration Options for the keys named here.