Geometric SGS variance: isentropic form of the horizontal term

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Author: szy21

NEWS.md

@@ -3,6 +3,16 @@ ClimaAtmos.jl Release Notes
 
 main
 ----
+- ![][badge-🔥behavioralΔ] The horizontal resolved-gradient (geometric) SGS variance
+  term is now the isentropic form: `q′q′` carries `|∇_h q_tot − r ∇_h θ_li|²`, the
+  gradient of `q_tot` along the `θ_li` surface with `r = (∂q_tot/∂z)/(∂θ_li/∂z)`, the
+  stratification floored at the provisional `sgs_variance_isentropic_min_dtheta_dz`,
+  there is no horizontal `θ′θ′` (or `θ′q′`) term,
+  and the term carries a stability weight on the saturated moist `N²` — the hard cut
+  `1[N²_sat > 0]` by default, smoothed over the provisional
+  `sgs_variance_isentropic_stability_width` when that is positive — so it vanishes where
+  the layer is moist-neutral or unstable. Changes every run with a nonzero
+  `sgs_variance_horizontal_scale_factor`.
 - ![][badge-🔥behavioralΔ] Add a resolved vertical-gradient block
   `c_g (c_Δz Δz)² (∂_z ψ)²` to `σ_T²`, `σ_q²` (uniform-box exact with the
   defaults `c_g = 1/12`, `c_Δz = 1`; independent of `l_mix`), so the

config/default_configs/default_config.yml

@@ -262,7 +262,7 @@ use_sgs_quadrature:
   help: "Enable subgrid-scale quadrature integration for microphysics tendencies. When true, microphysics tendencies are integrated over the joint SGS PDF of (T, q_tot). [`true` (default), `false`]"
   value: true
 tq_correlation_model:
-  help: "Closure for the SGS T-q correlation that couples σ_T and σ_q in the sampled joint PDF. `constant` uses the prescribed `Tq_correlation_coefficient`. `diagnosed` builds T′q′ as the sum of three blocks consistent with the variance decomposition — turbulent (prescribed ρ_turb · σ_T,turb · σ_q,turb), horizontal geometric (c_g (c_Δx Δx_h)² ∇_h θ_li · ∇_h q_tot), and vertical uniform-box (Δz² ∂_z θ_li · ∂_z q_tot / 12) — then forms ρ = T′q′ / √(T′T′ q′q′), falling back to the prescribed value where the variances vanish. [`constant` (default), `diagnosed`]"
+  help: "Closure for the SGS T-q correlation that couples σ_T and σ_q in the sampled joint PDF. `constant` uses the prescribed `Tq_correlation_coefficient`. `diagnosed` builds T′q′ as the sum of the covariance blocks consistent with the variance decomposition — turbulent (prescribed ρ_turb · σ_T,turb · σ_q,turb) and vertical uniform-box (Δz² ∂_z θ_li · ∂_z q_tot / 12); the isentropic horizontal block carries no θ′ and contributes none — then forms ρ = T′q′ / √(T′T′ q′q′), falling back to the prescribed value where the variances vanish. [`constant` (default), `diagnosed`]"
   value: "constant"
 sgs_distribution:
   help: "Distribution type for SGS quadrature over microphysics. [`gaussian` (default), `lognormal`, `mean` (grid-mean only, no SGS sampling)]"

docs/src/prophet_closures.md

@@ -542,21 +542,31 @@ a single isotropic intra-subdomain variance shared across subdomains, driven by
 the vertical grid-mean gradients with the grid-mean mixing length.
 
 An optional second term adds the subgrid variance implied by the resolved
