Merge pull request #224 from CliMA/aj/internal_energy_definition

Update dry air internal energy and enthalpy definitions (eq 8a and 12a)

NEWS.md

@@ -2,8 +2,12 @@ Thermodynamics.jl Release Notes
 ========================
 
 main
+----
+
+v0.12.8
 -------
 - ![][badge-🤖precisionΔ] Change the tolerance of PhaseEquil constructor to 1e-4
+- ![][badge-🔥behavioralΔ] Change the definition of dry air internal energy and enthalpy
 
 v0.12.7
 -------

Project.toml

@@ -1,7 +1,7 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.12.7"
+version = "0.12.8"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

docs/src/Formulation.md

@@ -113,7 +113,7 @@ The specific internal energies of the constituents of moist air can be written a
 ```math
 \begin{equation}
 \begin{aligned}
-I_d(T) & = c_{vd} (T - T_0),  \\
+I_d(T) & = c_{vd} (T - T_0) - R_d T_0,  \\
 I_v(T) & = c_{vv} (T - T_0) + I_{v,0},\\
 I_l(T) & = c_{vl} (T - T_0), \\
 I_i(T) & = c_{vi} (T - T_0) - I_{i,0}.
@@ -126,15 +126,15 @@ Here, the reference specific internal energy ``I_{v,0}`` is the difference in sp
 \begin{equation}
 \begin{aligned}
      I(T, q) & = (1-q_t) I_d(T) + q_v I_v(T) + q_l I_l(T) + q_i I_i(T)\\
-          & = c_{vm}(q) (T - T_0)  + q_v I_{v,0} - q_i I_{i,0}.
+          & = c_{vm}(q) (T - T_0)  + q_v I_{v,0} - q_i I_{i,0} - (1 - q_t) R_d T_0.
 \end{aligned}
 \label{e:totalInternalEnergy}
 \end{equation}
 ```
 The internal energy can be inverted to obtain the temperature given ``I`` and the specific humidities,
 ```math
 \begin{equation}
-    T = T_0 + \frac{I - (q_t - q_l) I_{v,0} + q_i (I_{i,0} + I_{v,0})}{c_{vm}(q)},
+    T = T_0 + \frac{I - (q_t - q_l - q_i) I_{v,0} + q_i I_{i,0} + (1 - q_t) R_d T_0}{c_{vm}(q)},
     \label{e:temperature}
 \end{equation}
 ```
@@ -158,7 +158,7 @@ The specific enthalpies of the constituents of moist air are obtained by adding
 \begin{equation}
 \label{e:Enthalpies}
 \begin{aligned}
-    h_d(T) = I_d(T) + R_d T &= c_{pd}(T-T_0) + R_d T_0, \\
+    h_d(T) = I_d(T) + R_d T &= c_{pd}(T-T_0), \\
     h_v(T) = I_v(T) + R_v T &= c_{pv}(T-T_0) + L_{v,0}, \\
     h_l(T) = I_l(T) &= c_{pl}(T-T_0), \\
     h_i(T) = I_i(T) &= c_{pi}(T-T_0) - L_{f,0}.
@@ -170,7 +170,7 @@ The enthalpy of moist air is the weighted sum of the constituent enthalpies:
 \begin{equation}
 \begin{split}\label{e:enthalpy_definition}
     h(T, q)  &= (1-q_t) h_d + q_v h_v + q_l h_l + q_i h_i \\
-        &= c_{pm}(q) (T-T_0) + q_v L_{v,0} - q_i L_{f,0} + (1-q_t) R_d T_0\\
+        &= c_{pm}(q) (T-T_0) + q_v L_{v,0} - q_i L_{f,0}\\
         &= I(q, T) + R_m T,
 \end{split}
 \end{equation}
@@ -241,7 +241,7 @@ where ``λ_i`` interpolates between 0 at the temperature of homogeneous ice nucl
 
 ## Saturation Specific Humidity
 
