feat(edmfx): add the momentum component of the horizontal SGS diffusive flux

Apply the horizontal weak divergence of the SGS stress `τ = -2 K_u S`, with
the full three-dimensional strain rate of the grid-mean velocity, to the
grid-mean momentum: 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 eddy viscosity is the horizontal `K_u` from the
TKE-based closure with the mixing length limited by the horizontal node
spacing. Enabled under the existing `edmfx_sgs_horizontal_diffusive_flux`
option; the vertical stress divergence remains with the vertical diffusion
pathway.

The shear production of TKE from horizontal gradients is not included here.

Extend the tests: horizontally uniform velocity produces no momentum
tendency, perturbed velocity produces one, and the single-column tendency is
identically zero.

Author: haakon-e

NEWS.md

@@ -4,6 +4,10 @@ ClimaAtmos.jl Release Notes
 main
 ----
 
+- ![][badge-✨feature/enhancement] Add the momentum component of the horizontal
+  EDMFX SGS diffusive flux: the horizontal weak divergence of the stress
+  `τ = -2 K_u S` with the full strain rate, applied under the same
+  `edmfx_sgs_horizontal_diffusive_flux` option.
 - [#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

docs/src/edmf_horizontal_diffusion.md

@@ -35,7 +35,22 @@ tracers — i.e. the same set the vertical flux diffuses. The microphysics and p
 tracers are scaled by the tracer diffusion factor ``\alpha`` (`α_vert_diff_tracer`,
 shared with the vertical flux). The total-specific-humidity flux additionally enters the
 air-mass tendency ``\partial_t \rho``; condensate and precipitation fluxes do not.
-Horizontal momentum diffusion is not included.
+
+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 shear production of TKE from horizontal gradients is not
+included.
 
 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.

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -359,8 +359,8 @@ function edmfx_sgs_horizontal_diffusive_flux_tendency!(
         buoyancy_gradient = ᶜlinear_buoygrad,
     )
 
-    ᶜK_u_h =
-        @. lazy(eddy_viscosity(turbconv_params, ᶜtke, ᶜmixing_length_h))
+    ᶜK_u_h = p.scratch.ᶜtemp_scalar_4
+    @. ᶜK_u_h = eddy_viscosity(turbconv_params, ᶜtke, ᶜmixing_length_h)
     ᶜprandtl_nvec = @. lazy(
         turbulent_prandtl_number(params, ᶜlinear_buoygrad, ᶜstrain_rate_norm),
     )
@@ -395,5 +395,17 @@ function edmfx_sgs_horizontal_diffusive_flux_tendency!(
         @. Yₜ.c.ρtke += wdivₕ(ᶜρ * ᶜK_u_h * gradₕ(ᶜtke))
     end
 
+    # Momentum: horizontal weak divergence of the SGS stress `τ = -2 K_u S`
+    # with the full strain rate. The vertical stress divergence is handled by
+    # the vertical diffusion pathway.
+    (; ᶜu, ᶠu) = p.precomputed
+    ᶠρ = @. p.scratch.ᶠtemp_scalar = ᶠinterp(ᶜρ)
+    ᶜτ_h = compute_strain_rate_center_full!(p.scratch.ᶜtemp_UVWxUVW, ᶜu, ᶠu)
+    @. ᶜτ_h = -2 * ᶜK_u_h * ᶜτ_h
+    @. Yₜ.c.uₕ -= C12(wdivₕ(ᶜρ * ᶜτ_h) / ᶜρ)
+    ᶠτ_h = compute_strain_rate_face_full!(p.scratch.ᶠtemp_UVWxUVW, ᶜu, ᶠu)
+    @. ᶠτ_h = -2 * ᶠinterp(ᶜK_u_h) * ᶠτ_h
+    @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_h) / ᶠρ)
+
     return nothing
 end

test/prognostic_equations/edmfx_horizontal_diffusion_tests.jl

@@ -74,16 +74,21 @@ box_config_dict(; extra...) = bomex_edmfx_config_dict(;
     uniform_max = Dict(
         name => maximum(abs, parent(getproperty(Yₜ.c, name))) for name in rows
     )
+    uniform_max_uₕ = maximum(abs, parent(Yₜ.c.uₕ))
+    uniform_max_u₃ = maximum(abs, parent(Yₜ.f.u₃))
 
     # Perturb the state horizontally to activate every flux pathway. The
-    # enthalpy flux reads the precomputed `ᶜh_tot`, which is perturbed
-    # directly.
+    # enthalpy flux reads the precomputed `ᶜh_tot` and the momentum stress
+    # reads the precomputed `ᶜu`/`ᶠu`, which are perturbed directly.
     ᶜx = Fields.coordinate_field(Y.c).x
+    ᶠx = Fields.coordinate_field(Y.f).x
     ᶜpert = @. 1 + FT(0.1) * sin(FT(2π) * ᶜx / FT(6400))
     @. Y.c.ρq_tot *= ᶜpert
     @. Y.c.ρq_rai = FT(1e-5) * Y.c.ρ * ᶜpert
     @. Y.c.ρtke *= ᶜpert
     @. p.precomputed.ᶜh_tot *= 1 + FT(1e-3) * sin(FT(2π) * ᶜx / FT(6400))
+    @. p.precomputed.ᶜu *= ᶜpert
+    @. p.precomputed.ᶠu *= 1 + FT(0.1) * sin(FT(2π) * ᶠx / FT(6400))
 
     Yₜ .= zero(eltype(Yₜ))
     CA.edmfx_sgs_horizontal_diffusive_flux_tendency!(
@@ -93,11 +98,16 @@ box_config_dict(; extra...) = bomex_edmfx_config_dict(;
         name => maximum(abs, parent(getproperty(Yₜ.c, name))) for name in rows
     )
 
+    perturbed_max_uₕ = maximum(abs, parent(Yₜ.c.uₕ))
+    perturbed_max_u₃ = maximum(abs, parent(Yₜ.f.u₃))
+
     # Non-vacuous fluxes for every perturbed pathway.
     @test perturbed_max[:ρq_tot] > 0
     @test perturbed_max[:ρe_tot] > 0
     @test perturbed_max[:ρq_rai] > 0
     @test perturbed_max[:ρtke] > 0
+    @test perturbed_max_uₕ > 0
+    @test perturbed_max_u₃ > 0
 
     # The moist air mass tendency matches the ρq_tot tendency bitwise.
     @test parent(Yₜ.c.ρ) == parent(Yₜ.c.ρq_tot)
@@ -108,6 +118,8 @@ box_config_dict(; extra...) = bomex_edmfx_config_dict(;
         reference = max(perturbed_max[name], perturbed_max[:ρq_tot])
         @test uniform_max[name] <= FT(1e-10) * reference
     end
+    @test uniform_max_uₕ <= FT(1e-10) * perturbed_max_uₕ
+    @test uniform_max_u₃ <= FT(1e-10) * perturbed_max_u₃
     # Identically zero tracers stay identically zero.
     for name in (:ρq_lcl, :ρq_icl, :ρq_sno)
         @test uniform_max[name] == 0
@@ -165,6 +177,8 @@ end
     for name in (:ρ, :ρe_tot, :ρq_tot, :ρq_lcl, :ρq_icl, :ρq_rai, :ρq_sno, :ρtke)
         @test maximum(abs, parent(getproperty(Yₜ.c, name))) == 0
     end
+    @test maximum(abs, parent(Yₜ.c.uₕ)) == 0
+    @test maximum(abs, parent(Yₜ.f.u₃)) == 0
 end
 
 @testset "Incompatible horizontal SGS closures are rejected" begin