feat(edmfx): add horizontal component to the EDMFX SGS diffusive flux

The EDMFX environment SGS diffusive flux was vertical-only. Add its horizontal
component, enabled by the opt-in option `edmfx_sgs_horizontal_diffusive_flux`
(default false), for high-resolution box runs. On single-column geometry the
tendency returns immediately, so column results are unchanged by the option.

The scalar fluxes mirror the vertical EDMFX diffusive flux: the same variables
(total enthalpy, total specific humidity with its air-mass coupling, the
environment SGS tracers scaled by `α_vert_diff_tracer`, and TKE) and the same
TKE-based eddy diffusivity, with horizontal rather than vertical operators.
The momentum tendency is the horizontal weak divergence of the SGS stress
`τ = -2 K_u S`, with the full three-dimensional strain rate, on cell centers
for the horizontal wind and on cell faces for the vertical wind, mirroring the
Smagorinsky-Lilly stress split by divergence direction. The corresponding TKE
shear production from horizontal gradients, `2 ρ K_u S_h : S_h` with `S_h` the
strain rate built from horizontal gradients only, is positive definite; the
production from vertical gradients and its stencil are unchanged. All terms
are always explicit and never enter the implicit Jacobian; the vertical flux
stays implicit when `implicit_diffusion` is set.

The anisotropy follows the Smagorinsky-Lilly length-scale split: the eddy
diffusivity is isotropic in TKE, so only the mixing length differs by direction
through its grid-scale limiter (the cell thickness for the vertical flux, the
horizontal node spacing for the horizontal flux). `mixing_length_lopez_gomez_2020`
generalizes its limiter argument (`ᶜdz` to `grid_scale`); `ᶜmixing_length`
gains `grid_scale` and `buoyancy_gradient` keywords (defaults: the vertical
cell thickness and the centered `ᶜlinear_buoygrad`), and the horizontal call
sites pass both explicitly.

Add diagnostics for the horizontal closure fields (`lmixh`, `edth`, `evuh`).
Combining the option with Smagorinsky-Lilly or anisotropic minimum dissipation,
which already apply horizontal SGS diffusion, is rejected at model construction.

Add tests: horizontally uniform states and single-column geometry produce no
tendency, the moist air mass tendency matches the `ρq_tot` tendency bitwise,
the grid-scale limit is directional and attains a binding cap, the TKE
production is positive definite under uniform TKE, and incompatible closure
combinations are rejected.

Add docs/src/edmf_horizontal_diffusion.md.

Author: haakon-e

NEWS.md

@@ -4,6 +4,7 @@ ClimaAtmos.jl Release Notes
 main
 ----
 
