diff --git a/NEWS.md b/NEWS.md index 265571b0..aa43ff15 100644 --- a/NEWS.md +++ b/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 ------- diff --git a/Project.toml b/Project.toml index 6f37b6ac..ce653ca9 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/Formulation.md b/docs/src/Formulation.md index a2cbe96f..241410b3 100644 --- a/docs/src/Formulation.md +++ b/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,7 +126,7 @@ 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} @@ -134,7 +134,7 @@ Here, the reference specific internal energy ``I_{v,0}`` is the difference in sp 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 diff --git a/src/relations.jl b/src/relations.jl index 721bf11b..e69ebef6 100644 --- a/src/relations.jl +++ b/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ⁿ))^α diff --git a/test/relations.jl b/test/relations.jl index 9bda6181..53e76650 100644 --- a/test/relations.jl +++ b/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))