A noniteratiave approach to modelling moist thermodynamics

Abstract. 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, re5 spectively. The proposed parameterizations offer high accuracy (mean absolute errors of 0.017°C 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, as well as 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 valuable noniterative alternative, Bakhshaii and Stull (2013) offer an approximate solution devised using gene-expression programming (GEP).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 sub-domains, 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 Bakhshaii and Stull (2013) is the only existing noniterative solution to approximate saturated pseudoadiabats.
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 Section 2, normalizes the raw data before fitting it with polynomials.The resultant approximation is evaluated against the "truth" (the iterated solution) and summarized in Sections 3 and 4, respectively.As Supplementary Material 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 bolded 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 thermo diagram in Figure 1) are taken as "truth", to which we fit our polynomial-based optimization.The non-iterative 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 Figure 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 non-linear 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 precision level.Since we are examining a fixed range of temperatures relevant to atmospheric applications, the potentially extreme oscillatory behavior 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 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 II in Supplementary Material.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 bi-exponential and arctan functions using as little as three variable parameters.While efficient, the results of such 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 well-behaved.Discontinuity and asymptomatic behavior 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 behavior, which in turn is likely to produce a complex and discontinuous error field, such as appeared in Bakhshaii and Stull (2013).
Polynomial fitting doesn't 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 I (Supplementary Material).The resulting modelled (non-iterative) moist adiabats can be seen in Figure 1, compared to the truth (iterated) values.

Approximating θ w (P, T )
A similar approach can be used to produce a non-iterative approximation for θ w (P, T ).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 IV in Supplementary Material).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 III in Supplementary Material).
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 Figure 2, indicating errors on the order of few hundredths of a degree throughout most of the domain.Warmer values near the top of the domain tend to be modelled least accurately.Mean absolute error (MAE) for the entire modelled thermodynamic region is 0.016°C.Error contours for θ w (P, T ) are shown in Figure 3, with errors on the order of 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 improvement 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 U.S.
Air Force (USAF), Environment Canada (EC) and Air Transport Association of America (ATAA), 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 effects the definition of "truth" used in this work, however, 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 re-applied 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 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.Deep vertical extent of tropical thunderstorms, 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.
Individual steps to directly compute T (P, θ w ) and θ w (P, T ) are summarized below.This sample procedure, along with the required coefficient tables are provided in a ready-to-use form in the attached spreadsheet (Supplementary Material).

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.1) Compute coefficients k n (θ w ) using polynomial coefficients a 20 , ..., a 0 in Table I in Supplementary Material: for θ w in °C.
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, mountainwave clouds, and many other phenomena where a saturated air parcel moves vertically.Cloud base of convective clouds marks the bottom of saturated ascent, and cloud top marks the top.
For example, suppose that the forecast at some tropical weather station is P = 100 kPa, T = 32°C with dewpoint T d = 21°C (corresponding to a water vapor mixing ratio of approximately r = 16 g kg −1 ).Further suppose that a force (e.g., buoyancy, frontal uplift, or orographic uplift) causes an air parcel with these initial conditions to rise.Initially this air parcel is unsaturated (not cloudy), so we don't need to use the polynomial or iterative equations.Instead, simpler non-iterative 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°Cadiabat) 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 base and cloud top pressures of the previous example, the tally of computer operations to find both the initial and the final temperature are: 230 additions and subtractions, 2365 multiplies (where rational numbers to integer powers are counted as sequential multiplies).Compare that to the computation tally for the "truth" iterative solution, requiring 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 numerical weather prediction 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 two 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.