A synthesis of previous studies suggests the need for new, more accurate approximations for ice–liquid water potential temperature (θil), a thermodynamic variable utilized in some numerical models. Starting from equations presented in a previous study, two new approximate formulations of θil are derived, along with their governing equations. The new formulations are significant improvements over previous ones because no terms are dropped during their derivation and no Taylor series approximations are utilized. The governing equations for the new formulations reveal conditions under which θil can be considered a conserved variable.
Potential temperature lapse rates determined from a reference thermodynamic equation are compared numerically against lapse rates determined from several approximations of θil. Many of the findings agree with previous studies. However, the results show that a commonly used formulation does not account for the specific heats of water, and thus has an inherent cold bias. When the tendency to θil is assumed to be zero (by approximation), one of the new formulations is superior to all other formulations that have been presented previously. When the full governing equation is used, the new θil formulations produce results identical to those from the reference equation.
The liquid water potential temperature (θl) was introduced by Betts (1973). It is the potential temperature (θ) an air parcel would have if all liquid water in the parcel were evaporated. There is an analogous thermodynamic variable if one considers the ice phase; that is, ice water potential temperature (θi) is the potential temperature an air parcel would have if all ice in the parcel were sublimated. Tripoli and Cotton (1981) introduced ice–liquid water potential temperature (θil), which, presumably, is the potential temperature an air parcel would have if all liquid were evaporated and all ice were sublimated.
The usefulness of using θl (for applications involving shallow convection) or θil (for applications involving deep convection) as a predictive variable in numerical models has been discussed by Betts (1973), Deardorff (1976), and Tripoli and Cotton (1981). Because θil is (or is assumed to be) conserved during phase changes of water, useful thermodynamic variables such as temperature, water vapor mixing ratio, and/or hydrometeor mixing ratio can be diagnosed from θil. Thus, the use of θil as a predictive variable in numerical models is the most common alternative to “saturation adjustment” techniques in which θ (or temperature) is integrated without accounting for phase changes in one step and then “adjusted” to account for phase changes in a second step. Another advantage of θil is its usefulness in subgrid-scale turbulence closures. Because θil is conserved in turbulent mixing conditions, it is often utilized in first-order “mixing length” turbulence closures.
There are several mathematical formulations of θil. Betts (1973) introduced the first approximate formulation of θl, which has been shown to be sufficiently accurate for studies of shallow convection. Tripoli and Cotton (1981, hereafter referred to as TC81) presented several approximate formulations of θil in order to extend the utility of θl to deep precipitating convection modeling. Through numerical tests, TC81 quantified the accuracy of these formulations under a number of realistic scenarios, including moist parcel ascent and subsaturated parcel descent. Based on their tests, TC81 advocated a semi-empirical formulation of θil that produced the most accurate results. Another formulation of θl was presented by Emanuel (1994), which is exact for reversible evaporation of liquid water in exactly saturated environments.
One purpose of this paper is to present new formulations of θil suitable for use in deep atmospheric numerical models. The most common argument against using θil in numerical models has been the approximations used in their derivations. Thus, we seek to derive new, more accurate formulations by minimizing the assumptions that are made and without invoking Taylor series approximations.
The second purpose of this paper is to evaluate numerically the accuracy of several formulations of θil. The numerical results confirm the unprecedented accuracy of the new formulations. In addition, these numerical tests provide new insight into commonly used formulations and resolve discrepancies presented in previously published results.
2. Previous formulations of θil
Three previously derived formulations of θil are evaluated in this study (a complete list of symbols is provided in the appendix):
where Lυ(T0) and Ls(T0) represent constant values of Lυ and Ls evaluated at a reference temperature T0 = 273.15 K. The subscript notation “,i” where i = 1, 2, 3, … is used here to distinguish the various mathematical formulations of θil. A derivation of (1)–(3) was presented by TC81.
The first formulation, θil,1, has the form of the liquid water potential temperature defined by Betts (1973). Although its problems with large liquid and frozen water contents have been well established (e.g., Wilhelmson 1977; TC81), it is retained in this study because of its common usage (although typically only in studies of shallow cumuli).
The second formulation, θil,2, was derived by TC81 using a different Taylor series approximation than was used by Betts (1973). Hauf and Höller (1987) analyzed this formulation and concluded that it implicitly assumes exact saturation (i.e., relative humidity must always be 100%). Thus, Hauf and Höller argue that the use of this formulation as a governing variable in a numerical model would be “inconsistent” with its derivation, because nonsaturated conditions will surely occur in most applications.
The third formulation, θil,3, is the form advocated by TC81 as the most accurate, based upon results from their numerical tests. The empirical correction “max(T,253)” in θil,3 is argued by TC81 to account adequately for approximations made during the derivation of θil,2, most notable being 1) the neglect of the specific heats of water (i.e., cpv = cl = ci = 0); 2) the neglect of certain tendencies to θil [specifically, the neglect of ɛ1 and ɛ2 by TC81 (p. 1096)]; and 3) the assumption of saturation (according to Hauf and Höller 1987).
Pointin (1984, hereafter referred to as P84) derived two nearly conservative thermodynamic variables—wet equivalent potential temperature and wet equivalent enthalpy—using a minimal number of assumptions. P84 also showed numerically, in a similar manner to TC81, that the accuracy of these two variables is generally acceptable for use in a numerical model. More relevant to the present study is the fact that P84 also included θil,3 in the numerical evaluation. Interestingly, P84 showed errors for θil,3 that are much larger than those shown by TC81. That is, errors in diagnosing θ from θil,3 were shown to be less than 1 K by TC81, but are about −7 K according to P84. It seems that this discrepancy between studies has not been resolved, leaving the accuracy of θil,3 in question.
Based on this analysis, we conclude that a more accurate θil formulation is needed. We derive two new formulations in the following section.
3. New formulations of θil
All conservative thermodynamic variables are derived from the first law of thermodynamics and the Gibbs function (e.g., Iribarne and Godson 1981). Furthermore, as shown by Hauf and Höller (1987), the derivation of most commonly used thermodynamic variables requires additional assumptions and/or specifications. The best conservative thermodynamic variable, Hauf and Höller (1987) argue, is one in which the fewest additional approximations have been made.
With these points in mind, our derivation of new θil formulations begins with Eq. (14) of TC81,
A detailed derivation of (4) has been presented in TC81 (p. 1095) and, in a different manner, in Hauf and Höller (1987, p. 2896). As shown in TC81, the derivation requires one to assume that water vapor moves with the air parcel; assume several physical variables are constant (such as R, cp, and cl); neglect the sensible heat transport by hydrometeor fallout; ignore the surface curvature effects of water drops and ice particles; and invoke a few commonly used equations, such as the Kirchoff and Clausius–Clapeyron equations.
In this study, we do not seek to verify the accuracy or universal applicability of (4). Rather, because (4) is the equation TC81 used to derive their θil formulations, we follow their lead and begin here. One of the main goals of this study is to show that new mathematical formulations of θil can be derived from (4) without neglecting any terms and without using Taylor series approximations, in contrast to previous studies.
The derivation begins by substituting the mathematical definition of θ in the left-hand side of (4), yielding
In this study, the conservation of total water is given by
wherein hydrometeor fallout has been excluded for simplicity. We have conducted a derivation in which fallout terms are included in (6); this derivation includes numerous fallout terms that need to be retained throughout the derivation, but add little insight to the exercise. Furthermore, the ultimate mathematical formulation of θil does not change when the fallout terms are included. Specifically, the inclusion of fallout adds several external flux terms [analogous to Eqs. (33) and (50) of TC81] to the governing equation for θil. For practical applications of θil where hydrometeor fallout is included, we recommend the “adjustment method” of TC81 (p. 1097).
Using (6) and Kirchoff's relations,
the second and third terms on the right-hand side of (5) can be rewritten, yielding
Next, utilizing (7) and the Clausius–Clapeyron equations,
it can be shown that
which can be used to rearrange (8), yielding
Next, we define the relative humidities with respect to liquid (Hl) and with respect to ice (Hi) as
It can be shown using the equation of state and the definition of rυ that
Utilizing (6), the third term on the right-hand side can be rewritten, yielding
Dividing by (cp + cpvrt), pulling terms that are constant into the derivatives [noting that rt is a constant here, due to (6)], and adding
to the left-hand side, we come to the equation
Equation (18) is equivalent to (4), assuming the validity of the relations given here [most notably (7) and (9)]. In contrast to previous derivations of θil, no terms have been neglected, nor have any Taylor series approximations been used.
The entire left-hand side of (18) can be combined into one derivative, that is,
which can be expressed as
If the ice phase is neglected (ri = 0) and if saturation is imposed (Hl = 1), then θil,4 reduces to the liquid water potential temperature presented by Emanuel (1994), and the right-hand side of (22) is zero.
Another formulation for θil can be derived from (21) by incorporating parts of the right-hand side into the formulation of θil, that is,
Equations (22) and (24) are equivalent, so either could be used in a numerical model. However, (24) would be more convenient for most model applications, since the right-hand side of (24) involves time derivatives that are readily available from microphysical schemes. In contrast, some awkward calculations would be necessary to compute the right-hand side of (22).
According to our analysis, θil is a conserved variable only under certain conditions (where, in this context, “conserved” means the time tendency of a variable is zero in the absence of external tendencies, hydrometeor fallout, and turbulent mixing). For the new formulations, the right-hand sides of (22) and (24) make it clear that θil,4 and θil,5 are not always constant when water changes phase, and therefore cannot generally be considered as conserved variables. If the right-hand sides of (22) and (24) are declared to be zero as an approximation, it is uncertain how accurate these two formulations would be under typical atmospheric conditions. To address this concern, numerical tests are presented in the next section.
4. Numerical evaluation
a. Formulation of the model
As in Wilhelmson (1977), TC81, and P84, a one-dimensional numerical model (i.e., parcel model) is used to evaluate numerically the accuracy of the θil formulations. The model used in this paper is simpler than the one used by TC81 and P84 because a detailed microphysics module is not included here for tests of saturated ascent. The exact specification of phase changes has a negligible effect on qualitative conclusions about the comparative accuracy of several θil formulations because realistic distributions of hydrometeors are produced by the model. Furthermore, the results presented here compare well with those presented in other studies, despite the relative simplicity of the model. The goal of simplicity is to make it easier for others to reproduce the results.
For tests based on ascent, a saturated parcel with specified potential temperature is assumed to begin ascent with no liquid water or ice. The parcel is raised from an initial specified pressure to the 100-mb level in 1-mb increments. For temperatures above 273.15 K, supersaturation with respect to liquid water is not allowed, hence
where rvsl is the saturation mixing ratio with respect to liquid water. Below 233.15 K, supersaturation with respect to ice is not allowed, hence
where rvsi is the saturation mixing ratio with respect to ice. At temperatures between 233.15 and 273.15 K, condensation is calculated using a weighted technique similar to that presented by Tao et al. (1989). Specifically, in this temperature range, rυ is diagnosed by
where Fl is a weighting factor for liquid water,
Fi is a weighting factor for ice,
and where 0 ≤ Fl ≤ 1 and 0 ≤ Fi ≤ 1 are enforced. The mixing ratios of ice and liquid water are then diagnosed from
This formulation for condensation ensures a reasonably realistic growth rate of ice between 0° and −40°C, along with a realistic freezing rate of liquid water, without the need for a complicated mixed-phase microphysics scheme.
The latent heats are linear functions of temperature, thereby satisfying Kirchoff's relations:
All specific heats are assumed to be constant.
The saturation mixing ratios for liquid and ice are calculated by
where esl and esi are the saturation vapor pressures with respect to liquid and frozen water, respectively. As in Lipps and Hemler (1980), since Lυ and Ls are linear, the saturation vapor pressures can be defined by simply integrating the Clausius–Clapeyron equation:
For the second set of cases—descent in a subsaturated environment—a saturated parcel begins at a specified pressure and potential temperature with a specified liquid water mixing ratio. The ice phase is neglected for this test. The parcel descends at 1-mb increments with the only governing process being evaporation. After the first pressure increment, the relative humidity is less than 100% during the entire descent. The evaporation formulation from Klemp and Wilhelmson (1978) governs the rate at which the liquid water is evaporated.
Admittedly, some of the choices made in the model are simplifications of reality. For example, there are more accurate formulations for saturation vapor pressure, latent heat, and specific heats based on experimental observations, and the consequences of assuming simple relations are not negligible (e.g., Tag 1980). On the other hand, the approximations made in this section are consistent with those made during derivations in the previous sections. Thus, the choices made here allow for a fair test of the mathematical formulations. However, one could argue that following through on the same approximations predetermines the results [following, e.g., the assertions of Hauf and Höller (1987)]. These are all points that need to be kept in mind when interpreting the numerical results that follow. They might also explain discrepancies between the studies of TC81 and P84.
A total of nine equations (labeled A–I) are evaluated with the numerical code (Table 1). Equation (A) is the benchmark equation utilized by TC81. This equation is also used as the benchmark here. Equation (B) neglects the specific heats of liquid and frozen water, which is commonly done in more complex three-dimensional numerical models. Equation (C) further neglects the specific heat of water vapor.
Equations (D)–(I) govern the various formulations of θil. For (D)–(H), the tendency of θil is explicitly assumed to be zero. Equation (I) is the exact formulation, (22), derived in the previous section. Because (22) and (24) are mathematically equivalent, it is unnecessary to evaluate both.
As in TC81, θ is either predicted from one of the thermodynamic equations [(A)–(C)] or diagnosed from θil [for (D)–(I)]. Since the condensation parameterization is temperature dependent, iteration is required to determine θ. We use a Newton–Raphson technique that converges to within machine accuracy in two–four steps. This efficient retrieval technique is an important component of the new θil formulations, given the numerically expensive exponents in their formulations; an efficient retrieval technique is necessary to make the new formulations computationally feasible in three-dimensional time-dependent numerical models.
Also following the methods of TC81, the equations are integrated in the model using average values for thermodynamic variables when necessary. For example, (A) is integrated by
where the subscripts 1 and 2 refer to values at the last pressure level and new pressure level, respectively, and overbars represent an average of values at levels 1 and 2. The right-hand sides of (B), (C), and (I) are obtained in a similar manner. Since (D)–(H) are diagnostic equations, no such averaging is necessary.
Hydrometeor fallout is not included in any test presented here. As discussed by TC81 and P84, results from tests without fallout overestimate errors that arise in typical conditions, since large hydrometeor mixing ratios are carried to unusual heights in this study. On the other hand, this model design helps reveal the inherent errors of approximate formulations of θil. Furthermore, equations that are shown to be accurate even under such extreme conditions should be preferred over formulations that are accurate only under certain conditions.
b. Results of saturated ascent
Results for parcels rising from 900, 700, and 500 mb with an initial potential temperature of 300 K are presented in Fig. 1. The results from all equations are similar for small displacements (i.e., for a few hundred millibars). However, after lifting to very low pressures, results from the nine equations begin to diverge. These differences are especially large above 300 mb, where nearly all water vapor has been converted to ice.
Also shown are potential temperature errors for (B)– (H) in Fig. 2, where results from (A) are assumed to be the benchmark (or control) solution. Results from (I) are not shown in Fig. 2 because these were found to be identical to those from (A). The total root-mean-square errors during parcel ascent are listed in Table 2, where (A) is assumed as the benchmark, and where errors are calculated every pressure increment during lifting.
As found by TC81, (D) and (E) (i.e., θil,1 and θil,2) produce a considerably higher temperature than does the benchmark equation, (A) (Fig. 2). The differences between θil,1 and θil,2 are lower with smaller liquid water contents and small parcel displacements. Therefore, any of the formulations would be adequate for studies of shallow convection, as argued previously by several studies (e.g., Betts 1973; Deardorff 1976; TC81).
In contrast to the findings of TC81, the results from (F) (θil,3, the formulation advocated by TC81) are several K too low, especially above 350 mb where the temperature falls below 253 K. This is the same qualitative result noted by P84. The two approximate formulations of the thermodynamic equation that ignore the specific heats of water, ice, and/or vapor [i.e., (B) and (C)] also produce temperatures much too low [compared to (A), and compared to the results of TC81]. TC81 indicated that (C) and (F) are accurate to within about 1 K, but our results show much larger departures.1
The discrepancy between these results and those by TC81—which is qualitatively similar to the discrepancy between those by P84 and TC81—might be due to an inadvertent omission of the specific heats of water (i.e., cpv = cl = ci = 0) in the code used by TC81. This conclusion is supported by plotting potential temperature differences when Eq. (C) (the formulation that ignores cpv, cl, and ci) is used as a benchmark. These results (Fig. 3) are quite similar to those reported in TC81. Specifically, the errors for ascent starting at 900 mb from (F) are slightly too high until about 300 mb, and then become too low with further decreases in pressure (Fig. 3a); this is the same qualitative pattern presented for this formulation by TC81. The differences might also be explained by different assumptions made in the numerical codes used by us and by TC81.
Inspection of the mathematical formulation of θil,3 also leads to the same conclusion. If all water vapor is condensed, and if θil, rl, and ri are constant, then (3) cannot predict any further change in θ (assuming T is less than 253 K). For an ascending parcel without hydrometeor fallout, the potential temperature should always continue to increase, owing to the exchange of energy from the hydrometeors to the air, even when all water vapor has been condensed (e.g., Emanuel 1994, p. 133). This behavior is shown in Fig. 1, wherein the potential temperature of the benchmark parcel, (A), continues to increase above 200 mb. In contrast, the potential temperature of the parcel from (F) becomes constant above 200 mb; this is a property of a moist adiabat that does not account for the specific heat of hydrometeors. For these reasons, we conclude that the formulation of θil that has been advocated by TC81 (θil,3) does not adequately account for the specific heats of water.
The results from the new formulations [(G)–(I)] are very encouraging. Results from (I) are identical (to within machine accuracy) to those from (A). Even when the tendency to θil is ignored [i.e., for (G) and (H)], the potential temperature errors reach a maximum of only −0.41 K for the case of a parcel starting at 900 mb. A detailed analysis shows that results from (G) and (H) are identical to results from (A) until the parcel temperature falls below freezing. Since supersaturation with respect to water and subsaturation with respect to ice exist when the temperature is between 273.15 and 233.15 K, errors accumulate for (G) and (H), owing to the neglect of the tendency terms. Nevertheless, despite the approximation that has been made, (G) and (H) are superior to all other formulations of θil that have been presented previously in published literature (Table 2).
c. Results of subsaturated descent
For this test, an initially saturated parcel with a specified amount of liquid water is moved downward while liquid water is evaporated. Relative humidity is less than 100% during the entire descent. Potential temperature errors evaluated against results from (A) are shown in Fig. 4. For descent from 500 and 700 mb (Figs. 4b and 4c), the temperature of the parcel is greater than 253 K at all levels, so the results from (E) and (F) are identical. As before, results from (I) are identical to results from (A), and so are not plotted in Fig. 4. Root-mean-square errors are presented in Table 3.
For this test, results from (F) are comparable to those reported by TC81. The root-mean-square error for the parcel starting at 700 mb is identical (0.18 K) to that by TC81, and is about double for the parcel starting at 500 mb (0.35 K here, compared to 0.15 K by TC81). Perhaps this result is related to differences in model design for this test (e.g., the evaporation formulation is different in this model compared to the formulation used by TC81). However, as TC81 stated, (F) is superior to all other previously analyzed formulations [i.e., compared to (B)–(E)].
One of the new approximate formulations, (G), is less accurate than both (E) and (F) for this test. However, the other new approximate formulation, (H), produces the most accurate results for all of the subsaturated descent tests (Table 3). Interestingly, all other approximate formulations have errors that grow rapidly for small displacements, while errors for (H) grow more gradually (Fig. 4). The main advantage of (H) now becomes clear: if one must assume that the tendency to θil is zero, then (H) (i.e., θil,5) emerges as the preferred formulation. Presumably, (H) tends to be more accurate than (G) because the right-hand side of (24) tends to be smaller in magnitude than the right-hand side of (22).
A review of previous studies suggests the need for new, more accurate formulations of θil. In this paper, we derived two new approximate formulations and their attendant governing equations. The derivation of the new formulations (22)–(25), from the reference equation, (4), was exact, assuming the validity (and applicability) of a minimal number of commonly used relations, such as Kirchoff's relations and the Clausius– Clapeyron equation. Since no terms were dropped from their derivations, and since no Taylor series approximations were invoked, these formulations represent a significant improvement over previous formulations of θil.
A one-dimensional numerical model was used to evaluate the accuracy of three previously presented formulations of θil and the two new formulations. Many of the findings agree with previous studies. However, the semi-empirical formulation of TC81 is found not to properly account for the specific heats of water, and thus tends to have a cold bias (when θ is retrieved from θil). For the new formulations, when the tendency to θil is assumed to be zero (by approximation), then one formulation, (25), produces the best results. When the complete governing equation is used, the numerical tests confirm the equivalency of the new equation and the benchmark equation, (4).
One potential drawback to the use of the new θil formulations is the computationally expensive exponential operators in their formulations. On the other hand, an efficient Newton–Raphson technique can be used to facilitate the retrieval of T (or θ), given a known value of θil. Furthermore, as mentioned by TC81, the percent volume of a numerical model domain that contains liquid or frozen water tends to be small, so no retrieval is needed over most of the domain. Since these computational issues can be overcome, the accuracy of the new formulations make θil an attractive thermodynamic variable to use for numerical applications.
We reiterate that the new formulations are still approximations, except under specific, limited conditions. The approximations and assumptions inherent in their derivations limit their universal applicability. Nevertheless, the new formulations are clearly superior to all previously presented definitions under all tests we have conducted.
We appreciate the insight provided by Dave Stauffer, Jerry Harrington, Richard James, Jason Knievel, Richard Rotunno, and Mark Kelly. We would also like to acknowledge the creative suggestions provided by anonymous reviewers of this and an earlier version of this paper. All figures were created using the Grid Analysis and Display System (GrADS). This work was supported by NSF Grants ATM 9806309 and ATM0133506, and by the Advanced Study Program at NCAR.
List of Symbols
ci Specific heat of ice water at constant pressure (=2106 J kg−1 K−1)
cl Specific heat of liquid water at constant pressure (=4186 J kg−1 K−1)
cp Specific heat of dry air at constant pressure (=1004 J kg−1 K−1)
cpm Specific heat of moist air at constant pressure (cpm ≡ cp + cpvrυ)
cpml Total specific heat at constant pressure (cpml ≡ cp + cpvrυ + clrl + ciri)
cpv Specific heat of water vapor at constant pressure (=1885 J kg−1 K−1)
e Partial pressure of water vapor (Pa)
es0 Reference value of saturation vapor pressure (=611.2 Pa)
esi Saturation vapor pressure with respect to ice (Pa)
esl Saturation vapor pressure with respect to liquid water (Pa)
Fi Weighting factor for saturation mixing ratio with respect to ice
Fl Weighting factor for saturation mixing ratio with respect to liquid
Hi Relative humidity with respect to ice
Hl Relative humidity with respect to liquid
ls1 ls1 ≡ Ls0 + (ci − cpv)T0
ls2 ls2 ≡ ci − cpv
lυ1 lυ1 ≡ Lυ0 + (cl − cpv)T0
lυ2 lυ2 ≡ cl − cpv
Lf Latent heat of freezing (J kg−1)
Ls Latent heat of sublimation (J kg−1)
Ls0 Reference value of Ls (=2.836 × 106 J kg−1)
Lυ Latent heat of vaporization (J kg−1)
Lυ0 Reference value of Lυ (=2.5 × 106 J kg−1)
p Pressure (Pa)
p0 Reference pressure (=1 × 105 Pa)
ri Ice water mixing ratio
rl Liquid water mixing ratio
rt Total water mixing ratio
rυ Water vapor mixing ratio
rvsi Saturation mixing ratio with respect to ice
rvsl Saturation mixing ratio with respect to liquid water
R Gas constant of dry air (=287 J kg−1 K−1)
Rm Gas constant of moist air (Rm ≡ R + Rυrυ)
Rυ Gas constant of water vapor (461 J kg−1 K−1)
Hi Relative humidity with respect to ice
Hl Relative humidity with respect to liquid water
T Temperature (K)
T0 Reference temperature (=273.15 K)
γ Ratio defined by (20)
ɛ Ratio of R to Rυ (ɛ ≡ R/Rυ)
ɛ1 Term ignored in the derivation of θil,2 by TC81
ɛ2 Term ignored in the derivation of θil,2 by TC81
θ Potential temperature (K)
θi Ice water potential temperature (K)
θil Ice–liquid water potential temperature (K)
θl Liquid water potential temperature (K)
χ Ratio defined by (19)
Corresponding author address: George H. Bryan, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307. Email: firstname.lastname@example.org
The National Center for Atmospheric Research is sponsored by the National Science Foundation.