## Abstract

This paper introduces the Mean Airflow as Lagrangian Dynamics Approximation (MAFALDA), a new method designed to extract thermodynamic cycles from numerical simulations of turbulent atmospheric flows. This approach relies on two key steps. First, mean trajectories are obtained by computing the mean circulation using height and equivalent potential temperature as coordinates. Second, thermodynamic properties along these trajectories are approximated by using their conditionally averaged values at the same height and *θ*_{e}. This yields a complete description of the properties of air parcels that undergo a set of idealized thermodynamic cycles.

MAFALDA is applied to analyze the behavior of an atmosphere in radiative–convective equilibrium. The convective overturning is decomposed into 20 thermodynamic cycles, each accounting for 5% of the total mass transport. The work done by each cycle can be expressed as the difference between the maximum work that would have been done by an equivalent Carnot cycle and a penalty that arises from the injection and removal of water at different values of its Gibbs free energy. The analysis indicates that the Gibbs penalty reduces the work done by all thermodynamic cycles by about 55%. The cycles are also compared with those obtained for doubling the atmospheric carbon dioxide, which in the model used here leads to an increase in surface temperature of about 3.4 K. It is shown that warming greatly increases both the energy transport and work done per unit mass of air circulated. As a result, the ratio of the kinetic energy generation to the convective mass flux increases by about 20% in the simulations.

## 1. Introduction

The atmospheric circulation is characterized by a complex interplay between dynamics and thermodynamics. As air parcels move around the atmosphere, they go through a wide variety of transformations, such as expansion, compression, phase transition, radiative cooling, and mixing. These thermodynamics transformations in return affect air density and henceforth the evolution of the flow. In particular, the atmosphere as a whole acts as a heat engine that produces kinetic energy by transporting energy from a warm source to a cold sink, and, consequently, thermodynamic processes determine the amount of kinetic energy available to maintain atmospheric motions.

However, not all heat engines are equivalent: the mechanical output of a heat engine depends not only on the amount of energy transported and on temperature difference between the energy sources and sinks, but also on the nature of the thermodynamics transformations involved. In atmospheric science, the concept of heat engine is too often equated with that of a Carnot cycle (e.g., Emanuel 1986; Emanuel and Bister 1996; Bohren and Albrecht 1998; Ambaum 2010). Unfortunately, the Carnot cycle is a poor model for how the atmosphere generates kinetic energy. The primary reason for this lies in the role played by Earth’s hydrological cycle. A large fraction of the radiative energy received by Earth is used to evaporate water vapor at the surface of the ocean, which is then transported upward and removed through condensation and precipitation.

This active hydrological cycle has a strong impact on the atmosphere’s ability to generate kinetic energy. First, falling hydrometeors are associated with microscale Stoke flows that slow their fall and dissipate kinetic energy (Pauluis et al. 2000). The amount of dissipation in falling precipitation is on the same order of magnitude as the dissipation of kinetic energy by the large-scale circulation (Pauluis and Dias 2013). Second, the upward transport of latent heat in the atmosphere is not a very efficient heat engine. Pauluis (2011) shows that the mechanical work produced by a steam cycle—an idealized heat engine similar to a Carnot cycle except that the energy source is due to evaporation—is always less than the work produced in a Carnot cycle. The difference between the two cycles can be directly explained by the fact that water vapor evaporates in unsaturated air parcels but is removed at saturation and can be quantified in terms of the Gibbs free energy of the water vapor when it is added and removed from the atmosphere.

These two effects—precipitation-induced dissipation and the Gibbs penalty due to the injection of water vapor at low relative humidity—greatly reduce the amount of work that can be produced by the atmosphere. Pauluis and Held (2002a,b) analyze the entropy budget in simulations of radiative–convective equilibrium and show that the generation of kinetic energy by resolved atmospheric motions is one order of magnitude smaller than the maximum work that would have been produced by a Carnot cycle for the same external heating. They argue that moistening of dry air must occur through diffusion of water vapor and irreversible phase transition. These irreversible processes are associated with an internal entropy production and reduce the overall efficiency of the atmospheric heat engine. In their numerical simulations, Pauluis and Held (2002a) found that this irreversible moistening of dry air accounted for two-thirds of the total entropy production by deep convection so that the generation of kinetic energy by moist convection was only a fraction of what would be predicted if one assumes that the atmosphere behaves as a Carnot cycle. The overwhelming implication is that the amount of work generated by atmospheric motion depends crucially on the behavior of the hydrological cycle, a result that has been recently confirmed by Laliberte et al. (2015) in global climate models.

The main purpose of this paper is to develop a systematic framework to analyze the thermodynamic behavior of atmospheric flows. Such flows are highly turbulent and involve multiple scales of motion and a wide range of thermodynamic processes. Furthermore, performing a thorough thermodynamic budget, such as that of Pauluis and Held (2002a,b) requires detailed diagnostics on diffusion, dissipation, and microphysics, which are rarely available as standard model output. To address these issues, the thermodynamic analysis developed in Pauluis (2011) will be adapted to study the thermodynamics behavior of moist convection in high-resolution numerical simulations. The central challenge here is to apply the concept of thermodynamic cycle, in which a single air parcel follows a well-defined set of thermodynamic transformations, to a turbulent flow in which each individual air parcel exhibits a distinct trajectory and evolution.

To do so, a new technique, the Mean Airflow as Lagrangian Dynamics Approximation (MAFALDA) is introduced in section 2. First, the mean parcel trajectories in height and equivalent potential temperature are computed using the isentropic analysis method developed by Pauluis and Mrowiec (2013). Second, the various properties of air parcels along these trajectories are obtained from their conditionally averaged values on isentropic surfaces. In section 3, we analyze the MAFALDA cycles obtained in an idealized radiative–convective equilibrium simulation. It is shown that the hydrological cycle greatly reduces the amount of kinetic energy that can be generated by atmospheric motions. Section 4 discusses the fact that deep convection is much more efficient at generating kinetic energy than shallow overturning, which is explained by scaling arguments for the upward energy and water transports. Section 5 investigates the impacts of a doubling of carbon dioxide concentration on the thermodynamic behavior of moist convection. Section 6 discusses the approximations underlying MAFALDA and evaluates the accuracy of the reconstructed cycles. The conclusions are presented in the last section.

## 2. The Mean Airflow as Lagrangian Dynamics Approximation

