fix(Rosenbrock): saturation adjustment on max(S_ice, S_liq)

docs/src/RosenbrockNumerics.md

@@ -34,8 +34,8 @@ substeps from the latent heat of the realized increment.
 
 - **`TendencyLimiter`** — a limiter applied to the realized substep increment.
   - `NoLimiter`.
-  - `EndStateSaturationAdjustment` — scale the increment so the latent-heated end state does not cross ice
-    saturation (see below).
+  - `EndStateSaturationAdjustment` — scale the increment so the latent-heated end state does not cross saturation
+    over its more-supersaturated phase, `max(S_ice, S_liq)` (see below).
 
 Three preset configurations are supported:
 
@@ -87,10 +87,11 @@ physical saturation limit:
   modes and is well-conditioned at any substep size. The exact off-diagonal couplings (which the donor-based
   matrix drops) are retained, so accuracy at cold, supersaturated cells is better than the donor scheme.
 - **`EndStateSaturationAdjustment`** scales the substep increment by the largest `s ∈ [0, 1]` for which the
-  latent-heated end state stays at or above ice saturation. It acts only on cells that begin at or above
-  saturation (a subsaturated, evaporating or sublimating cell cannot over-deposit, so its increment is returned
-  unchanged). It is a no-op at fine substeps and engages only when the full step would cross saturation. The
-  bisection count is set from the float precision.
+  latent-heated end state stays at or above saturation over its more-supersaturated phase (`max(S_ice, S_liq)`;
+  see the next section). It acts only on cells that begin at or above saturation (a subsaturated, evaporating or
+  sublimating cell cannot over-deposit, so its increment is returned unchanged). It is a no-op at fine substeps
+  and engages only when the full step would cross saturation. The bisection count is set from the float
+  precision.
 
 Both pieces are required: zeroing the growth diagonal alone leaves the explicit growth unbounded at the coarsest
 single-substep steps, and the saturation adjustment supplies the missing nonlinear bound. Together they make the
@@ -101,6 +102,44 @@ exact scheme robust across the resolved time-step envelope.
     precipitation at coarse time steps. Two or more substeps recover accurate precipitation; the saturation
     adjustment is then rarely active.
 
+## The mixed-phase saturation criterion
+
+A cell has two saturation thresholds, over liquid and over ice, and the condensation and deposition processes draw
+on a single shared vapor reservoir. `EndStateSaturationAdjustment` limits the increment on the more-supersaturated
+phase, `max(S_ice, S_liq)`: it keeps the latent-heated end state at or above the saturation of whichever phase
+carries the larger supersaturation. Equivalently the vapor floor is the lower of the two saturation specific
+humidities, since the smaller `q_sat` is the larger supersaturation.
+
+The two saturation curves cross at the freezing point (panel (c) below). Below freezing `q_sat_ice < q_sat_liq`, so
+ice carries the larger supersaturation and the criterion binds on ice — reducing exactly to the ice-saturation
+limit that bounds the ice-growth instability, so the cold behaviour is unchanged. Above freezing the curves swap
+and the criterion binds on liquid; an ice-only limit there would stop vapor depletion early and leave the cell
+supersaturated over liquid, suppressing warm cloud, which binding on the more-supersaturated phase avoids.
+
+A single 0-D parcel that exchanges vapor with cloud liquid and cloud ice only (no collection or precipitation)
+isolates the mixed-phase physics behind the two thresholds. Cloud liquid is kinetically fast and cloud ice slow, so
+while liquid is present it pins the vapor near water saturation (panel (a), `S_liq ≈ 0`); the parcel stays
+supersaturated over ice (`S_ice > 0`) and ice deposits, drawing mass from the evaporating liquid (panel (b)) — the
+Wegener–Bergeron–Findeisen transfer. Only once the liquid is exhausted does the vapor relax to ice saturation
+(`S_ice → 0`). The drawdown from water saturation to ice saturation is therefore inherently a multi-step process at
+the resolved time step.
+
+```@example
+include("plots/SatAdjustmentWBF_plots.jl")
+```
+![](SatAdjustmentWBF.svg)
+
+Binding on `max(S_ice, S_liq)` is a single shared scalar on the whole increment, which is what keeps the limiter
+stable: scaling processes independently breaks the coupling between paired transfers and the shared vapor draw. The
+criterion binds the correct phase wherever a single phase grows from vapor — ice only (the cold deposition cells,
+where it reduces to the ice limit) or liquid only (warm cells). When both phases are supersaturated below freezing,
+a vigorous mixed-phase updraft core, the criterion binds on the ice floor (the more-supersaturated phase there) and
+so condenses the co-present, still-growing liquid past its own saturation in a single step rather than transferring
+it gradually. A fully per-phase floor — binding on the first growing phase to saturate and letting the slower phase
+relax over subsequent steps — would represent the gradual transfer more faithfully and is a candidate refinement.
+The distinction is inactive at fine substeps, where the limiter rarely engages, and does not affect the cold
+instability resolution, which is set by the ice floor.
+
 ## Approaches that did not resolve the instability
 
 These were tried and are not part of the supported framework.

