## Abstract

The zonal and meridional components of the atmospheric general circulation are used to define a global thermodynamic streamfunction in dry static energy versus latent heat coordinates. Diabatic motions in the tropical circulations and fluxes driven by midlatitude eddies are found to form a single, global thermodynamic cycle. Calculations based on the Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) dataset indicate that the cycle has a peak transport of 428 Sv (Sv ≡ 10^{9} kg s^{−1}). The thermodynamic cycle encapsulates a globally interconnected heat and water cycle comprising ascent of moist air where latent heat is converted into dry static energy, radiative cooling where dry air loses dry static energy, and a moistening branch where air is warmed and moistened. It approximately follows a tropical moist adiabat and is bounded by the Clausius–Clapeyron relationship for near-surface air. The variability of the atmospheric general circulation is related to ENSO events using reanalysis data from recent years (1979–2009) and historical simulations from the EC-Earth Consortium (EC-Earth) coupled climate model (1850–2005). The thermodynamic cycle in both EC-Earth and ERA-Interim widens and weakens with positive ENSO phases and narrows and strengthens during negative ENSO phases with a high correlation coefficient. Weakening in amplitude suggests a weakening of the large-scale circulation, while widening suggests an increase in mean tropical near-surface moist static energy.

## 1. Introduction

The atmospheric general circulation forms as a response to differential solar heating and oceanic heat transport, acting to redistribute energy from warm regions to cold regions, mainly through the exchange of dry static energy (DSE) and latent heat (LH) (Trenberth and Stepaniak 2003a). The geographical distribution of these exchanges dictates, among other climatic features, global patterns of precipitation and evaporation (Muller and O’Gorman 2011). The atmospheric circulation is known to vary as the driving forces (solar heating and ocean heat fluxes) vary.

El Niño–Southern Oscillation (ENSO) is one of the dominant modes of interannual variability in the climate system (Lu et al. 2008). It is characterized by weakened (strengthened) winds and zonal gradients of sea surface temperature (SST) over the Pacific during positive (negative) ENSO phases (Bjerknes 1969; Oort and Yienger 1996). It has been shown that ENSO events also affect tropical-mean surface temperatures (Seager et al. 2003) and the meridional overturning circulation (Lu et al. 2008). The Walker circulation is a zonally asymmetric feature of the total circulation that plays a key role in ENSO events (Bjerknes 1969) and is closely tied to the zonal SST gradient in the tropical Pacific. It redistributes DSE and LH zonally (Trenberth and Stepaniak 2003b) and is driven by convection over the western tropical Pacific and by subsidence over the eastern Pacific and the Atlantic (Peixoto and Oort 1992). Because of its east–west orientation, the Walker circulation is not captured by zonal-mean diagnostics of the meridional circulation. It is often viewed schematically and seldom studied using closed streamfunctions, since its poleward extent is not well defined. This makes most definitions of the Walker circulation strength ambiguous [attempts include Trenberth et al. (2000), Vecchi et al. (2006), Tokinaga et al. (2012), and L’Heureux et al. (2013)].

The meridional transport of dry and moist static energy (MSE, equal to the sum of DSE and LH) can be represented by the zonal-mean meridional overturning circulations in dry and moist isentropic coordinates (Fig. 1) (Townsend and Johnson 1985; Karoly et al. 1997; Held and Schneider 1999; Pauluis et al. 2008; Döös and Nilsson 2011). DSE is approximately conserved for dry isentropic motions (adiabatic motions, where potential temperature *θ* is constant). In the same way, MSE is almost conserved for moist isentropic motions [moist adiabatic motions, where the equivalent (or “moist”) potential temperature *θ*_{e} is constant] (Emanuel 1994). Following Döös and Nilsson (2011), we compute the meridional overturning circulations in isobaric, dry isentropic and moist isentropic coordinates [Fig. 1, and Eq. (A1)] using 6-hourly Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) data. Details of the calculations are presented in appendix A. In the midlatitudes, mass fluxes by synoptic eddies are lost when zonally averaging in isobaric coordinates but included when using dry or moist isentropic coordinates (Peixoto and Oort 1992; McIntosh and McDougall 1996; Held and Schneider 1999; Pauluis et al. 2008; Döös and Nilsson 2011; Laliberté et al. 2012; Kjellsson and Döös 2012). Furthermore, moist isentropic coordinates include fluxes of moist air that cancel in dry isentropic coordinates (Pauluis et al. 2010). As a result, in the midlatitudes the dry and moist isentropic circulations are much stronger than the isobaric circulation, and the moist isentropic circulation is stronger than the dry isentropic circulation, which can be seen in Fig. 1. This indicates the importance of considering both DSE and LH when studying the global atmospheric circulation.

