From 5a9a3233c4b079362f6805687c2287b20d8e6182 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 11 Nov 2025 10:09:35 -0700 Subject: [PATCH 1/7] Build abstraction for mixed phase surfaces --- src/Thermodynamics/vapor_saturation.jl | 83 ++++++++++++++++++++++---- 1 file changed, 70 insertions(+), 13 deletions(-) diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index 7e94aa1c..a166c488 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -1,41 +1,60 @@ """ $(TYPEDSIGNATURES) -Compute the saturation vapor pressure ``pᵛ⁺`` over a planar surface -composed of the "``β``-phase" (liquid, or ice) +Compute the [saturation vapor pressure](https://en.wikipedia.org/wiki/Vapor_pressure) +``pᵛ⁺`` over a surface labeled ``β`` (for example, a planar liquid surface, or curved ice surface) using the Clausius-Clapeyron relation, ```math dpᵛ⁺ / dT = pᵛ⁺ ℒᵝ(T) / (Rᵛ T^2) , ``` -where the latent heat is ``ℒᵝ(T) = ℒᵝ(T=0) + Δcᵝ T``, with ``Δcᵝ ≡ (cᵖᵛ - cᵝ)`` . +where the temperature-dependent latent heat of the surfaceis ``ℒᵝ(T)``. -The saturation vapor pressure ``pᵛ⁺`` is obtained after integrating the above from -the triple point, i.e., ``p(Tᵗʳ) = pᵗʳ`` to get +Using a model for the latent heat that is linear in temperature, eg ```math -pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵝ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒᵝ(T=0) / Rᵛ \\right ] . +ℒᵝ = ℒᵝ₀ + Δcᵝ T, ``` -Note that latent heat ``ℒᵝ(T=0)`` is the difference between the enthalpy of water vapor -and the phase ``β`` when the temperature is ``T = 0``K, or absolute zero. +where ``ℒᵝ₀ ≡ ℒᵝ(T=0)`` is the latent heat at absolute zero and +``Δcᵝ ≡ (cᵖᵛ - cᵝ)`` is the constant difference between the vapor specific heat +and the specific heat of phase ``β``. -We define the latent heat using its value ``ℒᵝᵣ = ℒᵝ(T=Tᵣ)`` at the "energy reference temperature" -``T = Tᵣ`` (usually ``Tᵣ ≡ 273.15``K ``= 0^∘``C), such that +Note that we typically parameterize the latent heat interms of a reference +temperature ``T = Tᵣ`` that is well above absolute zero. In that case, +the latent heat is written + +```math +ℒᵝ = ℒᵝᵣ + Δcᵝ (T - Tᵣ), +\qquad \text{and} \qquad +ℒᵝ₀ = ℒᵝᵣ - Δcᵝ Tᵣ . +``` + +Integrating the Clausius-Clapeyron relation with a temperature-linear latent heat model, +from the triple point pressure and temperature ``(pᵗʳ, Tᵗʳ)`` to pressure ``pᵛ⁺`` +and temperature ``T``, we obtain + +```math +log(pᵛ⁺ / pᵗʳ) = - ℒᵝ₀ / (Rᵛ T) + ℒᵝ₀ / (Rᵛ Tᵗʳ) + log(Δcᵝ / Rᵛ * T / Tᵗʳ) +``` + +Which then becomes ```math -ℒᵝ(T) = ℒᵝᵣ + Δcᵝ (T - Tᵣ), \\quad \\text{and} \\quad ℒᵝ(T=0) = ℒᵝᵣ - Δcᵝ Tᵣ``. +pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵝ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒᵝ₀ / Rᵛ \\right ] . ``` + """ -@inline function saturation_vapor_pressure(T, thermo, phase::CondensedPhase) - ℒᵣ = phase.reference_latent_heat # at thermo.energy_reference_temperature +@inline function saturation_vapor_pressure(T, thermo, surface) + ℒ₀ = absolute_zero_latent_heat(thermo, surface) Tᵣ = thermo.energy_reference_temperature Tᵗʳ = thermo.triple_point_temperature pᵗʳ = thermo.triple_point_pressure Rᵛ = vapor_gas_constant(thermo) + cᵖᵛ = thermo.vapor.heat_capacity cᵝ = phase.heat_capacity Δcᵝ = cᵖᵛ - cᵝ @@ -43,6 +62,44 @@ We define the latent heat using its value ``ℒᵝᵣ = ℒᵝ(T=Tᵣ)`` at the return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * (ℒᵣ - Δcᵝ * Tᵣ) / Rᵛ) end +@inline function specific_heat_difference(thermo, phase::CondensedPhase) + cᵖᵛ = thermo.vapor.heat_capacity + cᵝ = phase.heat_capacity + return cᵖᵛ - cᵝ +end + +@inline function absolute_zero_latent_heat(thermo, phase::CondensedPhase) + ℒᵣ = phase.reference_latent_heat # at thermo.energy_reference_temperature + Δcᵝ = specific_heat_difference(thermo, phase) + return ℒᵣ - Δcᵝ * Tᵣ +end + +struct PlanarLiquidSurface end +struct PlanarIceSurface end + +struct PlanarMixedPhaseSurface{FT} + liquid_fraction :: FT +end + +@inline specific_heat_difference(thermo, ::PlanarIceSurface) = specific_heat_difference(thermo, thermo.ice) +@inline absolute_zero_latent_heat(thermo, ::PlanarIceSurface) = absolute_zero_latent_heat(thermo, thermo.ice) +@inline absolute_zero_latent_heat(thermo, ::PlanarLiquidSurface) = absolute_zero_latent_heat(thermo, thermo.liquid) +@inline specific_heat_difference(thermo, ::PlanarLiquidSurface) = specific_heat_difference(thermo, thermo.liquid) + +@inline function specific_heat_difference(thermo, surf::PlanarMixedPhaseSurface) + Δcˡ = specific_heat_difference(thermo, thermo.liquid) + Δcⁱ = specific_heat_difference(thermo, thermo.ice) + λ = surf.liquid_fraction + return λ * Δcˡ + (1 - λ) * Δcⁱ +end + +@inline function absolute_zero_latent_heat(thermo, surf::PlanarMixedPhaseSurface) + ℒˡ₀ = absolute_zero_latent_heat(thermo, thermo.liquid) + ℒⁱ₀ = absolute_zero_latent_heat(thermo, thermo.ice) + λ = surf.liquid_fraction + return λ * ℒˡ₀ + (1 - λ) * ℒⁱ₀ +end + """ $(TYPEDSIGNATURES) From fc17007d988e33213bbc9ce3a4596154b1e1898f Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 11 Nov 2025 13:43:10 -0700 Subject: [PATCH 2/7] add tests --- src/Thermodynamics/vapor_saturation.jl | 12 ++---- test/saturation_vapor_pressure.jl | 60 ++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 8 deletions(-) create mode 100644 test/saturation_vapor_pressure.jl diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index a166c488..c39e1e2e 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -27,7 +27,7 @@ the latent heat is written ```math ℒᵝ = ℒᵝᵣ + Δcᵝ (T - Tᵣ), -\qquad \text{and} \qquad +\\qquad \\text{and} \\qquad ℒᵝ₀ = ℒᵝᵣ - Δcᵝ Tᵣ . ``` @@ -48,18 +48,13 @@ pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵝ / Rᵛ} \\exp \ """ @inline function saturation_vapor_pressure(T, thermo, surface) ℒ₀ = absolute_zero_latent_heat(thermo, surface) - Tᵣ = thermo.energy_reference_temperature + Δcᵝ = specific_heat_difference(thermo, surface) Tᵗʳ = thermo.triple_point_temperature pᵗʳ = thermo.triple_point_pressure Rᵛ = vapor_gas_constant(thermo) - - cᵖᵛ = thermo.vapor.heat_capacity - cᵝ = phase.heat_capacity - Δcᵝ = cᵖᵛ - cᵝ - - return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * (ℒᵣ - Δcᵝ * Tᵣ) / Rᵛ) + return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * ℒ₀ / Rᵛ) end @inline function specific_heat_difference(thermo, phase::CondensedPhase) @@ -71,6 +66,7 @@ end @inline function absolute_zero_latent_heat(thermo, phase::CondensedPhase) ℒᵣ = phase.reference_latent_heat # at thermo.energy_reference_temperature Δcᵝ = specific_heat_difference(thermo, phase) + Tᵣ = thermo.energy_reference_temperature return ℒᵣ - Δcᵝ * Tᵣ end diff --git a/test/saturation_vapor_pressure.jl b/test/saturation_vapor_pressure.jl new file mode 100644 index 00000000..bcc8e885 --- /dev/null +++ b/test/saturation_vapor_pressure.jl @@ -0,0 +1,60 @@ +using Breeze +using Test + +using Breeze.Thermodynamics: + ThermodynamicConstants, + saturation_vapor_pressure, + PlanarLiquidSurface, + PlanarIceSurface, + PlanarMixedPhaseSurface, + absolute_zero_latent_heat, + specific_heat_difference, + vapor_gas_constant + +function reference_mixed_surface_pressure(T, thermo, λ) + ℒˡ₀ = absolute_zero_latent_heat(thermo, thermo.liquid) + ℒⁱ₀ = absolute_zero_latent_heat(thermo, thermo.ice) + Δcˡ = specific_heat_difference(thermo, thermo.liquid) + Δcⁱ = specific_heat_difference(thermo, thermo.ice) + + ℒ₀ = λ * ℒˡ₀ + (one(λ) - λ) * ℒⁱ₀ + Δcᵝ = λ * Δcˡ + (one(λ) - λ) * Δcⁱ + + Tᵗʳ = thermo.triple_point_temperature + pᵗʳ = thermo.triple_point_pressure + Rᵛ = vapor_gas_constant(thermo) + + return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((one(T) / Tᵗʳ - one(T) / T) * ℒ₀ / Rᵛ) +end + +@testset "Saturation vapor pressure surfaces [$FT]" for FT in (Float32, Float64) + thermo = ThermodynamicConstants(FT) + Tᵗʳ = thermo.triple_point_temperature + temperatures = FT.((Tᵗʳ * FT(0.9), Tᵗʳ, Tᵗʳ * FT(1.1))) + + liquid_surface = PlanarLiquidSurface() + ice_surface = PlanarIceSurface() + rtol = FT === Float64 ? 1e-12 : FT(1e-5) + + @testset "Planar homogeneous surfaces" begin + for T in temperatures + pˡ = saturation_vapor_pressure(T, thermo, thermo.liquid) + pⁱ = saturation_vapor_pressure(T, thermo, thermo.ice) + + @test saturation_vapor_pressure(T, thermo, liquid_surface) ≈ pˡ rtol=rtol + @test saturation_vapor_pressure(T, thermo, ice_surface) ≈ pⁱ rtol=rtol + end + end + + @testset "Planar mixed-phase surfaces" begin + for λ in (zero(FT), FT(0.25), FT(0.5), FT(0.75), one(FT)) + surface = PlanarMixedPhaseSurface(λ) + for T in temperatures + p_surface = saturation_vapor_pressure(T, thermo, surface) + p_reference = reference_mixed_surface_pressure(T, thermo, λ) + + @test p_surface ≈ p_reference rtol=rtol + end + end + end +end From eed15866b3d6616db01ea90179cd51d987d8f9b7 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 11 Nov 2025 15:25:41 -0700 Subject: [PATCH 3/7] doctest --- src/Thermodynamics/vapor_saturation.jl | 49 +++++++++++++++++++++++--- 1 file changed, 44 insertions(+), 5 deletions(-) diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index c39e1e2e..7c535131 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -44,7 +44,6 @@ Which then becomes ```math pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵝ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒᵝ₀ / Rᵛ \\right ] . ``` - """ @inline function saturation_vapor_pressure(T, thermo, surface) ℒ₀ = absolute_zero_latent_heat(thermo, surface) @@ -100,17 +99,57 @@ end $(TYPEDSIGNATURES) Compute the saturation specific humidity for a gas at temperature `T`, total -density `ρ`, `thermo`dynamics, and `phase` via: +density `ρ`, `thermo`dynamics, and over `surface` via: ```math qᵛ⁺ = pᵛ⁺ / (ρ Rᵛ T) , ``` -where ``pᵛ⁺`` is the [`saturation_vapor_pressure`](@ref), ``ρ`` is total density, +where ``pᵛ⁺`` is the [`saturation_vapor_pressure`](@ref) over `surface`, ``ρ`` is total density, and ``Rᵛ`` is the specific gas constant for water vapor. + +# Examples + +First we compute the saturation specific humidity over a liquid surface: + +```jldoctest saturation +using Breeze +using Breeze.Thermodynamics: PlanarLiquidSurface, PlanarIceSurface, PlanarMixedPhaseSurface + +thermo = ThermodynamicConstants() +T = 288.0 # Room temperature (K) +p = 101325.0 # Mean sea-level pressure +Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo) +q = zero(Breeze.Thermodynamics.MoistureMassFractions{Float64}) +ρ = Breeze.Thermodynamics.density(p, T, q, thermo) +qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, PlanarLiquidSurface()) + +# output +0.010359995391195264 +``` + +Note, this is slightly smaller than the saturation specific humidity over an ice surface: + +```jldoctest saturation +qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, PlanarIceSurface()) +0.011945100768555072 +``` + +If a medium contains a mixture of 40% water and 60% ice that has (somehow) acquired +thermodynamic equilibrium, we can compute the saturation specific humidity +over the mixed phase surface, + +```jldoctest saturation +mixed_surface = PlanarMixedPhaseSurface(0.4) +qᵛ⁺ᵐ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, mixed_surface) + +# output +0.01128386068542303 +``` + """ -@inline function saturation_specific_humidity(T, ρ, thermo::ThermodynamicConstants, phase::CondensedPhase) - pᵛ⁺ = saturation_vapor_pressure(T, thermo, phase) +@inline function saturation_specific_humidity(T, ρ, thermo::ThermodynamicConstants, surface) + pᵛ⁺ = saturation_vapor_pressure(T, thermo, surface) Rᵛ = vapor_gas_constant(thermo) return pᵛ⁺ / (ρ * Rᵛ * T) end From ba86b0377448416f29d922f489b789c9e8e6b6f5 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 11 Nov 2025 15:25:50 -0700 Subject: [PATCH 4/7] add zero --- src/Thermodynamics/thermodynamics_constants.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Thermodynamics/thermodynamics_constants.jl b/src/Thermodynamics/thermodynamics_constants.jl index 3a685fef..4e43517a 100644 --- a/src/Thermodynamics/thermodynamics_constants.jl +++ b/src/Thermodynamics/thermodynamics_constants.jl @@ -251,6 +251,8 @@ struct MoistureMassFractions{FT} ice :: FT end +Base.zero(::Type{MoistureMassFractions{FT}}) where FT = MoistureMassFractions(zero(FT), zero(FT), zero(FT)) + function Base.summary(q::MoistureMassFractions{FT}) where FT return string("MoistureMassFractions{$FT}(vapor=", prettysummary(q.vapor), ", liquid=", prettysummary(q.liquid), ", ice=", prettysummary(q.ice), ")") From dc0dc66dfa691b2af34234adf2a2e6e46d3f904e Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 11 Nov 2025 15:26:04 -0700 Subject: [PATCH 5/7] update clausuis claperyon docs --- docs/src/thermodynamics.md | 195 +++++++++++++++++++++++++++++++------ 1 file changed, 164 insertions(+), 31 deletions(-) diff --git a/docs/src/thermodynamics.md b/docs/src/thermodynamics.md index a0ea714a..00ba732d 100644 --- a/docs/src/thermodynamics.md +++ b/docs/src/thermodynamics.md @@ -392,34 +392,70 @@ cᵖᵐ - cᵖᵈ ``` -## The Clausius--Clapeyron relation and saturation specific humidity +## The Clausius--Clapeyron relation and saturation vapor pressure The [Clausius--Clapeyron relation](https://en.wikipedia.org/wiki/Clausius%E2%80%93Clapeyron_relation) -for an ideal gas +for an ideal gas describes how saturation vapor pressure changes with temperature: ```math \frac{\mathrm{d} pᵛ⁺}{\mathrm{d} T} = \frac{pᵛ⁺ ℒ^β(T)}{Rᵛ T^2} , ``` -where ``pᵛ⁺`` is saturation vapor pressure, ``T`` is temperature, ``Rᵛ`` is the specific -gas constant for vapor, ``ℒ^β(T)`` is the latent heat of the transition from vapor to the -``β`` phase (e.g., ``β = l`` for vapor → liquid and ``β = i`` for vapor to ice). +where ``pᵛ⁺`` is saturation vapor pressure over a surface of condensed phase ``β``, +``T`` is temperature, ``Rᵛ`` is the specific gas constant for vapor, and +``ℒ^β(T)`` is the latent heat of the phase transition from vapor to phase ``β``. +For atmospheric moist air, the relevant condensed phases are liquid water (``β = l``) +and ice (``β = i``). + +### Temperature-dependent latent heat For a thermodynamic formulation that uses constant (i.e. temperature-independent) specific -heats, the latent heat of a phase transition is linear in temperature. -For example, for phase change from vapor to liquid, +heats, the latent heat of a phase transition is linear in temperature: + +```math +ℒ^β(T) = ℒ^β_0 + \Delta c^β \, T , +``` + +where ``ℒ^β_0 ≡ ℒ^β(T=0)`` is the latent heat at absolute zero and +``\Delta c^β ≡ c_p^v - c^β`` is the constant difference between the vapor specific heat +capacity at constant pressure and the specific heat capacity of the condensed phase ``β``. + +Note that we typically parameterize the latent heat in terms of a reference +temperature ``T_r`` that is well above absolute zero. In that case, +the latent heat is written + +```math +ℒ^β(T) = ℒ^β_r + \Delta c^β (T - T_r), \qquad \text{and} \qquad +ℒ^β_0 = ℒ^β_r - \Delta c^β T_r , +``` + +where ``ℒ^β_r`` is the latent heat at the reference temperature ``T_r``. + +### Integration of the Clausius-Clapeyron relation + +To find the saturation vapor pressure as a function of temperature, we integrate +the Clausius-Clapeyron relation with the temperature-linear latent heat model +from the triple point pressure and temperature ``(p^{tr}, T^{tr})`` to a generic +pressure ``pᵛ⁺`` and temperature ``T``: + +```math +\int_{p^{tr}}^{pᵛ⁺} \frac{\mathrm{d} p}{p} = \int_{T^{tr}}^{T} \frac{ℒ^β_0 + \Delta c^β T'}{Rᵛ T'^2} \, \mathrm{d} T' . +``` + +Evaluating the integrals yields ```math -ℒˡ(T) = ℒˡ(T=0) + \big ( \underbrace{cᵖᵛ - cˡ}_{≡Δcˡ} \big ) T , +\log\left(\frac{pᵛ⁺}{p^{tr}}\right) = -\frac{ℒ^β_0}{Rᵛ T} + \frac{ℒ^β_0}{Rᵛ T^{tr}} + \frac{\Delta c^β}{Rᵛ} \log\left(\frac{T}{T^{tr}}\right) . ``` -where ``ℒˡ(T=0)`` is the latent heat at absolute zero, ``T = 0 \; \mathrm{K}``. -By integrating from the triple-point temperature ``Tᵗʳ`` for which ``p(Tᵗʳ) = pᵗʳ``, we get +Exponentiating both sides gives the closed-form solution: ```math -pᵛ⁺(T) = pᵗʳ \left ( \frac{T}{Tᵗʳ} \right )^{Δcˡ / Rᵛ} \exp \left [ \frac{ℒˡ(T=0)}{Rᵛ} \left (\frac{1}{Tᵗʳ} - \frac{1}{T} \right ) \right ] . +pᵛ⁺(T) = p^{tr} \left ( \frac{T}{T^{tr}} \right )^{\Delta c^β / Rᵛ} \exp \left [ \frac{ℒ^β_0}{Rᵛ} \left (\frac{1}{T^{tr}} - \frac{1}{T} \right ) \right ] . ``` +### Example: liquid water and ice parameters + Consider parameters for liquid water, ```@example thermo @@ -427,17 +463,63 @@ using Breeze.Thermodynamics: CondensedPhase liquid_water = CondensedPhase(reference_latent_heat=2500800, heat_capacity=4181) ``` -or water ice, +and water ice, ```@example thermo water_ice = CondensedPhase(reference_latent_heat=2834000, heat_capacity=2108) ``` -The saturation vapor pressure is +These represent the latent heat of vaporization at the reference temperature and +the specific heat capacity of each condensed phase. We can compute the specific heat +difference ``\Delta c^β`` for liquid water: + +```@example thermo +using Breeze.Thermodynamics: vapor_gas_constant +cᵖᵛ = thermo.vapor.heat_capacity +cˡ = thermo.liquid.heat_capacity +Δcˡ = cᵖᵛ - cˡ +``` + +This difference ``\Delta c^l ≈`` $(round(1885 - 4181, digits=1)) J/(kg⋅K) is negative because +water vapor has a lower heat capacity than liquid water. + +### Mixed-phase saturation vapor pressure + +In atmospheric conditions near the freezing point, condensate may exist as a mixture of +liquid and ice. Following [Pressel2015](@citet), we model the saturation vapor pressure +over a mixed-phase surface using a liquid fraction ``λ`` that varies smoothly between +0 (pure ice) and 1 (pure liquid). The effective latent heat and specific heat difference +for the mixture are computed as weighted averages: + +```math +ℒ^{li}_0 = λ \, ℒ^l_0 + (1 - λ) \, ℒ^i_0 , +``` + +```math +\Delta c^{li} = λ \, \Delta c^l + (1 - λ) \, \Delta c^i . +``` + +These effective properties are then used in the Clausius-Clapeyron formula to compute +the saturation vapor pressure over the mixed-phase surface. This approach ensures +thermodynamic consistency and smooth transitions between pure liquid and pure ice states. + +We can illustrate this by computing the mixed-phase specific heat difference for a +50/50 mixture: + +```@example thermo +Δcⁱ = thermo.vapor.heat_capacity - thermo.ice.heat_capacity +λ = 0.5 +Δcˡⁱ = λ * Δcˡ + (1 - λ) * Δcⁱ +``` + +### Visualizing saturation vapor pressure + +The saturation vapor pressure over liquid, ice, and mixed-phase surfaces can be computed +and visualized: ```@example using Breeze -using Breeze.Thermodynamics: saturation_vapor_pressure +using Breeze.Thermodynamics: saturation_vapor_pressure, PlanarMixedPhaseSurface thermo = ThermodynamicConstants() @@ -446,48 +528,99 @@ pᵛˡ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.liquid) for Tⁱ in pᵛⁱ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.ice) for Tⁱ in T] pᵛⁱ⁺[T .> thermo.triple_point_temperature] .= NaN +# Mixed-phase surface with 50% liquid, 50% ice +mixed_surface = PlanarMixedPhaseSurface(0.5) +pᵛᵐ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, mixed_surface) for Tⁱ in T] + using CairoMakie fig = Figure() -ax = Axis(fig[1, 1], xlabel="Temperature (ᵒK)", ylabel="Saturation vapor pressure pᵛ⁺ (Pa)", yscale = log10, xticks=200:20:320) -lines!(ax, T, pᵛˡ⁺, label="vapor pressure over liquid") -lines!(ax, T, pᵛⁱ⁺, linestyle=:dash, label="vapor pressure over ice") +ax = Axis(fig[1, 1], xlabel="Temperature (ᵒK)", ylabel="Saturation vapor pressure pᵛ⁺ (Pa)", + yscale = log10, xticks=200:20:320) +lines!(ax, T, pᵛˡ⁺, label="liquid", linewidth=2) +lines!(ax, T, pᵛⁱ⁺, label="ice", linestyle=:dash, linewidth=2) +lines!(ax, T, pᵛᵐ⁺, label="mixed (λ=0.5)", linestyle=:dot, linewidth=2, color=:purple) axislegend(ax, position=:rb) fig ``` -The saturation specific humidity is +The mixed-phase saturation vapor pressure lies between the liquid and ice curves, +providing a smooth interpolation between the two pure phases. + +## Saturation specific humidity + +The saturation specific humidity ``qᵛ⁺`` is the maximum amount of water vapor that +can exist in equilibrium with a condensed phase at a given temperature and density. +It is related to the saturation vapor pressure by: ```math -qᵛ⁺ ≡ \frac{ρᵛ⁺}{ρ} = \frac{pᵛ⁺}{ρ Rᵐ T} . +qᵛ⁺ ≡ \frac{ρᵛ⁺}{ρ} = \frac{pᵛ⁺}{ρ Rᵛ T} , ``` -and this is what it looks like: +where ``ρᵛ⁺`` is the vapor density at saturation, ``ρ`` is the total air density, +and ``Rᵛ`` is the specific gas constant for water vapor. + +### Visualizing saturation vapor pressure and specific humidity + +We can visualize how both saturation vapor pressure and saturation specific humidity +vary with temperature for different liquid fractions, demonstrating the smooth +interpolation provided by the mixed-phase model: ```@example using Breeze -using Breeze.Thermodynamics: saturation_specific_humidity +using Breeze.Thermodynamics: saturation_vapor_pressure, saturation_specific_humidity, PlanarMixedPhaseSurface thermo = ThermodynamicConstants() -p₀ = 101325 +# Temperature range covering typical atmospheric conditions +T = collect(250:0.1:320) +p₀ = 101325 # Surface pressure (Pa) Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo) -T = collect(273.2:0.1:313.2) -qᵛ⁺ = zeros(length(T)) -for i = 1:length(T) - ρ = p₀ / (Rᵈ * T[i]) - qᵛ⁺[i] = saturation_specific_humidity(T[i], ρ, thermo, thermo.liquid) -end +# Liquid fractions to visualize +λ_values = [0.0, 0.25, 0.5, 0.75, 1.0] +labels = ["ice (λ=0)", "λ=0.25", "λ=0.5", "λ=0.75", "liquid (λ=1)"] +colors = [:blue, :cyan, :purple, :orange, :red] +linestyles = [:solid, :dash, :dot, :dashdot, :solid] using CairoMakie -fig = Figure() -ax = Axis(fig[1, 1], xlabel="Temperature (ᵒK)", ylabel="Saturation specific humidity qᵛ⁺ (kg kg⁻¹)") -lines!(ax, T, qᵛ⁺) +fig = Figure(size=(1000, 400)) + +# Panel 1: Saturation vapor pressure +ax1 = Axis(fig[1, 1], xlabel="Temperature (K)", ylabel="Saturation vapor pressure (Pa)", + yscale=log10, title="Saturation vapor pressure") + +for (i, λ) in enumerate(λ_values) + surface = PlanarMixedPhaseSurface(λ) + pᵛ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, surface) for Tⁱ in T] + lines!(ax1, T, pᵛ⁺, label=labels[i], color=colors[i], linestyle=linestyles[i], linewidth=2) +end + +axislegend(ax1, position=:lt) + +# Panel 2: Saturation specific humidity +ax2 = Axis(fig[1, 2], xlabel="Temperature (K)", ylabel="Saturation specific humidity (kg/kg)", + title="Saturation specific humidity") + +for (i, λ) in enumerate(λ_values) + surface = PlanarMixedPhaseSurface(λ) + qᵛ⁺ = zeros(length(T)) + for (j, Tⁱ) in enumerate(T) + ρ = p₀ / (Rᵈ * Tⁱ) # Approximate density using dry air + qᵛ⁺[j] = saturation_specific_humidity(Tⁱ, ρ, thermo, surface) + end + lines!(ax2, T, qᵛ⁺, label=labels[i], color=colors[i], linestyle=linestyles[i], linewidth=2) +end + fig ``` +This figure shows how the liquid fraction ``λ`` smoothly interpolates between pure ice +(``λ = 0``) and pure liquid (``λ = 1``). At lower temperatures, the differences between +phases are more pronounced. The mixed-phase model allows for realistic representation of +conditions near the freezing point where both liquid and ice may coexist. + ## Moist static energy For moist air, a convenient thermodynamic invariant that couples temperature, composition, and height is the moist static energy (MSE), From f822fb4b4e7411744f090cb47c619cb30583e58c Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Tue, 11 Nov 2025 16:25:27 -0700 Subject: [PATCH 6/7] Update src/Thermodynamics/vapor_saturation.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Mosè Giordano <765740+giordano@users.noreply.github.com> --- src/Thermodynamics/vapor_saturation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index 7c535131..50912863 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -131,7 +131,7 @@ qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, Pl Note, this is slightly smaller than the saturation specific humidity over an ice surface: ```jldoctest saturation -qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, PlanarIceSurface()) +julia> qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, thermo, PlanarIceSurface()) 0.011945100768555072 ``` From f78b4331d1db0bc52bbeca38a7a9300cd90df60d Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Tue, 11 Nov 2025 16:26:28 -0700 Subject: [PATCH 7/7] Update src/Thermodynamics/vapor_saturation.jl --- src/Thermodynamics/vapor_saturation.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index 50912863..ec60f018 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -76,10 +76,11 @@ struct PlanarMixedPhaseSurface{FT} liquid_fraction :: FT end +@inline specific_heat_difference(thermo, ::PlanarLiquidSurface) = specific_heat_difference(thermo, thermo.liquid) @inline specific_heat_difference(thermo, ::PlanarIceSurface) = specific_heat_difference(thermo, thermo.ice) -@inline absolute_zero_latent_heat(thermo, ::PlanarIceSurface) = absolute_zero_latent_heat(thermo, thermo.ice) @inline absolute_zero_latent_heat(thermo, ::PlanarLiquidSurface) = absolute_zero_latent_heat(thermo, thermo.liquid) -@inline specific_heat_difference(thermo, ::PlanarLiquidSurface) = specific_heat_difference(thermo, thermo.liquid) +@inline absolute_zero_latent_heat(thermo, ::PlanarIceSurface) = absolute_zero_latent_heat(thermo, thermo.ice) + @inline function specific_heat_difference(thermo, surf::PlanarMixedPhaseSurface) Δcˡ = specific_heat_difference(thermo, thermo.liquid)