-From the saturation vapor pressure ``p_v^*``, the saturation specific humidity can be computed using the ideal gas law \eqref{e:eos}, giving the density of water vapor at saturation ``ρ_v^* = p_v^*(T)/(R_v T)``, and hence the saturation specific humidity 
+From the saturation vapor pressure ``p_v^*``, the saturation specific humidity can be computed using the ideal gas law \eqref{e:eos}, giving the density of water vapor at saturation ``ρ_v^* = p_v^*(T)/(R_v T)``, and hence the saturation specific humidity
 ```math
 \begin{equation}
      q_v^* = \frac{ρ_v^*}{ρ} = \frac{p_v^*(T)}{ρ R_v T}.
@@ -281,7 +281,7 @@ A zeroth-order approximation of the temperature ``T`` satisfying the saturation
     T_1 = T_0 + \frac{I - q_t I_{v,0}}{c_{vm}^*}.
 \end{equation}
 ```
-Here, the isochoric specific heat capacity in equilibrium, ``c_{vm}^* = c_{vm}(q^*)``, is the specific heat capacity under equilibrium partitioning ``q^*`` of the phases, which here, for unsaturated conditions, means ``q^*=(q_t; q_l=0, q_i=0)``. If the total specific humidity ``q_t`` is less than the saturation specific humidity at ``T_1`` (``q_t \le q_v^*(T_1, ρ)``), the air is indeed unsaturated, and ``T=T_1`` is the exact temperature consistent given ``I``, ``ρ``, and ``q_t``. 
+Here, the isochoric specific heat capacity in equilibrium, ``c_{vm}^* = c_{vm}(q^*)``, is the specific heat capacity under equilibrium partitioning ``q^*`` of the phases, which here, for unsaturated conditions, means ``q^*=(q_t; q_l=0, q_i=0)``. If the total specific humidity ``q_t`` is less than the saturation specific humidity at ``T_1`` (``q_t \le q_v^*(T_1, ρ)``), the air is indeed unsaturated, and ``T=T_1`` is the exact temperature consistent given ``I``, ``ρ``, and ``q_t``.
 
 If the air is saturated (``q_t > q_v^*(T_1, ρ)``), successively improved temperature estimates ``T_{n+1}`` can be obtained from the temperature ``T_n`` (``n=1,\dots``) by Newton's method, with analytical gradients. Linearizing the saturation internal energy ``I^*(T; ρ, q_t)`` around the temperature ``T_n`` gives
 ```math
@@ -298,7 +298,7 @@ and solving for the temperature ``T`` gives the first-order Newton update
 The derivative ``\partial I^*/\partial T|_{T_n}`` is obtained by differentiation of the internal energy \eqref{e:total_internal_energy}, \hl{[add derivatives of phase partitioning function and make that function smooth]}
 ```math
 \begin{multline}
-     \left.\frac{\partial I^*(T; ρ, q_t)}{\partial T}\right|_{T_n} 
+     \left.\frac{\partial I^*(T; ρ, q_t)}{\partial T}\right|_{T_n}
      = c_{vm}^*(q_t, T_n) \\
      +  \left( I_{v,0} + [1-λ_p(T_n)]I_{i,0} + (T_n - T_0) \left. \frac{dc_{vm}^*}{dq_v^*}\right|_{T_n} \right) \left. \frac{\partial q_v^*(T; ρ, q_t)}{\partial T}\right|_{T_n},
 \end{multline}