This paper analyzes a set of simulations of convection using the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall 2003). This widely used cloud-resolving model combines an anelastic dynamical core with a variety of physical parameterizations for cloud microphysics, turbulence, and radiation. The configuration used in our simulations includes a single-moment microphysics in which the condensed water is split between cloud water, cloud ice, rain, snow and graupel, a Smagorinsky closure for turbulent diffusion, and an interactive radiative transfer. Convection is simulated on a doubly periodic domain with a horizontal resolution of 1 km. The vertical resolution is staggered from 75 m near the surface to 500 m in the upper troposphere. For the lower boundary, we use a 2-m-thick slab ocean of uniform temperature. The model is used to simulate radiative convective equilibrium for a 380-ppm concentration of carbon dioxide. The radiation algorithm assumes a solar constant of 685 W m^{−2} and a zenith angle of 51.7°, with no annual or diurnal cycle. This configuration yields the same average incoming solar radiation as that at the equator. A uniform energy sink of 100 W m^{−2} is imposed in the ocean slab, which mimics the energy export out of the equatorial regions by the global oceanic and atmospheric circulations, albeit its primary purpose here is to obtain a reasonable equilibrium temperature. Because of the interactive ocean model, the convergence toward radiative convective equilibrium is rather slow (hundreds of days), and the last 60 days of the simulations are used for the analysis presented in this section. The average ocean temperature in this reference experiment is 300.5 K. With the exception of the interactive slab ocean, the model is identical to the one used in Pauluis and Mrowiec (2013).

MAFALDA is a systematic procedure meant to extract the thermodynamic cycles from the standard output of a numerical simulations. In effect, it requires a sequence of three-dimensional snapshots for vertical velocity *w*; temperature *T*; mixing ratio for water vapor , liquid water , and ice ; and as total pressure. It proceeds in the five following steps:

Compute the isentropic streamfunction in height and equivalent potential temperature coordinates to determine the mean parcel trajectories (section 2a).

Compute the conditional average values of various state variables as function of

*z*and (section 2b).Define a set of mean trajectories in

*z*–*θ*_{e}space based on the streamfunction (section 2c).Evaluate the value of various thermodynamic variables along these trajectories (section 2d).

Analyze the Gibbs relationship along these trajectories (section 3).

A schematic representation of this process is shown in Fig. 1. The first four steps aim at extracting a set of thermodynamic cycles that capture the various transformations implicit in the simulations and are discussed in sections 2a–d. The last step analyzes the behavior of these cycles through the lens of the Gibbs relationship and is detailed in section 3. The model here is treated as a black box in that MAFALDA only requires a set of three-dimensional outputs of standard model variables (temperature, water content, and vertical velocity).

### a. Isentropic streamfunction

The methodology introduced by Pauluis and Mrowiec (2013) is used to analyze a convective overturning in isentropic coordinates , with the equivalent potential temperature with respect to ice. At the core of the isentropic analysis lies the concept of an integral on an isentropic surface:

Here, *T* is the period used for time averaging, and are the horizontal extent of the domain, respectively, and *δ* is the Dirac function. In practice, the integral in (1) is computed by sorting the function *f* on finite bins, which amounts to approximating the Dirac function as a constant function of value on an interval of width . The isentropic integral is expressed in units of *f* per kelvin.

In this study, isentropes are defined as surfaces of constant , which is given by

Here, , , , and are the mixing ratios for water vapor, liquid water, ice, and total water, respectively; , , , and are the specific heat capacities at constant pressure of dry air, ice, liquid water, and water vapor, respectively; and are the specific gas constants for dry air and water vapor, respectively; and are the latent heat of vaporization and freezing, respectively; is the partial pressure; *T* is the temperature; is the relative humidity with respect to liquid water; is the reference pressure; and is the freezing temperature for water under atmospheric pressure. From a physical point of view, corresponds to the temperature that an air parcel would have after first being lifted so that all the water content is in the ice phase, then being compressed back to the reference pressure without allowing any phase transition (including melting). It is larger than the equivalent potential temperature over liquid water (Emanuel 1994) because of the inclusion of the latent heat of freezing. Its use here is motivated by the desire to better capture deep convective motions above the freezing level. For simplicity of notation, the subscript *i* will be dropped for the rest of the paper, but it should be understood that the equivalent potential temperature here refers to its value over ice.

The isentropic distribution of vertical mass flux is computed as the isentropic integral of the vertical mass transport following (1). An isentropic streamfunction in the *z*–*θ*_{e} space is then defined as

with *H* the Heaviside function. The isentropic streamfunction is equal to the net vertical mass flux at height *z* of all air parcels with an equivalent potential temperature less than . Figure 2a shows the isentropic distribution of vertical mass , and Fig. 2b shows the corresponding streamfunction . The isentropic flow here is very similar to the one discussed in Pauluis and Mrowiec (2013). It is characterized by the ascent of warm moist air at a high value of the equivalent potential temperature and subsidence of air with lower . The mass transport peaks in the lower troposphere, a signature of the vigorous mixing associated with shallow convection. The equivalent potential temperature of the rising air parcels decreases with height, which indicates entrainment of environmental air into the updrafts. Downward motions occur at values of that are very close to the horizontal mean (solid black line in Fig. 2), which confirms that most of the downward mass flux occurs through slow subsidence in the environment. The interested reader should refer to Pauluis and Mrowiec (2013) for a more detailed discussion of the physical interpretation of the isentropic streamfunction and its use for analyzing convective overturning.

### b. Conditional averaging

Finally, for any variable *f*, we define its conditional mean for all air parcels at a given height and equivalent potential temperature as