Variations of the atmospheric general circulation can occur in both the meridional and zonal components. There are, however, few diagnostics that include both zonal and meridional overturning circulations. A similar problem exists in the ocean. Recent studies (Döös et al. 2012; Zika et al. 2012) have made use of purely thermodynamic (temperature–salinity) coordinates to represent the circulation in a single picture that describes global exchanges of heat and freshwater as well as water mass transformations. The unified picture then allows for analysis of the interconnected nature and variations of the global ocean circulation (Ballarotta 2013; Zika et al. 2013). In this study we use a similar technique to explore the redistribution of DSE and LH in the atmospheric general circulation. Previous studies have investigated the atmospheric circulation in *θ*–*θ*_{e} (potential and equivalent potential temperature, similar to DSE–MSE) coordinates (Pauluis et al. 2008; Laliberté et al. 2012, 2013) but only included meridional mass fluxes. In a recent study, Pauluis and Mrowiec (2013) examined the circulation in *z*–*θ*_{e} coordinates using a cloud-resolving model, but this only included vertical mass fluxes. We include mass fluxes in the zonal, meridional, and vertical directions to obtain a new description of the global atmospheric circulation as a cycle in purely thermodynamic coordinates (LH–DSE space). This new description gives a unified picture of the zonal and meridional circulations and allows us to estimate the strength of its zonally asymmetric component. We diagnose the variability of the global thermodynamic cycle with a focus on annual and interannual variations and the impacts of ENSO, and also compare the results from ERA-Interim and the EC-Earth Consortium (EC-Earth) climate model.

## 2. Data and method

We use 6-hourly surface pressure *p*_{s}, temperature *T*, specific humidity *q*, and horizontal winds *u*, *υ* on model levels from ERA-Interim (Berrisford et al. 2009) for the period 1979–2009. We use data with 1.25° horizontal resolution and 60 model levels, of which roughly 30 are in the troposphere. Using surface pressure, the pressure at each model level interface *k* is *p*_{k} = *A*_{k} + *B*_{k}*p*_{s}, where *A*_{k} and *B*_{k} are constants (Simmons and Burridge 1981). The geopotential height *z* (m) is calculated by integrating the hydrostatic equation vertically using temperature, specific humidity, and pressure. Vertical velocity *ω* (Pa s^{−1}) is calculated using the model continuity equation following Simmons and Burridge (1981) [see Kjellsson and Döös (2012) for a similar derivation]. We wish to stress that the calculated vertical velocity does not include any explicit parameterization of convection. Rather, it is a diagnostic of horizontal motions at each level. The vertical mass flux is therefore realistic in spatial coordinates. However, the LH and DSE values in a grid box may be different from those that would have been in a convective updraft if convection were fully resolved. In such a case the mass fluxes in *l*–*s* space may be positioned differently and the integral *ψ*(*l*, *s*) could be different. The implications of this are discussed in section 6 of this paper.

We also make use of atmospheric data from historical simulations of years 1850–2005 from the EC-Earth coupled climate model (Hazeleger et al. 2011). The atmospheric component of the EC-Earth model, Integrated Forecasting System (IFS), is also used by ECMWF when generating the ERA-Interim. After more than 1000 years of coupled spinup integration in which radiative forcing and atmospheric composition (including aerosols) are held constant at preindustrial conditions, the model is integrated with observed or reconstructed radiative forcing and atmospheric composition (Meinshausen et al. 2011) from the years 1850–2005. The data has 1.125° horizontal resolution, 62 model levels, and is stored every 6 h. As for the ERA-Interim data, we use *u*, *υ*, *T*, *q*, and *p*_{s} to calculate *z* and *ω*.

Using the six variables *u*, *υ*, *T*, *q*, *z*, and *p*, we calculate LH, *l* = *L*_{υ}*q* (J kg^{−1}); DSE, *s* = *c*_{p}*T* + *gz* (J kg^{−1}); and MSE, *h* = *l* + *s* (J kg^{−1}), where *L*_{υ} = 2.5 × 10^{6} J kg^{−1} is the latent heat of vaporization, *c*_{p} = 1004 J kg^{−1} K^{−1} is the specific heat capacity of dry air, and *g* = 9.81 m s^{−2} is the gravitational acceleration at Earth’s surface (Emanuel 1994; Wallace and Hobbs 2006).

We begin by binning the material derivative of DSE, , from the Eulerian space into an LH–DSE space:

where *dM* = *g*^{−1}*dxdydp* is a mass element, *δ*(*ϕ*) is the Dirac delta function [*δ*(*ϕ*) = 1 when *ϕ* = 0 and 0 otherwise], and primed variables are in *x*–*y*–*p* space. The result is shown in Fig. 2a and the method to obtain it is described in appendix B. We then apply the same binning procedure to the material derivative of LH, :

In Eqs. (1) and (2), and are mass integrals of and , respectively. Both have the unit of kg s^{−1}(kJ kg^{−1})^{−1}. It can be shown that with these definitions and therefore the binning of DSE tendencies is associated with a closed and unique streamfunction; that is,

As global atmospheric circulations tend to have magnitudes of approximately 10^{10} kg s^{−1} (cf. Fig. 1), we will present streamfunctions in units of Sverdrups (Sv), where l Sv ≡ 10^{9} kg s^{−1}.

