Assessment of algorithms for computing moist available potential energy

Atmospheric moist available potential energy (MAPE) has been traditionally defined as the potential energy of a moist atmosphere relative to that of the adiabatically sorted reference state defining a global potential energy minimum. Finding such a reference state was recently shown to be a linear assignment problem, and therefore exactly solvable. However, this is computationally extremely expensive, so there has been much interest in developing heuristic methods for computing MAPE in practice. Comparisons of the accuracy of such approximate algorithms have so far been limited to a small number of test cases; this work provides an assessment of the algorithms' performance across a wide range of atmospheric soundings, in two different locations. We determine that the divide-and-conquer algorithm is the best suited to practical application, but suffers from the previously overlooked shortcoming that it can produce a reference state with higher potential energy than the actual state, resulting in a negative value of MAPE. Additionally, we show that it is possible to construct an algorithm exploiting a theoretical expression linking MAPE to Convective Available Potential Energy (CAPE) previously derived by Kerry Emanuel. This approach has a similar accuracy to existing approximate sorting algorithms, whilst providing greater insight into the physical source of MAPE. In light of these results, we discuss how to make progress towards constructing a satisfactory moist APE theory for the atmosphere. We also outline a method for vectorising the adiabatic lifting of moist air parcels, which increases the computational efficiency of algorithms for calculating MAPE, and could be used for other applications such as convection schemes.