+- [#4657](https://github.com/CliMA/ClimaAtmos.jl/pull/4657) ![][badge-✨feature/enhancement] Add a horizontal component to the EDMFX SGS diffusive flux, enabled by the opt-in `edmfx_sgs_horizontal_diffusive_flux` config option (default `false`), for high-resolution box configurations: scalar and TKE fluxes, the momentum stress `τ = -2 K_u S` with the full strain rate, the corresponding TKE shear production from horizontal gradients, and diagnostics `lmixh`, `edth`, and `evuh`.
 - [#4648](https://github.com/CliMA/ClimaAtmos.jl/pull/4648) ![][badge-🔥behavioralΔ] Carry the energy of sedimenting water at each subdomain's own value under `PrognosticEDMFX`. Each sedimenting species transports its specific internal, potential, and kinetic energy `e_int(T) + Φ + Kin(w, u)`; the grid-mean sedimentation energy flux is set equal to the sum of the updraft and environment fluxes, each the subdomain sedimentation mass flux times its specific energy (with `Φ` common to all subdomains and the environment mass flux taken as the grid-mean residual `ρqw - ρaʲqʲwʲ`). The updraft kinetic energy uses the updraft velocities; the environment kinetic energy uses the grid-mean terminal velocity, since the environment terminal velocity is not stored separately.
 - [#4653](https://github.com/CliMA/ClimaAtmos.jl/pull/4653) ![][badge-🔥behavioralΔ] Improve the manual implicit Jacobian (the tendencies are unchanged; behavior changes only through the single-Newton-iteration implicit solve).
   - Document that the existing height-form pressure-gradient-force Jacobian blocks are the exact linearization of the Exner-form tendency `cp_d θ_v ∇Π` (stencils unchanged): the thermal-buoyancy and pressure-buoyancy terms cancel via the equation-of-state identity `δρ/ρ = (1-κ_d) δp/p - δθ_v/θ_v`, leaving the acoustic operator `∇δp/ρ` in every column perturbed at fixed density plus a single buoyancy term in the `ρ` column, so sound and gravity waves are both treated fully implicitly.

config/default_configs/default_config.yml

@@ -455,6 +455,9 @@ edmfx_sgs_mass_flux:
 edmfx_sgs_diffusive_flux:
   help: "If set to true, it switches on EDMFX SGS diffusive flux.  [`true`, `false` (default)]"
   value: false
+edmfx_sgs_horizontal_diffusive_flux:
+  help: "If set to true, it switches on the horizontal component of the EDMFX SGS diffusive flux. Intended for high-resolution box runs.  [`true`, `false` (default)]"
+  value: false
 updraft_number:
   help: "Sets the number of updrafts for the EDMF scheme"
   value: 1

docs/make.jl

@@ -56,6 +56,7 @@ makedocs(;
             "Physics & Parameterizations" => [
                 "Microphysics" => "microphysics.md",
                 "EDMF: Prognostic Equations" => "edmf_equations.md",
+                "EDMF: Horizontal Diffusion" => "edmf_horizontal_diffusion.md",
                 "Gravity Wave Drag" => "gravity_wave.md",
                 "Ocean Surface Albedo" => "surface_albedo.md",
                 "Topography Representation" => "topography.md",

docs/src/edmf_horizontal_diffusion.md

@@ -0,0 +1,94 @@
+# Horizontal EDMF diffusion
+
+The EDMFX environment sub-grid scale (SGS) diffusive flux is, by default, vertical only: it parameterizes turbulent transport along the column with an eddy-diffusivity closure and applies the resulting tendency to the grid mean. The horizontal component adds the analogous horizontal down-gradient flux. It is enabled with the configuration option `edmfx_sgs_horizontal_diffusive_flux` (default `false`).
+
+## Formulation
+
+The horizontal flux uses the same TKE-based eddy diffusivity as the vertical flux,
+
+```math
+K = c_m \, l \, \sqrt{e} \,/\, \mathrm{Pr}_t ,
+```
+
+where ``l`` is the environment mixing length, ``e`` is the environment turbulent kinetic energy (TKE), ``c_m`` is the eddy-diffusivity coefficient, and ``\mathrm{Pr}_t`` the turbulent Prandtl number. The grid-mean tendencies are, on cell centers,
+
+```math
+\partial_t (\rho e_\text{tot}) \mathrel{+}= \nabla_h \cdot \left[ \rho \, K_h \left( \nabla_h s_\text{d} + \sum_{\mu} (h_\mu + \Phi) \, \nabla_h q_\mu \right) \right],
+```
+
+```math
+\partial_t (\rho \chi) \mathrel{+}= \nabla_h \cdot (\rho \, K_h \, \nabla_h \chi),
+\qquad
+\partial_t (\rho e) \mathrel{+}= \nabla_h \cdot (\rho \, K_u \, \nabla_h e),
+```
+
+where ``\nabla_h`` and ``\nabla_h\cdot`` are the horizontal gradient and weak divergence, ``s_\text{d}`` is the dry static energy, ``h_\mu`` and ``q_\mu`` are the specific enthalpy and specific humidity of the water species ``\mu \in \{\text{vap}, \text{liq}, \text{ice}\}``, ``\Phi`` is the geopotential, ``K_u`` is the eddy viscosity and ``K_h = K_u/\mathrm{Pr}_t``. The total enthalpy uses the dry-static-energy + water-enthalpy decomposition of the vertical flux, which avoids the spurious enthalpy flux that diffusing ``h_\text{tot}`` directly would carry with the dry-air mass. The tracer ``\chi`` ranges over the total specific humidity and the environment SGS tracers — the microphysics species (cloud, precipitation, and, for two-moment schemes, number concentrations) and any passive tracers — i.e. the same set the vertical flux diffuses. Sedimenting microphysics species are scaled by the tracer diffusion factor ``\alpha`` (`α_vert_diff_tracer`); other tracers use the unscaled ``K_h``, matching the vertical flux. The total-specific-humidity flux additionally enters the air-mass tendency ``\partial_t \rho``; condensate and precipitation fluxes do not.
+
+The momentum tendency is the horizontal weak divergence of the subgrid-scale stress ``\tau = -2 K_u \mathcal{S}``, with ``\mathcal{S}`` the full three-dimensional strain rate of the grid-mean velocity,
+
+```math
+\partial_t u_h \mathrel{-}= \frac{1}{\rho} \nabla_h \cdot (\rho \, \tau),
+\qquad
+\partial_t u_3 \mathrel{-}= \frac{1}{\rho} \nabla_h \cdot (\rho \, \tau),
+```
+
+evaluated on cell centers for the horizontal wind and on cell faces for the vertical wind, mirroring the Smagorinsky-Lilly stress split by divergence direction; the vertical stress divergence is handled by the vertical diffusion pathway.
+
+The corresponding shear production of TKE uses the strain rate built from horizontal gradients only, ``\mathcal{S}_h``,
+
+```math
+\partial_t (\rho e) \mathrel{+}= 2 \rho \, K_u \, \mathcal{S}_h : \mathcal{S}_h ,
+```
+
+which is positive definite; the production from vertical gradients and its stencil are unchanged. This mirrors the decoupled Smagorinsky-Lilly split, in which each directional norm is built from that direction's gradients.
+
+The horizontal flux mirrors the vertical EDMFX diffusive flux: the same variables, the same tracer set, and the same scaling, with horizontal rather than vertical operators.
+
+## Anisotropic length scale
+
+The eddy diffusivity depends on direction only through the mixing length ``l``: the TKE ``e`` is isotropic. In the Lopez-Gomez et al. (2020) closure the mixing length is the smooth minimum of physical scales (wall distance, TKE production–dissipation balance, static stability), limited from above by a grid scale,
+
+```math
+l = \min(l_\text{phys}, \, \Delta).
+```
+
+Unlike the Smagorinsky-Lilly closure, whose length scale is purely geometric (``l = c_s \Delta``), the EDMF length is physically based; only its grid-scale limiter is geometric. The vertical and horizontal fluxes therefore share the same physical length but use different limiters — the cell thickness ``\Delta z`` for the vertical flux and the horizontal node spacing ``\Delta x`` for the horizontal flux:
+
+```math
+l_v = \min(l_\text{phys}, \, \Delta z),
+\qquad
+l_h = \min(l_\text{phys}, \, \Delta x).
+```
+
+This mirrors the Smagorinsky-Lilly horizontal/vertical split, in which the two length scales are ``c_s \Delta x`` and ``c_s \Delta z``. Two limits follow:
+
+  - where ``l_\text{phys} < \min(\Delta x, \Delta z)`` (for example small near-surface
+    eddies), ``l_h = l_v = l_\text{phys}`` and the diffusivity is isotropic;
+  - where ``l_\text{phys} > \max(\Delta x, \Delta z)``, the ratio ``l_h / l_v = \Delta x / \Delta z`` recovers the Smagorinsky grid anisotropy.
+
+## When to enable it
+
+Horizontal EDMF diffusion is intended for box configurations whose horizontal grid spacing is comparable to or smaller than the physical mixing length, so that the horizontal limiter ``\min(l_\text{phys}, \Delta x)`` binds. This is the high-resolution and gray-zone regime, where the horizontal grid is fine enough for an anisotropic SGS length scale to matter.
+
+At coarse horizontal resolution (for example global runs with ``\Delta x`` of tens of kilometers) the horizontal limiter rarely binds, the horizontal diffusivity reduces to the isotropic environment value, and the term is typically negligible next to resolved horizontal transport. The option is off by default for that reason.
+
+The horizontal term is a grid-mean closure applied independently of the vertical flux and of the diffusion mode: it is always explicit, since the horizontal operators are not part of the column implicit solve, while the vertical flux remains implicit when `implicit_diffusion` is set. The Smagorinsky-Lilly and anisotropic-minimum-dissipation closures already supply horizontal SGS diffusion of the same fields, so combining either with this option is rejected at model construction.
+
+Because the term is explicit, it adds a horizontal diffusive stability limit on the timestep, ``\Delta t \lesssim \Delta x^2 / (2 K_h)``. At the fine horizontal resolutions where this closure is intended, the timestep is in practice already set by the explicit horizontal acoustic limit ``\Delta t \lesssim \Delta x / c_s``, which is the more restrictive of the two.
+
+## Configuration
+
+```yaml
+edmfx_sgs_horizontal_diffusive_flux: true   # default: false
+```
+
+## Diagnostics
+
+Three diagnostic variables expose the horizontal closure fields:
+
+  - `lmixh`: the mixing length with the grid-scale limit set by the horizontal node
+    spacing, ``l_h`` [m];
+  - `edth`: the horizontal eddy diffusivity for scalars, ``K_h`` [m² s⁻¹];
+  - `evuh`: the horizontal eddy viscosity, ``K_u`` [m² s⁻¹].
+
+The vertical counterparts are `lmix`, `edt`, and `evu`.

src/config/model_getters.jl

@@ -664,6 +664,16 @@ function check_case_consistency(parsed_args)
         "Unknown `config = $(repr(config))`. Valid options are: $(join(valid_configs, ", "))."
     )
 
+    if parsed_args["edmfx_sgs_horizontal_diffusive_flux"] && (
+        !isnothing(parsed_args["smagorinsky_lilly"]) || parsed_args["amd_les"]
+    )
+        error(
+            "`edmfx_sgs_horizontal_diffusive_flux` cannot be combined with \
+             `smagorinsky_lilly` or `amd_les`, which already apply horizontal \
+             SGS diffusion to the same fields",
+        )
+    end
+
     # ISDAC consistency: when initial_condition is ISDAC, surface/rad/external
     # forcing must all be set to the matching ISDAC variants. Subsidence,
     # scm_coriolis, and ls_adv are owned by the setup, not the YAML schema.
@@ -769,6 +779,7 @@ function AtmosTurbconv(config::AtmosConfig, params, ::Type{FT}) where {FT}
         detr_model = get_detrainment_model(pa),
         sgs_mass_flux = pa["edmfx_sgs_mass_flux"],
         sgs_diffusive_flux = pa["edmfx_sgs_diffusive_flux"],
+        sgs_diffusive_flux_horizontal = pa["edmfx_sgs_horizontal_diffusive_flux"],
         nh_pressure = pa["edmfx_nh_pressure"],
         vertical_diffusion = pa["edmfx_vertical_diffusion"],
         filter = pa["edmfx_filter"],

src/diagnostics/edmfx_diagnostics.jl

@@ -851,6 +851,87 @@ add_diagnostic_variable!(short_name = "evu", units = "m^2 s^-1",
     compute = compute_evu,
 )
 
+###
+# Horizontal mixing length (3d)
+###
+compute_lmixh(state, cache, time) =
+    compute_lmixh(state, cache, time, cache.atmos.turbconv_model)
+compute_lmixh(_, _, _, turbconv_model) =
+    error_diagnostic_variable("lmixh", turbconv_model)
+
+function compute_lmixh(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
+    Δx = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(state.c)))
+    return ᶜmixing_length(
+        state, cache;
+        grid_scale = Δx, buoyancy_gradient = cache.precomputed.ᶜlinear_buoygrad,
+    )
+end
+
+add_diagnostic_variable!(short_name = "lmixh", units = "m",
+    long_name = "Environment Horizontal Mixing Length",
+    comments = "Mixing length with the grid-scale limit set by the horizontal \
+                node spacing rather than the vertical cell thickness",
+    compute = compute_lmixh,
+)
+
+###
+# Horizontal diffusivity of heat (3d)
+###
+compute_edth(state, cache, time) =
+    compute_edth(state, cache, time, cache.atmos.turbconv_model)
+compute_edth(_, _, _, turbconv_model) =
+    error_diagnostic_variable("edth", turbconv_model)
+
+function compute_edth(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
+    turbconv_params = CAP.turbconv_params(cache.params)
+    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = cache.precomputed
+    (; params) = cache
+
+    ᶜtke = @. lazy(specific(state.c.ρtke, state.c.ρ))
+    Δx = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(state.c)))
+    ᶜmixing_length_field = ᶜmixing_length(
+        state, cache; grid_scale = Δx, buoyancy_gradient = ᶜlinear_buoygrad,
+    )
+    ᶜK_u = @. lazy(eddy_viscosity(turbconv_params, ᶜtke, ᶜmixing_length_field))
+    ᶜprandtl_nvec =
+        @. lazy(turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm))
+    return @. lazy(eddy_diffusivity(ᶜK_u, ᶜprandtl_nvec))
+end
+
+add_diagnostic_variable!(short_name = "edth", units = "m^2 s^-1",
+    long_name = "Horizontal Eddy Diffusivity Coefficient for Temperature",
+    comments = "Horizontal diffusion coefficient for scalars in the EDMFX \
+                horizontal SGS diffusive flux",
+    compute = compute_edth,
+)
+
+###
+# Horizontal diffusivity of momentum (3d)
+###
+compute_evuh(state, cache, time) =
+    compute_evuh(state, cache, time, cache.atmos.turbconv_model)
+compute_evuh(_, _, _, turbconv_model) =
+    error_diagnostic_variable("evuh", turbconv_model)
+
+function compute_evuh(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
+    turbconv_params = CAP.turbconv_params(cache.params)
+    ᶜtke = @. lazy(specific(state.c.ρtke, state.c.ρ))
+    Δx = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(state.c)))
+    ᶜmixing_length_field = ᶜmixing_length(
+        state, cache;
+        grid_scale = Δx, buoyancy_gradient = cache.precomputed.ᶜlinear_buoygrad,
+    )
+    return @. lazy(eddy_viscosity(turbconv_params, ᶜtke, ᶜmixing_length_field))
+end
+
+add_diagnostic_variable!(short_name = "evuh", units = "m^2 s^-1",
+    long_name = "Horizontal Eddy Viscosity Coefficient for Momentum",
+    comments = "Horizontal eddy viscosity from the TKE-based closure with the \
+                mixing length limited by the horizontal node spacing; used to \
+                diffuse TKE in the EDMFX horizontal SGS diffusive flux",
+    compute = compute_evuh,
+)
+
 ###
 # Updraft passive gas tracer A (3d)
 ###

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -256,7 +256,7 @@ end
     ustar,
     ᶜz,
     z_sfc,