@@ -315,11 +315,11 @@ obtained from the definition \eqref{e:SpecificHeat} of the specific heat of mois
     \left. \frac{\partial q_v^*(T; ρ, q_t)}{\partial T}\right|_{T_n} = q_v^*(T_n) \frac{L}{R_v T_n^2} \quad \text{with} \quad L = λ_p(T_n) L_v + [1-λ_p(T_n)] L_s,
 \end{equation}
 ```
-obtained from the Clausius-Clapeyron relation \eqref{e:Clausius_Clapeyron} and the relation \eqref{e:sat_shum} between specific humidity and vapor pressure. 
+obtained from the Clausius-Clapeyron relation \eqref{e:Clausius_Clapeyron} and the relation \eqref{e:sat_shum} between specific humidity and vapor pressure.
 
-The resulting successive Newton approximations ``T_n`` generally converge quadratically. Because condensate specific humidities are usually small, ``T_1`` provides a close initial estimate, and few iterations are needed. Even the first-order approximation ``T\approx T_2`` often suffices. However, convergence may not be achieved near the phase transition at the freezing temperature ``T_{\mathrm{freeze}}`` because the derivative of ``I^*`` with respect to temperature is discontinuous there. In that case, the number of iterations needs to be limited (2--3 iterations generally suffice). 
+The resulting successive Newton approximations ``T_n`` generally converge quadratically. Because condensate specific humidities are usually small, ``T_1`` provides a close initial estimate, and few iterations are needed. Even the first-order approximation ``T\approx T_2`` often suffices. However, convergence may not be achieved near the phase transition at the freezing temperature ``T_{\mathrm{freeze}}`` because the derivative of ``I^*`` with respect to temperature is discontinuous there. In that case, the number of iterations needs to be limited (2--3 iterations generally suffice).
 
-Using saturation adjustment makes it possible to construct a moist dynamical core that has the total specific humidity ``q_t`` as the only prognostic moisture variable. The price for this simplicity is the necessity to solve a nonlinear problem iteratively (or approximately) at each time step, and being confined to an equilibrium thermodynamics framework which cannot adequately account for non-equilibrium processes. Using explicit tracers for the condensates ``q_l`` and ``q_i`` in addition to ``q_t`` avoids iterations at each time step and allows the inclusion of explicit non-equilibrium processes, such as those leading to the formation of supercooled liquid in mixed-phase clouds. 
+Using saturation adjustment makes it possible to construct a moist dynamical core that has the total specific humidity ``q_t`` as the only prognostic moisture variable. The price for this simplicity is the necessity to solve a nonlinear problem iteratively (or approximately) at each time step, and being confined to an equilibrium thermodynamics framework which cannot adequately account for non-equilibrium processes. Using explicit tracers for the condensates ``q_l`` and ``q_i`` in addition to ``q_t`` avoids iterations at each time step and allows the inclusion of explicit non-equilibrium processes, such as those leading to the formation of supercooled liquid in mixed-phase clouds.
 
 ## Auxiliary Thermodynamic Functions
 

src/relations.jl

@@ -382,11 +382,14 @@ and, optionally,
     cvm = cv_m(param_set, q),
 ) where {FT <: Real}
     T_0 = TP.T_0(param_set)
+    R_d = TP.R_d(param_set)
     e_int_v0 = TP.e_int_v0(param_set)
     e_int_i0 = TP.e_int_i0(param_set)
     return T_0 +
            (
-        e_int - (q.tot - q.liq) * e_int_v0 + q.ice * (e_int_v0 + e_int_i0)
+        e_int - (q.tot - q.liq - q.ice) * e_int_v0 +
+        q.ice * e_int_i0 +
+        (1 - q.tot) * R_d * T_0
     ) / cvm
 end
 
@@ -407,15 +410,9 @@ and, optionally,
 ) where {FT <: Real}
     cp_m_ = cp_m(param_set, q)
     T_0 = TP.T_0(param_set)
-    R_d = TP.R_d(param_set)
     LH_v0 = TP.LH_v0(param_set)
     LH_f0 = TP.LH_f0(param_set)
-    e_int_i0 = TP.e_int_i0(param_set)
-    q_vap = vapor_specific_humidity(q)
-    return (
-        h + cp_m_ * T_0 - q_vap * LH_v0 + q.ice * LH_f0 -
-        (1 - q.tot) * R_d * T_0
-    ) / cp_m_
+    return T_0 + (h - (q.tot - q.liq - q.ice) * LH_v0 + q.ice * LH_f0) / cp_m_
 end
 
 """
@@ -469,10 +466,11 @@ and, optionally,
     cvm = cv_m(param_set, q),
 ) where {FT <: Real}
     T_0 = TP.T_0(param_set)
+    R_d = TP.R_d(param_set)
     e_int_v0 = TP.e_int_v0(param_set)
     e_int_i0 = TP.e_int_i0(param_set)
-    return cvm * (T - T_0) + (q.tot - q.liq) * e_int_v0 -
-           q.ice * (e_int_v0 + e_int_i0)
+    return cvm * (T - T_0) + (q.tot - q.liq - q.ice) * e_int_v0 -
+           q.ice * e_int_i0 - (1 - q.tot) * R_d * T_0
 end
 
 """
@@ -516,8 +514,8 @@ The dry air internal energy
 @inline function internal_energy_dry(param_set::APS, T::FT) where {FT <: Real}
     T_0 = TP.T_0(param_set)
     cv_d = TP.cv_d(param_set)
-
-    return cv_d * (T - T_0)
+    R_d = TP.R_d(param_set)
+    return cv_d * (T - T_0) - R_d * T_0
 end
 
 """
@@ -1331,7 +1329,7 @@ is a function that is 1 above `T_freeze` and goes to zero below `T_icenuc`.
     # Model for cloud water condensate probabilty when temperature is below
     # freezing (all liquid condensate), but above the temperature of homogeneous
     # ice nucleation (all ice condensate). In it's simplest form, this is simply
-    # a number that varies between 0 and 1. For example, see figure 6 in 
+    # a number that varies between 0 and 1. For example, see figure 6 in
     # Hu et al., https://doi.org/10.1029/2009JD012384, JGR (2010).
     λᵖ = ((T - Tⁿ) / (Tᶠ - Tⁿ))^α
 

test/relations.jl

@@ -221,7 +221,7 @@ end
     # test the wrapper for q_vap_saturation over liquid water and ice
     ρ = FT(1)
     ρu = FT[1, 2, 3]
-    ρe = FT(1100)
+    ρe = FT(1100) - _R_d * _T_0
     e_pot = FT(93)
     e_int = internal_energy(ρ, ρe, ρu, e_pot)
     q_pt = PhasePartition(FT(0.02), FT(0.002), FT(0.002))
@@ -291,31 +291,31 @@ end
     T = FT(300)
     e_kin = FT(11)
     e_pot = FT(13)
-    @test air_temperature(param_set, _cv_d * (T - _T_0)) === FT(T)
+    @test air_temperature(param_set, _cv_d * (T - _T_0) - _R_d * _T_0) === FT(T)
     @test air_temperature(
         param_set,
-        _cv_d * (T - _T_0),
+        _cv_d * (T - _T_0) - _R_d * _T_0,
         PhasePartition(FT(0)),
     ) === FT(T)
 
     @test air_temperature(
         param_set,
-        cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0),
+        cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0) - _R_d * _T_0,
         PhasePartition(FT(0)),
     ) === FT(T)
     @test air_temperature(
         param_set,
-        cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) +
-        q_tot * _e_int_v0,
+        cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) -
+        (1 - q_tot) * _R_d * _T_0 + q_tot * _e_int_v0,
         PhasePartition(q_tot),
     ) ≈ FT(T)
 
     @test total_energy(param_set, FT(e_kin), FT(e_pot), _T_0) ===
-          FT(e_kin) + FT(e_pot)
+          FT(e_kin) + FT(e_pot) - _R_d * _T_0
     @test total_energy(param_set, FT(e_kin), FT(e_pot), FT(T)) ≈
-          FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0)
+          FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0) - _R_d * _T_0
     @test total_energy(param_set, FT(0), FT(0), _T_0, PhasePartition(q_tot)) ≈
-          q_tot * _e_int_v0
+          q_tot * _e_int_v0 - (1 - q_tot) * _R_d * _T_0
 
     # phase partitioning in equilibrium
     q_liq = FT(0.1)
@@ -595,7 +595,13 @@ end
             )
         _e_int = (e_int_upper .+ e_int_lower) / 2
         ts = PhaseEquil_ρeq.(param_set, ρ, _e_int, q_tot)
-        @test all(air_temperature.(param_set, ts) .== Ref(_T_freeze))
+        @test all(
+            isapprox.(
+                air_temperature.(param_set, ts),
+                Ref(_T_freeze),
+                rtol = rtol_temperature,
+            ),
+        )
 
         # Args needs to be in sync with PhaseEquil:
         ts =
@@ -608,7 +614,13 @@ end
                 FT(1e-1),
                 RS.SecantMethod,
             )
-        @test all(air_temperature.(param_set, ts) .== Ref(_T_freeze))
+        @test all(
+            isapprox.(
+                air_temperature.(param_set, ts),
+                Ref(_T_freeze),
+                rtol = rtol_temperature,
+            ),
+        )
 
         # PhaseEquil
         ts_exact = PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot, 100, FT(1e-6))