add Clausius_Clapeyron_relation

docs/Project.toml

@@ -3,6 +3,7 @@ CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"

docs/make.jl

@@ -14,6 +14,7 @@ pages = Any[
     "Tested profiles" => "TestedProfiles.md",
     "Temperature profiles" => "TemperatureProfiles.md",
     "Developer docs" => "DevDocs.md",
+    "Clausius Clapeyron relation" => "Clausius_Clapeyron.md",
     "Thermodynamics overview" => "Formulation.md",
     "References" => "References.md",
 ]

docs/src/Clausius_Clapeyron.jl

@@ -0,0 +1,115 @@
+import ForwardDiff
+import Thermodynamics as TD
+import CLIMAParameters as CP
+using Thermodynamics.TestedProfiles
+import Thermodynamics.Parameters as TP
+
+function get_parameter_set(::Type{FT}) where {FT}
+    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
+    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
+    CP.log_parameter_information(toml_dict, logfilepath)
+    return param_set
+end
+const param_set_Float64 = get_parameter_set(Float64)
+const param_set_Float32 = get_parameter_set(Float32)
+parameter_set(::Type{Float64}) = param_set_Float64
+parameter_set(::Type{Float32}) = param_set_Float32
+ArrayType = Array{Float64}
+
+FT = eltype(ArrayType)
+param_set = parameter_set(FT)
+profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
+(; T, p, e_int, ρ, θ_liq_ice, phase_type) = profiles
+(; q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot) = profiles
+
+k = findfirst(q -> q > 0.01, q_tot) # test for one value with q_tot above some threshhold
+ts_sol = TD.PhaseEquil_ρTq(param_set, ρ[k], T[k], q_tot[k])
+
+function q_vap_sat(_T::FT) where {FT}
+    _ρ = TD.air_density(param_set, ts_sol)
+    _q_tot = TD.total_specific_humidity(param_set, ts_sol)
+    _phase_type = TD.PhaseEquil{FT}
+    _q_pt = TD.PhasePartition_equil(
+        param_set,
+        _T,
+        oftype(_T, _ρ),
+        oftype(_T, _q_tot),
+        _phase_type,
+    )
+    return TD.q_vap_saturation(
+        param_set,
+        _T,
+        oftype(_T, _ρ),
+        _phase_type,
+        _q_pt,
+    )
+end
+
+function ∂q_vap_sat_∂T_vs_T(_T::FT) where {FT}
+    _ρ = TD.air_density(param_set, ts_sol)
+    _q_tot = TD.total_specific_humidity(param_set, ts_sol)
+    _λ = TD.liquid_fraction(param_set, ts_sol)
+    _phase_type = TD.PhaseEquil{FT}
+    _q_pt = TD.PhasePartition_equil(
+        param_set,
+        _T,
+        oftype(_T, _ρ),
+        oftype(_T, _q_tot),
+        _phase_type,
+    )
+    _q_vap_sat = TD.q_vap_saturation(param_set, _T, _ρ, _phase_type, _q_pt)
+    return TD.∂q_vap_sat_∂T(
+        param_set,
+        oftype(_T, _λ),
+        _T,
+        oftype(_T, _q_vap_sat),
+    )
+end
+
+∂q_vap_sat_∂T_fd = _T -> ForwardDiff.derivative(q_vap_sat, _T)
+
+_ρ = TD.air_density(param_set, ts_sol)
+_q_tot = TD.total_specific_humidity(param_set, ts_sol)
+
+T_sorted = sort(T)
+import Plots
+p1 = Plots.plot()
+Plots.plot!(
+    T_sorted,
+    ∂q_vap_sat_∂T_fd.(T_sorted);
+    label = "∂qvsat_∂T ForwardDiff",
+    yaxis = :log,
+)
+Plots.plot!(
+    T_sorted,
+    ∂q_vap_sat_∂T_vs_T.(T_sorted);
+    label = "∂qvsat_∂T Analytic",
+    yaxis = :log,
+)
+Plots.plot!(; xlabel = "T [K]", legend = :topleft)
+
+p2 = Plots.plot()
+Plots.plot!(
+    T_sorted,
+    ∂q_vap_sat_∂T_fd.(T_sorted) .- ∂q_vap_sat_∂T_vs_T.(T_sorted);
+    label = "error",
+)
+Plots.plot!(; xlabel = "T [K]", legend = :topleft)
+
+p3 = Plots.plot()
+Plots.plot!(T_sorted, q_vap_sat.(T_sorted); label = "q_vap_sat")
+Plots.plot!(; xlabel = "T [K]")
+ρq_vals = "ρ=$_ρ, q_tot=$_q_tot"
+Plots.plot(
+    p1,
+    p2,
+    p3;
+    layout = Plots.grid(3, 1),
+    plot_title = "$ρq_vals",
+    titlefontsizes = 5,
+)
+
+Plots.savefig("Clausius_Clapeyron.svg")