docs/src/plots/SatAdjustmentWBF_plots.jl

@@ -0,0 +1,85 @@
+import CloudMicrophysics.ThermodynamicsInterface as TDI
+using CairoMakie
+
+# Thermodynamics handles (a representative air density is fixed throughout).
+const FT = Float64
+const tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
+const Tf = TDI.T_freeze(tps)
+const cp_d = TDI.TD.Parameters.cp_d(tps)
+const Lv = TDI.TD.Parameters.LH_v0(tps)
+const Ls = TDI.TD.Parameters.LH_s0(tps)
+const RHO = FT(0.8)
+
+# Saturation specific humidity over each phase, recovered from the supersaturation interface
+# (S = q_vap/q_sat - 1 evaluated at zero condensate, so q_sat = q_ref / (1 + S)).
+const QREF = FT(1e-2)
+q_sat_liq(T) = QREF / (1 + TDI.supersaturation_over_liquid(tps, QREF, FT(0), FT(0), RHO, T))
+q_sat_ice(T) = QREF / (1 + TDI.supersaturation_over_ice(tps, QREF, FT(0), FT(0), RHO, T))
+
+# Simplified 0-D vapor-exchange-only model: cloud liquid and cloud ice exchange mass with the shared vapor by
+# relaxation toward their own saturation, with liquid kinetically fast and ice slow. Collection and precipitation
+# are excluded so the Wegener-Bergeron-Findeisen transfer is isolated.
+function wbf_box(; T0, q_liq0, q_ice0, tau_liq, tau_ice, t_end, dt)
+    q_tot = q_sat_liq(T0) + q_liq0 + q_ice0           # start at liquid saturation
+    q_liq, q_ice, T = q_liq0, q_ice0, T0
+    n = Int(round(t_end / dt))
+    ts = zeros(n + 1)
+    QL = zeros(n + 1)
+    QI = zeros(n + 1)
+    SL = zeros(n + 1)
+    SI = zeros(n + 1)
+    for k in 0:n
+        q_vap = q_tot - q_liq - q_ice
+        ts[k + 1] = k * dt
+        QL[k + 1] = q_liq
+        QI[k + 1] = q_ice
+        SL[k + 1] = q_vap / q_sat_liq(T) - 1
+        SI[k + 1] = q_vap / q_sat_ice(T) - 1
+        cond = (q_liq > 0 || q_vap > q_sat_liq(T)) ? (q_vap - q_sat_liq(T)) / tau_liq : 0.0
+        dep = (q_ice > 0 || q_vap > q_sat_ice(T)) ? (q_vap - q_sat_ice(T)) / tau_ice : 0.0
+        dq_liq = max(cond * dt, -q_liq)               # cannot evaporate more liquid than present
+        dq_ice = max(dep * dt, -q_ice)
+        q_liq += dq_liq
+        q_ice += dq_ice
+        T += (Lv * dq_liq + Ls * dq_ice) / cp_d
+    end
+    (; ts, QL, QI, SL, SI)
+end
+
+colors = Makie.wong_colors()
+fig = Figure(size = (1200, 360))
+
+r = wbf_box(; T0 = Tf - 15, q_liq0 = 1.5e-3, q_ice0 = 1e-5, tau_liq = 3.0, tau_ice = 50.0, t_end = 700.0, dt = 0.5)
+
+# Panel (a): supersaturation -- the WBF signature. While liquid is present the fast liquid pins the vapor at water
+# saturation (S_liq ~ 0), so the parcel stays supersaturated over ice (S_ice > 0) and ice deposits; once the liquid
+# is gone the vapor relaxes to ice saturation (S_ice -> 0, S_liq < 0).
+axa = Axis(fig[1, 1], xlabel = "time [s]", ylabel = "supersaturation",
+    title = "(a) WBF supersaturation, T = T_f − 15 K")
+lines!(axa, r.ts, r.SL, color = colors[1], linewidth = 2, label = "S_liq")
+lines!(axa, r.ts, r.SI, color = colors[2], linewidth = 2, label = "S_ice")
+hlines!(axa, [0.0], color = :gray, linestyle = :dot)
+axislegend(axa, position = :rt, framevisible = false)
+
+# Panel (b): the mass transfer -- liquid evaporates and ice grows at sustained S_ice > 0 (the WBF transfer).
+axb = Axis(fig[1, 2], xlabel = "time [s]", ylabel = "specific content [g/kg]",
+    title = "(b) WBF mass transfer")
+lines!(axb, r.ts, r.QL .* 1e3, color = colors[1], linewidth = 2, label = "cloud liquid")
+lines!(axb, r.ts, r.QI .* 1e3, color = colors[2], linewidth = 2, label = "cloud ice")
+axislegend(axb, position = :rc, framevisible = false)
+
+# Panel (c): the thermodynamic basis of the max(S_ice, S_liq) criterion. The two saturation curves cross at the
+# freezing point; the limiter holds the vapor at the saturation of the more-supersaturated phase = the lower of the
+# two q_sat curves (smaller q_sat is larger S). The binding phase changes from ice to liquid at T_f.
+Ts = collect(range(Tf - 25, Tf + 8, length = 200))
+qsl = q_sat_liq.(Ts) .* 1e3
+qsi = q_sat_ice.(Ts) .* 1e3
+axc = Axis(fig[1, 3], xlabel = "temperature [K]", ylabel = "saturation specific humidity [g/kg]",
+    title = "(c) max-criterion floor")
+lines!(axc, Ts, qsl, color = colors[1], linewidth = 2, label = "q_sat over liquid")
+lines!(axc, Ts, qsi, color = colors[2], linewidth = 2, label = "q_sat over ice")
+lines!(axc, Ts, min.(qsl, qsi), color = colors[6], linewidth = 4, label = "vapor floor = min(q_sat)")
+vlines!(axc, [Tf], color = :gray, linestyle = :dash)
+axislegend(axc, position = :lt, framevisible = false)
+
+save("SatAdjustmentWBF.svg", fig)