Introduction
Available Potential Energy (APE) theory, as originally outlined by Lorenz (1955), provides a framework to study the energy available to atmospheric motions. The theory is underpinned by the concept of an atmospheric background or reference state.
Such a state has been traditionally envisioned as being obtained through an adiabatic mass rearrangement such that the sum of the internal and potential energies of the atmosphere (total potential energy) is minimised. The APE is then found as the difference between the total potential energy of the atmosphere and the total potential energy of the reference state. In its reference state, the atmosphere is at rest and in hydrostatic equilibrium; its density stratification is therefore statically stable and horizontally uniform, and no further conversion with kinetic energy can take place. The APE thus gives the total potential energy that is available for reversible conversions into kinetic energy. Assuming hydrostatic balance, minimisation of the total potential energy is equivalent to minimisation of the enthalpy H, so that where h is the specific enthalpy, and the integral is over all the mass in the considered atmospheric domain.
For a moist atmosphere, the rearrangements are made via reversible adiabatic processes conserving total water content (Lorenz 1978). In this work, we refer to the APE of a moist atmosphere as MAPE (Moist Available Potential Energy), following the terminology of Stansifer et al. (2017), and we focus only on the vertical component of MAPE. Unlike the dry case, for which reference pressure is uniquely determined by sorting potential temperature, there is no known analytical solution for obtaining the moist reference state from the distribution of entropy and specific humidity. As a result, previous methods of calculating MAPE have relied on heuristic approaches involving discretising atmospheric domains into parcels of equal mass and sorting them according to density at differing pressure levels to obtain a reference state. From a computational viewpoint, the discretised approach to computing MAPE is equivalent to finding the permutation of the actual state with the lowest total potential energy. Tailleux and Grandpeix (2004) characterised such a problem as an asymmetric travelling salesman problem, but recently, it was realised by Hieronymus and Nycander (2015) that the computation of such a reference state was in fact a linear assignment problem that can be solved by using the Munkres algorithm (Munkres 1957). Whilst the Munkres algorithm is exact, it is also computationally expensive, and therefore it is still desirable to use approximate algorithms for speed. The time taken for the algorithms to compute MAPE is explored in more detail in Appendix A. The Munkres algorithm can be used when the considered atmospheric domain comprises a small number of parcels n, but the runtime of the algorithm increases as n 3 (Stansifer et al. 2017), so it quickly becomes infeasible for large domains.
Approximate sorting algorithms have been employed to investigate the intensity of extratropical storm tracks (O'Gorman 2010), using Lorenz's algorithm (Lorenz 1979) to calculate MAPE, and the energetics of tropical cyclones (Wong et al. 2016

Algorithms for Computing MAPE
In this section we describe the algorithms that may be used to compute the reference state, and hence the MAPE, of an atmospheric sounding. We assume here that the sounding has been discretised into parcels of equal mass. To begin, we outline the Munkres algorithm, which finds the reference state corresponding to the exact minimum enthalpy rearrangement of the parcels.
We then describe the parcel-sorting algorithms that have been designed to find approximations to the reference state. Due to their approximate nature, these methods are less computationally expensive than the Munkres algorithm, but their typical accuracy compared to the Munkres algorithm is unknown; this will be investigated in Section 4. Finally we describe a method for calculating MAPE that does not rely on a sorting procedure, but instead makes use of the relationship between MAPE and CAPE, which was suggested by Emanuel (1994).

Munkres algorithm
The Munkres algorithm (Munkres 1957) may be used to obtain the exact minimum enthalpy rearrangement of a set of air parcels, by treating the computation of the parcels' reference pressures as a linear assignment problem (Hieronymus and Nycander 2015;Stansifer et al. 2017). This method first calculates a cost matrix C, in which the entry c ij is the enthalpy of the i th parcel at the j th pressure level. Using this cost matrix, the algorithm allocates parcels to the pressure levels resulting in a minimised total enthalpy. This is done by using the linear algebra procedure described by Munkres (1957), which tracks how difficult it is to find a low-enthalpy position for each parcel during the rearrangement process.

Lorenz's algorithm
The first algorithm for approximating the minimised enthalpy reference state of a moist sounding was developed by Lorenz (1979). For a set of n parcels at pressures p 1 < p 2 < . . . < pn, this algorithm begins by calculating the virtual temperatures of all parcels as if they were lifted reversibly and adiabatically to p 1 , denoted T v1 , and if they were similarly lifted to pn, denoted Tvn.
The algorithm first finds a parcel to assign to pressure level p 1 , and then moves to progressively higher pressures. This assignment is determined as follows: at each level p j , the unassigned parcels

Top-down algorithm
The top-down algorithm was used to compute reference states in the study of APE in tropical cyclones by Wong et al. (2016).
The performance of the top-down algorithm was also analysed by Stansifer et al. (2017), who referred to it as the "greedy algorithm". The top-down algorithm for n parcels proceeds as follows: all n air parcels are moved reversibly adiabatically to p 1 , the lowest pressure in the sounding. Their densities at this pressure are calculated, and the parcel with the lowest density is assigned to have pn as its reference pressure. This parcel is then eliminated from sorting. The remaining n − 1 parcels are moved to p 2 , and again their densities are calculated, and the least dense parcel assigned to p 2 . The algorithm continues in this way until all parcels have been assigned to a reference pressure level.

Bottom-up algorithm
Bottom-up sorting works similarly to top-down sorting, but the parcels are first moved to the highest pressure pn, assigning the parcel with the highest density to this level, and proceeding to lower pressure levels p n−1 , p n−2 . . . . Bottom-up sorting was suggested by Wong et al. (2016) to limit the inclusion of CAPE in the definition of MAPE. This may be desirable in practice since not all the CAPE present in the atmosphere will be released, for example due to the presence of Convective Inhibition (CIN).

Divide-and-conquer algorithm
The divide-and-conquer algorithm was introduced by  Emanuel (1994), which we will compare (in Section 4) to the MAPE computed by the sorting algorithms described above. Emanuel (1994) supposes that MAPE is due solely to the presence of CAPE in a thin boundary layer of depth ∆p b . In this case an approximation to the MAPE is given by where CAPE b is the mean CAPE in the boundary layer, κ = R d c pd , and ∆θv is the change in the virtual potential temperature between the top of the boundary layer at p b,top and the boundary layer's level of neutral buoyancy, p LNB . The overbar denotes a θv-weighted average from p b,top to p LNB . The first term of Eq.
(2) corresponds to the release of CAPE when the boundary layer rises upwards to its LNB. The second term accounts for the energy change that occurs as a result of the remaining air parcels descending by ∆p b . We will henceforth refer to MAPE calculated using Eq.
We compute the Emanuel MAPE by calculating the value of Eq.
(2) for ∆p b depths ranging from 0 mb to 150 mb, and selecting the maximum value of MAPE returned by any of these ∆p b values. We increment ∆p b simply by including the next lowest parcel in the sounding. Theoretically, it would be possible to use smaller increments in ∆p b , and include fractions of parcels in the boundary layer. We have not done this because the sorting algorithms discussed earlier in this section are only able to rearrange whole parcels, so allowing this CAPE-based algorithm to only lift whole parcels provides a fairer comparison of the MAPE.
To compute the boundary layer CAPE, CAPE b , we use a parcel with a value of θ given by the pressure-weighted mean of θ in the boundary layer, and q given by the mean q in the boundary layer.
The CAPE is then where αp is the specific volume of the parcel when it is lifted reversibly adiabatically, and αe is the environmental specific volume. The parcel is lifted from its initial position p i , which we take to be the bottom of the boundary layer (i.e. the surface), to its highest level of neutral buoyancy.

Data
To calculate the MAPE of a sounding, the sorting algorithms outlined in Section 2 require the input of the temperature, pressure and total specific humidity profiles. The atmospheric profiles used to compare the algorithms are data obtained through soundings from the Atmospheric Radiation Measurement (ARM) Program (Stokes and Schwartz 1994). We assume that the total specific humidity q T in the soundings is equal to the specific humidity q, i.e. that no liquid water is present in the atmosphere. This widens our choice of data since we do not require liquid water measurements, and is justified since we do not expect large quantities of liquid water to be residing in the atmosphere for long periods of time.
We have used soundings from Nauru dating from 1 April 2001 to 16 August 2006. These soundings contain data that have been interpolated onto 5 mb pressure levels and quality controlled as described by Holloway and Neelin (2009). We take all soundings with at least 150 valid measurements of temperature and specific humidity, for which the valid measurements span at least the interval from 1000 mb to 100 mb. Any missing temperature or humidity measurements are filled in by linear interpolation.
This results in 3130 soundings for which we can use the sorting algorithms to compute the MAPE in the 1000 mb to 100 mb layer using 181 parcels of 5 mb depth.
To verify whether the performance of the algorithms is significantly affected if the soundings are from a different location, we have also used soundings from the ARM Southern Great Plains (SGP) sites during the Intensive Observation Period from 4 June 1997 to 7 July 1997; this dataset is the one used by Tailleux and Grandpeix (2004  The median profiles of temperature T and specific humidity q are shown for each location in Figure 1, along with the 25 th to 75 th percentiles (dark shading) and 10 th to 90 th percentiles (light shading). The profiles are similar in the two locations, with Nauru soundings exhibiting higher moisture at lower levels (this is reasonable because we have kept Nauru data at higher pressure levels, whereas there were insufficient measurements to do so for the SGP data). The Nauru soundings also show colder temperatures at high altitude. It is notable that there is very little variation about the median Nauru temperature profile, and therefore differences in the ability of the algorithms to accurately calculate MAPE here will be mostly due to the differences in humidity profiles between the soundings.

Comparison of Algorithms
All the sorting algorithms discussed in Section 2 were used to calculate the MAPE of each of the 3714 ARM soundings described in Section 3. To summarise, these algorithms are: To quantify the accuracy of each algorithm we define the percentage relative difference in MAPE as    than relying on the physically unconstrained rearrangements of a sorting algorithm, as will be discussed further in Section 5.

Discussion
The results presented in Section 4 allow us to make a more informed assessment of which algorithms are most suitable for the computation of MAPE, based on an analysis of a wider range of soundings than in previous studies. The key challenge for MAPE algorithms stems from the fact that MAPE is ultimately a residual between the positive work due to the release of CAPE minus the negative work due to compensating subsidence. As a result, sorting a vertical sounding according to decreasing density, which is the approach underlying the majority of algorithms, may occasionally result in a reference state with a larger potential energy than the actual state, if the negative work exceeds the positive work. This is in contrast to the case of a dry atmosphere, for which sorting the actual state according to potential temperature always returns the state of minimum potential energy. Without an explicit procedure to forbid it, most heuristics for computing MAPE are bound to return a negative value in some cases. That such situations may occur in practice appears to have been overlooked in previous studies, with the exception of Randall and Wang (1992), but is clearly established for the particular soundings analysed here.
In terms of performance, we have found that the Lorenz and top-down algorithms have nearly identical levels of accuracy, which we did not anticipate. However, our results also indicate that both algorithms are so prone to returning a negative MAPE that they are not suited to practical application. Such an issue was not mentioned in previous studies using these algorithms So far, the implicit assumption of the present study and others has been that it is legitimate or most useful to define the APE of a moist atmosphere in terms of the reference state that defines the absolute minimum in potential energy, but this is not necessarily the case. For a moist atmosphere, it is a priori possible to construct alternative sorted reference states that define only a local minimum in potential energy. Although such reference states would result in a lesser global value of APE, it is unclear why this would necessarily invalidate their use. In tropical cyclones, for instance, numerical simulations reveal that boundary layer parcels away from the eyewall may have CAPE whose release is suppressed by the subsidence in that region, as pointed out by Wong et al. (2016). Since the CAPE of such parcels can rarely if ever be released, it is unclear why it should be included in the definition of a tropical cyclone APE, as will normally be the case if the reference state defining a global potential energy minimum is selected. From a practical viewpoint, it is important to remark that the choice of reference state affects the overall value of APE as well as its diabatic generation rate G(AP E), but neither affects the energy conversion between APE and kinetic energy, nor the general form of the APE evolution equation, given by: where the index i is used to indicate dependence on the reference state chosen. Eq. (5) states that the conversion C(KE, AP E) between kinetic energy and APE always appears as a residual between the APE storage term dAP E i /dt and the APE generation rate G(AP E) i (see Pauluis (2007) for a discussion of how moist processes may affect the latter). From a theoretical viewpoint, Eq. (5) represents a balance between three terms, of which the storage term is the least interesting or meaningful. For this reason, Wong et al. (2016) argued that the reference state should be chosen so as to minimise the storage term, in order to potentially make it possible to predict the APE/KE conversion from the knowledge of the APE generation rate. In this regard, Wong et al. (2016) found the use of the bottom-up sorted reference state to yield a lower storage term that the top-down sorted reference state, but more research is required to establish whether this can be regarded as a general result.
Given the computational and conceptual difficulties entailing their use, it is important to question whether sorting algorithms are really needed to study the energetics of a moist atmosphere.
The idea that simpler alternatives might exist is indeed justified by the fact that some recent APE studies successfully moved away from the use of sorting algorithms by resorting to p.d.f. approaches instead, as in the case of Saenz et al. (2015), itself an extension of Tseng and Ferziger (2001), although it is unclear how such a method could be applied to a moist atmosphere.
Also, it has long been known from the works of Andrews (1981) and Holliday and Mcintyre (1981) that it is possible to construct a local theory of APE based on an arbitrary reference state defined by a reference pressure p 0 (z, t) and specific volume α 0 (z, t) in hydrostatic equilibrium. Based on Tailleux (2013) and Novak and Tailleux (2017), this would lead one to define the APE density for a moist atmosphere as the work that a fluid parcel needs to perform to move from its reference pressure pr to its actual pressure p, viz., where θ l is liquid potential temperature and q T is total water content. An alternative formulation for APE density in a compressible atmosphere, based on a modified potential temperature, has been proposed by Peng et al. (2015). In contrast to what is often assumed, a sorting algorithm is not required to calculate the reference pressure pr. Indeed, as shown by Tailleux (2013), if α 0 (p, t) is known at all times as a function of pressure, pr can be simply estimated by solving the so-called Level of Neutral Buoyancy (LNB) equation: α(θ l , q T , pr) = α 0 (pr, t).
This corresponds to the use of a level of neutral buoyancy in the

A. Increasing Algorithm Efficiency
All the algorithms used to compute MAPE in this work require the repeated calculation of the temperature of an air parcel when it moves reversibly adiabatically to a given pressure level. For the approximate sorting algorithms, this calculation is required in order to sort the parcels by density at each pressure level of the sorting process. The Emanuel algorithm requires the calculation in order to find the specific volume of lifted parcels for the calculation of CAPE (see Eq. (3)). For the Munkres algorithm, it is required in order to compute the cost matrix C, where c ij is the enthalpy of the i th parcel if it were moved to the j th pressure level.
This cost matrix is then manipulated to find the minimum enthalpy rearrangement of the parcels. Such temperature calculations are a time-consuming stage of the algorithms because they require the use of a root-finding procedure. In this appendix we show how the bisection method can be employed to enable vectorised calculation of the new temperatures, so that the runtime of the algorithms is reduced by moving many parcels to new pressure levels at once.
The problem can be formulated as follows: we know the initial temperature T ( • C), pressure p (mb) and total specific humidity q T (kgkg -1 ) of the parcel, and the target pressure p ′ . We must find T * , the temperature such that ∆s = s(T, p, q T ) − s(T * , p ′ , q T ) = 0, where s (Jkg -1 K -1 ) is the specific entropy of the parcel (so that the movement of the parcel is adiabatic). Using a built-in root-finding function in Python (such as Brent's method) to obtain the value of T * , as has been done previously, necessitates looping through all parcels individually to calculate their temperatures, since built-in root-finding functions can compute the root of only one function at a time.
Alternatively, it is possible to use a simple bisection method to compute T * . This converges less quickly than other rootfinding procedures (its convergence is linear), but makes it easy to compute T * for many parcels at once. We initially take a wide temperature interval from 0 K to 1000 K. This interval is successively halved, at each stage identifying which half T * lies in by computing the sign of ∆s at the midpoint of the interval.
The bisection method is guaranteed to converge to the root T * provided that the function ∆s(T ) is monotonic, ∆s (0 K) < 0, and ∆s (1000 K) > 0. The number of iterations can be adapted depending on the required accuracy of T * ; we have ensured convergence to within 10 −6 K.
To assess the impact of using the bisection method in this way, we have created two versions of the Munkres algorithm provided by Stansifer et al. (2017): one that calculates the cost matrix by looping through all parcels and using Brent's method to compute temperature changes (similar to Stansifer et al.'s original algorithm); and one that uses the bisection method to compute one row of the cost matrix at a time. We use both these algorithms to compute the MAPE for soundings with a varying number of parcels. For each number of parcels n we create the desired sounding by taking the ARM Nauru data from 1200 UTC on 1 April 2001 and linearly interpolating it to n pressure levels.
The time taken to compute the MAPE on a personal computer is displayed in Figure 5; in all cases it was verified that the two algorithms computed the same MAPE. It is clear that use of the bisection method results in a much faster algorithm for large numbers of parcels.
The bisection method can be similarly employed for the approximate sorting algorithms. Figure 6 shows the time taken to compute MAPE for our implementations of each of the algorithms described in the main body of the paper, for a varying number of parcels. As expected, the exact Munkres procedure is slowest for large numbers of parcels. The divide-and-conquer algorithm of Stansifer et al. (2017) is easily the fastest, making it a good compromise between speed and accuracy.