feat(edmfx): add horizontal diffusion of prognostic EDMFX updrafts

Diffuse the prognostic updraft variables horizontally, mirroring
`edmfx_vertical_diffusion` with the mixing length limited by the horizontal
node spacing: moist static energy and total specific humidity at `K_h`, and
the updraft SGS tracers at `α_vert_diff_tracer · K_h`, all weighted by the
updraft density. The updraft dry-air mass is unchanged by the water flux, so
`ρa` receives the counterpart tendency `ρa/(1 - q_tot)` times the `q_tot`
tendency, mirroring the hyperdiffusion treatment. Always explicit; enabled by
the opt-in `edmfx_horizontal_diffusion` config option (default false).

Extend the tests: horizontally uniform updraft fields produce no tendency,
perturbed fields produce one with the stated `ρa`/`q_tot` pairing, and the
single-column tendency is identically zero.

Author: haakon-e

NEWS.md

@@ -4,6 +4,7 @@ ClimaAtmos.jl Release Notes
 main
 ----
 
+- [#4659](https://github.com/CliMA/ClimaAtmos.jl/pull/4659) ![][badge-✨feature/enhancement] Add horizontal diffusion of prognostic EDMFX updrafts (moist static energy, total specific humidity with its area-weighted density counterpart, and SGS tracers), enabled by the opt-in `edmfx_horizontal_diffusion` config option (default `false`).
 - [#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).

config/default_configs/default_config.yml

@@ -434,6 +434,9 @@ edmfx_nh_pressure:
 edmfx_vertical_diffusion:
   help: "If set to true, it switches on vertical diffusion of prognostic EDMFX updrafts. [`true`, `false` (default)]"
   value: false
+edmfx_horizontal_diffusion:
+  help: "If set to true, it switches on horizontal diffusion of prognostic EDMFX updrafts, with the mixing length limited by the horizontal node spacing. Intended for high-resolution box runs.  [`true`, `false` (default)]"
+  value: false
 edmfx_entr_model:
   help: "EDMFX entrainment closure.  [`nothing` (default), `PiGroups`, `Generalized`]"
   value: ~

docs/src/edmf_horizontal_diffusion.md

@@ -42,6 +42,16 @@ The corresponding shear production of TKE uses the strain rate built from horizo
 
 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.
 
+## Updraft horizontal diffusion
+
+The separate option `edmfx_horizontal_diffusion` (default `false`) switches on horizontal diffusion of the prognostic updraft variables, mirroring the vertical updraft diffusion of `edmfx_vertical_diffusion` with the horizontal mixing length: for each updraft,
+
+```math
+\partial_t \chi^j \mathrel{+}= \frac{1}{\rho^j} \nabla_h \cdot (\rho^j \, K_h \, \nabla_h \chi^j),
+```
+
+for the moist static energy, the total specific humidity, and the updraft SGS tracers (the latter scaled by ``\alpha``). The updraft dry-air mass is unchanged by the water flux, so the area-weighted density ``\rho a^j`` receives the counterpart tendency ``\rho a^j/(1 - q_\text{tot}^j)`` times the ``q_\text{tot}^j`` tendency, mirroring the hyperdiffusion treatment.
+
 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

src/config/model_getters.jl

@@ -782,6 +782,7 @@ function AtmosTurbconv(config::AtmosConfig, params, ::Type{FT}) where {FT}
         sgs_diffusive_flux_horizontal = pa["edmfx_sgs_horizontal_diffusive_flux"],
         nh_pressure = pa["edmfx_nh_pressure"],
         vertical_diffusion = pa["edmfx_vertical_diffusion"],
+        horizontal_diffusion = pa["edmfx_horizontal_diffusion"],
         filter = pa["edmfx_filter"],
         scale_blending_method,
     )

src/prognostic_equations/mass_flux_closures.jl

@@ -236,6 +236,65 @@ function edmfx_vertical_diffusion_tendency!(
     end
 end
 
+edmfx_horizontal_diffusion_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
+
+function edmfx_horizontal_diffusion_tendency!(
+    Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX,
+)
+    p.atmos.edmfx_model.horizontal_diffusion isa Val{true} || return nothing
+    iscolumn(axes(Y.c)) && return nothing
+    (; params) = p
+    (; ᶜρʲs, ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
+    turbconv_params = CAP.turbconv_params(params)
+    n = n_mass_flux_subdomains(turbconv_model)
+
+    ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ))
+    Δx = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(Y.c)))
+    ᶜmixing_length_h = p.scratch.ᶜtemp_scalar
+    ᶜmixing_length_h .= ᶜmixing_length(
+        Y, p; grid_scale = Δx, buoyancy_gradient = ᶜlinear_buoygrad,
+    )
+    ᶜK_u_h = @. lazy(eddy_viscosity(turbconv_params, ᶜtke, ᶜmixing_length_h))
+    ᶜprandtl_nvec =
+        @. lazy(turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm))
+    ᶜK_h_h = p.scratch.ᶜtemp_scalar_2
+    @. ᶜK_h_h = eddy_diffusivity(ᶜK_u_h, ᶜprandtl_nvec)
+
+    ᶜq_totʲₜ_diffusion = p.scratch.ᶜtemp_scalar_3
+    for j in 1:n
+        ᶜρʲ = ᶜρʲs.:($j)
+        ᶜmseʲ = Y.c.sgsʲs.:($j).mse
+        ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
+        @. Yₜ.c.sgsʲs.:($$j).mse +=
+            wdivₕ(ᶜρʲ * ᶜK_h_h * gradₕ(ᶜmseʲ)) / ᶜρʲ
+        @. ᶜq_totʲₜ_diffusion =
+            wdivₕ(ᶜρʲ * ᶜK_h_h * gradₕ(ᶜq_totʲ)) / ᶜρʲ
+        @. Yₜ.c.sgsʲs.:($$j).q_tot += ᶜq_totʲₜ_diffusion
+        # The updraft dry-air mass is unchanged by the water flux, so the
+        # area-weighted density responds to the change in total specific
+        # humidity, mirroring the hyperdiffusion treatment.
+        @. Yₜ.c.sgsʲs.:($$j).ρa +=
+            Y.c.sgsʲs.:($$j).ρa / (1 - ᶜq_totʲ) * ᶜq_totʲₜ_diffusion
+    end
+
+    if !isempty(sgs_tracer_names(Y))
+        α_diff_microphysics = CAP.α_vert_diff_tracer(params)
+        ᶜρʲ = ᶜρʲs.:(1)
+        # Sedimenting microphysics species are diffused with
+        # α_vert_diff_tracer * K_h, passive tracers with the unscaled K_h,
+        # matching the vertical updraft diffusion.
+        for χ_name in sgs_tracer_names(Y)
+            α =
+                χ_name in sgs_sedimenting_tracer_candidates ?
+                α_diff_microphysics : one(α_diff_microphysics)
+            ᶜχʲ = MatrixFields.get_field(Y.c.sgsʲs.:(1), χ_name)
+            ᶜχʲₜ = MatrixFields.get_field(Yₜ.c.sgsʲs.:(1), χ_name)
+            @. ᶜχʲₜ += wdivₕ(ᶜρʲ * ᶜK_h_h * α * gradₕ(ᶜχʲ)) / ᶜρʲ
+        end
+    end
+    return nothing
+end
+
 # Private helper: clips grid-mean condensate tracers to non-negative values and
 # rescales the condensate sum so it cannot exceed the available total moisture.
 function enforce_grid_mean_microphysics_constraints!(Y, p, t)

