pyPI (v1.3): Tropical Cyclone Potential Intensity Calculations in Python

Potential intensity (PI) is the maximum speed limit of a tropical cyclone found by modeling the storm as a thermal heat engine. Because there are significant correlations between PI and actual storm wind speeds, PI is a useful diagnostic for evaluating or predicting tropical cyclone intensity climatology and variability. Previous studies have calculated PI given a set of atmospheric and oceanographic conditions, but although a PI algorithm—originally developed by Kerry Emanuel—is in widespread use, it remains under-documented. The Tropical Cyclone Potential Intensity Calculations in Python (pyPI, v1.3) 5 package develops the PI algorithm in Python, and for the first time details the full background and algorithm (line-by-line) used to compute tropical cyclone potential intensity constrained by thermodynamics. The pyPI package (1) provides a freely available, flexible, validated Python PI algorithm, (2) carefully documents the PI algorithm and its Python implementation, and (3) demonstrates and encourages the use of PI theory in tropical cyclone analyses. Validation shows pyPI output is nearly identical to the previous potential intensity computation, but is an improvement on the algorithm’s consistency and handling of 10 missing data. Example calculations with reanalyses data demonstrate pyPI’s usefulness in climatological and meteorological research. Planned future improvements will improve on pyPI’s assumptions, flexibility, and range of applications and tropical cyclone thermodynamic calculations.