As indicated by the definition (4), these are functions of both height *z* and equivalent potential temperature and are defined as the mass-weighted average of the quantity *f* for all air parcels at a given height and equivalent potential temperature.^{1 }Figure 3 shows the conditional average of various fields that will be used in the thermodynamic diagnostics in section 3. Temperature is shown in Fig. 3a. It decreases systematically with height and is close to a moist adiabatic profile through the troposphere. At a fixed height, temperature increases with equivalent potential temperature. Moist entropy (Fig. 3b) here is defined as the moist entropy per unit mass of dry air and is defined in the appendix. The close relationship between moist entropy and potential temperature is quite apparent, with the entropy increasing systematically with , albeit it is not a one-to-one relationship. The mixing ratio, shown in Fig. 3c, exhibits also a systematic decrease with height that is consistent with a Clausius–Clapeyron relationship. At a given height, the mixing ratio increases with as well. The specific Gibbs free energy for water vapor is defined in the appendix. It can be approximated by . The conditional mean distribution , shown in Fig. 3d, reflects primarily the variations in relative humidity: is close to 0 at low altitude and high , which corresponds to cloudy air. Above the freezing level, cloudy air is saturated with respect to ice but is unsaturated with respect to liquid water, corresponding to a negative value for the Gibbs free energy. Finally, unsaturated air is characterized by large negative values for . The buoyancy shown in Fig. 3e is small for unsaturated air but can reach values as high as 0.1 m s^{−2} for saturated air parcels with high . The vertical velocity distribution in Fig. 3f indicates the presence of strong updraft, with vertical velocity up to 35 m s^{−2} in the upper troposphere for air panels with of about 340 K.

### c. MAFALDA cycles

As discussed in Pauluis and Mrowiec (2013), the isolines of the isentropic streamfunction correspond to the mean trajectories in the *z*–*θ*_{e} space in the sense that the conditional averaged vertical velocity and diabatic tendency are parallel to the isolines of *ψ*. In MAFALDA, we use the isopleths of the streamfunction to represent “parcel” trajectories. In practice, we select a number of trajectories with :

The parameter *λ* is used here as a free parameter along a given trajectory. Large values of *k* correspond to a more negative value of the streamfunction and are typically associated with shallow overturning, while lower values of *k* are associated with deeper overturning. The positive values of the streamfunction are omitted from this analysis. These are associated with the overshoot of convective towers near the tropopause and the breaking gravity waves in the stratosphere. While present in the isentropic analysis, these account for a small fraction (less than 5%) of the total mass transport.

### d. Thermodynamic cycles

The various properties of the flow along the streamlines are defined as the isentropic average of the property at the same value of *z* and . For example, the moist entropy per unit mass of dry air *s* along the *k*th trajectory is obtained as follows:

This approach makes it possible to define a set of trajectories based on the isentropic streamfunction in *z*–*θ*_{e} space and by using the conditional average to determine the thermodynamic properties of the air parcels along these trajectories.

Figure 4 shows the MAFALDA trajectories from our simulations in four different coordinates: moist entropy and temperature (*s*–*T*), buoyancy and height (*b*–*z*), water content and height (*r*_{T}–*z*), and water content and Gibbs free energy (*r*_{T}–*g*_{υ}). The coordinate pairs were each chosen for their significance for the thermodynamical behavior of convection, which will be discussed in section 3. The definitions of the thermodynamic variables can be found in the appendix. The choice of axes in Fig. 4 is such that all trajectories are counterclockwise. Three points labeled 1–3 have been added to the figures: point 1 corresponds to cold dry air near the surface, point 2 corresponds to warm moist air near the surface, and point 3 corresponds to upper-tropospheric conditions. The thermodynamic trajectories can be divided into three steps: isobaric heating and moistening , warm and moist expansion , and cold, dry compression .

The cycles in *s*–*T* coordinates shown in Fig. 4 are quite similar to the isentropic streamfunction in *z*–*θ*_{e} coordinates. This is not surprising, as temperature can be thought of as a proxy for height and entropy as one for the equivalent potential temperature. Surface heating and moistening appears as an increase in entropy at high temperature. It is followed by the expansion characterized by a decrease in temperature. This expansion is associated with a decrease in entropy resulting from the diffusive loss of water vapor due to entrainment of dry air in the updrafts. For the deeper cycles, there is a slight increase in entropy at temperature below freezing. The entropy *s* used here is the moist entropy with respect to liquid water. As ice here has a negative specific entropy, a loss of ice corresponds to a gain of entropy. Finally, during compression , temperature increases. In the first part of the compression , entropy decreases because of the radiative cooling, but in the second portion , entropy increases because of the humidity gained from detrainment. Overall, the thermodynamic cycles in *s*–*T* coordinates exhibit a structure similar to that of a heat engine, with the parcel’s entropy increasing at high temperature and decreasing at lower temperature, but do not correspond to the theoretical Carnot cycle, as entropy changes do not occur at constant temperature but rather are spread through the entire atmosphere.

In the buoyancy–height coordinates (*b*–*z* in Fig. 4), the warming and moistening step appears as a slight increase in buoyancy. During the ascent , the buoyancy is at first small, but it increases markedly above the cloud base (2a in the figure). Buoyancy is highest for the deeper cycles, as these are associated with higher values of during the ascent. There is also a noticeable increase in buoyancy near the freezing level (2b in the figure). During compression , buoyancy is close to zero or slightly negative for all cycles, which confirms that compression occurs primarily through the slow subsidence in the environment. Above the cloud base, the buoyancy of rising air parcels is higher than that of the descending air, and the cycles are associated with a net generation of kinetic energy. Near the surface, however, several cycles are characterized by the ascent of heavier air and subsidence of lighter air. This indicates the presence of a convective inhibition near the surface, in which air parcels must first be mechanically forced before becoming unstable once they reach their level of neutral buoyancy.

In water content and height coordinates (*r*_{T}–*z* in Fig. 4), the moistening step shows an increase in water content due to evaporation at the surface. The expansion shows the decrease of water content due to precipitation and entrainment. The compression step shows a continuous increase in humidity due to detrainment and reevaporation of precipitation in the free troposphere. For all cycles, the water content is higher for ascending air than for descending air, corresponding to a net upward transport of water.

Figure 4 shows the thermodynamic cycles in the mixing ratio and Gibbs free energy (*r*_{T}–*g*_{υ}) coordinates. The specific Gibbs free energy of water vapor is given here by

The Gibbs free energy depends primarily on relative humidity, with the impacts of the temperature being comparatively small. Surface warming and moistening appears here as an increase in water content occurring for and . The negative value of the Gibbs free energy confirms that water is evaporating in unsaturated air. The ascent can be subdivided into three substeps. First, between the surface and the cloud base, , unsaturated parcels rise with approximately constant mixing ratio, but their relative humidity increases, and therefore so does their Gibbs free energy. Between the cloud base and freezing level , the Gibbs free energy of water vapor is equal to that of liquid water, which is small because of the choice of the reference state (see the appendix). Water gradually condenses and precipitates. Above the freezing point , the water content decreases at negative value of the Gibbs free energy. The negative Gibbs free energy here arises from the fact that water vapor is in equilibrium with ice but unsaturated for liquid water. Subsidence is characterized by an increase of water content in unsaturated air (i.e., corresponding to negative values of the Gibbs free energy). Note that the Gibbs free energy of water vapor is systematically higher when water is removed than when it is being added. The importance of this will be discussed in the next section.

