Technical note: A noniterative approach to modelling moist thermodynamics

Formulation of noniterative mathematical expressions for moist thermodynamics presents a challenge for both numerical and theoretical modellers. This technical note offers a simple and efficient tool for approximating two common thermodynamic relationships: temperature, T , at a given pressure, P , along a saturated adiabat, T (P,θw), as well as its corresponding inverse form θw(P,T ), where θw is wet-bulb potential temperature. Our method allows direct calculation of T (P,θw) and θw(P,T ) on a thermodynamic domain bounded by −70≤ θw < 40 C, P > 1 kPa and −100≤ T < 40 C, P > 1 kPa, respectively. The proposed parameterizations offer high accuracy (mean absolute errors of 0.016 and 0.002 C for T (P,θw) and θw(P,T ), respectively) on a notably larger thermodynamic region than previously studied. The paper includes a method summary and a ready-to-use tool to aid atmospheric physicists in their practical applications.


Introduction
Saturated thermodynamics commonly present a challenge for theoretical studies because moist convective condensation, such as deep cumulus precipitation, often involves pseudoadiabtic (irreversible) processes.The latent heat released during water vapour condensation is important for estimating thunderstorm intensity and thickness, precipitation amount and phase, global climate and atmospheric general circulation (Stull, 2017).These processes are governed by nonlinear equations that require iteration to solve.Numerical weather prediction (NWP) models, hence, suffer from the added computational cost to their cloud, precipitation, convection and turbulence schemes and parameterizations because of the iterations required during each timestep of the NWP integration.
A common iterative approach, such as described by Caballero (2014), uses step-wise numerical integration along a saturated adiabat for any constant wet-bulb potential temperature, θ w .The moist adiabatic lapse rate is derived from conservation of moist entropy as a function of temperature, T , and saturated mixing ratio, r s , which itself is a nonlinear function of T and pressure P .To improve efficiency Davies-Jones (2008) proposed a different iterative method, based on inverting Bolton's formula (Bolton, 1980) for equivalent potential temperature valid for the pressure range 10 ≤ P ≤ 105 kPa and wet-bulb potential temperatures −20 < θ w < 40 • C. As a noniterative alternative, Bakhshaii and Stull (2013) offer an approximate solution devised using gene-expression programming.They provide two separate sets of equations for determining T (P , θ w ) and θ w (P , T ) for the domain bounded by −30 < θ w < 40 • C, P > 20 kPa and −60 < T < 40 • C. The complex nature of the problem required their splitting of the modelled region into subdomains, resulting in error discontinuity.The method also produced fairly large errors (on the order of 1-2 • C) in the upper atmosphere.Despite the limitations, to our knowledge the only existing noniterative solution to approximate saturated pseudoadiabats is that by Bakhshaii and Stull (2013).
Our current study presents a different approach for directly calculating T (P , θ w ) and θ w (P , T ) offering improved accuracy for a larger thermodynamic domain.The method, described in Sect.2, normalizes the raw data before fitting it with polynomials.The resultant approximation is evaluated against the "truth" (the iterated solution) and summa- Emagram plot showing select "true" (solid black) and modelled (dashed red) moist adiabats θ w (difference not apparent at this scale).Temperature and pressure domains are restricted for clarity.An emagram (energy mass diagram) is a thermodynamic diagram with the log of pressure on the vertical axis, plotted with maximum and minimum values reversed so that higher in the diagram corresponds to higher in the atmosphere, where pressures are lower.The noniterative results presented in this paper can be plotted on any thermodynamic diagram, including tephigrams and skew-T diagrams.rized in Sects.3 and 4, respectively.As a Supplement we offer the readers a ready-to-use spreadsheet implementing our methodology.
The goal of this paper is to provide a simple tool that can aid analytical modellers in their theoretical work as well as numerical modellers in reducing the computational cost of their simulations.

Data
In order to obtain a set of truth curves for T (P , θ w ) we have used an iterative approach to numerically integrate the equation for dT dP (Tables 1 and 2) for values in the range of −100 ≤ θ w < 100 • C between 105 ≥ P > 1 kPa.The constants used to devise this solution are consistent with Bolton (1980), unless otherwise indicated in Table 1.Note that throughout this technical note we will rely on a common meteorological convention, by which wet-bulb potential temperature at standard pressure of 100 kPa is used to label moist adiabats.Such references, hence, represent curves, rather than constants, and are written in bold for clarity.We found that numerical integration along a saturated adiabat θ w from the bottom to the top of the domain required an increasingly refined pressure step, as all adiabats tend to absolute zero near the top of the atmosphere, and each consecutive pressure step corresponds to a larger temperature jump.For our numerical integration we used 10 −4 kPa step for 105 ≥ P > 10 kPa, 10 −5 kPa step for 10 ≥ P > 2 kPa and 10 −6 kPa step for 2 ≥ P > 1 kPa.The resulting curves (shown on the thermodynamic dia-gram in Fig. 1) are taken as truth, to which we fit our polynomial-based optimization.The noniterative approximations for T (P , θ w ) and θ w (P , T ) described below are valid for thermodynamic ranges bounded by −70 ≤ θ w < 40 • C and −100 ≤ T < 40 • C, respectively.