The streamfunction *ψ*^{tot}(*l*, *s*) represents mass fluxes across a DSE surface and is due to both the movement of the surface (a local Eulerian change ∂_{t}*s*) and the transport of air from one DSE value to another (advective change ). In Eq. (3), there are thus two effects that are binned onto the *l*–*s* plane. The total streamfunction can therefore be seen as the sum of two streamfunctions: *ψ*^{tot} = *ψ*^{loc} + *ψ*^{adv}. The former streamfunction, *ψ*^{loc}, is a streamfunction based only on ∂_{t}*s*. It represents air at a certain point in *x*–*y*–*p* coordinates changing its LH and/or DSE from one time to another. In particular, it will be nondivergent if all *x*, *y*, *p* points have the same initial and final positions in *l*–*s* coordinates. The latter streamfunction, *ψ*^{adv}, is based on advection of DSE. At any instant, it represents global mass fluxes projected onto the LH–DSE plane. For *ψ*^{adv}, nonclosed streamlines arise only in the presence of sources and sinks of mass. Sources and sinks of mass are present in our data since neither reanalysis nor climate models are perfectly mass conserving (Berrisford et al. 2011). We also make the important note that for a purely adiabatic flow, *ψ*^{tot} = 0, in which case *ψ*^{adv} and *ψ*^{loc} are either both zero or equal and of opposite sign.

We show *ψ*^{loc} and *ψ*^{adv} calculated from ERA-Interim data during 1979–2009 in Figs. 2b and 2c, respectively. In Fig. 2c, not all streamlines close. Unclosed streamlines mean that there are net mass fluxes across some DSE surfaces. The maximum mass flux outside closed streamlines is small and corresponds to a maximum error in mass conservation of less than 5% of the amplitude. In Fig. 2b, there are also errors in *ψ*^{loc} corresponding to about 13% of the amplitude. These errors are mostly due to the fact that a closed circulation in *ψ*^{loc} can only be obtained if a DSE surface returns to its initial position at the end. We find that *ψ*^{adv} contributes much more to *ψ*^{tot} than *ψ*^{loc} does. For this reason, we will assume that for the data that we are analyzing

However, we do not argue that this approximation holds generally (e.g., for different datasets).

We denote *ψ*^{adv} as the hydrothermal streamfunction to be consistent with terminology used for the oceanic thermohaline streamfunction in temperature–salinity coordinates in previous studies (Zika et al. 2012; Döös et al. 2012; a quantification of *ψ*^{adv} is presented in Groeskamp et al. 2013, manuscript submitted to *J. Phys. Oceanogr.*). We will drop the superscript and denote the hydrothermal streamfunction as *ψ*(*l*, *s*). The hydrothermal streamfunction can be expressed succinctly as

where *μ*[*l* − *l*′] is the Heaviside function, which is the integral of the Dirac delta function, *δ*[*l* − *l*′]. The possible effects of including *ψ*^{loc} in the analysis will be discussed in section 6.

## 3. The hydrothermal cycle

We show *ψ*(*l*, *s*) calculated from ERA-Interim data 1979–2009 (Fig. 2c). The result is a single anticlockwise cycle with an amplitude of −428 Sv. We define the hydrothermal streamfunction maximum as the difference between the maximum and minimum value along the LH axis for each value of DSE (Fig. 2d) and note that the maximum value occurs at DSE ≈ 320 kJ kg^{−1}.

The outer streamlines (e.g., −30 Sv in Fig. 2c) of *ψ*(*l*, *s*) have three branches that share characteristics with different atmospheric processes. One branch extends from *l* = 45 kJ kg^{−1}, *s* = 300 kJ kg^{−1} to *l* = 0 kJ kg^{−1}, *s* = 345 kJ kg^{−1} and is bounded by a moist adiabat (constant MSE, *h* ≈ 345 kJ kg^{−1}), shown in Figs. 2a–c as a dashed line. It exhibits the behavior of ascent with precipitation where LH is converted into DSE such that MSE is approximately conserved (Emanuel 1994). It is worth noting that this branch has a minimum of MSE at about *l* = 15 kJ kg^{−1}, *s* = 320 kJ kg^{−1}. We speculate that convective detrainment of convecting near-surface tropical air and its subsequent mixing with drier air in the tropical midtroposphere contribute to this minimum. A similar minimum was also found by Pauluis and Mrowiec (2013). We also identify a branch where dry air with *l* ~ 0 kJ kg^{−1} loses DSE through radiative cooling (Peixoto and Oort 1992) and we speculate that it is associated with large-scale subsiding motions. This branch is bounded to the left by the limit of completely dry air (*l* = 0 kJ kg^{−1}). The third branch extends from *l* = 0 kJ kg^{−1}, *s* = 260 kJ kg^{−1} to *l* = 45 kJ kg^{−1}, *s* = 300 kJ kg^{−1} and is bounded by a curve representing saturated air at the surface. This curve shows the LH and DSE of saturated surface air, which can be derived from the Clausius–Clapeyron relationship as discussed in appendix C. In this branch, cold and dry lower-tropospheric air masses are warmed and moistened as they are mixed with near-surface air [a phenomenon similar to that described by Laliberté et al. (2012)]. However, following a streamline at about 400 Sv, we find that moistening of dry air also occurs at relatively high DSE values (~315 kJ kg^{−1}), which may not be near-surface air but rather mixing with air moistened by shallow convection. Given these processes and for the remainder of this work, we will refer to the first, second, and third branches as the precipitating, cooling, and moistening branches, respectively.