-horizontal field — the leading-order scale-similarity estimate from the
-horizontal gradient and the horizontal grid scale alone,
-
-```math
-\sigma_{\psi,\mathrm{geo}}^2 = c_g \, (c_{\Delta x} \Delta x_h)^2 \, |\nabla_h \psi|^2 ,
-```
-
-added to both the ``\theta_{li}`` and ``q_t`` variances. The term is off by
-default (``c_{\Delta x} = 0``). It is not Richardson-weighted: this is a
-physical feature of the resolved mean field and is independent of whether the
-turbulent closure (above) is active, so the two blocks contribute additively
-without double-counting (they are not two estimates of the same quantity).
-The total-water standard deviation is then bounded by ``\sigma_q \le r_{\max} q_t`` (``r_{\max} = 0.5``), because a wider Gaussian places a quadrature node
-at negative total water, and condensing the clamped excess can drive the
-grid-mean vapor negative.
+horizontal field. Large-scale stirring of the free troposphere is
+quasi-adiabatic, so parcels move along surfaces of ``\theta_{li}``; the
+moisture contrast a parcel brings to the cell's level is therefore the
+contrast *along the isentrope*, and after the adiabatic displacement to that
+level only ``q_t`` differs. The term is the leading-order scale-similarity
+estimate built on that gradient,
+
+```math
+\sigma_{q,\mathrm{geo}}^2 = c_g \, (c_{\Delta x} \Delta x_h)^2 \,
+  |\nabla_h q_t - r \nabla_h \theta_{li}|^2 \, w(N^2_{\mathrm{sat}}) ,
+\qquad r = \frac{\partial_z q_t}{\max(\partial_z \theta_{li}, \theta_{z,\min})} ,
+```
+
+added to the ``q_t`` variance only; there is no horizontal ``\theta_{li}``
+variance or ``\theta_{li}``–``q_t`` covariance term. ``r`` is the slope ratio
+``(\partial_z q_t)/(\partial_z \theta_{li})`` with the stratification floored at
+``\theta_{z,\min}`` (`sgs_variance_isentropic_min_dtheta_dz`), which bounds ``|r|``
+by ``|\partial_z q_t| / \theta_{z,\min}`` and treats a neutral or overturned layer
+as one of minimal stable stratification rather than as a surface of opposite tilt.
+The stability weight ``w`` is the hard cut
+``1[N^2_{\mathrm{sat}} > 0]`` on the saturated moist buoyancy frequency, or
+``(1 + \tanh(N^2_{\mathrm{sat}}/N^2_w))/2`` when the provisional width
+`sgs_variance_isentropic_stability_width` ``N^2_w`` is positive: where
+saturated air is buoyant the condensation of the subgrid moist tail is
+convective and belongs to the updrafts, so the term is zeroed.
 
 A resolved vertical-gradient block is also added,
 
@@ -571,9 +581,11 @@ uniform position in a box of size ``\Delta z``. The scale factor
 switch (``c_{\Delta z} = 0`` disables the term). This block keeps
 ``\sigma_T^2``, ``\sigma_q^2`` non-zero where the turbulent closure collapses,
 so the quadrature integrand retains in-cell variability in laminar cells and
-single-column configurations. It is not Richardson-weighted — the resolved
-vertical gradient is a physical feature of the mean profile, not a turbulent
-property.
+single-column configurations.
+
+The total-water standard deviation is then bounded by ``\sigma_q \le r_{\max} q_t`` (``r_{\max} = 0.5``),
+because a wider Gaussian places a quadrature node at negative total water, and
+condensing the clamped excess can drive the grid-mean vapor negative.
 
 The cross-covariance is set either as a prescribed correlation,
 
@@ -582,17 +594,19 @@ The cross-covariance is set either as a prescribed correlation,
 \qquad (\texttt{tq\_correlation\_model: constant}, \text{default})
 ```
 
-or as a diagnosed covariance that combines the three variance blocks with the
+or as a diagnosed covariance that combines the variance blocks with the
 correlation structure appropriate to each
 (``\texttt{tq\_correlation\_model: diagnosed}``):
 
 ```math
 \langle T' q_t' \rangle
   = r_{\mathrm{turb}} \sigma_{T,\mathrm{turb}} \sigma_{q,\mathrm{turb}}