## 3. Thermodynamic analysis of moist convection

The MAFALDA trajectories defined above are cyclical trajectories, along which one can compute the value of any thermodynamic variables. These trajectories can thus be viewed as thermodynamic cycles and analyzed under the framework laid out in Pauluis (2011). Changes of entropy can be related to changes in enthalpy, pressure, and air composition by the Gibbs relationship:

Here, *s* is the moist entropy per unit mass of dry air, *h* is the total enthalpy per unit mass of dry air, is the specific volume per unit mass of dry air, and , , and are the Gibbs free energy of water vapor, liquid water, and ice, respectively. The definitions of the different thermodynamic variables are provided in the appendix. Equation (7) is written assuming a fixed mass of dry air, but water content is allowed to vary freely.

For a cycle, the work per unit mass of dry air is given by the integral

This can be rewritten using the relationship (7) as

Here, is the work that would be done by a Carnot cycle transporting the same amount of entropy. We refer to the term as the Gibbs penalty. It accounts for the thermodynamic consequences of adding and removing water under different thermodynamic conditions. As discussed in Pauluis (2011), water is typically added as unsaturated water vapor at low Gibbs free energy and is removed as liquid or ice with a comparatively higher value of the Gibbs free energy. Hence, the Gibbs free energy increases through the thermodynamic cycle, which corresponds to a reduction of the mechanical output [ in (9)]. Note that, in his analysis of the idealized steam cycle, Pauluis (2011) only considers changes in water vapor mixing ratio and assumes that, when present, liquid water is in thermodynamic equilibrium with water vapor. The expression for the Gibbs penalty used here includes changes of water in all three phases and can take into account nonequilibrium phase transitions, such as the freezing of supercooled water or evaporation of water vapor in unsaturated air.

In an anelastic model, such as SAM, the pressure term in the momentum equation is split between a vertical acceleration due to buoyancy and a nonhydrostatic pressure term. In such cases, the mechanical work is related to the buoyancy flux and vertical transport of water vapor by

Here, is the specific volume per unit mass of moist air; *G* is the gravitational acceleration. Under the anelastic approximation, and are the specific volume and density profiles of the reference atmosphere. The buoyancy *b* is defined as . The higher-order term (HOT) in (10) is associated with the cross product and is significantly smaller than the buoyancy terms, as total water vapor content is small [ within Earth’s atmosphere]. The first term on the right-hand side cancels out when integrated over a cyclical trajectory so that the work in (8) is given by

The term corresponds to the upward transport of buoyancy by convective motions and is directly related to the generation of kinetic energy by the resolved motions. The second term is the work necessary to lift the water in all its phase and is dissipated as precipitation falls down [see Pauluis et al. (2000) for further discussion].

The units of the various integrals in (12) are in joules per kilogram. These correspond to a work per unit of mass of dry air circulating through the cycle. The left-hand side corresponds to the maximum work that could be produced by a Carnot cycle. This maximum work can be decomposed—on the right-hand side—as the sum of the generation of kinetic energy by the resolved motions , the increase of the geopotential energy of the condensed water , and the Gibbs penalty (i.e., the reduction in mechanical work required to balance the increase of the Gibbs free energy of the water). Equation (12) thus makes it possible to compare the mechanical work produced by a thermodynamic cycle to that of a Carnot cycle and to identify why the generation of kinetic energy may be less than the theoretical maximum.

The terms in (12) can be directly computed for the mean trajectories obtained under MAFALDA:

where the suffix *k* refers to the cycle index. Graphically, the terms , , and are equal to the areas contained within the *k*th cycles in Fig. 4. Similarly, is equal to the area within the cycle in Fig. 4 multiplied by the gravitational acceleration *G*. Figure 5 shows the decomposition of (12) for 20 cycles. Figure 5 shows the contribution of the last three terms normalized by the maximum work.

The maximum work for the deepest cycle is about . If all this work were produced and converted into kinetic energy, this would yield a velocity of approximately . However, the actual generation of kinetic energy is much lower, with . The maximum vertical velocity along cycle 1 is only , corresponding to a kinetic energy of , which indicates that most of the buoyancy work is transferred to the rest of the fluid through the nonhydrostatic pressure field. The lifting of water by the same cycle is associated with an increase of geopotential energy of about . This corresponds to lifting of water—the maximum amount of water in cycle 1—by about 3200 m. The total mechanical work in cycle 1 is thus about . In contrast, the Gibbs penalty amounts to about , which accounts for 42% of the maximum work. This Gibbs penalty corresponds to adding of water vapor at a relative humidity of 75% and removing it at saturation.

Cycles 1–5 correspond to deep convection with cloud tops between 10.5 and 8 km, respectively. The maximum work decreases from for cycle 1 to for cycle 5. These large values of for deep convection can be attributed both to the large difference of entropy between the expansion and compression portions of the cycle and to the large temperature difference between the energy source at the surface and its sink in the upper troposphere, as seen from Fig. 4. The buoyancy flux drops more sharply than : while the buoyancy flux accounts for 32% of the maximum work in cycle 1, its contribution drops to of the maximum work in cycle 5. Conversely, the relative contribution of Gibbs penalty increases from cycle 1 to 5, going from 42% of the maximum work in cycle 1 to 62% in the fifth cycle. The relative contribution associated with the lifting of water appears to be relatively constant, at about one-third of the maximum work. Thus, shallower cycles not only have a smaller maximum work but also exhibit a comparatively larger Gibbs penalty. Both these factors lead to a significant reduction of the production of kinetic energy in shallow cycles.

