Skip to content

Commit e2bcdba

Browse files
committed
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.
1 parent 3060c01 commit e2bcdba

16 files changed

Lines changed: 652 additions & 3 deletions

File tree

NEWS.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ ClimaAtmos.jl Release Notes
33

44
main
55
----
6+
- [#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`.
67
- [#4699](https://github.com/CliMA/ClimaAtmos.jl/pull/4699) ![][badge-🔥behavioralΔ] Select the tracers that receive the `α_vert_diff_tracer` eddy-diffusivity scaling in the boundary-layer vertical diffusion from the shared `gs_sedimenting_tracer_candidates` list instead of a hardcoded tuple of species. The tracer diffusivity scaling is now consistent across the boundary-layer diffusion, the EDMFX SGS flux, the EDMFX updraft vertical diffusion, and the implicit Jacobian.
78

89
0.42.0

config/default_configs/default_config.yml

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -455,6 +455,9 @@ edmfx_sgs_mass_flux:
455455
edmfx_sgs_diffusive_flux:
456456
help: "If set to true, it switches on EDMFX SGS diffusive flux. [`true`, `false` (default)]"
457457
value: false
458+
edmfx_sgs_horizontal_diffusive_flux:
459+
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)]"
460+
value: false
458461
updraft_number:
459462
help: "Sets the number of updrafts for the EDMF scheme"
460463
value: 1

docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ makedocs(;
5656
"Physics & Parameterizations" => [
5757
"Microphysics" => "microphysics.md",
5858
"EDMF: Prognostic Equations" => "edmf_equations.md",
59+
"EDMF: Horizontal Diffusion" => "edmf_horizontal_diffusion.md",
5960
"Gravity Wave Drag" => "gravity_wave.md",
6061
"Ocean Surface Albedo" => "surface_albedo.md",
6162
"Topography Representation" => "topography.md",
Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
# Horizontal EDMF diffusion
2+
3+
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`).
4+
5+
## Formulation
6+
7+
The horizontal flux uses the same TKE-based eddy diffusivity as the vertical flux,
8+
9+
```math
10+
K = c_m \, l \, \sqrt{e} \,/\, \mathrm{Pr}_t ,
11+
```
12+
13+
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,
14+
15+
```math
16+
\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],
17+
```
18+
19+
```math
20+
\partial_t (\rho \chi) \mathrel{+}= \nabla_h \cdot (\rho \, K_h \, \nabla_h \chi),
21+
\qquad
22+
\partial_t (\rho e) \mathrel{+}= \nabla_h \cdot (\rho \, K_u \, \nabla_h e),
23+
```
24+
25+
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.
26+
27+
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,
28+
29+
```math
30+
\partial_t u_h \mathrel{-}= \frac{1}{\rho} \nabla_h \cdot (\rho \, \tau),
31+
\qquad
32+
\partial_t u_3 \mathrel{-}= \frac{1}{\rho} \nabla_h \cdot (\rho \, \tau),
33+
```
34+
35+
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.
36+
37+
The corresponding shear production of TKE uses the strain rate built from horizontal gradients only, ``\mathcal{S}_h``,
38+
39+
```math
40+
\partial_t (\rho e) \mathrel{+}= 2 \rho \, K_u \, \mathcal{S}_h : \mathcal{S}_h ,
41+
```
42+
43+
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.
44+
45+
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.
46+
47+
## Anisotropic length scale
48+
49+
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,
50+
51+
```math
52+
l = \min(l_\text{phys}, \, \Delta).
53+
```
54+
55+
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:
56+
57+
```math
58+
l_v = \min(l_\text{phys}, \, \Delta z),
59+
\qquad
60+
l_h = \min(l_\text{phys}, \, \Delta x).
61+
```
62+
63+
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:
64+
65+
- where ``l_\text{phys} < \min(\Delta x, \Delta z)`` (for example small near-surface
66+
eddies), ``l_h = l_v = l_\text{phys}`` and the diffusivity is isotropic;
67+
- where ``l_\text{phys} > \max(\Delta x, \Delta z)``, the ratio ``l_h / l_v = \Delta x / \Delta z`` recovers the Smagorinsky grid anisotropy.
68+
69+
## When to enable it
70+
71+
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.
72+
73+
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.
74+
75+
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.
76+
77+
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.
78+
79+
## Configuration
80+
81+
```yaml
82+
edmfx_sgs_horizontal_diffusive_flux: true # default: false
83+
```
84+
85+
## Diagnostics
86+
87+
Three diagnostic variables expose the horizontal closure fields:
88+
89+
- `lmixh`: the mixing length with the grid-scale limit set by the horizontal node
90+
spacing, ``l_h`` [m];
91+
- `edth`: the horizontal eddy diffusivity for scalars, ``K_h`` [m² s⁻¹];
92+
- `evuh`: the horizontal eddy viscosity, ``K_u`` [m² s⁻¹].
93+
94+
The vertical counterparts are `lmix`, `edt`, and `evu`.

src/cache/precomputed_quantities.jl

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,17 @@ function precomputed_quantities(Y, atmos)
197197
ᶠK_u = similar(Y.f, FT),
198198
ᶠK_entr = similar(Y.f, FT),
199199
ᶜl_mix = similar(Y.c, FT),
200+
# Horizontal eddy viscosity/diffusivity of the TKE-based
201+
# closure, with the mixing length capped by the horizontal
202+
# node spacing; written by `set_horizontal_diffusivities!`.
203+
(
204+
atmos.edmfx_model.sgs_diffusive_flux_horizontal isa
205+
Val{true} ?
206+
(;
207+
ᶜK_u_h = similar(Y.c, FT),
208+
ᶜK_h_h = similar(Y.c, FT),
209+
) : (;)
210+
)...,
200211
) : (;)
201212
)...,
202213
ᶜstrain_rate_norm = similar(Y.c, FT),
@@ -722,6 +733,11 @@ NVTX.@annotate function set_explicit_precomputed_quantities!(Y, p, t)
722733
# from the covariance/cloud-fraction update above.
723734
set_face_diffusivities!(Y, p)
724735

736+
if turbconv_model isa AbstractEDMF &&
737+
p.atmos.edmfx_model.sgs_diffusive_flux_horizontal isa Val{true}
738+
set_horizontal_diffusivities!(Y, p)
739+
end
740+
725741
# Master mixing length at centers for consumers that live at centers
726742
# (TKE dissipation, covariance closure, updraft internal diffusion,
727743
# diagnostics). When the configuration uses (co)variances, ᶜl_mix is

src/config/model_getters.jl

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -664,6 +664,16 @@ function check_case_consistency(parsed_args)
664664
"Unknown `config = $(repr(config))`. Valid options are: $(join(valid_configs, ", "))."
665665
)
666666

667+
if parsed_args["edmfx_sgs_horizontal_diffusive_flux"] && (
668+
!isnothing(parsed_args["smagorinsky_lilly"]) || parsed_args["amd_les"]
669+
)
670+
error(
671+
"`edmfx_sgs_horizontal_diffusive_flux` cannot be combined with \
672+
`smagorinsky_lilly` or `amd_les`, which already apply horizontal \
673+
SGS diffusion to the same fields",
674+
)
675+
end
676+
667677
# ISDAC consistency: when initial_condition is ISDAC, surface/rad/external
668678
# forcing must all be set to the matching ISDAC variants. Subsidence,
669679
# 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}
769779
detr_model = get_detrainment_model(pa),
770780
sgs_mass_flux = pa["edmfx_sgs_mass_flux"],
771781
sgs_diffusive_flux = pa["edmfx_sgs_diffusive_flux"],
782+
sgs_diffusive_flux_horizontal = pa["edmfx_sgs_horizontal_diffusive_flux"],
772783
nh_pressure = pa["edmfx_nh_pressure"],
773784
vertical_diffusion = pa["edmfx_vertical_diffusion"],
774785
filter = pa["edmfx_filter"],

src/diagnostics/Diagnostics.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,7 @@ import ..get_physical_w
6666
import ..compute_gm_mixing_length
6767

6868
import ..horizontal_integral_at_boundary
69+
import ..horizontal_filter_scale
6970
import ..ᶜmixing_length
7071
import ..eddy_diffusivity
7172
import ..eddy_viscosity

src/diagnostics/edmfx_diagnostics.jl

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -874,6 +874,74 @@ add_diagnostic_variable!(short_name = "kentr", units = "m^2 s^-1",
874874
compute = compute_kentr,
875875
)
876876

877+
###
878+
# Horizontal mixing length (3d)
879+
###
880+
compute_lmixh(state, cache, time) =
881+
compute_lmixh(state, cache, time, cache.atmos.turbconv_model)
882+
compute_lmixh(_, _, _, turbconv_model) =
883+
error_diagnostic_variable("lmixh", turbconv_model)
884+
885+
function compute_lmixh(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
886+
Δx_h = horizontal_filter_scale(axes(state.c))
887+
return ᶜmixing_length(state, cache; grid_scale = Δx_h)
888+
end
889+
890+
add_diagnostic_variable!(short_name = "lmixh", units = "m",
891+
long_name = "Environment Horizontal Mixing Length",
892+
comments = "Mixing length with the grid-scale limit set by the horizontal \
893+
node spacing rather than the resolvability filter scale",
894+
compute = compute_lmixh,
895+
)
896+
897+
###
898+
# Horizontal diffusivity of heat (3d)
899+
###
900+
compute_edth(state, cache, time) =
901+
compute_edth(state, cache, time, cache.atmos.turbconv_model)
902+
compute_edth(_, _, _, turbconv_model) =
903+
error_diagnostic_variable("edth", turbconv_model)
904+
905+
function compute_edth(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
906+
hasproperty(cache.precomputed, :ᶜK_h_h) || error_diagnostic_variable(
907+
"Can only compute the horizontal eddy diffusivity when \
908+
`edmfx_sgs_horizontal_diffusive_flux` is enabled",
909+
)
910+
return cache.precomputed.ᶜK_h_h
911+
end
912+
913+
add_diagnostic_variable!(short_name = "edth", units = "m^2 s^-1",
914+
long_name = "Horizontal Eddy Diffusivity Coefficient for Temperature",
915+
comments = "Horizontal diffusion coefficient for scalars in the EDMFX \
916+
horizontal SGS diffusive flux, reported from the cached \
917+
`ᶜK_h_h` applied by the tendency",
918+
compute = compute_edth,
919+
)
920+
921+
###
922+
# Horizontal diffusivity of momentum (3d)
923+
###
924+
compute_evuh(state, cache, time) =
925+
compute_evuh(state, cache, time, cache.atmos.turbconv_model)
926+
compute_evuh(_, _, _, turbconv_model) =
927+
error_diagnostic_variable("evuh", turbconv_model)
928+
929+
function compute_evuh(state, cache, _, ::Union{EDOnlyEDMFX, PrognosticEDMFX})
930+
hasproperty(cache.precomputed, :ᶜK_u_h) || error_diagnostic_variable(
931+
"Can only compute the horizontal eddy viscosity when \
932+
`edmfx_sgs_horizontal_diffusive_flux` is enabled",
933+
)
934+
return cache.precomputed.ᶜK_u_h
935+
end
936+
937+
add_diagnostic_variable!(short_name = "evuh", units = "m^2 s^-1",
938+
long_name = "Horizontal Eddy Viscosity Coefficient for Momentum",
939+
comments = "Horizontal eddy viscosity from the TKE-based closure with the \
940+
mixing length limited by the horizontal node spacing, reported \
941+
from the cached `ᶜK_u_h` applied by the tendency",
942+
compute = compute_evuh,
943+
)
944+
877945
###
878946
# Updraft passive gas tracer A (3d)
879947
###

src/prognostic_equations/eddy_diffusion_closures.jl

Lines changed: 38 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -797,7 +797,22 @@ end
797797
return min(ℓ_phys, ml.l_grid)
798798
end
799799

800-
function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
800+
"""
801+
ᶜmixing_length(Y, p, property = Val(:master); grid_scale)
802+
803+
Lazy center-space mixing length of the TKE-based closure, with `property`
804+
selecting a `MixingLength` component (see `get_mixing_length_field`).
805+
806+
# Keyword Arguments
807+
808+
- `grid_scale`: upper bound on the mixing length. By default, the
809+
resolvability filter scale `max(Δx_h, Δz)` (see
810+
[`resolvability_filter_scale`](@ref)).
811+
"""
812+
function ᶜmixing_length(
813+
Y, p, property::Val{P} = Val{:master}();
814+
grid_scale = resolvability_filter_scale(axes(Y.c)),
815+
) where {P}
801816
(; params) = p
802817
(; ustar, obukhov_length) = p.precomputed.sfc_conditions
803818
# Stability-biased buoyancy gradient: registers unresolved inversions
@@ -807,7 +822,7 @@ function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
807822
(; ᶜN²_eff, ᶜbuoygrad, ᶜstrain_rate_norm) = p.precomputed
808823
ᶜz = Fields.coordinate_field(Y.c).z
809824
z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
810-
ᶜΔ_f = resolvability_filter_scale(axes(Y.c))
825+
ᶜΔ_f = grid_scale
811826

812827
# ᶜmixing_length is only evaluated for AbstractEDMF, which always carries
813828
# Y.c.ρtke.
@@ -846,6 +861,27 @@ function ᶜmixing_length(Y, p, property::Val{P} = Val{:master}()) where {P}
846861
return @. lazy(get_mixing_length_field(ᶜmixing_length_tuple, property))
847862
end
848863

864+
"""
865+
set_horizontal_diffusivities!(Y, p)
866+
867+
Compute and cache the horizontal eddy viscosity `ᶜK_u_h` and eddy diffusivity
868+
`ᶜK_h_h` of the TKE-based closure, with the mixing length capped by the
869+
horizontal node spacing, `l_h = min(l_phys, Δx_h)`.
870+
"""
871+
function set_horizontal_diffusivities!(Y, p)
872+
(; params) = p
873+
(; ᶜK_u_h, ᶜK_h_h, ᶜN²_eff, ᶜstrain_rate_norm) = p.precomputed
874+
turbconv_params = CAP.turbconv_params(params)
875+
Δx_h = horizontal_filter_scale(axes(Y.c))
876+
ᶜl_h = ᶜmixing_length(Y, p; grid_scale = Δx_h)
877+
ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ))
878+
@. ᶜK_u_h = eddy_viscosity(turbconv_params, ᶜtke, ᶜl_h)
879+
ᶜprandtl_nvec =
880+
@. lazy(turbulent_prandtl_number(params, ᶜN²_eff, ᶜstrain_rate_norm))
881+
@. ᶜK_h_h = eddy_diffusivity(ᶜK_u_h, ᶜprandtl_nvec)
882+
return nothing
883+
end
884+
849885
"""
850886
ᶜdiffusive_flux_divergenceᵥ(ᶠcoef, ᶜχ)
851887

0 commit comments

Comments
 (0)