We now estimate the spatial structure of the hydrothermal cycle by projecting the streamfunction onto the time-averaged location of air masses. We select air masses between the −100 and −400 Sv streamlines of the hydrothermal cycle and divide them into sections similar to the hours on a clock (Döös et al. 2012) centered at the minimum value of the streamfunction min_{l,s}(*ψ*). Each section is colored such that 12 o’clock (green) is high DSE, 6 o’clock (dark blue) is low DSE, 3 o’clock (purple) is high LH, and 9 o’clock (light blue) is low LH. We then project these sections onto the time-averaged three-dimensional (longitude–latitude–pressure) LH and DSE values for ERA-Interim data (Fig. 3). The zonal mean of the projection (Fig. 3a) shows that the precipitating branch of the hydrothermal cycle occurs in the tropics, the cooling branch occurs at high altitudes, and the moistening branch occurs at low altitudes in midlatitudes and the subtropics. The hydrothermal cycle thus includes a meridional overturning cell in each hemisphere with ascent in the tropical regions and subsidence almost everywhere poleward of 30°. The reader should note that Fig. 3 is a projection of a thermodynamic streamfunction, which represents both mean and time-varying processes, into a time-averaged atmospheric state that does not take temporal variations into account. In particular, the projection neglects seasonal changes and cannot represent convection associated with cyclones in the midlatitude storm tracks (Pauluis et al. 2010), which play an important role in the meridional overturning circulation. Given this important caveat, these results are in qualitative agreement with the circulations that have been described using isentropic-mean or residual-mean frameworks (Karoly et al. 1997; Held and Schneider 1999; Juckes 2001).

A slice through the midtroposphere (the 570-hPa level is shown in Fig. 3b) shows that the precipitating branch is located over the equatorial western Pacific warm pool and the intertropical convergence zone (ITCZ), where we generally find precipitation and heating across the column from LH release (Kållberg et al. 2005). In the tropics, the cooling branch is found over the equatorial Atlantic and eastern Pacific—regions that are associated with radiative cooling. The zonal asymmetries in ascending and subsiding motions are associated with the Walker circulation and give a sense of its poleward extent. Note that, in midlatitudes and polar regions, the cooling branch is highly zonally symmetric.

As previously noted, the maximum mass flux of the meridional overturning circulation in moist isentropic coordinates is stronger than in that in dry isentropic coordinates (Fig. 1). We calculate the combined maximum mass flux of the meridional overturning circulation as the difference between the maximum and minimum value of the meridional overturning streamfunction marked as black dots in Fig. 1. The resulting mass fluxes are 191 Sv when using dry isentropic coordinates and 282 Sv when using moist isentropic coordinates, which is in close agreement with previous studies (Pauluis et al. 2008; Döös and Nilsson 2011). While moist isentropic coordinates take both DSE and LH variations into account, the maximum mass flux is still less than the hydrothermal cycle maxima of 428 Sv. Based on these results, we conclude that the meridional overturning component can at most account for 282 Sv (66%) of the total hydrothermal streamfunction, suggesting that zonally oriented flow features account for at least the remaining 146 Sv (34%). The projection onto spatial coordinates (Fig. 3b) shows stationary zonal asymmetries mainly in the tropics with ascent over the western tropical Pacific and subsidence to the east—both associated with the Walker circulation. We therefore suggest that the Walker circulation is a dominant feature of the zonally asymmetric component of the hydrothermal cycle.

## 4. Annual and interannual variability

The hydrothermal streamfunction projects the general circulation of the atmosphere onto a plane in thermodynamical (LH–DSE) coordinates. The resulting hydrothermal cycle thus includes the Hadley and Walker circulations, as well as midlatitude eddies, and presents an approximation of their combined mass fluxes across DSE surfaces. We study the temporal variability of the hydrothermal streamfunction by calculating it for every month for the period 1979–2009 using ERA-Interim data. We define the strength of the hydrothermal cycle as the average of the hydrothermal streamfunction maximum between 310 and 325 kJ kg^{−1} and denote this as Δ*ψ*. The values 310 and 325 kJ kg^{−1} were chosen since the amplitude of the hydrothermal streamfunction is somewhat constant over those values and because the maximum mass flux occurs somewhere in that range.

While Δ*ψ* is a measure of the maximum mass flux in the hydrothermal streamfunction, we also present a measure of the range in LH over which mass fluxes occur. We define the width of the hydrothermal cycle as the difference between the maximum and minimum LH value within the −50-Sv streamline at DSE *s* = 310 kJ kg^{−1}. The width is then a measure of the difference in LH between the precipitating branch and the cooling branch, and we denote it as ΔLH. The precipitating branch approximately follows a moist adiabat, and thus ΔLH is in fact a measure of changes in near-surface MSE. Furthermore, as DSE and LH of near-surface air are connected by the Clausius–Clapeyron relationship (Fig. 2c), ΔLH can be seen as a measure of changes in near-surface DSE and thus near-surface temperature. For ascending motions over the ocean, we can therefore expect a correlation between ΔLH and SST. We also note that ΔLH is fairly insensitive to the choice of DSE value at which it is evaluated.

We plot the resulting time series of Δ*ψ* as both seasonal [December–February (DJF), March–May (MAM), June–August (JJA), and September–November (SON)] and annual averages in Fig. 4. We find that the amplitude of the hydrothermal streamfunction undergoes large seasonal variations, but also that there is a clear interannual variability. By calculating a climatological annual cycle from the monthly data, we find that both Δ*ψ* and ΔLH vary semiannually with peaks in the boreal winter and austral winter (Fig. 5). The strength of the meridional overturning streamfunction in dry isentropic coordinates previously calculated (Fig. 1b) has a semiannual variation of similar amplitude (not shown) since the meridional overturning circulation is stronger in the winter hemispheres. We therefore suggest that a large part of the seasonal variations of the hydrothermal streamfunction is due to the variations in the strength of the meridional overturning. We further speculate that the semiannual variability in ΔLH is due to the fact that it is connected to the surface air temperature in regions of moist convection that are mostly in the tropics where surface air temperatures are known to have semiannual peaks.