Cycle 6 reaches an intermediate height of about 5.5 km, similar to that of congestus clouds. While the corresponding maximum work is still substantial, with , the corresponding buoyancy flux is weak, about . Cycles 7–10 correspond to shallow overturning, with cloud top between 2 and 4 km. They are associated with very small generation of kinetic energy. In fact, the buoyancy fluxes for cycles 8 and 9 are negative, indicating that these cycles are mechanically forced. The Gibbs penalty in these cycles accounts for the bulk of the maximum work. These cycles still perform some mechanical work, between 10% and 28% of the maximum work, but they are almost entirely used to lift water. Cycles 11–20 never reach saturation: they correspond to shallow dry convection confined to the subcloud layer. They are associated with very low values of . It should be noted that the horizontal resolution of the model (1 km) is insufficient to accurately resolve shallow overturning cells. While cycles 8–20 are thermodynamically consistent, they correspond to the behavior of underresolved shallow overturning in a numerical model and may not be representative of the behavior of shallow convection.

The integral in (12) yields the amount of work generated per unit mass of air circulated in each cycle. The total work produced by atmospheric motions can be obtained by multiplying the work by the mass transport in each cycle:

where is the mass transported in cycle *k*. In MAFALDA [see (5)], the mass transport is the same across all cycles, with .

The results of the global analysis are shown in Table 1. The maximum work is . The generation of kinetic energy accounts for , which amounts to 18% of the total. The lifting of water consumes kinetic energy at a rate of . The Gibbs penalty for the whole simulation is . The total budget is strongly dominated by the deepest cycles, as these are the ones that are associated with the largest mechanical output per unit of air circulated. It confirms that convection in our simulation is greatly affected by moist processes. The Gibbs penalty and lifting of condensed water greatly reduce the amount of work that is available to sustain convective motions, and kinetic energy is generated at a rate that is one order of magnitude smaller than what would be expected from an equivalent Carnot cycle.

## 4. Depth of convection and mechanical efficiency

The behavior of these thermodynamic cycles indicates that the partitioning of the maximum work between kinetic energy generation, water lifting, and Gibbs penalty is significantly affected by the depth of the overturning. In this discussion, we focus primarily on the partitioning of the maximum work between Gibbs penalty, water loading, and buoyancy flux for deep convection (e.g., cycles 1–6). In particular, Fig. 5 shows that

the relative contribution of the water lifting remains roughly unchanged;

the relative contribution of the Gibbs penalty decreases for deeper cycles;

the relative contribution of the buoyancy flux increases for deeper cycles.

This behavior can be understood in light of the differences between the cycles, as shown in Fig. 4. A deeper cycle corresponds not only to an increase in the temperature difference at which the cycle is operated but also to an increase in the entropy difference between the updraft and downdraft. This means that the maximum work, which is proportional to the product of both, will increase rapidly. Similarly, the contribution of the lifting of water increases not only because water is transported to a deeper level, but also because more water is being transported. The entropy difference between the updrafts and downdrafts is approximately proportional to the difference in water vapor content: that is,

where the subscripts “up” and “down” indicate that the value of the variable is taken during the ascending and subsiding portions of the cycle. Similarly, the temperature difference between the surface and the top is proportional to the height difference

Here, the subscripts “sfc” and “top” refer to the value of the variables at the surface and at the top of the cycle. Hence, the contribution from the water lifting is proportional to the maximum work:

The temperature scale height is defined as and has a value of about 50 km for a lapse rate of 6 K km^{−1}, which is typical of the lower troposphere. The scaling (17) indicates that the increase in geopotential energy of the water accounts for about 20% of the maximum work (i.e., ). This corresponds to the lower end of the relative contribution for in Fig. 5. The underlying argument in (17) ignores the lifting of condensed water. Including the water loading would increase the ratio . The value of observed in cycles 1–5 corresponds to a situation where the condensed water mixing ratio in the updraft is equal to half the difference in water vapor between updraft and downdraft [i.e., ].

The Gibbs penalty shows a smaller sensitivity to the depth of the cycle. This can be again understood by looking at Fig. 4. While the difference in Gibbs free energy between the ascending and descending branch increases for deeper cycle, the total change in humidity remains approximately constant. Once convection reaches a depth of 5 km or more (corresponding approximately to twice the scale height for saturation water vapor), most of the water vapor present in the updrafts condenses during its ascent. Further increasing the depth of the convection does not substantially increase condensation. The difference in Gibbs free energy between the ascending and subsiding air parcels is primarily determined by the difference in water content:

It follows that the Gibbs penalty can be approximated by

The relative contribution of the Gibbs penalty to the cycle then can be approximated as follows:

For deep convection, one can assume that all the water condenses and . The fraction of the maximum work that is lost because of the Gibbs penalty is approximately given by

This implies that the relative importance of the Gibbs penalty is inversely proportional to the depth of the cycle, measured here in terms of the temperature difference between the heat source and heat sink. Within Earth’s atmosphere, the deepest cycle would be associated with a temperature difference *T*_{sfc} minus *T*_{top} corresponding to the temperature difference between Earth’s surface and the tropopause (i.e., ). For such a deep cycle, the Gibbs penalty would account for about one-third of the maximum work. The relative contribution of the Gibbs penalty decreases with the depth convection, and it is similar to the maximum work for . This would correspond to tropical convection reaching slightly above the freezing level. The scaling in (21) explains the increase of the relative contribution of the Gibbs penalty from cycle 1 to cycle 5 in Fig. 5. As the relative contribution of the water lifting is relatively constant, as implied by the scaling in (17), the relative increase in must be compensated by a decrease of the contribution from the buoyancy flux . As convection becomes deeper, the impact of the Gibbs penalty is reduced, and convection can generate more kinetic energy per unit mass of air transported.

For shallow convection, the scaling implies that the Gibbs penalty should be on the same order of magnitude as the maximum work. Pauluis (2011) shows that, in an unsaturated cycle, the Gibbs penalty accounts for approximately of the maximum work, which matches the behavior of the shallow cycles 10–20. Finally, it should be stressed that the scaling in (4) assumes that the difference of entropy between ascending and descending parcels is mostly due to the difference in mixing ratio [see (15)]. If this is not the case, the arguments above break down. In our simulations of oceanic convection, the Bowen ratio (i.e., the ratio of the sensible heat flux to the latent heat flux) is small, and this assumption is well justified. However, as discussed in Pauluis (2011), an increase in the Bowen ratio would lead to an increase in the mechanical efficiency of the heat engine.

## 5. Impacts of temperature increase on thermodynamic cycles