-  + c_g (c_{\Delta x} \Delta x_h)^2 \, \nabla_h \theta_{li} \cdot \nabla_h q_t
   + c_g (c_{\Delta z} \Delta z)^2 \, \partial_z \theta_{li} \, \partial_z q_t ,
 ```
 
+the isentropic horizontal block contributing nothing (no ``\theta_{li}'`` along
+the surface), so its ``q_t`` variance only dilutes the correlation,
+
 followed by the thermodynamic Jacobian transform to ``T`` basis and
 ``r = \langle T' q_t' \rangle / \sqrt{\sigma_T^2 \sigma_q^2}``, which is
 naturally in ``[-1, 1]`` because each variance block satisfies Cauchy-Schwarz

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-425
+426
 
 # **README**
 #
@@ -32,6 +32,10 @@
 # 3) (optional) leave a link to the buildkite run that prompted this ref counter bump.
 
 #=
+426
+- The horizontal geometric SGS variance term is the isentropic form (gradient of q_tot
+  along the θ_li surface in q′q′ only, with a saturated-N² stability weight).
+
 425
 - Add SGS T–q correlation closure and resolved vertical-gradient variance block.
 

src/config/model_getters.jl

@@ -1142,10 +1142,6 @@ function AtmosWater(config::AtmosConfig, params, ::Type{FT}) where {FT}
     microphysics_model = get_microphysics_model(pa)
     sgs_quadrature = get_sgs_quadrature(pa, params)
     tq_correlation_model = get_tq_correlation_model(pa)
-    # No horizontal-term requirement for `DiagnosedTqCorrelation`: the vertical
-    # resolved-gradient block (uniform-box, always on) and the prescribed-ρ_turb
-    # turbulent cross together give a well-defined T′q′ even in single-column
-    # configurations where `sgs_variance_horizontal_scale_factor = 0`.
 
     if microphysics_model isa DryModel
         @warn "Running simulations without any moisture present."

src/parameterized_tendencies/microphysics/cloud_fraction.jl

@@ -223,6 +223,73 @@ function materialized_mixing_length!(Y, p)
     end
 end
 
+"""
+    isentropic_validity_weight(N²_sat, width)
+
+Validity weight `w ∈ [0, 1]` of the isentropic geometric variance term from the
+saturated buoyancy frequency squared `N²_sat` [s⁻²]: the sub-grid moist tail the
+term feeds into the quadrature condenses in place only where saturated air is
+stably stratified; where it is buoyant its condensation is convective and
+belongs to the EDMF plumes, so the term is zeroed there.
+With `width = 0` (`sgs_variance_isentropic_stability_width`) the
+weight is the hard cut `1[N²_sat > 0]`; a positive width smooths it, with
+w = (1 + tanh(N²_sat / width)) / 2
+"""
+@inline function isentropic_validity_weight(N²_sat, width)
+    FT = typeof(N²_sat)
+    smooth = (one(FT) + tanh(N²_sat / max(width, eps(FT)))) / 2
+    hard = ifelse(N²_sat > 0, one(FT), zero(FT))
+    return ifelse(width > zero(FT), smooth, hard)
+end
+
+"""
+    ᶜsgs_isentropic_weight(Y, p)
+
+Lazy broadcast of `isentropic_validity_weight` with the width
+`sgs_variance_isentropic_stability_width`, evaluated on the squared buoyancy
+frequency of saturated air (`blended_N²` with cloud fraction 1).
+"""
+function ᶜsgs_isentropic_weight(Y, p)
+    (; ᶜbg_coeffs, ᶜgradᵥ_θ_liq_ice, ᶜgradᵥ_q_tot) = p.precomputed
+    FT = eltype(Y.c.ρ)
+    width = FT(CAP.sgs_variance_isentropic_stability_width(p.params))
+    ᶜlg = Fields.local_geometry_field(Y.c)
+    return @. lazy(
+        isentropic_validity_weight(
+            blended_N²(
+                ᶜbg_coeffs,
+                one(FT),
+                projected_vector_data(C3, ᶜgradᵥ_θ_liq_ice, ᶜlg),
+                projected_vector_data(C3, ᶜgradᵥ_q_tot, ᶜlg),
+            ),
+            width,
+        ),
+    )
+end
+
+"""
+    isentropic_slope_ratio(∂q∂z, ∂θ∂z, θz_min)
+
+Ratio `r = (∂q_tot/∂z) / max(∂θ_li/∂z, θz_min)` that converts the horizontal
+gradient of `θ_li` into the vertical displacement a horizontal excursion along a
+sloping `θ_li` surface implies, times the vertical `q_tot` gradient (the isentropic
+geometric variance term). The floor `θz_min` (`sgs_variance_isentropic_min_dtheta_dz`)
+on the stratification keeps the ratio finite and bounds it by `|∂q/∂z| / θz_min`: a
+neutral or overturned layer (`∂θ_li/∂z ≤ θz_min`) is treated as one of minimal stable
+stratification.
+"""
+@inline isentropic_slope_ratio(∂q∂z, ∂θ∂z, θz_min) = ∂q∂z / max(∂θ∂z, θz_min)
+
+"""
+    isentropic_gradient_invariant(inv_q, inv_θq, inv_θ, r)
+
+`|∇_h q_tot − r ∇_h θ_li|² = |∇_h q|² − 2 r (∇_h θ_li · ∇_h q) + r² |∇_h θ_li|²`, the
+squared horizontal gradient of `q_tot` along the `θ_li` surface; floored
+at zero against round-off.
+"""
+@inline isentropic_gradient_invariant(inv_q, inv_θq, inv_θ, r) =
+    max(inv_q - 2 * r * inv_θq + r^2 * inv_θ, zero(inv_q))
+
 """
     hgrad_invariant!(ᶜinv, ᶜψ, p)
 
@@ -264,8 +331,9 @@ Diagnosed SGS T–q correlation `clamp(T′q′ / √(T′T′ q′q′), -1, 1)
 `fallback` where the product of the variances is at the round-off floor.
 
 In exact arithmetic each variance block satisfies its own Cauchy-Schwarz