Approximating T (P , θ w )
While the moist adiabiatic curves θ w in Fig. 1 look smooth and fairly similar, it is challenging for most common optimization routines to capture all of them with one analytical expression with high accuracy.Due to the inherently nonlinear nature of the process, there is no simple way to collapse the curves into a single shape.However, to aid fitting, we can normalize our curves by modelling θ w as a function of a reference moist adiabat θ ref .That allows us to model only the deviations from a reference curve.For our example we used θ ref = −70 • C.This particular choice of θ ref implies no theoretical importance.It is possible to choose any of the directly calculated normalized adiabats to represent θ ref .Depending on the choice, the resulting transformed adiabats shift around the θ ref unity line.The single consideration for choosing a particular θ ref is the ease and accuracy with which it can be fit by a particular optimization tool.
We use polynomial fitting to describe T (P ) for θ ref .This is convenient, since polynomials are generally well-behaved and are computationally easy to use.In particular, they are both continuous and smooth, while being able to capture a wide variety of curve shapes.Moreover, they have well understood properties and a simple form, allowing the model to be easily implemented in a basic spreadsheet.The choice of the degree of polynomial depends on the desired preci- gas constant for water vapour (J K −1 kg −1 ) C pd = 1005.7 specific heat of dry air at constant pressure (J K −1 kg −1 ) T 0 = 273.15reference temperature (K) P 0 = 100 reference pressure (kPa) e 0 = 0.611657 Clausius-Clapeyron constant (kPa) (Koutsoyiannis, 2012) = 0.6220 ratio of gas constants (kg kg −1 ) Table 2. Variable definitions.

Variable
Description (unit) θ w saturated adiabat where the value of T defined at P = P 0 is defined as wet-bulb potential temperature (K) saturation vapour pressure (kPa) saturation mixing ratio (kg kg −1 ) ) change of temperature with pressure along a saturated adiabat, which can be iterated to find T vs. P (K kPa −1 ) sion level.Since we are examining a fixed range of temperatures relevant to atmospheric applications, the potentially extreme oscillatory behaviour of high-degree polynomials outside of the modelled domain is not a primary concern.The fitted polynomials have no predictive value outside of the modelled range and serve purely as an interpolation function.While the large number of possible inflection points associated with high-degree polynomials may be of a concern near the edges of the fitting interval, a problem known as Runge's phenomenon (Epperson, 1987), the current algorithm relies on a least-squares method to minimize the effect and achieve a high-quality fit.For this example, the aim was to ensure that the mean absolute error (MAE) is on the order of 10 −2 • C, requiring a 20th degree polynomial to achieve such fit.The coefficients for this polynomial are provided in Table S2 in the Supplement.The next step is to choose a single functional form to represent the entire family of the transformed curves (i.e., the moist adiabat deviations from θ ref ).Each given shape of a particular curve is then controlled by variable parameters of the same function.A number of simple functions exists that are able to model the above relationship.For this work we tested biexponential, arctan, rational and polynomial functions.Generally, a reasonable (on the order of 1-2 • C) fit can be achieved with both biexponential and arctan functions using as few as three variable parameters.While efficient, results of such a fit are unlikely to be sufficiently accurate to be useful for real-life modelling applications and, more importantly, only constitute a part of the solution.The bigger concern with these choices is that, unlike polynomials, they produce variable parameters that do not appear wellbehaved.Discontinuity and asymptomatic behaviour arising from error minimization for all transformed adiabats renders the parameter curves very difficult to model.A variety of functions would be necessary to capture the parameter behaviour, which in turn is likely to produce a complex and discontinuous error field, such as appeared in Bakhshaii and Stull (2013).
Polynomial fitting does not appear to suffer from such issues.Moreover, the accuracy can be controlled by changing the degree of the polynomial and, hence, allowing a higher number of variable parameters.In this example, the curves were modelled using 10th degree polynomials, resulting in 11 variable parameters.Conveniently, and unlike other functional forms mentioned above, these parameters are also well-behaved.They can, again, be modelled using high-degree polynomials to the desired level of accuracy.Results of parameter fitting for this given example were again produced using 20th degree polynomials, with fit coefficients provided in Table S1 (Supplement).The resulting modelled (noniterative) moist adiabats can be seen in Fig. 1, compared to the truth (iterated) values.Figure 2. Approximation error between iterated ("truth") and modelled T along moist adiabats θ w .