docs/src/Clausius_Clapeyron.md

@@ -0,0 +1,15 @@
+# Clausius Clapeyron relation
+
+This script plots the Clausius Clapeyron relation for a range of temperatures. The analytically derived expression is compared with a solution computed using ForwardDiff.jl.
+
+
+!!! warn
+
+    This script is decoupled from the implementation in the test suite,
+    and should be unified to ensure that tests and plots stay.
+
+```@example
+include("Clausius_Clapeyron.jl")
+```
+
+![](Clausius_Clapeyron.svg)

src/relations.jl

@@ -222,6 +222,31 @@ vapor_specific_humidity(q::PhasePartition) = max(0, q.tot - q.liq - q.ice)
 vapor_specific_humidity(param_set::APS, ts::ThermodynamicState) =
     vapor_specific_humidity(PhasePartition(param_set, ts))
 
+"""
+    ∂q_vap_sat_∂T(param_set, ts)
+
+The change in saturation vapor specific humidity with temperature given by the Clausius Clapeyron relation and the definition of specific humidity `q = p_v/(ρ R_v T)` in terms of vapor pressure `p_v`.
+ - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
+ - `ts` ThermodynamicState
+"""
+function ∂q_vap_sat_∂T(
+    param_set::APS,
+    λ::FT,
+    T::FT,
+    q_vap_sat::FT,
+) where {FT <: Real}
+    R_v::FT = TP.R_v(param_set)
+    L = weighted_latent_heat(param_set, T, λ)
+    return q_vap_sat * (L / (R_v * T^2) - 1 / T)
+end
+
+function ∂q_vap_sat_∂T(param_set::APS, ts::ThermodynamicState)
+    λ = liquid_fraction(param_set, ts)
+    T = air_temperature(param_set, ts)
+    q_vap_sat = vapor_specific_humidity(param_set, ts)
+    return ∂q_vap_sat_∂T(param_set, λ, T, q_vap_sat)
+end
+
 """
     cp_m(param_set, q::PhasePartition)
 
@@ -781,6 +806,19 @@ function latent_heat_generic(
     return LH_0 + Δcp * (T - T_0)
 end
 
+"""
+    weighted_latent_heat(param_set, T, λ)
+
+Weighted latent heats, computed from
+ - `param_set` - a `ThermodynamicsParameters` struct
+ - `T` air temperature
+ - `λ` liquid fraction
+"""
+function weighted_latent_heat(param_set::APS, T::FT, λ::FT) where {FT <: Real}
+    L_v = latent_heat_vapor(param_set, T)
+    L_s = latent_heat_sublim(param_set, T)
+    return λ * L_v + (1 - λ) * L_s
+end
 
 """
     Phase