Radiative–convective equilibrium is simulated a second time after increasing the CO_{2} concentration from 380 to 760 ppm, while keeping everything else unchanged. The SST equilibrates to about 304.3 K, and both temperature and humidity increase throughout the troposphere. The climate sensitivity here is high when compared to global climate models, with a temperature increase of about 3.8 K for a doubling of CO_{2} here. Several aspects of the simulations might contribute to such high climate sensitivity: the use of a “tropical only” domain can overestimate the water vapor and cloud feedbacks; the resolution is insufficient to properly capture the dynamics of shallow convection, so the low-cloud feedback is likely inadequate; and the ocean cooling of does not account for possible changes in the global circulation. The purpose of this experiment here is not to be an accurate representation of atmospheric response to a CO_{2} doubling, but rather to illustrate how MAFALDA can be used to diagnose a change in the thermodynamic behavior of convection.

The isentropic mass flux and streamfunction for the 2×CO_{2} experiment are shown in Fig. 6 and can be directly compared to the reference case in Fig. 2. The warming is evident through the shift of the circulation toward higher values of . The increase in of about 15 K is much larger than the temperature increase and reflects the large increase in water vapor content. Convection also deepens slightly, reaching 12 km in the 2×CO_{2} case instead of 10.5 km in the control run. The streamfunction broadens in the sense that the difference in equivalent potential temperature between ascending and descending air increases. If one assumes that the relative humidity of the descending air remains unchanged in the 2×CO_{2} experiment, the difference of water content between updrafts and downdrafts increases with temperature following the Clausius–Clapeyron relationship, which is captured here by a broadening of the isentropic streamfunction. Finally, there is a slight reduction in the overall mass transport, with the absolute minimum of the streamfunction decreasing from 0.028 to 0.026 kg m^{−2} s^{−1}.

The MAFALDA trajectories for the 2×CO2 experiment are shown in Fig. 7. In *T*–*s* coordinates (Fig. 7a), heating at Earth’s surface corresponds to the proportion of the cycle characterized by entropy increase at almost constant temperature. The warming is clearly evident through shift of the warming toward higher temperatures. It is also visible in the shift of all cycles toward higher entropy values. The entropy difference between rising and subsiding air also increases, a reflection of the similar increase in the difference in . The trajectories in *b*–*z* coordinates show mostly the deepening of convection. Updrafts exhibit a slight increase in their buoyancy. In both *r*_{T}–*g*_{υ} and *r*_{T}–*z* coordinates, the moistening of the atmosphere is evident, with the water content in ascending air increasing from about 16 to 20 g kg^{−1}.

The thermodynamic integrals in the 2×CO_{2} experiment are shown in Fig. 8. The maximum work increases for all cycles. For the deepest cycle, the maximum work increases from 1570 to 2020 J kg^{−1}. This large increase is due primarily to the larger entropy transport by the cycle and, to a lesser extent, to the deepening of convection. The increase in also translates into systematic increases for all three terms on the right-hand side of (12). As temperature increases, the increase in water content also leads to an increase in the energy that can be transported per unit of mass of air. The work performed per unit mass of air increases significantly, and individual convective events become more intense.

The relative contributions of the buoyancy flux , water loading , and Gibbs penalty are only slightly affected by the warming. The Gibbs penalty remains the dominant term for all cycles, followed by the water loading and the buoyancy flux. A closer look indicates that the relative contribution of the water loading increases slightly. The scaling in (17) indicates that the relative contribution of is proportional to the temperature scale height . In an atmosphere in radiative–convective equilibrium, the temperature scale height is close to that of moist adiabatic ascent and increases substantially with temperature. Form a physical point of view, a warmer and moister atmosphere exhibits a smaller lapse rate, which implies that condensation within updrafts occurs on average at a higher level, and thus more work is done to lift water.

The total work done by the atmosphere is obtained by summing the contributions from all cycles [see (14)]. The results for the 2×CO_{2} experiment are shown in Table 1. The maximum work for the 2×CO_{2} experiment is . Generation of kinetic energy accounts for , the lifting of water for , and the Gibbs penalty for about . The increase in maximum work is about 20%, which is less than the 30% increase observed in the maximum work per unit mass of air. The difference is due to a 10% reduction in the convective mass transport in the 2×CO_{2} experiment when compared to the reference simulations. While individual convective events are more intense in a warmer atmosphere, they are also less frequent, as indicated by the reduction of the convective mass transport.

Figure 9 shows that, in the 2×CO_{2} experiment, the energy transport per unit mass increases by about 25%, while the mass transport decreases by about 10% when compared to the reference simulation. As a result, the upward energy transport increases by about 15%. Our results differ from the analysis of Held and Soden (2006), who argue that, in a warming atmosphere, radiative constraints prevent the upward energy transport from increasing significantly and that, as a result, the overall convective mass transport should be reduced in a warming atmosphere. They speculate that the convective mass transport should be reduced by about 5% per 1 K of warming. In contrast, the convective overturning in our simulations weakens only by 2.5% per 1 K of warming, in large part because there is a substantial increase in the upward energy transport. Held and Soden (2006) note a similar discrepancy between global climate models and one-dimensional radiative–convective equilibrium, but the cause of this discrepancy is unclear. Had the mass transport decreased by 20%, consistent with Held and Soden (2006), we would expect the maximum work to increase by about 10% rather than the 20% increase reported in Table 1.

## 6. MAFALDA cycles versus Lagrangian trajectories

Ideally, one would rather perform the thermodynamic analysis on the actual trajectories of air parcels. However, we have to confront the fact that there are an infinite number of parcel trajectories, that they are all distinct, and that these trajectories are not periodic (i.e., they cannot be described as a thermodynamic cycle). A practical alternative, which is the basis for MAFALDA, is to study instead some idealized trajectories that are representative of the flow. In doing so, we have made two fundamental approximations. First, the mean flow is interpreted as parcel trajectories. Second, the thermodynamic properties along these trajectories are given by their conditionally averaged value. While both of these approximations should be regarded as ad hoc solutions, we offer here some basis for their use in MAFALDA and discuss the inherent uncertainties they create.

It should first be recognized that “mean flow” is not as simple a concept as it may seem, Indeed, it has been long known in atmospheric and oceanic sciences that the mean flow depends critically on choice of the coordinate system in which the averaging is done. Differences between the meridional circulation in Eulerian and isentropic coordinates (Held and Schneider 1999; Pauluis et al. 2008, 2010) or the disappearance of the Deacon cell when the oceanic circulation is averaged on isopycnal surfaces (Karoly et al. 1997; McDougall and McIntosh 2001) are well-known examples of how the coordinate system in which the flow is averaged affects the mean circulation.