Increased SSTs in the eastern equatorial Pacific that define positive ENSO phases are associated with a weakening of the zonal surface pressure gradient over the Pacific Ocean giving weaker zonal winds (Bjerknes 1969). We determine the monthly ENSO phase for the period 1979–2009 using the Niño-3.4 index (Barnston et al. 1997) from the extended reconstructed SST, version 3b (ERSST.v3b) dataset (Smith et al. 2008). Following Trenberth (1997), the Niño-3.4 index is defined as the SST averaged over the region 5°S–5°N, 170°–120°E. We then define seasonal positive or negative ENSO phases when the seasonally averaged Niño-3.4 index is ≥ 0.4 K or ≤ −0.4 K, respectively. In the same way, we denote a year as a positive or negative ENSO year when the annual-mean Niño-3.4 is ≥ 0.4 or ≤ −0.4 K, respectively.

We now correlate the seasonal and annual time series of Δ*ψ* and ΔLH with the Niño-3.4 index. We detrend the data by fitting a linear curve *y* = *at* + *b*, where *y* is either Δ*ψ* or ΔLH and *t* is time, and subtracting this curve from all time series. In ERA-Interim, Δ*ψ* has a linear trend of 0.46 Sv yr^{−1}, while ΔLH has no significant trend. We note that a positive trend in the Walker circulation strength was also indicated by L’Heureux et al. (2013) and Tokinaga et al. (2012) using reanalysis of recent years. The annual detrended time series are shown in Fig. 6. In the annual data we find correlations of *R* = −0.42 and *R* = 0.63 for Δ*ψ* and ΔLH, respectively, both significant above the 95% level (*p* = 0.02 and *p* < 0.01, respectively). Note that Δ*ψ* has been multiplied by −1 in Fig. 6 to show the anticorrelation to Niño-3.4 more clearly. From the correlation coefficients, we find that a linear fit to the Niño-3.4 time series explains *R*^{2} = (−0.42)^{2} = 18% and *R*^{2} = 0.66^{2} = 44% of the variance in Δ*ψ* and ΔLH in Fig. 6. For the seasonal data (not shown), the correlations are somewhat lower: *R* = −0.36 and *R* = 0.41 for Δ*ψ* and ΔLH, respectively. From Fig. 6, we see that the correlation with Δ*ψ* and ΔLH to Niño-3.4 is higher during strong El Niño and La Niña events, and that it is lower during weaker events or neutral ENSO phases. If all neutral ENSO phases are removed from the data (not shown), the correlations for the annual time series improve to *R* = −0.51 and *R* = 0.76. Furthermore, we note that the correlations are especially low in the period 1991–93, which could be because of changes in the radiative balance and atmospheric circulation due to the eruption of Mount Pinatubo and the decrease in global-mean surface air temperature that it caused (Bender et al. 2010). We also note a low correlation for Δ*ψ* around 1982 at the time of the eruption of El Chichón. As volcanic eruptions have been shown to sometimes mask the effects of ENSO events (Angell 1988) it is conceivable that they can also affect the hydrothermal cycle, albeit differently depending on the magnitude and location.

The negative correlation for Δ*ψ* and positive correlation for ΔLH with Niño-3.4 imply that the hydrothermal cycle weakens and widens during El Niño and strengthens and narrows during La Niña years. We show this by calculating composites of the hydrothermal streamfunctions during El Niño events and La Niña events in Fig. 7. The widening during El Niño implies that the precipitating branch of the cycle shifts toward a moist adiabat of higher MSE than during neutral ENSO phases. This is consistent with previous findings that El Niño events are associated with an overall increase in tropical surface air temperatures (Seager et al. 2003), which results in increased tropical DSE and, by Clausius–Clapeyron scaling, increased LH. Because the precipitating branch approximately follows moist adiabats (MSE conserved), an increase in surface MSE results in an increase in MSE throughout the branch. Our analysis also shows that the hydrothermal cycle weakens during El Niño events (Fig. 6). This implies less mass flux across DSE surfaces, and we suggest that this is mainly because of the weakened Walker circulation. Previous studies using spatial coordinates have found that El Niño events are also associated with stronger Hadley circulation (Oort and Yienger 1996; Lu et al. 2008). We also find correlations (not shown) indicating that this is true in our analysis. Since the hydrothermal cycle weakens, it suggests that the strengthening of the Hadley circulation is not enough to compensate for the weakening of the Walker circulation, leading to less mass fluxes across DSE surfaces globally, but distributed over a wider range of LH.

## 5. EC-Earth historical simulation