src/BMT_rosenbrock.jl

@@ -672,10 +672,10 @@ end
     _apply_limiter(limiter, x, d, ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps)
 
 Limit the substep increment `d` at state `x`. [`NoLimiter`](@ref) returns `d`.
-[`EndStateSaturationAdjustment`](@ref), for a cell at or above ice saturation
-whose full-increment end state would drop below it, scales `d` by `s ∈ [0, 1]`
-to keep the latent-heated end state at or above ice saturation; otherwise it
-returns `d`.
+[`EndStateSaturationAdjustment`](@ref), for a cell at or above saturation over its
+more-supersaturated phase whose full-increment end state would drop below it, scales
+`d` by `s ∈ [0, 1]` to keep that latent-heated end state at or above saturation over
+that phase; otherwise it returns `d`.
 """
 @inline _apply_limiter(::NoLimiter, x, d, ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps) = d
 
@@ -688,24 +688,23 @@ precision of `FT`.
 @inline _saturation_bisection_count(::Type{FT}) where {FT} = ceil(Int, -log2(eps(FT)))
 
 """
-    _saturation_bisection(Sice, latent, x, d, Tsub)
+    _saturation_bisection(Ssat, latent, x, d, Tsub)
 
-Scale the increment `d` at state `x` so the latent-heated end state stays at or
-above ice saturation, for a state that begins at or above it; return `d`
-unchanged otherwise. `Sice(x, T)` and `latent(d)` close over the substep
-context.
+Scale the increment `d` at state `x` so the latent-heated end state keeps
+`Ssat >= 0`, for a state that begins with `Ssat >= 0`; return `d` unchanged
+otherwise. `Ssat(x, T)` and `latent(d)` close over the substep context.
 """
 @inline function _saturation_bisection(
-    Sice::FS, latent::FL, x::SA.StaticVector{N, FT}, d, Tsub,
+    Ssat::FS, latent::FL, x::SA.StaticVector{N, FT}, d, Tsub,
 ) where {FS, FL, N, FT}
     xf = max.(x .+ d, 0)
-    if Sice(x, Tsub) >= 0 && Sice(xf, Tsub + latent(xf .- x)) < 0
+    if Ssat(x, Tsub) >= 0 && Ssat(xf, Tsub + latent(xf .- x)) < 0
         lo = zero(FT)
         hi = one(FT)
         for _ in 1:_saturation_bisection_count(FT)
             s = (lo + hi) / 2
             xs = max.(x .+ s .* d, 0)
-            if Sice(xs, Tsub + latent(xs .- x)) >= 0
+            if Ssat(xs, Tsub + latent(xs .- x)) >= 0
                 lo = s
             else
                 hi = s
@@ -720,19 +719,24 @@ end
     x::MicroState1M{FT}, d::MicroState1M{FT},
     ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps,
 ) where {FT}
-    Sice(xx, TT) =
-        TDI.supersaturation_over_ice(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_icl + xx.q_sno, ρ, TT)
+    Ssat(xx, TT) = max(
+        TDI.supersaturation_over_ice(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_icl + xx.q_sno, ρ, TT),
+        TDI.supersaturation_over_liquid(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_icl + xx.q_sno, ρ, TT),
+    )
     latent(dd) = Lv_over_cp * (dd.q_lcl + dd.q_rai) + Ls_over_cp * (dd.q_icl + dd.q_sno)
-    return _saturation_bisection(Sice, latent, x, d, Tsub)
+    return _saturation_bisection(Ssat, latent, x, d, Tsub)
 end
 
 @inline function _apply_limiter(::EndStateSaturationAdjustment,
     x::MicroState2MP3{FT}, d::MicroState2MP3{FT},
     ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps,
 ) where {FT}
-    Sice(xx, TT) = TDI.supersaturation_over_ice(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_ice, ρ, TT)
+    Ssat(xx, TT) = max(
+        TDI.supersaturation_over_ice(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_ice, ρ, TT),
+        TDI.supersaturation_over_liquid(tps, q_tot, xx.q_lcl + xx.q_rai, xx.q_ice, ρ, TT),
+    )
     latent(dd) = Lv_over_cp * (dd.q_lcl + dd.q_rai) + Ls_over_cp * dd.q_ice
-    return _saturation_bisection(Sice, latent, x, d, Tsub)
+    return _saturation_bisection(Ssat, latent, x, d, Tsub)
 end
 
 "Exact ForwardDiff Jacobian provider for [`_rosenbrock_average_1m`](@ref)."

src/BulkMicrophysicsTendencies.jl

@@ -206,9 +206,10 @@ struct NoLimiter <: TendencyLimiter end
 """
     EndStateSaturationAdjustment <: TendencyLimiter
 