-    ᶜdz,
+    grid_scale,
     sfc_tke,
     ᶜN²_eff,
     ᶜtke,
@@ -273,7 +273,9 @@ where:
 - `ustar`: Friction velocity [m/s].
 - `ᶜz`: Cell center height [m].
 - `z_sfc`: Surface elevation [m].
-- `ᶜdz`: Cell vertical thickness [m].
+- `grid_scale`: Grid length scale used as the upper bound on the mixing length [m].
+  The vertical cell thickness for vertical diffusion, the horizontal node spacing for
+  horizontal diffusion.
 - `sfc_tke`: TKE near the surface (e.g., first cell center) [m^2/s^2].
 - `ᶜN²_eff`: Effective squared Brunt-Väisälä frequency [1/s^2].
 - `ᶜtke`: Turbulent kinetic energy at cell center [m^2/s^2].
@@ -299,7 +301,7 @@ function mixing_length_lopez_gomez_2020(
     ustar,
     ᶜz,
     z_sfc,
-    ᶜdz,
+    grid_scale,
     sfc_tke,
     ᶜN²_eff,
     ᶜtke,
@@ -400,7 +402,7 @@ function mixing_length_lopez_gomez_2020(
     l_limited_phys_wall = min(l_smin, l_z)
 
     # 2. Impose the grid-scale limit (TODO: replace by volumetric grid scale)
-    l_grid = ᶜdz   # TODO include costant rescaling factor
+    l_grid = grid_scale   # TODO include constant rescaling factor
     l_final = min(l_limited_phys_wall, l_grid)
 
     # Final check: guarantee that the mixing length is at least a small positive
@@ -421,20 +423,23 @@ end
 @inline get_mixing_length_field(ml::MixingLength, ::Val{:buoy}) = ml.buoy
 @inline get_mixing_length_field(ml::MixingLength, ::Val{:l_grid}) = ml.l_grid
 
-function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
+function ᶜmixing_length(
+    Y, p, property::Val{P} = Val{:master}();
+    grid_scale = Fields.Δz_field(axes(Y.c)),
+    buoyancy_gradient = p.precomputed.ᶜlinear_buoygrad,
+) where {P}
     (; params) = p
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
-    (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
+    (; ᶜstrain_rate_norm) = p.precomputed
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
-    ᶜdz = Fields.Δz_field(axes(Y.c))
 
     ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ))
     sfc_tke = Fields.level(ᶜtke, 1)
 
     ᶜprandtl_nvec = p.scratch.ᶜtemp_scalar_5
     @. ᶜprandtl_nvec =
-        turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm)
+        turbulent_prandtl_number(params, buoyancy_gradient, ᶜstrain_rate_norm)
 
     # Extract sub-parameters before the lazy broadcast to avoid capturing
     # the full ClimaAtmosParameters struct (~4 KiB) in GPU kernel parameters.
@@ -450,9 +455,9 @@ function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
             ustar,
             ᶜz,
             z_sfc,
-            ᶜdz,
+            grid_scale,
             sfc_tke,
-            ᶜlinear_buoygrad,
+            buoyancy_gradient,
             ᶜtke,
             obukhov_length,
             ᶜstrain_rate_norm,