We now calculate the hydrothermal streamfunction from a historical simulation with the EC-Earth model and show the result in Fig. 8. Comparing to the results from ERA-Interim data (Fig. 2c), we find that the precipitating branch in EC-Earth has lower DSE and LH (and thus MSE) than ERA-Interim. As MSE is approximately conserved in this branch and LH is connected to DSE by the Clausius–Clapeyron relation, this indicates that EC-Earth has lower surface air temperatures in moist convective regions, which is consistent with findings of Hazeleger et al. (2011). Less DSE and LH would imply less DSE and LH fluxes in the precipitating branch of the hydrothermal cycle, but EC-Earth compensates this by having more ascent of moist air and thus a stronger hydrothermal cycle. The time-averaged hydrothermal streamfunction in EC-Earth has an amplitude of −514 Sv. As the EC-Earth historical simulation uses observed or reconstructed radiative forcing and atmospheric composition, the incoming and outgoing radiation at the top of the atmosphere should agree between ERA-Interim and EC-Earth. The upward flux of MSE should therefore be similar, suggesting that a hydrothermal cycle with low MSE in the precipitating branch will have a high amplitude. This argument is also consistent with the thermodynamic reasoning of Held and Soden (2006), where higher surface air temperatures lead to less convective mass flux. Furthermore, we find that the meridional overturning in EC-Earth (not shown) is stronger than in ERA-Interim and can account for 303 Sv (59%) of the amplitude of the hydrothermal streamfunction, suggesting that the zonal flows account for 211 Sv (41%).

The Niño-3.4 index for 1979–2009 indicates that only 17 of the 31 years of ERA-Interim data can be associated with ENSO events. To test our statistic on a longer time series, we repeat our analysis using the 155-yr historical simulation with EC-Earth. Sterl et al. (2012) found that EC-Earth can reproduce the ENSO phenomenon well, although the ENSO phases are somewhat weaker and less frequent than observations. The period 1850–2005 results in 60 years that can be associated with ENSO events. In the EC-Earth historical simulation, Δ*ψ* has a trend of −0.05 Sv yr^{−1} and ΔLH has a trend of 0.02 kJ kg^{−1} yr^{−1}, both statistically significant (*p* < 0.01) and indicating a weakening and widening of the cycle consistent with the thermodynamic scaling arguments of Held and Soden (2006) and Vecchi et al. (2006).

We also calculate the standard deviations of detrended annual-mean Niño-3.4, Δ*ψ*, and ΔLH using both ERA-Interim and EC-Earth. We find that standard deviations of all three metrics are lower in EC-Earth than in ERA-Interim, indicating that ENSO events are weaker in EC-Earth. Moreover, we identify 60 ENSO years of the total 155 years in the EC-Earth historical simulation, while there are 17 ENSO years in 31 years of ERA-Interim data. This indicates that ENSO events are also less frequent in EC-Earth than recent observations.

As done for the ERA-Interim data, we calculate a mean annual cycle for Δ*ψ* and ΔLH in EC-Earth. The results, shown in Fig. 9, are similar to the results for ERA-Interim with a semiannual variation with both Δ*ψ* and ΔLH peaking in January and July. The amplitudes of the annual cycles of Δ*ψ* and ΔLH are slightly higher in EC-Earth than in ERA-Interim, which is likely due to the meridional overturning (not shown) being slightly stronger, as mentioned earlier.

We show the detrended, annual-mean Niño-3.4 index together with −Δ*ψ* and ΔLH from the EC-Earth data in Fig. 10. Note that, just as in Fig. 6, Δ*ψ* has been multiplied by −1. We find that Δ*ψ* and ΔLH are correlated with the Niño-3.4 index by *R* = −0.81 and *R* = 0.45, respectively, both significant at the 99% level (both *p* < 0.01). Hence, just as in the 30 years of ERA-Interim data, the hydrothermal cycle widens and weakens during El Niño events and narrows and intensifies during La Niña events. Composites of the hydrothermal streamfunction for El Niño and La Ninã events in EC-Earth (not shown) are similar to those shown for ERA-Interim in Fig. 7. Compared to the correlations obtained in ERA-Interim the correlation is higher for Δ*ψ* but lower for ΔLH. Hazeleger et al. (2011) showed that the SST variations in EC-Earth associated with ENSO events are more confined to the equator than in observations, which we speculate could partially explain the lower correlation between ΔLH and Niño-3.4. Both datasets include observed volcanic gases, and aerosols therefore include the effects of eruptions. This is particularly discernible in the ERA-Interim time series (Fig. 6) in which both Δ*ψ* and ΔLH decrease in the year 1992 when the Niño-3.4 index increases. Such events have a larger impact on the correlations using ERA-Interim data since that dataset is much shorter than the EC-Earth simulation.

## 6. Discussion

We have shown (Figs. 2c and 8) that global mass fluxes form a single hydrothermal cycle in latent heat (LH) versus dry static energy (DSE) coordinates. Results from ERA-Interim data during 1979–2009 show that the cycle has an amplitude of 428 Sv. This thermodynamic cycle describes three processes: heating and moistening of near-surface air masses following the Clausius–Clapeyron relationship, heating and drying of air masses along moist adiabats in the tropics, and radiative cooling of air masses at the top of the atmosphere and at high latitudes. A projection onto spatial coordinates (Fig. 3) shows that the hydrothermal streamfunction captures both the Walker circulation and the meridional overturning circulation. The hydrothermal streamfunction thus includes both the zonal and meridional overturning circulations owing to their exchanges of LH and DSE and does not differentiate between them because of their orientation in spatial coordinates. Less than 70% of the total mass flux in the hydrothermal streamfunction can be explained by the meridional circulation in moist isentropic coordinates (Pauluis et al. 2008; Döös and Nilsson 2011). Because they play a central role in the global heat budget, we suggested that the unexplained 30% of the total mass fluxes is associated with zonally oriented flows such as the Walker circulation. From our analysis based on ERA-Interim data, this would lead to about 148 Sv attributable to the Walker circulation, which is in line with an estimate of more than 100 Sv from other studies (Trenberth et al. 2000).