-inequality and the Cauchy-Schwarz bound on the sum of the three (turbulent,
-horizontal geometric, vertical uniform-box) covariances also holds, so
+inequality and the Cauchy-Schwarz bound on the sum of the two (turbulent,
+vertical uniform-box) covariances also holds — the isentropic horizontal block
+contributes variance to q′q′ only and no covariance — so
 `|ρ| ≤ 1` by construction. The explicit clamp is a cheap safeguard against
 (a) Float32 roundoff at the Cauchy-Schwarz boundary (which can produce
 `|ρ| = 1 + O(ε)` and feed NaN into downstream `√(1 − ρ²)` in the bivariate
@@ -322,10 +390,10 @@ Pipeline:
     `ρ_turb` scalar as `ρ_turb · σ_T,turb · σ_q,turb` (not the gradient-direction
     product), so the turbulent block carries a physical correlation independent of
     the sign of `∂_z θ · ∂_z q`.
- 3. Add the horizontal resolved-gradient (geometric) term to θ′θ′ (skipped when
-    `sgs_variance_horizontal_scale_factor` is 0). For the diagnosed correlation,
-    add the geometric cross `∇_h θ_li · ∇_h q_tot` with the same element filter.
-    Not Richardson-weighted (consistent with the vertical uniform-box block).
+ 3. Add the horizontal resolved-gradient (geometric) term, the isentropic form
+    `|∇_h q_tot − r ∇_h θ_li|²` times the stability weight `w(N²_sat)`, to q′q′
+    only — no horizontal θ′θ′ and no horizontal cross term (skipped when
+    `sgs_variance_horizontal_scale_factor` is 0).
  4. Add the resolved vertical-gradient (uniform-box) term
     `c_g (c_Δz Δz)² (∂_z ψ)²` to both variances, and the exact cross
     `c_g (c_Δz Δz)² ∂_z θ · ∂_z q` to θ′q′ (skipped when
@@ -414,27 +482,52 @@ function set_covariance_cache!(Y, p, thermo_params)
     # ----------------------------------------------------------------------
     # 2. Horizontal resolved-gradient (geometric) block.
     #    `geo_h |∇_h ψ|²`: leading-order scale-similarity estimate from the
-    #    resolved horizontal field. Not Richardson-weighted — like the
-    #    vertical uniform-box block below, the resolved horizontal gradient
-    #    is a physical feature of the mean profile, independent of whether
-    #    the turbulent closure is active. q′q′ is in native q basis so adding
+    #    resolved horizontal field. q′q′ is in native q basis so adding
     #    here (before the θ→T transform) is equivalent to adding after.
+    #    The block is the isentropic form: no horizontal θ′θ′ term, because
+    #    θ_li is constant along the surface the excursion follows, so the whole
+    #    geometric variance is carried by q′q′ as the gradient of q_tot ALONG
+    #    the θ_li surface,
+    #    |∇_θ q|² = |∇_h q − r ∇_h θ_li|², r = (∂q/∂z)/max(∂θ_li/∂z, θz_min)
+    #    (`isentropic_slope_ratio`): a horizontal excursion δ on a sloping
+    #    isentrope carries the parcel down/up by r-weighted amounts, so aligned
+    #    (warm = moist) gradients amplify and anti-aligned ones cancel the
+    #    saturation variance. The validity weight 1[N²_sat > 0]
+    #    (`ᶜsgs_isentropic_weight`) zeroes the WHOLE term where the layer is
+    #    moist-neutral or unstable (no isentrope to displace along; the SGS
+    #    variability there is convective, carried by the EDMF plumes). The
+    #    weight is optionally smoothed over `sgs_variance_isentropic_stability_width`.
+    #    With the diagnosed T–q correlation this block adds nothing to T′q′ 
+    #    (no θ′ along the surface), so the correlation is carried by the
+    #    turbulent and vertical blocks. Scratch: `ᶜtemp_scalar_2` is free until
+    #    the θ→T transform below, `_3` (θ_li), `_5` (|∇q|²), `_6` (cross term).
     # ----------------------------------------------------------------------
     if use_geometric
         (; ᶜT, ᶜq_liq, ᶜq_ice) = p.precomputed
         ᶜθ_li = p.scratch.ᶜtemp_scalar_3
         @. ᶜθ_li = TD.liquid_ice_pottemp(
             thermo_params, ᶜT, Y.c.ρ, ᶜq_tot_nonneg, ᶜq_liq, ᶜq_ice,
         )
-        ᶜinv = p.scratch.ᶜtemp_scalar_5                # reused for 2–3 invariants
-        hgrad_invariant!(ᶜinv, ᶜθ_li, p)
-        @. ᶜT′T′ += geo_h * ᶜinv
-        hgrad_invariant!(ᶜinv, ᶜq_tot_nonneg, p)
-        @. ᶜq′q′ += geo_h * ᶜinv
-        if diag_corr
-            hgrad_cross_invariant!(ᶜinv, ᶜθ_li, ᶜq_tot_nonneg, p)
-            @. p.precomputed.ᶜT′q′ += geo_h * ᶜinv
-        end
+        ᶜinv_q = p.scratch.ᶜtemp_scalar_5
+        ᶜinv_θ = p.scratch.ᶜtemp_scalar_2
+        ᶜinv_θq = p.scratch.ᶜtemp_scalar_6
+        hgrad_invariant!(ᶜinv_q, ᶜq_tot_nonneg, p)
+        hgrad_invariant!(ᶜinv_θ, ᶜθ_li, p)
+        hgrad_cross_invariant!(ᶜinv_θq, ᶜθ_li, ᶜq_tot_nonneg, p)
+        θz_min = CAP.sgs_variance_isentropic_min_dtheta_dz(p.params)
+        ᶜlg = Fields.local_geometry_field(Y.c)
+        @. ᶜinv_q = isentropic_gradient_invariant(
+            ᶜinv_q,
+            ᶜinv_θq,
+            ᶜinv_θ,
+            isentropic_slope_ratio(
+                projected_vector_data(C3, ᶜgradᵥ_q_tot, ᶜlg),
+                projected_vector_data(C3, ᶜgradᵥ_θ_liq_ice, ᶜlg),
+                θz_min,
+            ),
+        )
+        ᶜisen_weight = ᶜsgs_isentropic_weight(Y, p)
+        @. ᶜq′q′ += ᶜisen_weight * (geo_h * ᶜinv_q)
     end
 
     # ----------------------------------------------------------------------

src/parameters/Parameters.jl

@@ -270,6 +270,14 @@ not yet define.
     SGS quadrature, in `[-1, 1]`. Used verbatim by `tq_correlation_model: constant`;
     used as the turbulent-block cross `ρ_turb` (and the fallback where all variance
     blocks vanish) by `tq_correlation_model: diagnosed` [-].
+  - `sgs_variance_isentropic_min_dtheta_dz`: Floor `θz_min` on the stratification in
+    the slope ratio `r = (∂q/∂z) / max(∂θ_li/∂z, θz_min)` of the isentropic geometric
+    variance term; keeps `r` finite and bounded by `|∂q/∂z| / θz_min` in neutral and
+    overturned layers [K m^-1].
+  - `sgs_variance_isentropic_stability_width`: Width of the stability weight
+    `w = (1 + tanh(N²_sat / width)) / 2` of the isentropic term
+    (`isentropic_validity_weight`); `0` (the default) is the hard cut
+    `1[N²_sat > 0]` [s^-2].
   - `cloud_fraction_steepness_scale`: Steepness scale `α` of the
     cloud-fraction/condensate relationship; 1 for exact Gaussian or lognormal SGS
     distributions [-].
@@ -302,6 +310,8 @@ Base.@kwdef struct SGSQuadratureParameters{FT} <: ASQP
     cloud_fraction_floor_release_abs_margin::FT
     cloud_fraction_floor_release_sharpness::FT
     cloud_fraction_floor_residual::FT
+    sgs_variance_isentropic_min_dtheta_dz::FT
+    sgs_variance_isentropic_stability_width::FT
 end
 
 """

src/parameters/create_parameters.jl

@@ -594,6 +594,11 @@ function SGSQuadratureParameters(
         cloud_fraction_floor_release_abs_margin = FT(1),
         cloud_fraction_floor_release_sharpness = FT(1),
         cloud_fraction_floor_residual = FT(1),
+        # Isentropic geometric variance term:
+        # floor on ∂θ_li/∂z [K/m] in the slope ratio.
+        sgs_variance_isentropic_min_dtheta_dz = FT(1e-3),
+        # Smoothing width of the stability weight; 0 = hard cut 1[N²_sat > 0].
+        sgs_variance_isentropic_stability_width = FT(0),
     )
     provisional_present = filter(collect(keys(provisional_defaults))) do name
         haskey(toml_dict.data, string(name))

test/parameter_tests.jl

@@ -121,6 +121,8 @@ end
         @test CAP.sgs_variance_max_rel_std(sq) == FT(0.5)
         @test CAP.sgs_liquid_uniform_fraction(sq) == FT(1)
         @test CAP.sgs_ice_uniform_fraction(sq) == FT(1)
+        @test CAP.sgs_variance_isentropic_min_dtheta_dz(sq) == FT(1e-3)
+        @test CAP.sgs_variance_isentropic_stability_width(sq) == FT(0)
     end
     # A run toml overrides the defaults.
     mktemp() do path, io