@@ -911,13 +949,12 @@ function saturation_vapor_pressure(
     cp_l::FT = TP.cp_l(param_set)
     cp_i::FT = TP.cp_i(param_set)
     # get phase partitioning
-    liquid_frac = liquid_fraction(param_set, T, phase_type, q)
-    ice_frac = 1 - liquid_frac
+    λ = liquid_fraction(param_set, T, phase_type, q)
 
     # effective latent heat at T_0 and effective difference in isobaric specific
     # heats of the mixture
-    LH_0 = liquid_frac * LH_v0 + ice_frac * LH_s0
-    Δcp = liquid_frac * (cp_v - cp_l) + ice_frac * (cp_v - cp_i)
+    LH_0 = λ * LH_v0 + (1 - λ) * LH_s0
+    Δcp = λ * (cp_v - cp_l) + (1 - λ) * (cp_v - cp_i)
 
     # saturation vapor pressure over possible mixture of liquid and ice
     return saturation_vapor_pressure(param_set, T, LH_0, Δcp)
@@ -1277,7 +1314,7 @@ liquid_fraction(param_set::APS, ts::ThermodynamicState) = liquid_fraction(
 
 """
     PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
-    PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, liquid_frac)
+    PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, λ)
 
 Partition the phases in equilibrium, returning a [`PhasePartition`](@ref) object using the
 [`liquid_fraction`](@ref) function where
@@ -1288,7 +1325,7 @@ Partition the phases in equilibrium, returning a [`PhasePartition`](@ref) object
  - `q_tot` total specific humidity
  - `phase_type` a thermodynamic state type
  - `p_vap_sat` saturation vapor pressure
- - `liquid_frac` liquid fraction
+ - `λ` liquid fraction
 
 The residual `q.tot - q.liq - q.ice` is the vapor specific humidity.
 """
@@ -1298,11 +1335,11 @@ function PhasePartition_equil(
     ρ::FT,
     q_tot::FT,
     p_vap_sat::FT,
-    liquid_frac::FT,
+    λ::FT,
 ) where {FT <: Real}
     q_c = saturation_excess(param_set, T, ρ, p_vap_sat, PhasePartition(q_tot)) # condensate specific humidity
-    q_liq = liquid_frac * q_c                                                  # liquid specific humidity
-    q_ice = (1 - liquid_frac) * q_c                                            # ice specific humidity
+    q_liq = λ * q_c                                                            # liquid specific humidity
+    q_ice = (1 - λ) * q_c                                                      # ice specific humidity
     return PhasePartition(q_tot, q_liq, q_ice)
 end
 
@@ -1314,8 +1351,8 @@ function PhasePartition_equil(
     ::Type{phase_type},
 ) where {FT <: Real, phase_type <: ThermodynamicState}
     p_vap_sat = saturation_vapor_pressure(param_set, phase_type, T)
-    liquid_frac = liquid_fraction(param_set, T, phase_type) # fraction of condensate that is liquid
-    return PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, liquid_frac)
+    λ = liquid_fraction(param_set, T, phase_type) # fraction of condensate that is liquid
+    return PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, λ)
 end
 
 PhasePartition_equil(param_set::APS, ts::AbstractPhaseNonEquil) =
@@ -1363,9 +1400,9 @@ function PhasePartition(param_set::APS, ts::AbstractPhaseEquil)
     q_tot = total_specific_humidity(param_set, ts)
     phase_type = typeof(ts)
     p_vap_sat = saturation_vapor_pressure(param_set, phase_type, T)
-    liquid_frac = liquid_fraction(param_set, T, phase_type) # fraction of condensate that is liquid
+    λ = liquid_fraction(param_set, T, phase_type) # fraction of condensate that is liquid
 
-    return PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, liquid_frac)
+    return PhasePartition_equil(param_set, T, ρ, q_tot, p_vap_sat, λ)
 end
 PhasePartition(param_set::APS, ts::AbstractPhaseNonEquil) = ts.q
 
@@ -1377,9 +1414,6 @@ function ∂e_int_∂T(
     q_tot::FT,
     ::Type{phase_type},
 ) where {FT <: Real, phase_type <: PhaseEquil}
-    LH_v0::FT = TP.LH_v0(param_set)
-    LH_s0::FT = TP.LH_s0(param_set)
-    R_v::FT = TP.R_v(param_set)
     T_0::FT = TP.T_0(param_set)
     cv_v::FT = TP.cv_v(param_set)
     cv_l::FT = TP.cv_l(param_set)
