Build abstraction for computing saturation vapor pressure over mixed phase surfaces

src/Thermodynamics/vapor_saturation.jl

@@ -1,46 +1,99 @@
 """
 $(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
-    Tᵣ = thermo.energy_reference_temperature
+@inline function saturation_vapor_pressure(T, thermo, surface)
+    ℒ₀ = absolute_zero_latent_heat(thermo, surface)
+    Δcᵝ = specific_heat_difference(thermo, surface)
 
     Tᵗʳ = thermo.triple_point_temperature
     pᵗʳ = thermo.triple_point_pressure
     Rᵛ = vapor_gas_constant(thermo)
 
+    return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * ℒ₀ / Rᵛ)
+end
+
+@inline function specific_heat_difference(thermo, phase::CondensedPhase)
     cᵖᵛ = thermo.vapor.heat_capacity
     cᵝ = phase.heat_capacity
-    Δcᵝ = cᵖᵛ - cᵝ
+    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)
+    Tᵣ = thermo.energy_reference_temperature
+    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
 
-    return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * (ℒᵣ - Δcᵝ * Tᵣ) / Rᵛ)
+@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
 
 """

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