The hydrothermal cycle gives a unified description of the heat fluxes within the Hadley and Walker circulations and within synoptic eddies in midlatitudes. By calculating it for both ERA-Interim and historical simulations from the EC-Earth climate model, we find that the hydrothermal cycle is narrower and stronger in EC-Earth (Fig. 8). We attribute this to the colder global-mean surface air temperature in EC-Earth leading to less surface DSE and LH (and thus MSE), which is compensated by more ascending motion such that the upward MSE flux is similar to that in ERA-Interim. We speculate that examining several climate models would show that those with colder surface air temperatures than observations would have a stronger and narrower hydrothermal cycle than those with warmer temperatures. The dependence on the surface temperatures can be tested by calculating the hydrothermal streamfunction for different ENSO phases. We have shown in both ERA-Interim and EC-Earth historical simulations that the hydrothermal cycle widens and weakens with increased SSTs and vice versa. On interannual time scales, variations in the Niño-3.4 index from the ERSST.v3b dataset (Smith et al. 2008) explain 18% and 44% of the variations in strength and width of the hydrothermal streamfunction, respectively. From 155-yr historical simulations with the EC-Earth climate model, the numbers are 66% and 20%. This shows that the results obtained with ERA-Interim are verified in longer simulated datasets. It also gives some evidence that EC-Earth is able to simulate both ENSO and its effects on the atmospheric general circulation, although ENSO events are somewhat weaker and less frequent than in observations.

We have defined the hydrothermal streamfunction using globally integrated advection of DSE and shown that it is related to advection of LH. The total mass flux across DSE surfaces is given by the globally integrated tendency . Calculating a streamfunction, *ψ*^{tot}, based on this tendency yields results very similar to those obtained for the hydrothermal streamfunction *ψ*. Repeating the analysis in this paper with *ψ*^{tot} instead of *ψ*, the correlations of Δ*ψ* and ΔLH to ENSO become slightly smaller in both ERA-Interim and EC-Earth but are still robust and statistically significant. The composites (Fig. 7) do not change. Hence, the hydrothermal streamfunction is a good approximation to the total mass fluxes across DSE surfaces in the data.

A caveat of the analysis in this paper is the absence of resolved convection in ERA-Interim and EC-Earth, as moist convection is an important process for the conversion of LH to DSE. Following Simmons and Burridge (1981), the vertical velocity *ω* is calculated from the horizontal mass fluxes and changes in surface pressure. As the vertical mass fluxes must balance the divergence of the horizontal mass fluxes, which in turn are calculated partly using assimilated observations, we argue that the vertical mass flux amount in spatial coordinates is realistic. Convection occurs mostly on subgrid scales so we expect both ERA-Interim and EC-Earth to underestimate the correlations between *ω* and DSE and LH. As a consequence, it is thus possible that some moist convection occurs at somewhat higher LH and DSE than the precipitating branch in Figs. 2c and 8. In a recent study Pauluis and Mrowiec (2013) show convection occurring at *θ*_{e} values approximately 5 K higher than the mean *θ*_{e} profile. This could have the effect that in our analysis some ascent occurs at the same LH and DSE as nearby descent, which would make us underestimate the mass flux to some extent.

Future applications of the analysis presented here should include a study of the hydrothermal cycle response to anthropogenic forcing in climate models. Results from an intercomparison of climate models (Meehl et al. 2007) suggest that the pattern of global warming is similar to that of an El Niño event, in which case we would expect a widening and weakening of the hydrothermal cycle. The hydrothermal cycle would reveal global changes in mass fluxes as well as exchanges of MSE, DSE, and LH, and thereby changes in precipitation and evaporation.

In concurrent work (F. Laliberté et al. 2013, unpublished manuscript), we have applied a similar methodology to study how irreversible entropy production associated with moist processes constrains the maximum work output of the atmospheric heat engine in a climate change scenario. The results indicate that a wider hydrothermal cycle size is associated with a reduction in maximum work output for the global atmospheric heat engine.

## Acknowledgments