@@ -1396,14 +1430,14 @@ function ∂e_int_∂T(
     q_c = condensate(q)
     cvm = cv_m(param_set, q)
     q_vap_sat = q_vap_saturation_from_density(param_set, T, ρ, p_vap_sat)
-    L = λ * LH_v0 + (1 - λ) * LH_s0
+    L = weighted_latent_heat(param_set, T, λ)
 
     ∂λ_∂T = (T_i < T < T_f) ? (1 / (T_f - T_i))^n_i * n_i * T^(n_i - 1) : FT(0)
-    ∂q_vap_sat_∂T = q_vap_sat * L / (R_v * T^2)
+    _∂q_vap_sat_∂T = ∂q_vap_sat_∂T(param_set, λ, T, q_vap_sat)
     dcvm_dq_vap = cv_v - λ * cv_l - (1 - λ) * cv_i
     return cvm +
            (e_int_v0 + (1 - λ) * e_int_i0 + (T - T_0) * dcvm_dq_vap) *
-           ∂q_vap_sat_∂T +
+           _∂q_vap_sat_∂T +
            q_c * e_int_i0 * ∂λ_∂T
 end
 
@@ -2843,7 +2877,7 @@ end
 """
     partial_pressure_dry(param_set, p, q)
 
-The specific entropy of water vapor, given
+The partial pressure of water vapor, given
 
  - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
  - `p` air pressure
@@ -2862,7 +2896,7 @@ end
 """
     partial_pressure_vapor(param_set, p, q)
 
-The specific entropy of water vapor, given
+The partial pressure of water vapor, given
 
  - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
  - `p` air pressure

test/relations.jl

@@ -13,6 +13,7 @@ import RootSolvers
 const RS = RootSolvers
 
 using LinearAlgebra
+import ForwardDiff
 
 const TP = TD.Parameters
 
@@ -516,6 +517,61 @@ end
         @test all(saturated.(param_set, ts[RH_sat_mask]))
         @test !any(saturated.(param_set, ts[RH_unsat_mask]))
 
+        # test Clausius Clapeyron relation
+        k = findfirst(q -> q > 0.01, q_tot) # test for one value with q_tot above some threshhold
+        ts_sol = TD.PhaseEquil_ρTq(param_set, ρ[k], T[k], q_tot[k])
+
+        function q_vap_sat(_T::FT) where {FT}
+            _ρ = TD.air_density(param_set, ts_sol)
+            _q_tot = TD.total_specific_humidity(param_set, ts_sol)
+            _phase_type = PhaseEquil{FT}
+            _q_pt = PhasePartition_equil(
+                param_set,
+                _T,
+                oftype(_T, _ρ),
+                oftype(_T, _q_tot),
+                _phase_type,
+            )
+            return TD.q_vap_saturation(
+                param_set,
+                _T,
+                oftype(_T, _ρ),
+                _phase_type,
+                _q_pt,
+            )
+        end
+
+        function ∂q_vap_sat_∂T_vs_T(_T::FT) where {FT}
+            _ρ = TD.air_density(param_set, ts_sol)
+            _q_tot = TD.total_specific_humidity(param_set, ts_sol)
+            _phase_type = PhaseEquil{FT}
+            _λ = TD.liquid_fraction(param_set, ts_sol)
+            _q_pt = TD.PhasePartition_equil(
+                param_set,
+                _T,
+                oftype(_T, _ρ),
+                oftype(_T, _q_tot),
+                _phase_type,
+            )
+            _q_vap_sat =
+                TD.q_vap_saturation(param_set, _T, _ρ, _phase_type, _q_pt)
+            return TD.∂q_vap_sat_∂T(
+                param_set,
+                oftype(_T, _λ),
+                _T,
+                oftype(_T, _q_vap_sat),
+            )
+        end
+
+        ∂q_vap_sat_∂T_fd = _T -> ForwardDiff.derivative(q_vap_sat, _T)
+        @test all(
+            isapprox.(
+                log.(∂q_vap_sat_∂T_fd.(T)),
+                log.(∂q_vap_sat_∂T_vs_T.(T));
+                rtol = 2e-2,
+            ),
+        )
+
         # PhaseEquil (freezing)
         _T_freeze = FT(TP.T_freeze(param_set))
         e_int_upper =