src/prognostic_equations/remaining_tendency.jl

@@ -337,6 +337,7 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     horizontal_constant_diffusion_tendency!(Yₜ, Y, p, t, chd)
 
     edmfx_sgs_horizontal_diffusive_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    edmfx_horizontal_diffusion_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 
     # Optional tendency to bring negative small tracers back from negative
     # at the cost of water vapor.

src/types.jl

@@ -811,7 +811,8 @@ const ValTF = Union{Val{true}, Val{false}}
 
 struct EDMFXModel{
     EEM, EDM,
-    ESMF <: ValTF, ESDF <: ValTF, ESDFH <: ValTF, ENP <: ValTF, EVD <: ValTF, EF <: ValTF,
+    ESMF <: ValTF, ESDF <: ValTF, ESDFH <: ValTF, ENP <: ValTF, EVD <: ValTF,
+    EHD <: ValTF, EF <: ValTF,
     SBM <: AbstractScaleBlendingMethod,
 }
     entr_model::EEM
@@ -821,6 +822,7 @@ struct EDMFXModel{
     sgs_diffusive_flux_horizontal::ESDFH
     nh_pressure::ENP
     vertical_diffusion::EVD
+    horizontal_diffusion::EHD
     filter::EF
     scale_blending_method::SBM
 end
@@ -836,6 +838,7 @@ function EDMFXModel(;
     sgs_diffusive_flux_horizontal::Union{Bool, ValTF} = false,
     nh_pressure::Union{Bool, ValTF} = false,
     vertical_diffusion::Union{Bool, ValTF} = false,
+    horizontal_diffusion::Union{Bool, ValTF} = false,
     filter::Union{Bool, ValTF} = false,
     scale_blending_method,
     kwargs...,
@@ -851,6 +854,7 @@ function EDMFXModel(;
         parse_val_tf(sgs_diffusive_flux_horizontal),
         parse_val_tf(nh_pressure),
         parse_val_tf(vertical_diffusion),
+        parse_val_tf(horizontal_diffusion),
         parse_val_tf(filter),
         scale_blending_method,
     )

test/prognostic_equations/edmfx_horizontal_diffusion_tests.jl

@@ -163,6 +163,52 @@ box_config_dict(; extra...) = bomex_edmfx_config_dict(;
     @test parent(ᶜl_buoygrad) == parent(ᶜl_default)
 end
 
+@testset "EDMFX updraft horizontal diffusion (Bomex box)" begin
+    config = CA.AtmosConfig(
+        box_config_dict(; edmfx_horizontal_diffusion = true);
+        job_id = "edmfx_updraft_horizontal_diffusion_box_test",
+    )
+    (; Y, p, simulation) = generate_test_simulation(config)
+    t = simulation.integrator.t
+    FT = eltype(Y)
+
+    @. Y.c.ρtke = FT(0.5) * Y.c.ρ
+    # A nonzero updraft area makes the ρa counterpart nonvacuous.
+    @. Y.c.sgsʲs.:(1).ρa = FT(0.1) * Y.c.ρ
+
+    up_row_max(Yₜ) = (
+        maximum(abs, parent(Yₜ.c.sgsʲs.:(1).mse)),
+        maximum(abs, parent(Yₜ.c.sgsʲs.:(1).q_tot)),
+        maximum(abs, parent(Yₜ.c.sgsʲs.:(1).ρa)),
+    )
+    Yₜ = similar(Y)
+    Yₜ .= zero(eltype(Yₜ))
+    CA.edmfx_horizontal_diffusion_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    uniform_mse, uniform_q_tot, uniform_ρa = up_row_max(Yₜ)
+
+    ᶜx = Fields.coordinate_field(Y.c).x
+    ᶜpert = @. 1 + FT(0.1) * sin(FT(2π) * ᶜx / FT(6400))
+    @. Y.c.sgsʲs.:(1).q_tot *= ᶜpert
+    @. Y.c.sgsʲs.:(1).mse *= 1 + FT(1e-3) * sin(FT(2π) * ᶜx / FT(6400))
+
+    Yₜ .= zero(eltype(Yₜ))
+    CA.edmfx_horizontal_diffusion_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    perturbed_mse, perturbed_q_tot, perturbed_ρa = up_row_max(Yₜ)
+    @test perturbed_mse > 0
+    @test perturbed_q_tot > 0
+    @test perturbed_ρa > 0
+    @test uniform_mse <= FT(1e-10) * perturbed_mse
+    @test uniform_q_tot <= FT(1e-10) * perturbed_q_tot
+    @test uniform_ρa <= FT(1e-10) * perturbed_ρa
+
+    # The ρa tendency is the q_tot tendency scaled by ρa / (1 - q_tot).
+    ᶜρa_expected = similar(Y.c.ρ)
+    @. ᶜρa_expected =
+        Y.c.sgsʲs.:(1).ρa / (1 - Y.c.sgsʲs.:(1).q_tot) *
+        Yₜ.c.sgsʲs.:(1).q_tot
+    @test parent(Yₜ.c.sgsʲs.:(1).ρa) ≈ parent(ᶜρa_expected) rtol = FT(1e-12)
+end
+
 @testset "EDMFX horizontal diffusive flux vanishes on a column" begin
     config = CA.AtmosConfig(
         bomex_edmfx_config_dict(;
@@ -172,6 +218,7 @@ end
             dt = "120secs",
             t_end = "30mins",
             z_elem = 30,
+            edmfx_horizontal_diffusion = true,
         );
         job_id = "edmfx_horizontal_diffusion_column_test",
     )
@@ -182,7 +229,7 @@ end
     @. Y.c.ρtke = FT(0.5) * Y.c.ρ
 
     # Spectral operators return exact zeros on a column, so all existing
-    # single-column results are unaffected by the option.
+    # single-column results are unaffected by the options.
     Yₜ = similar(Y)
     Yₜ .= zero(eltype(Yₜ))
     CA.edmfx_sgs_horizontal_diffusive_flux_tendency!(
@@ -193,6 +240,12 @@ end
     end
     @test maximum(abs, parent(Yₜ.c.uₕ)) == 0
     @test maximum(abs, parent(Yₜ.f.u₃)) == 0
+
+    Yₜ .= zero(eltype(Yₜ))
+    CA.edmfx_horizontal_diffusion_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    @test maximum(abs, parent(Yₜ.c.sgsʲs.:(1).mse)) == 0
+    @test maximum(abs, parent(Yₜ.c.sgsʲs.:(1).q_tot)) == 0
+    @test maximum(abs, parent(Yₜ.c.sgsʲs.:(1).ρa)) == 0
 end
 
 @testset "Incompatible horizontal SGS closures are rejected" begin