As moist convection is highly turbulent, the mean flow depends strongly on the coordinate system used for averaging. For example, in the simulations discussed above, the Eulerian mean circulation [i.e., the flow averaged in (*x*, *y*, *z*) coordinates] vanishes. Trajectories based on the Eulerian mean flow would not provide any useful information. This problem is circumvented in MAFALDA by averaging the flow in *z*–*θ*_{e} coordinates, which takes advantage of the fact that the equivalent potential temperature is an adiabatic invariant and thus conserved for reversible adiabatic ascent. Pauluis and Mrowiec (2013) show that the isentropic averaging indeed distinguishes between rising air parcels in high-*θ*_{e} updrafts and subsiding air with lower equivalent potential temperature. Furthermore, Pauluis and Mrowiec (2013) also find that the bulk of the ascent occurs in the region of the *z*–*θ*_{e} space, where the flow is mostly upward [see Fig. 2c in Pauluis and Mrowiec (2013)]. In other words, equivalent potential temperature is a very good choice to distinguish between rising and descending parcels.^{2} These results give us some confidence that the mean flow in MAFALDA does capture some key aspects of the parcel trajectories.

While the isentropic analysis does a good job at separating ascending and subsiding air parcels in the boundary layer and free troposphere, it is less accurate near the tropopause. As pointed out in Pauluis and Mrowiec (2013), convective overshoot near the tropopause corresponds to positive values of the streamfunction, which is visible between 12 and 14 km in Fig. 2. Unfortunately, this means that the streamfunction goes from negative through most of the troposphere to slightly positive near the tropopause. This makes it difficult to associate parcel trajectories with a unique value of the streamfunction near the tropopause. This could potentially be improved in future work by using a different thermodynamic coordinate, such as total water content, to better capture the overturning near the tropopause.

A second key approximation lies in the use of conditional averaging to estimate the thermodynamic state variables in MAFALDA. This assumes that an arbitrary thermodynamic variable, such as entropy or temperature, can be expressed as a function of height and equivalent potential temperature alone. This can be partially justified by Dalton’s law, which states that, under thermodynamic equilibrium, all properties of moist air can be inferred from three state variables. The pressure of an air parcel can be approximated by the horizontal-mean hydrostatic pressure .^{3} If, in addition to the height and , one were to know the total water content, then any thermodynamic property would be uniquely determined. For example, one can write the entropy of a parcel as

Instead, MAFALDA uses the conditionally averaged entropy (i.e., the average entropy of all air parcels at a given value of the height and but varying total water content). The parcels within this subset do not have the same entropy. Variations of entropy at a fixed *z* and depend both on variations of total water content among all parcels with the same pressure and and on the partial derivative of entropy with respect to water content at fixed values of pressure and . To quantify this, we compute the standard deviation of a quantity *f* at a given height and equivalent potential temperature as

Note that this standard deviation corresponds to the variations of *f* among all air parcels with the same equivalent potential temperature and at the same height. Figure 10 shows the standard deviation for the total water content. The fluctuations are small, less than 0.0005 kg kg^{−1} for most of the atmosphere, but can reach up to about 0.002 kg kg^{−1} at high altitude and , which corresponds to strong updrafts. As these air parcels are saturated, the large value of corresponds primarily to fluctuation in the amount of condensed water present in the parcels. Figure 10 shows the standard deviation for temperature and moist entropy. The largest value of the standard deviation for temperature is about 0.5 K, while the entropy fluctuations are less than 1.5 J K^{−1}. The low values of the standard deviation for , *s*, *T*, and *b* (not shown) in comparison to the variations of these quantities through the various MAFALDA cycles strongly reinforces our confidence in the quantitative estimates for thermodynamic integral , , and . The standard deviation for the Gibbs free energy shows regions of relatively large value, with up to in the upper troposphere above 8 km. These are associated with fluctuations of relative humidity in the unsaturated environment. As the absolute amount of water vapor is small, these fluctuations have little impact on but can greatly affect the Gibbs free energy. Fortunately, air parcels with large uncertainty in the Gibbs free energy correspond to those with little change in water content, and their contribution to the Gibbs penalty in (12) is small. Overall, using the conditionally averaged thermodynamic variables to compute for the various thermodynamic integrals in (12) yields very accurate estimates for , , and , and slightly larger error for the Gibbs penalty .

Recently, Laliberte et al. (2015) have proposed a new framework to diagnose thermodynamic cycles in climate models. Their approach relies on computing the mean circulation in a three-dimensional thermodynamic space instead of the two-dimensional *z*–*θ*_{e} coordinates used in MAFALDA. This method is slightly more computationally expensive than MAFALDA because it requires the explicit computation of Lagrangian tendencies for different thermodynamic variables at all grid points. It offers, however, the additional advantage that all thermodynamic integrals can be computed exactly, without having to approximate the state variables by their conditional average. Future work should compare the idealized cycles obtained under MAFALDA with both Lagrangian trajectories and the thermodynamic reconstruction using the methodology of Laliberte et al. (2015). Understanding the differences between these approaches would shed some light on how turbulent motions can affect the range of thermodynamic transformations that take place in various atmospheric flows.

## 7. Conclusions

The paper proposes a new approach, the Mean Airflow as Lagrangian Dynamics Approximation, to diagnose the thermodynamic transformations associated with complex atmospheric motion. To do so, one first computes the isentropic streamfunction following the approach of Pauluis and Mrowiec (2013). The isolines of the streamfunction correspond to the mean flow in *z*–*θ*_{e} coordinates and are used to define a set of cyclical trajectories. One can then compute the value of any thermodynamic variables along these main trajectories through conditional averaging. Altogether, MAFALDA makes it possible to extract a set of idealized trajectories from a complex numerical simulation and to analyze the thermodynamic transformations that are taking place along these trajectories.