Approximating θ w (P , T )
A similar approach can be used to produce a noniterative approximation for θ w (P , T ).S4 in the Supplement).We then produce a set of transformed curves by plotting the isotherms as a function of T ref .We fit the transformed curves with 10th degree polynomials, obtaining a dataset for 11 variable parameters.Finally, we use polynomials to model the variable parameters (Table S3 in the Supplement).The following section discusses the results and accuracy of our optimization procedure.

Evaluation
To test the accuracy of the proposed method, we compared our modelled curves for T (P , θ w ) and θ w (P , T ) with those obtained through direct calculation (the truth iterative solution).The results of the evaluation for T (P , θ w ) are shown in Fig. 2, indicating errors on the order of a few hundredths of a degree throughout most of the domain.Warmer values near the top of the domain tend to be modelled least accurately.MAE for the entire modelled thermodynamic region is 0.016 • C. Error contours for θ w (P , T ) are shown in Fig. 3, with errors on the order of a few thousandths of a degree throughout most of the domain and overall MAE = 0.002 • C. Once again, values near the low-pressure limit tend to be least accurate.Notably, applying the above optimization on a slightly shallower pressure domain of P > 2 kPa allows im- provement of the overall MAE for both approximations by an additional order of magnitude.As mentioned earlier, improved accuracy may also be achieved with the use of even higher degrees of polynomials for parameter fits.However, such precision is unlikely to be necessary, as some of the thermodynamic relationships used in the truth iterative computations contain substantially larger errors than those introduced by the above optimization procedure (Davies-Jones, 2009;Koutsoyiannis, 2012).Moreover, conventional pseudoadiabatic diagrams, such as those used by US Air Force, Environment Canada and Air Transport Association of America, differ from each other by nearly 1 • C at the 20 kPa pressure level (Bakhshaii and Stull, 2013).The specific choice of thermodynamic constants and relationships undoubtedly affects the definition of truth used in this work; however, this has little effect on the overall validity of the approach.Should more precise constant values and/or thermodynamic relationships become established in the future, the proposed method can be reapplied to generate updated fitting coefficients without loss of accuracy (within the limits of the specific polynomial optimization routine used).
Though the upper 10 kPa of the atmosphere contains the largest errors with our proposed approach, this vertical subrange also presents the most significant challenge for direct (iterative) numerical modelling.Accurate numerical computation requires an increasingly refined vertical step for the top part of the atmosphere.Hence, despite the errors, the proposed approximation offers a more accurate solution than one would obtain with a direct iterative approach using a somewhat coarse yet computationally demanding 0.001 kPa pressure step.
While common weather phenomena generally remain in the troposphere, the validity of the current method on a notably larger vertical domain is particularly useful in the lower latitudes.The deep vertical extent of tropical thun-derstorms, hurricanes and typhoons in combination with the high tropopause altitude (18 km or 8 kPa; WMO, 2017) in the tropics contribute to large computational costs of modelling these potentially destructive events.

Summary of approach
Individual steps to directly compute T (P , θ w ) and θ w (P , T ) are summarized below.This sample procedure, along with the required coefficient tables, is provided in a ready-to-use form in the attached spreadsheet (Supplement).

Computing T (P , θ w )
Let n = 0, . .., 10 correspond to the index of individual polynomial coefficients and m = 20 be the degree of polynomial fits for θ ref (P ) and k n (θ w ), respectively.
3. Compute T (θ ref ): where T and θ ref are in Kelvins, and values of k 0,...,n correspond to polynomial coefficients calculated in Step 1.

Computing θ w (P , T )
Let n = 0, . .., 10 correspond to the index of individual polynomial coefficients and m = 20 be the degree of polynomial fits for T ref (P ) and κ n (T ), respectively.
3. Compute θ w (T ref ): where θ w and T ref are in • C, and values of κ 0,...,n correspond to polynomial coefficients calculated in Step 1.