15
Tropical cyclones pose significant risks to coastal societies, being among the costliest and deadliest of global natural hazards (e.g. Pielke et al., 2008;Rappaport, 2014;Hsiang and Jina, 2014). Damages increase exponentially with tropical cyclone intensity (∼5% per ms −1 Murnane and Elsner, 2012), so it is crucial to understand and accurately bound tropical cyclone maximum wind speeds. Theoretical and numerical models (e.g. Emanuel, 1987;Tsuboki et al., 2015;Sobel et al., 2016;Wehner et al., 2018) along with recent observations (Kossin et al., 2020) indicate that climate change has already increased storm 20 intensities-a trend expected to continue as the earth system warms. Emanuel (2005) showed the total destructive potential of tropical cyclones (derived from time-integrated maximum intensity) has increased since the 1970s.  Figure 1. The cross-section (along radius, R, and pressure, p) of an idealized and mature tropical cyclone and its thermodynamic cycle. pyPI inputs and outputs are in blue and red text, respectively, and are defined in Table 1. Bold blue lines and black letters indicate the four branches of the Carnot cycle, A through D (see text). The tropical cyclone maximum potential intensity (Vmax) is found at the radius of maximum winds (RMW) and is directed into-the-page in the northern hemisphere. The minimum central pressure (pmin) is found in the storm's eye.
Based on Carnot cycle illustrations in Emanuel (1987Emanuel ( , 2006. of the engine ( (Ts−T0) T0 ), and in terms of the heat source itself (i.e. thermodynamic disequilibrium). These quantities may be 90 estimated with atmospheric and oceanic observations (Sect. 5.1).
As derived in Bister and Emanuel (1998) and Emanuel (2003), the maximum (near-surface) potential intensity of a tropical cyclone, V max , may be approximated by: where C k and C D are the enthalpy and momentum surface exchange coefficients, respectively, h * o is the saturation moist static 95 energy at the sea surface, and h * is the saturation moist static energy of the air above the boundary layer (often evaluated at ∼500-600 hPa, cf. Wing et al. 2015). Tropical cyclone thermodynamic disequilibrium and efficiency are represented by the terms (h * o −h * ) and (Ts−T0) T0 , respectively; the ratio C k C D is a defined constant that may be estimated from theory or observations (Sect. 4.1).
In physically-based axisymmetric models, mature tropical cyclone wind speeds tend to reach their PI (e.g. Rotunno and 100 Emanuel, 1987). In contrast, observed storms rarely attain their thermodynamically-constrained potential intensities (see lim-itations discussed in Sect. 7). However, Emanuel (2000) combined climatologically-derived PI with observed tracks and intensities of real-world storms to show that any observed storm statistically has an equal likelihood of attaining any lifetime maximum speed between some lower bound (symbolized as "α" in Gilford et al., 2019) and the PI along its track. This is a powerful statistical property, because it implies that any shift in the PI distribution-either on short timescales such as during 105 an anomalously cold summer or on long timescales such as a response to a warming climate-will be accompanied by similar shifts in the observed intensity distribution (cf. Wing et al., 2007). Gilford et al. (2019) showed that randomly sampled observed intensity distributions which have at least 25 tropical cyclones (of hurricane strength or greater) will robustly follow their associated along-track potential intensity distributions.
Links between observed and potential intensities have been shown on seasonal (Tonkin et al., 2000;Gilford et al., 2019), 110 interannual (Wing et al., 2007;Shields et al., 2019), and climatological (e.g. Emanuel, 2000) timescales. The relationship is more robust when PI is evaluated along the track of a storm rather than as a basin-wide average (Wing et al., 2007;Gilford et al., 2019;Shields et al., 2019). Other studies have examined the roles of volcanic eruptions/lower stratospheric variability (e.g. Emanuel et al., 2013), the Montreal protocol (Polvani et al., 2016), or climate change on potential intensity (Emanuel, 2005;Vecchi et al., 2014;Sobel et al., 2016). Any oceanic or atmospheric variability or trend which alters the thermodynamic 115 environments of tropical cyclones could have some effect on PI (though the relative importance and/or statistical significance of these effects will vary). The connection between tropical cyclone PI and climate change will likely remain a critical topic to understand: the troposphere continues to warm and moisten in response to anthropogenic emissions of greenhouse gases, while simultaneously the lower stratosphere is cooling (e.g. IPCC, 2013;Santer et al., 2013).

BE02 PI Formulation
This section provides a detailed description and record (including relevant citations) of the PI algorithm with its thermodynamic/meteorological origins, assumptions, and computations. In addition to outlining the numerical computation here, pyPI's algorithm module has been commented for reference (Gilford, 2020).
Potential intensity may be derived following Bister and Emanuel (2002) (which is largely based on the formulation of 125 Emanuel 1995) idealizing a tropical cyclone as a Carnot heat engine (e.g. Emanuel, 1987), and assuming (1) the work done against friction by the outflow is ignored, (2) when the storm intensity reaches its maximum the anticyclone at the top of the storm is fully developed, and (3) the gradient wind may be approximated by cyclostrophic wind at the RMW. Under these conditions, the Carnot cycle formulation yields an expression for the maximum potential intensity (which roughly scales with the approximated PI expression, Eq. 1; Gilford et al. 2017;Wing et al. 2015): Algorithm 1 pyPI's iterative procedure to calculate tropical cyclone potential intensity.
input: Ts, p msl , p, T (p), r(p) where CAP E * is the convective available potential energy of saturated air lifted from sea level to the outflow level referencing the environmental profile, and CAP E env is the convective available potential energy of the environment. Because the final term is evaluated at the RMW and CAP E * is pressure dependant, an expression for the surface pressure at the RMW is needed.
Following Bister and Emanuel (2002) (cf. also Garner 2015, their Eq. 6), the minimum pressure of the tropical cyclone at the 135 RMW, p m is found with 1 : where T υ is the surface environmental virtual temperature, and CAP E| RM W is the environmental convective available potential energy evaluated at the RMW. Because the boundary layer water vapor mixing ratio is higher in the tropical cyclone eyewall than the storm's outer region (assuming a constant relative humidity in the boundary layer across the storm's radius),

140
CAP E| RM W is slightly larger than CAP E env (discussed more below).
The pressure dependence of CAP E requires solving Eq. 2-3 with numerical iteration, which pyPI performs with individual PI and CAP E modules. Algorithm 1 summarizes how pyPI computes maximum potential intensity by modeling a tropical cyclone as a Carnot heat engine. Algorithm inputs/outputs are provided in Table 1 and described in Sect. 4; meteorological constants are provided in the Appendix (Table A1). I begin by describing the CAP E calculation, which is used throughout the 145 pyPI algorithm.

CAP E Module
CAP E is defined as the sum of positive and negative areas of buoyancy energy of a lifted parcel on a sounding (e.g. Emanuel 1994, their Eq. 6.3.6, discussed more below) and is calculated by pyPI with the procedure in Algorithm 2.
1 Eq. 4 in Bister and Emanuel (2002) mistakenly replaces R d with cp. pyPI includes the correct factor of R d . Given an initial surface parcel temperature, pressure, and mixing ratio, the procedure begins by finding the parcel's reversible 150 entropy, which is conserved as it is lifted on the sounding. The parcel's water vapor pressure is found via the Ideal Gas Law (e.g. Bolton 1980, their Eq. 16): Saturation vapor pressure (in hPa) is given empirically as a function of T in • C by Bolton (1980), their Eq. 10, following the Clausius-Clapeyron relation: Then fractional relative humidity is defined as RH ≡ e es ≤ 1.0. Assuming the temperature dependence of specific heats is negligible over the range of temperatures in the tropical atmosphere, and integrating Kirchoff's equation (e.g. Emanuel 1994, Eq. 4.4.3-4.4.4) then the temperature dependence of the latent heat of vaporization is: with T in • C. Finally, we are equipped to calculate the parcel's reversible total specific entropy (per unit mass of dry air), s, which is conserved as the parcel is lifted along the sounding (Emanuel 1994, their Eq. 4.5.9): where r T is the total water content mixing ratio, which is identical to the parcel mixing ratio at the surface.
Having determined the parcel's initial moisture properties, we next find the lifting condensation level (LCL) of the parcel, 165 in order to partition the upcoming buoyancy calculation between saturated and unsaturated regions of the profile. The pressure of the LCL, p LCL , is found empirically with 2 : where A = 1669 and B = 122. Note that the LCL of a lifted parcel that is already saturated is identical to its original pressure level, i.e. p LCL ≡ p. Likewise, parcels at levels below the LCL are (by definition) not saturated. After finding p LCL the CAPE Starting with calculations at levels below the LCL (p j > p LCL ), at each j th level the algorithm calculates the unsaturated parcel temperature by following a dry adiabat with the same temperature as the surface parcel. Applying Poisson's equation:

175
Because CAP E is proportional to the positive and negative areas enclosed by the environmental and lifted parcel density temperatures (T ρ,e and T ρ , respectively), we calculate the density temperature as (Emanuel 1994, their Eq. 4.3.6): where the net water mixing ratio (r T ) is the same as the parcel water mixing ratio at the surface, and in the environment below the LCL before condensation has occurred (i.e. r T = r and r T,j = r j<LCL ).

180
Next the algorithm finds the density temperature differences for all levels above the LCL (p j < p LCL ). Because the parcel is saturated above the LCL, its moisture characteristics and temperature must be found iteratively at each j th level. First the algorithm solves for r by rearranging Eq. 4, then until the numerical iteration converges (with objective for the parcel temperature: ∆T < 0.001 K) it solves for parcel moisture characteristics which conserve the parcel's specific entropy s (Eq. 7) following a moist adiabat; finding these permits an estimation of the density temperature differences at each level. At the beginning of the loop, T and r are set equal to the previous iteration's findings, then the loop steps forward updating the parcel's temperature (and the dependant L v and water vapor mixing ratios) assuming s is conserved following saturated reversible adiabatic displacement. Following Newton's method (T n+1 = T n + s(Tn) ds(Tn)/dT ; Wallis 1685), when the difference between s and this iteration's entropy, s k , scaled by the rate of change of entropy with temperature, s , is small (i.e. s−s k s < AP * 0.001), then the algorithm will converge to estimate the parcel temperature at this level, T j . Here AP is a numerical step size employed 190 to speed convergence, which changes dynamically depending upon the number of iterations that have taken place. If at any given level the total number of iterations exceeds 500 (an excessive number of iterations), or if the water vapor pressure becomes unrealistically close to the level pressure, then the algorithm fails to converge and returns zero CAP E.
When the algorithm converges for a level, the final parcel mixing ratio is set depending on the ascent type chosen by the user (Sect. 4.1). For pseudoadiabatic ascent (ascent_flag= 1), liquid water condensed in the parcel during its ascent is 195 assumed to drop out of the parcel, such that the heat capacity of liquid water is neglected and the mixing ratio is a function of the final level temperature (i.e. r j = r(T j )). For reversible ascent (ascent_flag= 0) the total water (and its heat capacity) is retained following the parcel (i.e. r j = r j=0 = r T ). For intermediate fractions of ascent_flag, the mixing ratio scales Note that the density temperature difference (and hence a parcel's buoyancy) with height is not strictly higher under either 200 ascent assumption. Parcels lifted reversibly are always warmer than those lifted psuedoadiabatically, but the weight of the carried condensate also means these parcels are more dense until they reach the upper troposphere (Emanuel 1994, their Table   4.2). Accordingly, Gilford et al. (2017) found that psuedoadiabatic (typically more buoyant) PI calculations generally have higher altitude OTLs than reversible (typically less buoyant) PI calculations on monthly timescales.
Having determined r j , the algorithm computes the density temperature for the parcel and the environment (Eq. 10), and 205 calculates each level's density temperature differences, T ρ − T ρ,env . We are now equipped to calculate the lifted parcel's convective available potential energy. CAP E is given by the vertically integrated buoyant energy between the level from which the parcel is initially lifted (j = 0) and the level of neutral buoyancy (LNB; j = LN B). Following Emanuel (1994), their Eq. 6.3.6: areas (P A) are vertical regions of positive buoyancy which cause the parcel to rise assuming an initial upward displacement.
Note that CAP E is not defined for parcels without positive areas. By definition, the level of free convection (LFC) separates regions that are negatively buoyant (below) from regions that are positively buoyant (above). When (j = LF C) > (j = LCL), then regions of the profile above the LCL and below the LFC may still be negatively buoyant.
The CAP E algorithm numerically solves Eqs. 11-13 in five steps: Second, noting that dp/(p mean,j:j+1 ) ≈ dlog(p) j:j+1 at each layer over the levels j and j + 1-where the average pressure of each layer is p mean,j:j+1 = 1 2 (p j + p j+1 )-we find the positive and negative areas between the second-highest altitude level (i.e. j = 1) and the the maximum level of positive buoyancy.  Finally, the negative and positive areas are added together with the residuals following Eq. 11. After this last step, the algorithm flag is set to indicate the algorithm has successfully computed CAP E. Then the values of CAP E, T LN B , the LNB, and the flag are returned to the PI module.

PI Module
The PI module begins by checking to ensure that the input atmospheric profile is appropriate for the PI calculation. If not, 235 missing values are returned by the algorithm (see Sect. 4.2). Water vapor mixing ratios above the boundary layer do not influence the PI calculation (they are redundant in CAP E * − CAP E env ), so any missing r values above the surface are replaced with 0 g/g.
Following Algorithm 2 described above, pyPI computes CAP E env , assuming the environmental air parcel is lifted from the lowermost input level in p.

240
Next, pyPI iteratively solves Eqs. 2 and 3 (with objective for the minimum pressure at the RMW: ∆ i,i−1 p m < 0.5 hPa). The algorithm begins by calculating the convective available potential energy (computing Algorithm 2 at each i th iteration) iterating from the initial lowest-level environment inward toward the radius of maximum winds, CAP E| RM W . At each iteration the mixing ratio is updated to account for pressure dependence (r increases slightly as p → p m approaching the RMW; Bister and Emanuel 2002). 245 Next, we calculate the saturation convective available potential energy at the radius of maximum winds, CAP E * . This calculation assumes the parcel is lifted directly from the sea surface, such that T = T s and r = r s (T s ) (where r s is found given T s via Eqs. 5 and 4). pyPI defines the outflow temperature level (OTL) and T 0 as the LNB and T LN B found during the final iteration of CAP E * computation, respectively Note that the OTL and T 0 and could instead be defined from the CAP E env computation. Flexibility in the outflow definition is a planned improvement for pyPI. The choice to use CAP E * follows from 250 defining the outflow level with a fully saturated parcel lifted directly from the sea surface (see Sect. 6.2).
The ratio of sea-surface and outflow temperatures in Eq. 2 represents the scaling of PI by dissipative heating, which increases PI when Ts T0 > 1 (Bister and Emanuel, 1998). At each iteration this ratio is set with the fixed input T s and the current T 0 ≈ T LN B . The relevance of this ratio for the PI calculation is set by the user with the adjustable parameter diss_flag (Sect. 4.1).
If dissipative heating is permitted to impact the tropical cyclone potential intensity (diss_flag=1) then the ratio remains as 255 defined above. If dissipative heating is not considered (diss_flag=0) then the algorithm assumes : 1 T s /T 0 in the following calculations.
Next, p m is estimated in each iteration following Eq. 3. The surface environmental virtual temperature, T υ , is found as the average of virtual temperatures over the mean layer composed of the parcel (with temperature, T s ) and the lowest level, i.e.
T υ = 1 2 (T υ,s + T υ,j=0 ). The virtual temperature is identical to the density temperature (Eq. 10) at the surface (as r T = r(T s )), 260 and may be approximated at the lowest level with r T ≈ r j=0 . Combining Eqs. 2-3 to solve for p m , the algorithm iterates towards a new pressure estimate. If the number of iterations exceeds 200 (an excessive number of iterations), or if the estimated pressure drops below an unphysical 400 hPa, then the PI algorithm fails to converge and returns missing outputs.
When the algorithm has successfully converged on a stable value of p m , then the final central minimum pressure, p min , is set. Assuming cyclostrophic balance and that the azimuthal velocity in the eye is given by power law scaling with exponent, b (see also Emanuel 1995, their Eqs. 25-26): where pyPI assumes following Bister and Emanuel (2002)   The ratio C k C D an uncertain constant which depends on the sea state and linearly scales potential intensity; its value is an ongoing area of field and theoretical research (Emanuel, 2003). Table 1 includes the 1σ range of the ratio found with energy

ascent_flag (default=0)
The ascent process proportion determines whether the air parcels displaced in each CAP E calculation (cf. Sect. 3.2) follow re-295 versible adiabatic ascent (ascent_flag=0) or pseudoadiabatic ascent (ascent_flag=1). In the case of reversible ascent, the full moist entropy of the buoyant parcel is conserved along its displacement following a moist adiabat. In pseudoadiabatic ascent the heat capacity of liquid water is neglected. Liquid water is assumed to fall out of the parcel as it condenses, while the parcel ascends following the pseudoadiabatic moist adiabat; for more details see Emanuel (1994), their Sect. 4.7. For practical applications of pyPI, ascent_flag may be set to any value between 0.0 and 1.0, such that the proportion of ascent is any 300 fraction intermediate to fully reversible and fully pseudoadiabatic ascent.

diss_flag (default=1)
The dissipative heating flag determines whether dissipative heating is accounted for (diss_flag=1) or ignored (diss_flag=0) in potential intensity theory (see Bister and Emanuel 1998, their Eq. 22). When dissipative heating is included in the PI calculation, the leading factor in the BE02 algorithm (Eq. 2) is ( C k C D * Ts T0 ), where Ts T0 > 1. In the absence of dissipative heating, 305 the leading factor is ( C k C D * 1) following the original findings of Emanuel (1986Emanuel ( , 1995. PI is often considerably lower when dissipative heating is neglected.

V_reduc (default=0.8)
Raw potential intensities are maximum gradient wind speeds (Emanuel, 2000). Therefore, gradient winds calculated with the BE02 algorithm are not directly comparable with observed intensities at the near-surface without applying an approximate 310 scaling between gradient and 10 meter winds. Following Powell (1980), a crude reduction of 20% (V_reduc=0.8) is typically applied to scale PI for comparison with near-surface winds. The percent reduction in the gradient wind in terms of V_reduc is P reduc = 100% * (1 − V_reduc). Note that for some applications of PI, such as using it as thermodynamic parameter in climate science-e.g. incorporation into the Genesis Potential Index, Camargo et al. (2007)-V_reduc should be set to 1.0 (no reduction).

ptop (default=50)
The upper level pressure bound is the minimum pressure below which the input profile is ignored during PI computation.
Theoretically modeled tropical cyclone outflow can often exceed the tropical tropopause on climatological timescales (Gilford et al., 2017), so setting ptop> 100 hPa it is not advisable. Reducing the number of considered levels by increasing ptop may potentially increase the speed of calculations, at the risk of finding a too-low OTL and too-warm T 0 . Before altering ptop, 320 users should consider their particular application and the outflow levels they anticipate given the stability of their input profiles.

miss_handle (default=1)
The missing data flag prescribes how missing values are handled in the CAP E calculation (discussed below). Following the BE02 MATLAB code (miss_handle=0), if missing values are found in the input temperature profile, then the algorithm will attempt to calculate PI for all available levels above the missing values. However, the user may also conservatively choose that 325 any missing values in the input profile will immediately set the entire PI calculation output to missing (miss_handle=1).

Handling Missing Data
Mirroring the output flag convention of the BE02 MATLAB code: IFL=1 when the PI algorithm successfully returns valid potential intensity ouputs, IFL=0 when the algorithm fails because the input data is improper for a PI calculation (e.g. if T s < 5 • C), and IFL=2 when the algorithm fails to converge.

330
One major difference between pyPI and the BE02 MATLAB algorithm is the handling of missing data and the (related) flag provided in the output. By convention, missing input variables in pyPI are assigned Python's "Not-a-Number", NaN, to avoid errors. The BE02 MATLAB code default is that profiles may contain missing values (specifically temperatures on pressure levels), and the algorithm computes PI over the remaining valid levels.
Because missing values may sometimes be found at the surface-and the primary CAP E calculation (Sect. 3.2) relies 335 heavily on the assumption of lifting the parcel within the storm and environment from that level-errors could arise from estimating PI when ignoring near-surface buoyancy. In principle, PI should be calculated only over data points with existent sea-surface temperatures and lowest profile level temperatures. In practice, missing data may arise at the lowest profile level, which would lead to errant PI calculations if these profiles are input to the BE02 MATLAB code. pyPI addresses this challenge in three ways. First, an adjustable parameter (miss_handle) is implemented to allow the 340 user to specify how pyPI handles the missing values. If miss_handle=0, the code attempts to handle missing values akin to the way that the BE02 MATLAB code did, although there still remain some differences in the outputs between pyPI and the MATLAB algorithm. Specifically, pyPI's CAP E calculation proceeds as normal only as long as there are no missing values between the lowest valid (non-missing) level and the OTL; otherwise, CAP E module outputs (and hence PI module outputs) are returned as missing. Second, if miss_handle=1, then the CAP E function will automatically interpret temperature 345 profiles with missing data as invalid, and return missing values to the PI algorithm, resulting in the PI outputs being set to missing in the return. Third, a new output flag value (IFL=3) is introduced in pyPI which is returned when missing values in the temperature profile results in a missing output return from the PI module (i.e. in either of the two cases described above), which aids in data analyses. Figure 2 shows an example of the output algorithm status flags from pyPI calculations with a month of MERRA2 data (Sect.  In the example pyPI calculation (Fig. 2) there are no inputs for which the algorithm does not converge (cf. Bister and Emanuel 2002). One final complication is that BE02 MATLAB code occasionally returns V max = 0 as an output. In these cases, pyPI instead returns V max = NaN. pyPI outputs valid (non-missing) potential intensities over ∼65.6% of the ocean grid points in the global 2004 sample 360 dataset-compared with ∼64.41% returned by the BE02 MATLAB code. In addition to how missing data is handled, output differences may also arise from slight variations in numerical computation between Python and MATLAB. pyPI validation tests (section 5) are computed over all spatio-temporal locations for which both algorithms have non-missing/non-zero potential intensities (and with miss_handle=1).

Sample Reanalysis Data
The pyPI sample data are monthly means of state variables from the second Modern-Era Retrospective Analysis for Research and Applications (MERRA2, Gelaro et al. 2017) in 2004, interpolated onto a 2.5 • x 2.5 • global grid (Gilford et al., 2017).
Note that for these example pyPI calculations the water vapor mixing ratio, r, is approximated by substituting in the reanalysis specific humidity, q (as q ≡ r 1+r ≈ r, because r 1).
Potential intensity calculations are generally linear, i.e., mean potential intensities may be estimated as a function of mean environmental variables, where E[ · ] is the expected value of a function or variable. Using monthly mean environmental conditions to compute climatological monthly means of potential intensity (and the algorithm's other output variables) generates a small bias of <1 375 ms −1 globally and <0.5 ms −1 in the tropics . PI's linearity property is convenient, because it reduces the scale of data needed to compute PI: daily or hourly data are not needed for monthly or longer (climatological) applications.
Applications on shorter (e.g. operational or daily) timescales, however, should use appropriately shorter frequency inputs to the pyPI algorithm.
For consistency with Gilford et al. (2017), the sample data uses the land-sea mask from the European Centre for Medium- PI may be mistakenly calculated over land with the PI module. In these cases, users should assign all PI algorithm outputs over land to the missing value in post-processing. As an alternative, in this pyPI example input variables over land are set to missing in a pre-processing step. Note that the mismatch between using the ERA-I land-sea mask and MERRA2 data in this example 385 results in a set of minor output artifacts caused by missing (MERRA2 land grid points) input data over ERA-I defined ocean grid points. This artifact provides a useful demonstration of the missing data flag employed in pyPI (Sect. 4.2).

Validating Against the BE02 Implementation
Accompanying the environmental conditions in the sample data are outputs from the BE02 MATLAB code written by Kerry Emanuel (see ftp://texmex.mit.edu/pub/emanuel/TCMAX) and revised by Daniel Gilford for climatological research applica-390 tions (see Gilford et al. 2017Gilford et al. , 2019. Potential intensities calculated over September 2004 with pyPI and the extensively-used BE02 MATLAB code are compared in Figure 3 over the globe; their difference is computed and plotted in Figure 4. There is excellent agreement between the two algorithms; 98.5% of output potential intensity have absolute differences < 0.01 ms −1 . Potential intensities calculated with the Python algorithm exhibit a slightly negative bias relative to the MATLAB calculations, but these differences are negligible  Figure 5 shows the scatter between all potential intensity values calculated with the two algorithms over the sample data, plotted against a 1-to-1 line (values lying on this line exhibit perfect agreement). The R-squared of this comparison is R 2 ≈ 1.0 to seven significant digits, such that the calculations are nearly identical. All other output variables (cf. Table 1) from the two algorithms have similarly strong levels of agreement.

400
A minor PI difference between pyPI and the BE02 MATLAB algorithm arises when the pyPI-found outflow level has a higher pressure and warmer temperature than found by the BE02 code. The outflow property differences result from pyPI's correction of a minor error that was present in the BE02 algorithm, where was defined as 0.622 rather than directly calculated as R d Rv . The small rounding error results in a handful of profiles with lower altitude outflow and lower PI values calculated by pyPI. In the absence of correcting this error, the correlation between the two calculations is R 2 ≈ 1.0 to thirteen significant 405 digits, and the absolute maximum difference anywhere is 4.1 × 10 −5 ms −1 .
I conclude that the PI calculations made with the pyPI algorithm are adequately validated against the BE02 MATLAB code, and that pyPI is sufficiently accurate for use in research applications.

PI Seasonal Cycles
A slightly more sophisticated application of pyPI is the calculation of potential intensity seasonal cycles. Reproducing the methodology of Gilford et al. (2017) with pyPI calculations over 2004, Figure 7 shows the seasonal cycles of sea-surface 420 temperatures, and outflow temperatures, outflow temperature levels, and potential intensities in 2004 averaged over tropical cyclone main development regions (defined in Gilford et al. 2017, their Table 1).  Figure 8. Skew-T log-P thermodynamic diagram with isotherms (thin black curves), dry adiabats (green curves), and moist adiabats (blue curves). The bold black line is a mean environmental temperature profile from the North Atlantic region (Gilford et al. (2017), their Table 1), the magenta curve is the moist adiabat associated with a mean North Atlantic sea-surface temperature (SST), and the red curve is the moist adiabat associated with a sea-surface temperature 3 • C warmer than the mean.

J F M A M J J A S O N D Month
The seasonal cycles of each basin illustrate the complex relationship between sea-surface temperatures, OTLs, and outflow temperatures. Figure 8 diagrams this relationship in more detail, showing how one assumption in the pyPI algorithm impacts the output PI values. pyPI assumes that the outflow temperature and its level are derived by finding the LNB assuming a saturated parcel lifted from the sea-surface with temperature, T s . This implies that, following a moist adiabat, the level of the 435 neutral buoyancy is a function of only T s and the environmental temperature profile, T . Given a fixed temperature profile, then an increase in T s (e.g. from SST1 → SST2 in Fig. 8) requires that the associated OTL will be found at a higher altitude (OTL1 → OTL2 in Fig. 8) and the associated outflow temperature will likewise change. As the atmosphere's stratification increases into the lower stratosphere, increases in T s become less effective at changing the OTL and its temperature, with the effect nearly saturating when the outflow reaches the cold-point tropopause (e.g. 100 hPa in Fig. 8). At this point, T 0 variability is almost 440 completely decoupled from T s variability. Instead these T 0 values become influenced by tropopause region variability (e.g. Emanuel et al. 2013;Ramsay 2013;Wang et al. 2014;Wing et al. 2015;Gilford et al. 2017) which is controlled by radiation, dynamics, and deep convection (e.g. Fueglistaler et al. 2009;Randel and Wu 2014).
These properties are borne out in the example 2004 seasonal cycles computed in Fig. 7. In the North Atlantic basin seasurface temperature and OTL seasonal cycles are inversely proportional: colder sea-surface temperatures have higher pres-445 sure OTLs and warmer outflow temperatures found in the upper troposphere (where dT dz < 0) in all months except August-September. In these late summer months the OTL reaches near the cold-point tropopause, and the outflow temperature slightly increases, following the seasonal cycle of warmer tropopause temperatures (e.g. Yulaeva et al. 1994). A contrasting pattern is observed in the Western North Pacific, where OTLs have almost no seasonal cycle: in this basin the calculated outflow always reaches the lowermost stratosphere (OTL ≤ 90 hPa). Accordingly, the outflow temperature seasonal cycle perennially follows 450 the seasonality of lowermost stratospheric temperatures, which minimize in the boreal winter and maximize in the boreal summer (Yulaeva et al., 1994). Comparing with Eq. 1, this T 0 seasonality leads to relatively increased PI values in the boreal winter and relatively decreased PI values in the boreal summer. Overall, the PI seasonal cycle in the Western North Pacific is damped over the year, a pattern that is observed in real-world tropical cyclone intensities (cf. Gilford et al. 2019).

455
The relative contributions to potential intensity may be mathematically derived by decomposing Eq. 1. Taking the natural logarithm of both sides: then PI variability is related to variability in either tropical cyclone efficiency ( Ts−T0 T0 ) or thermodynamic disequilibrium (h * o − h * ); recall that C k C D is taken as a constant. As an example following Gilford et al. (2017), Eq. 16 is applied to pyPI calculated 460 2004 seasonal cycles of potential intensity (from Fig. 7). pyPI calculates PI directly, and efficiency may be directly computed from input T s and output T 0 ; following Wing et al. (2015) the disequilibrium term is taken as a residual from Eq. 16.
After finding each term in Eq. 16 over each basin and seasonal cycle, the amplitude (defined as the annual range evaluated with monthly observations) of each seasonal cycle is plotted in Figure 9. By convention (Gilford et al., 2017, a negative amplitude indicates the ∼sinusoidal seasonal cycle reached its maximum in the boreal winter and minimum in the boreal 465 summer. In all basins, the disequilibrium term drives the largest portion of the seasonal amplitude. This is consistent with sea-surface temperature seasonal cycles which dominate the disequilibrium variance (Fig. 7). The efficiency term is smaller, and in each basin it follows the same cycle as thermodynamic disequilibrium, with the exception of the Western North Pacific (where the efficiency seasonal cycle maximizes in the boreal winter and minimizes in the boreal summer). This opposite signed seasonality 2019). Notably, the Southern Hemisphere shares this outflow temperature seasonality, which actually amplifies the efficiency seasonal cycle through both sea-surface temperatures and outflow temperature intraseasonal variability. In all other basins, outflow temperature seasonality is offset by the sea-surface temperature seasonality, which acts to mute the efficiency term and 475 further contribute to disequilibrium dominating their seasonal cycles. The decomposition in Fig. 9 illustrates how the roles of environmental conditions in PI seasonality are basin dependant.
These simple examples show how pyPI may be used to study tropical cyclone intensities, and likewise demonstrate pyPI's ability to produce findings similar to those previously computed with the BE02 MATLAB code. pyPI is a Python package that models the maximum potential intensity (PI) of a tropical cyclone given its environmental conditions. pyPI(v1.3) is the first fully documented PI algorithm, advancing on the previous MATLAB code (Bister and Emanuel, 2002) which has been extensively used but under-documented in the literature. In addition to documenting PI computation, allowing dynamic parameter selection, and correcting minor errors in the previous algorithm, pyPI is also an open source, maintained, and archived project which permits reproducibility, continual updates and improvements, and accountability for 485 future PI calculations in the tropical meteorology community (Gilford, 2020).
pyPI calculations exactly reproduce outputs from the Bister and Emanuel (2002) algorithm, except in rare cases where the original algorithm's implementation was errant. Sample analyses show the flexibility and usefulness of PI calculations for understanding variability and thermodynamic contributions to climatological tropical cyclone maximum intensities.
Because of its statistical ties with observed storms (Emanuel, 2000;Wing et al., 2007;Gilford et al., 2019;Shields et al., 490 2019) PI is powerful tool for exploring past and future changes in real-world maximum intensities. pyPI computations have a broad range of possible applications, which could include operational meteorology (e.g. the PI maps produced by the Center for Land-Atmosphere Prediction, http://wxmaps.org/pix/hurpot, Emanuel et al. 2004) and climate change research (Sobel et al., 2016).
Potential intensity is a theoretical model with several notable limitations. Real-world tropical cyclones rarely are in quasi-495 steady state or meet the idealized conditions required for the Carnot cycle model. This makes PI less suitable for operational purposes, though it may still be incorporated into real-time genesis or intensification indices (see below). Furthermore, PI theory does not directly account for complicating factors such as vertical wind shear or large-scale subsidence, which are known to have important influences on tropical cyclone intensity. The ratio of surface exchange coefficients, C k C D , is also highly uncertain but important for PI magnitude. Previous studies have adapted PI to make it more suitable for various applications 500 (e.g Lin et al., 2013;Kieu and Zhang, 2018); pyPI users should carefully consider PI assumptions and applicability in their research problems (Gilford, 2018), adapting pyPI or suggesting package enhancements as appropriate.
Future planned improvements of pyPI include an expansion of the codebase to compute other tropical cyclone thermodynamic and statistical indices, including the Genesis Potential Index Zhang et al., 2016) and Ventilation Index (Tang and Emanuel, 2012). A direct disequilibrium calculation (e.g. Wing et al., 2015) module would permit comparisons 505 with the residual approach currently employed in pyPI (Sect. 6.3). Finally, further improvements in the algorithm's handling of missing data are warranted to reduce the algorithm run time and improve pyPI's applicability for a wider range of input profiles.

Appendix A: pyPI Constants
Constants used in to model potential intensity in the BE02 algorithm have been directly used in pyPI and are recorded in Table   A1.
Author contributions. DG conceived the study, wrote and tested pyPI software and its v1.3 modules, notebooks, and scripts in Python, and pyPI, and for his development of PI theory and the original PI algorithm. DG is supported by NSF Grant ICER-1663807 and NASA Grant 80NSSC17K0698. Thanks to Daniel Rothenberg for implementing a Numba optimization (Lam et al., 2015) of the primary pyPI module.
Thanks to Allison Wing, Dan Chavas, Jonathan Lin, and Raphael Rousseau-Rizzi for helpful comments and suggestions on pyPI and its documentation. This study is validated against a portion of previous research on PI seasonality, published in Gilford et al. (2017)