J. Kjellsson and K. Döös acknowledge the Bolin Centre for Climate Research and its associated Climate Research School, both at Stockholm University, Sweden. F. B. Laliberté was supported by the Natural Sciences and Engineering Research Council (NSERC) G8 Research Initiative Grant “ExArch: Climate analytics on distributed exascale data archives” as well as by an NSERC Strategic Project Grant. J. D. Zika would like to acknowledge support from the National Environment Research Council (United Kingdom) and the ARC Centre of Excellence for Climate System Science. All authors extend great thanks to SMHI and Laurent Brodeau for supplying the EC-Earth output and ECMWF for ERA-Interim data. The ERSST.v3b data (Niño-3.4 index) was downloaded from NOAA Climate Prediction Center (http://www.cpc.ncep.noaa.gov/data/indices/). All computations, analyses, and storage of data have been performed at the Vagn and Triolith supercomputer systems maintained by the National Supercomputer Centre (NSC) at Linköping University, Sweden. The authors wish to thank Olivier Pauluis and one anonymous reviewer whose comments greatly improved the paper. We also acknowledge Harry Bryden, Lawrence Mudryk, Trevor J. McDougall, Rune Grand Graversen, Johan Nilsson, and Maxime Ballarotta for reading early versions of the paper and for fruitful discussions.

### APPENDIX A

#### Meridional Overturning

We calculate meridional overturning streamfunctions from the ERA-Interim (Fig. 1) and EC-Earth (not shown) datasets using pressure *p*, dry static energy *s*, and moist static energy *h* as vertical coordinates, following Döös and Nilsson (2011):

where *χ* stands for vertical coordinates *p*, *s*, or *h*. The unit step function *μ* is defined as

As in the hydrothermal streamfunction, *t*_{0} is 1 January 1979 and *t*_{1} is 31 December 2009. We find that amplitudes and shapes of the streamfunctions were consistent with other studies (Pauluis et al. 2008; Döös and Nilsson 2011). We estimate the total mass flux by the meridional circulation in moist isentropic coordinates as max[*ψ*(*y*, *h*)] − min[*ψ*(*y*, *h*)].

### APPENDIX B

#### Discretization

We discretize Eq. (5) as follows and illustrate the setup in Fig. B1. We denote the zonal, meridional, and vertical mass fluxes as *U*, *V*, and *W* (kg s^{−1}), respectively. The zonal mass flux through the interface between grid box (*i*, *j*, *k*) and (*i* + l, *j*, *k*) is *U*_{i,j,k}. If *s*_{i+1,j,k} > *s*_{i,j,k}, then we find the smallest positive integer *n* such that *s*_{i,j,k} + *n*Δ*s* > *s*_{i+1,j,k} > *s*_{i,j,k}. This indicates a positive mass flux through the DSE isosurfaces *s* and *s* + Δ*s* at latent heat *l*_{i,j,k}. We then record a mass flux of size *U*_{i,j,k} at coordinates (*l*, *s*), (*l*, *s* + Δ*s*), …, (*l*, *s* + *n*Δ*s*). If *s*_{i+1,j,k} < *s*_{i,j,k}, then we find the smallest positive integer *n* such that *s*_{i,j,k} − *n*Δ*s* < *s*_{i+1,j,k} < *s*_{i,j,k}. We then record a mass flux of size −*U*_{i,j,k} at coordinates (*l*, *s*), (*l*, *s* − Δ*s*), …, (*l*, *s* − *n*Δ*s*). Contributions from meridional and vertical mass fluxes are added in a similar fashion. We then obtain the hydrothermal streamfunction by integrating the net mass fluxes *F* in hydrothermal coordinates:

where *m* and *r* are discrete coordinates in hydrothermal space such that *l*_{m} = *l*_{0} + *m*Δ*l* and *s*_{r} = *s*_{0} + *r*Δ*s*. The index *n* is the time step and *N* is the total number of time steps over which the streamfunction is integrated.

The local Eulerian part of Eq. (1) *ψ*^{loc} is discretized in a similar way. However, instead of using mass fluxes from one grid box to an adjacent one, we use the mass of a grid box and the local rate of change in DSE from time step *n* to *n* + 1. Hence, if *s*_{i,j,k,n+1} > *s*_{i,j,k,n}, then we find the smallest positive integer *a* such that *s*_{i,j,k} + *a*Δ*s* > *s*_{i,j,k,n+1} > *s*_{i,j,k,n}. The mass flux from one time step to another is *dM*/Δ*t* (i.e., the mass of the grid box and the length of the time step). As done for the advective part above, we then record a mass flux from *s*_{i,j,k,n} to *s*_{i,j,k,n+1} and integrate the mass fluxes.

### APPENDIX C

#### Clausius–Clapeyron

In Figs. 2a–c and 8, we show the saturated moist static energy (MSE) for surface air as a dashed line. We calculated it as follows. For a given dry static energy (DSE) *s*, we obtain the maximum temperature *T* by adiabatic descent to *z* = 0, in which case *T* = *s*/*c*_{p}. For a given temperature, we find the saturation vapor pressure *e*_{s} by integrating the Clausius–Clapeyron equation (Wallace and Hobbs 2006); *e*_{s} is related to the saturation mixing ratio *q*_{s} by the approximation *q*_{s} ≈ 0.622*e*_{s}/1013 hPa if we assume the ambient pressure is taken at the surface pressure *p*_{s} = 1013 hPa. Thus, for saturated air close to the surface, LH is a function of DSE. With the approximations of Wallace and Hobbs (2006), we thus obtain

where *e*_{0} = *e*_{s}(273 K) = 6.11 hPa. The surface saturated MSE is then *h* = *l* + *s*, shown as a dashed line in Figs. 2a–c.

## REFERENCES

*Climate Dyn.,*

**39,**2611–2629,

*Nat. Climatic Change,*

**3,**571–576.

*Climate Change 2007: The Physical Science Basis,*S. Solomon et al., Eds., Cambridge University Press, 747–845.

*J. Atmos. Sci.,*

**70,**3673–3688.

**43,**2095–2112.