Usage example
Meteorologists typically use both θ w (P , T ) and T (P , θ w ) for moist convection such as thunderstorms, frontal clouds, mountain-wave clouds and many other phenomena, where a saturated air parcel moves vertically.The cloud base of convective clouds marks the bottom of saturated ascent, and the cloud top marks the top.
For example, suppose that the forecast at some tropical weather station is P = 100 kPa and T = 32 • C with dew point T d = 21 • C (corresponding to a water vapour mixing ratio of approximately r = 16 g kg −1 ).Further suppose that a force (e.g., buoyancy, frontal uplift, orographic uplift) causes an air parcel with these initial conditions to rise.Initially this air parcel is unsaturated (not cloudy), so we do not need to use the polynomial or iterative equations.Instead, simpler noniterative equations apply for the thermodynamic state as the parcel rises dry adiabatically.Namely, its temperature cools at the dry adiabatic lapse rate (9.8 • C km −1 ), and the mixing ratio and potential temperature are constant.This air parcel will become saturated (i.e., cloud base) at the lifting condensation level (LCL).With this information, other thermodynamic equations (Stull, 2017) can be used to find conditions at the LCL: z LCL = 1.375 km, P LCL = 85.4 kPa and T LCL = 18.5 • C.
Given this initial P and T at the LCL, we can use the polynomial equations provided in this paper to compute which moist adiabat the cloudy air parcel will follow: θ w (P , T ) = 24.0• C. If this cloudy air parcel (still following the θ w (P , T ) = 24.0• C adiabat) rises to an altitude where the pressure is P = 24.0kPa, then we can use the second set of polynomial equations in this paper to find the final temperature of the air parcel at this new height: T (P , θ w ) = −39.8• C.

Discussion and conclusions
The polynomial method proposed here is accurate, smooth and computationally efficient.For example, given the cloud www.atmos-chem-phys.net/17/15037/2017/Atmos.Chem.Phys., 17, 15037-15043, 2017 base and cloud top pressures of the previous example, the tally of computer operations to find both the initial and the final temperature is 230 additions and subtractions and 2365 multiplications (where rational numbers to integer powers are counted as sequential multiplies).This can be compared to the computation tally for the truth iterative solution, which requires a total of 2 750 000 variable pressure steps, where each step has 8 additions and subtractions, 17 multiplies (where rational numbers to integer powers are counted as sequential multiplies), 9 divides and 2 math functions (e.g., log, exp, non-integer exponents), totalling to 988 200 000 operations from the bottom to the top of the domain.Also, for comparison, some NWP models use a look-up table to get the average saturated adiabatic lapse rate θ w / P as a function of P and T .While this method is fairly fast, it is also less accurate and approximates the saturated lapse rate as a series of short straight-line segments instead of a smooth curve.It also has discontinuous jumps of saturated lapse rate as T varies along an isobar.
While interpolating values from look-up tables generally results in random errors, iterative solutions with a coarse step could potentially suffer from a directional drift due to numerical integration errors, which may introduce a consistent bias into latent heating profiles.Moreover, near the top of the atmosphere, where each pressure step corresponds to a large temperature jump along the moist adiabats the numerical solutions tend to become unstable.Though both of these concerns are addressed with the proposed low-cost polynomial method, the broader challenge of our limited overall understanding of moist convection remains.Existing thermodynamic relationships are based on the assumption of either a reversible moist adiabatic or an irreversible pseudoadiabatic process.Real-world atmospheric processes are likely to be a combination of both (Iribarne and Godson, 1981).The uncertainty introduced by our limited knowledge of the true state of saturated air is likely to remain the central obstacle in capturing moist convection.
The polynomial method proposed here provides a computation of high accuracy and smooth variation across the whole thermodynamic diagram range, at intermediate computation speed compared to the other methods.Moreover, it helps to model moist thermodynamics on a wider temperature range with roughly 2 orders of magnitude MAE improvement over the existing solution.In addition to the reduced computational costs of obtaining solutions for T (P , θ w ) and θ w (P , T ) in numerical simulations and improving accuracy, we hope that our tool will aid analytical modellers in their theoretical work.
Data availability.Spreadsheet tool for calculating T (P , θ w ) and θ w (P , T ) is included as Supplement to this article.
Figure1.Emagram plot showing select "true" (solid black) and modelled (dashed red) moist adiabats θ w (difference not apparent at this scale).Temperature and pressure domains are restricted for clarity.An emagram (energy mass diagram) is a thermodynamic diagram with the log of pressure on the vertical axis, plotted with maximum and minimum values reversed so that higher in the diagram corresponds to higher in the atmosphere, where pressures are lower.The noniterative results presented in this paper can be plotted on any thermodynamic diagram, including tephigrams and skew-T diagrams.
To obtain a new set of curves representing lines of constant temperature T in θ w domain, we have used our existing dataset for −100 ≤ θ w < 100 • C to extract isotherms on a 0.5 • C and 0.1 kPa grid for −100 ≤ T < 40 • C and 105 ≥ P > 1 kPa.Similarly to our earlier approach, we select a single reference curve T ref = −100 • C and use a high-order polynomial to model it as a function of pressure (Table

Table 1 .
Table of constants.