The core approximation in MAFALDA lies in interpreting the mean trajectories in *z*–*θ*_{e} coordinates as parcel trajectories. As with any turbulent flow, mean velocity can differ significantly from any individual trajectory. This problem is partially alleviated in MAFALDA by the use of an adiabatic invariant *θ*_{e} as a coordinate. Pauluis and Mrowiec (2013) show that such isentropic averaging can be successfully used to quantify many features of convective motions. An important issue for further research is to understand the extent to which such isentropic averaging captures Lagrangian trajectories.

MAFALDA has been used here to reconstruct the thermodynamic cycles associated with convection in radiative–convective equilibrium in high-resolution simulations with SAM. The convective cycles vary between deep convection with cloud top above 5 km, shallow convection with cloud top between 1.5 and 5 km, and convection within the sublcoud layer for convective cycles that never reach saturation. Deep convection accounts for about 30% of the total mass flux, shallow convection for about 15%, and unsaturated cycles for the remaining 55%.

Following Pauluis (2011), the maximum work that could be performed by each individual cycle is decomposed into kinetic energy generation, increase in geopotential energy of water, and a Gibbs penalty that accounts for the reduction of the mechanical output of the atmospheric heat engine due to the hydrological cycle. The maximum work that could be performed by the MAFALDA cycles—defined as the mechanical output of a Carnot cycle with the same entropy sources and sinks—depends strongly on the depth of the cycle. Unsaturated cycles act within small temperature fluctuations and are associated with a very weak maximum work, typically less than . Shallow cycles exhibit intermediary values of the maximum work reaching up to . The maximum work is significantly larger for deeper cycles, reaching up to .

For all the cycles, the maximum work is much larger than the actual mechanical output. This reduction of the mechanical work is tied to the fact that water is injected as water vapor in unsaturated air but removed through condensation in saturated air. This implies that the Gibbs free energy of the water is lower when it is added to the atmosphere than when it is removed. As a result, there is a reduction of the mechanical work produced in the cycles that can be directly estimated in terms of a Gibbs penalty, as discussed in Pauluis (2011). For the shallow and unsaturated cycles, the Gibbs penalty accounts for up to 80% of the maximum work. The relative contribution of the Gibbs penalty is lower for deep cycles but still accounts for about half of the maximum work.

A large fraction of the work produced by the convective cycles is used to increase the geopotential energy of water rather than to produce kinetic energy. This geopotential energy is then dissipated in the microphysical shear zone surrounding falling hydrometeors (Pauluis et al. 2000; Pauluis and Dias 2013). In our simulations, most of the work is used for this purpose, and there is very little kinetic energy production. Even for the deepest convective cycle, lifting water still accounts for about half of the total mechanical work.

The thermodynamic behavior of the individual cycles can be aggregated together to obtain the mechanical energy production of the atmosphere as a whole. As the deeper cycles are also associated with the strongest energy transport, they dominate the aggregated energy budget. When considering the atmosphere as a whole, the Gibbs penalty accounts for a little more than half of the mechanical work, lifting of condensed water account for 30%, and the generation of kinetic energy is a little under 20% of the maximum work. These results are in good agreement with the analysis of the entropy budget by Pauluis and Held (2002a,b), who found that moist processes significantly reduce the production of kinetic energy by convection.

We also compared the thermodynamical behavior of convection after doubling the CO_{2}. We found that the increase in atmospheric temperature led to an increase of the upward energy transport per unit of mass of about 5% per 1 K of temperature increase. This leads to a similar increase in the maximum work per unit of mass associated with the convective cycles. The partitioning of the maximum work between kinetic energy generation, water lifting, and Gibbs penalty did not change much, albeit there is a small increase in the relative contribution of water loading. This increase is consistent with a reduction of the lapse rate and an increase of the level at which condensation occurs. Overall, our analysis confirms that a warmer and moister world should be associated with more intense but less frequent convective events.

The Mean Airflow as Lagrangian Dynamics Approximation and a similar approach recently introduced by Laliberte et al. (2015) offer new methods to systematically study the thermodynamic behavior of complex numerical simulations of the atmosphere. Previous studies, such as Pauluis and Held (2002a), rely on detailed diagnostics of the thermodynamic transformation as they occur in the model. While such approaches can be extremely informative, they also require an in-depth knowledge of many aspects of the model at hand, which makes them difficult to reproduce across a wide array of numerical models. In contrast, a key advantage of MAFALDA and the thermodynamic reconstruction of Laliberte et al. (2015) lies in that they rely solely on three-dimensional snapshots of the flow and can be implemented while treating the numerical model as a black box. This opens up the possibility for systematic investigations of the thermodynamic transformations that take place in atmospheric flows and should help us better understand how thermodynamics and dynamics interact with each other in the global climate system.

## Acknowledgments

Olivier Pauluis is supported by the New York University in Abu Dhabi Research Institute under Grant G1102.

### APPENDIX

#### Gibbs Free Energy of Water in Different Phases

The specific Gibbs free energy is defined as the difference between its enthalpy and its entropy multiplied by the temperature:

As the exact value of both the entropy and enthalpy of a system depends on the choice of the reference state, so does the Gibbs free energy. For the study of moist convection, a convenient choice for the reference state is liquid water at the freezing temperature . In this case, the specific enthalpies of water vapor , liquid water , and are

Similarly, the corresponding specific entropies , , and are

For this choice of the reference state, the specific Gibbs free energy is therefore:

In the Gibbs relationship in (7)

moist air is treated as a mixture of 1 kg of dry air, kg of water vapor, kg of liquid water, and kg of ice. The corresponding entropy and enthalpy per unit mass of dry air are

The specific volume is the specific volume per unit mass of dry air:

## REFERENCES

*Thermal Physics of the Atmosphere*. Wiley, 256 pp.

*Atmospheric Thermodynamics*. Oxford University Press, 416 pp.

*Atmospheric Convection*. Oxford University Press, 580 pp.

## Footnotes

^{1}

Note that, for an anelastic model, the density is a function of height alone and can be factored out of the right-hand side of (4). Nevertheless, the isentropic integral of density still depends on the equivalent potential temperature, with , where the isentropic integral is equal to the fractional area covered by air parcels with equivalent potential temperature at level *z*.

^{2}

Of course, vertical velocity does an even better job at distinguishing between rising and descending parcels, but, in contrast to , *w* is not an adiabatic invariant and evolves rapidly.

^{3}

In an anelastic model, such as the one used in this study, Pauluis (2008) shows that the Gibbs relationship is still valid as long as one uses the reference state instead of the total pressure.