-Scale a substep increment so the latent-heated end state stays at or above ice
-saturation, for cells that begin at or above saturation. Derived and analyzed in
-the P3 numerics documentation.
+Scale a substep increment so the latent-heated end state stays at or above
+saturation over its more-supersaturated phase (ice when ice dominates, liquid when
+liquid dominates), for cells that begin at or above saturation. Derived and analyzed
+in the Rosenbrock substepping documentation.
 """
 struct EndStateSaturationAdjustment <: TendencyLimiter end
 

test/rosenbrock_framework_tests.jl

@@ -77,6 +77,40 @@ function test_framework_1m(FT)
         end
     end
 
+    @testset "EndStateSaturationAdjustment limits the more-supersaturated phase ($FT)" begin
+        Lv = TDI.TD.Parameters.LH_v0(tps) / TDI.TD.Parameters.cp_d(tps)
+        Ls = TDI.TD.Parameters.LH_s0(tps) / TDI.TD.Parameters.cp_d(tps)
+        S_liq(x, T, ρ, qt) = TDI.supersaturation_over_liquid(tps, qt, x[1] + x[3], x[2] + x[4], ρ, T)
+        S_ice(x, T, ρ, qt) = TDI.supersaturation_over_ice(tps, qt, x[1] + x[3], x[2] + x[4], ρ, T)
+        function end_state_S(mode, ρ, T, qt, x, Δt, nsub)
+            d = net_vec_1m(
+                BMT.bulk_microphysics_tendencies(
+                    mode, BMT.Microphysics1Moment(), mp, tps, ρ, T, qt, x..., Δt, nsub,
+                ),
+            )
+            xe = x .+ Δt .* d
+            Te = T + Lv * ((xe[1] - x[1]) + (xe[3] - x[3])) + Ls * ((xe[2] - x[2]) + (xe[4] - x[4]))
+            (S_liq(xe, Te, ρ, qt), S_ice(xe, Te, ρ, qt))
+        end
+        exact = BMT.rosenbrock_exact()
+        unlimited = BMT.RosenbrockAverage(BMT.ExactJacobian(), BMT.ExplicitGrowthDiagonal(), BMT.NoLimiter())
+        tol = FT == Float64 ? FT(1e-3) : FT(3e-2)
+        # warm, liquid-supersaturated, coarse single step
+        let ρ = FT(1.0), T = T_frz + FT(8), qt = FT(0.022), x = FT[1e-4, 0, 1e-4, 0], Δt = FT(240)
+            Sl_exact, _ = end_state_S(exact, ρ, T, qt, x, Δt, 1)
+            Sl_unlim, _ = end_state_S(unlimited, ρ, T, qt, x, Δt, 1)
+            @test abs(Sl_exact) < tol
+            @test Sl_unlim < -tol
+        end
+        # cold, ice-supersaturated, coarse single step
+        let ρ = FT(0.9), T = T_frz - FT(20), qt = FT(2.0e-3), x = FT[0, 1e-4, 0, 1e-4], Δt = FT(240)
+            _, Si_exact = end_state_S(exact, ρ, T, qt, x, Δt, 1)
+            _, Si_unlim = end_state_S(unlimited, ρ, T, qt, x, Δt, 1)
+            @test abs(Si_exact) < tol
+            @test Si_unlim < -tol
+        end
+    end
+
     @testset "Verbose(rosenbrock_donor()) per-process sums to net ($FT)" begin
         rtol = FT == Float64 ? FT(1e-10) : FT(1e-4)
         for r in regimes, nsub in (1, 4, 16)