1 Introduction

The general circulation of the atmosphere is driven by spatial inhomogeneities in input of solar radiation and fluxes to/from the ocean. The atmosphere responds to these inhomogeneities by transporting latent heat (LH) and/or dry static energy (DSE) from the regions with an excess in energy (e.g. the tropics) to regions of deficit in energy (e.g. polar regions). If the atmosphere’s radiative properties or the fluxes to/from the ocean change the atmospheric circulation might change as well. For instance, the tropical overturning circulations, i.e. Hadley and Walker cells, are projected to weaken as a result of anthropogenic emissions and associated surface warming (Vecchi et al. 2006; Vecchi and Soden 2007; Lu et al. 2008; Seager et al. 2010). Changes have also been projected for the Indian Monsoon (Bollasina et al. 2011; Ming and Ramaswamy 2011) and midlatitude dry and moist isentropic circulations (Laliberté and Pauluis 2010; Wu and Pauluis 2013).

Man-made emissions of greenhouse gases (GHG) and aerosols result in positive trends in surface air temperature over most of the Earth in the twentieth and twenty-first centuries (IPCC 2013). Increased temperature leads to increased specific humidity following the Clausius–Clapeyron relation [cf. (Wallace and Hobbs 2006)],

$$\begin{aligned} \frac{d q_s}{d T}&= \frac{q_s L_v}{R_v T^2}, \end{aligned}$$
(1)
$$\begin{aligned} \frac{{\varDelta } q_s}{q_s}&= \alpha (T) {\varDelta } T = \frac{L_v}{R_v T^2} {\varDelta } T, \end{aligned}$$
(2)

where \(q_s\) is saturation specific humidity, \(T\) is temperature, \(R_v = 461.5\hbox { J}\hbox { kg}^{-1}\hbox { K}^{-1}\) is the gas constant for water vapour and \(L_v(T)\) is the latent heat of condensation at temperature \(T\). It is assumed throughout this paper that \(L_v = 2.5 \cdot 10^6\hbox { J}\hbox { kg}^{-1}\) and the total pressure, \(p\), is much larger than the water vapour pressure. Then, \(T = 273\hbox { K}\) yields \(\alpha \approx 0.07 \hbox { K}^{-1}\), i.e. saturation specific humidity increases by about 7 % per degree of warming (Eq. 2). Held and Soden (2006) found that the relative humidity is largely unchanged with global warming, and that specific humidity, \(q\), averaged over the tropics (\(30^{\circ } \hbox {N}\)\(30^{\circ }\hbox {S}\)) increases by about \(7.5\,\%\,\text {K}^{-1}\). They further proposed a thermodynamic scaling where precipitation, \(P\), can be approximated as a product of convective mass flux, \(M_c\), and specific humidity, \(q\), near the Earth’s surface,

$$\begin{aligned} P = M_c q \, \Rightarrow \, \frac{{\varDelta } P}{P} = \frac{{\varDelta } M_c}{M_c} + 0.07 {\varDelta } T, \end{aligned}$$
(3)

where \({\varDelta }\) represents the difference between 1980–2000 and 2080–2100 and \({\varDelta } q / q = 0.07 {\varDelta } T\) is used. Convective mass flux, \(M_c\), can however be impractical as a diagnostic since it is not a resolved variable in global climate models and is often stored at relatively low temporal resolution.

With an unchanged global atmospheric circulation, i.e. \({\varDelta } M_c = 0\), global warming would cause upward flux of latent heat and thus precipitation and latent heat release to all increase by 7 % K−1. However, by considering radiative balance Stephens and Hu (2010) argued that precipitation is constrained by radiative cooling and can thus not increase as rapidly as surface specific humidity. Several climate-model studies have found that precipitation increases by about \(1-3\,\%\,\text {K}^{-1}\) (Held and Soden 2006; Vecchi and Soden 2007; Stephens and Hu 2010) with differences mostly arising from aerosol treatments and cloud feedbacks. It thus follows from Eq. 3 that the upward mass flux and hence the global overturning circulation should weaken by 4–6 % K−1.

When studying the global atmospheric circulation it is important to consider both its zonal and meridional components. Some studies of the tropical circulations have hitherto diagnosed changes in its strength by studying e.g. sea-level pressure (SLP) and sea-surface temperature (SST) gradients (Vecchi et al. 2006; Tokinaga et al. 2012; Bayr and Dommenget 2013). However, SST gradients and SLP gradients in certain regions only approximate the strength of one part of a circulation (e.g. the surface branch of the Walker cell).

Another common metric is vertical mass flux at \(500 \hbox { hPa}\), \(\omega _{500}\) (Vecchi and Soden 2007; Merlis and Schneider 2011; Bony et al. 2013). However, upward mass fluxes include both adiabatic and diabatic motions, i.e. only a part is associated with air-mass transformations. Furthermore, data must have daily or sub-daily resolution to capture fluxes by midlatitude eddies correctly. Temporal resolution is also important in the tropics as there are significant differences between 6-h and monthly data when calculating mass and energy fluxes by the Hadley cells (not shown).

This study diagnoses future changes of the global atmospheric circulation in a set of coordinated climate-model simulations in the Coupled Model Inter-comparison Project, phase 5 (CMIP5). CMIP5 simulations were carried out for the United Nations Intergovernmental Panel on Climate Change (IPCC) 5th Assessment Report (AR5). The simulations are analysed partly using a recently developed thermodynamic stream function that captures both the zonal and meridional overturning circulations as a single cycle (Kjellsson et al. 2014; Laliberté et al. 2014). The stream function resides in LH–DSE space and therefore clearly highlights precipitation, radiative cooling, diabatic heating and moistening. Changes of the global circulation are discussed in light of these processes, mass and energy transports are calculated, and model results are compared to simple theoretical predictions. Changes in the meridional overturning circulation and the poleward energy transport are also studied separately to understand the relative changes in zonal and meridional overturning circulations.

2 Data

Table 1 The resolution of CMIP5 data output and the ensemble members used

Horizontal winds, \(u,v \,[\text {m}\hbox { s}^{-1}]\), temperature, \(T~[\text {K}]\), specific humidity, \(q \,[\text {kg}\hbox { kg}^{-1}]\), and surface pressure, \(p_s\,[\text {Pa}]\) are obtained from 11 coupled climate models and ERA-Interim reanalysis as listed in Table 1. All data are retrieved from the CMIP5 archive at the Earth System Grid FederationFootnote 1 except EC-Earth (downloaded directly from SMHI storage, Sweden), GFDL-CM3 (downloaded directly from GFDLFootnote 2), CCSM4 (downloaded from Earth System GridFootnote 3) and ERA-Interim (downloaded from ECMWFFootnote 4). All data have 6 h output frequency and are on model levels where the pressure at level \(k\) is \(p_k = a_k + b_k p_s\) or \(p_k = a_k p_0 + b_k p_s\). Pressure, \(p~[\text {Pa}]\), geopotential height, \(z\,[\text {m}]\), and vertical velocity, \(\omega \,[\text {Pa}\hbox { s}^{-1}]\) are calculated every 6 h on each model level using the hydrostatic approximation and mass continuity following Simmons and Burridge (1981) and Kjellsson and Döös (2012). Prior to analysis all data are re-gridded from their original resolution (Table 1) to a \(2.5^{\circ } \times 2.5^{\circ }\) C-grid (Mesinger and Arakawa 1976) using bilinear interpolation. For the BCC and CanESM models this means moving to a somewhat finer grid than the output, but throughout the analysis in this paper nothing suggests that this is a problem. See Table 2 for more details about the models.

Climate-model simulations in the CMIP5 models begin with a long “spin-up” period (typically \(>\)1,000 years) where radiative forcing and concentrations of GHG and aerosols are kept fixed at pre-industrial (i.e. mid-nineteenth century) levels, until the models are in a relatively steady climate. The models are then run with observed or reconstructed radiative forcing and GHG and aerosol concentrations typically starting in 1850, but for some models a few years later. These simulations of the past \(\sim \)150 years are named historical simulations. All historical simulations end in 2005. The RCP8.5 scenario is a scenario where anthropogenic GHG emissions are such as to cause a radiative forcing of \(8.5\hbox { W}\hbox { m}^{-2}\). It implies more than a tripling in equivalent \(\hbox {CO}_{2}\) concentration from 2000 to 2100. All RCP8.5 simulations start in 2006 and run to at least 2100. All CMIP5 simulations are described in detail by Taylor et al. (2012) and various RCP emission scenarios are described by Meinshausen et al. (2011). The effects of anthropogenic emissions in the RCP8.5 scenario are here assessed by calculating the differences between model-simulations of the late twentieth century (1980–2000) and the late twenty-first century in the RCP8.5 emission scenario (2080–2100). These simulations will be denoted as 20C and 21C respectively.

3 Methods

Convection is not resolved by any of the models in Table 1 and it is not parameterised when analysing the data since it is already parameterised by the models. Vertical velocity, \(\omega \), is calculated from mass continuity on grid-box scales and this analysis will therefore not capture sub-grid scale vertical fluxes of moisture and heat associated with convective plumes. Furthermore, vertical mass fluxes associated with up- and downdrafts may, to some extent, cancel when averaged over a grid box leading to an underestimation of the vertical mass fluxes in thermodynamic space.

LH, \(l = L_v q\), DSE, \(s = c_p T + g z\), and moist static energy (MSE), \(h = l + s\), are calculated using \(L_v = 2500\hbox { kJ}\hbox { kg}^{-1}\) as the latent heat of vaporisation (assumed constant), \(c_p = 1.004\hbox { kJ}\hbox { kg}^{-1}\hbox { K}^{-1}\) as the specific heat capacity of dry air at constant pressure and \(g = 9.81\hbox { m}\hbox { s}^{-2}\) as gravity. Note that DSE is similar to potential temperature, \(\theta \), which is conserved for adiabatic motions, and MSE is similar to equivalent potential temperature, \(\theta _e\), which is conserved for moist adiabatic motions.

Fig. 1
figure 1

Meridional overturning stream functions calculated from ERA-interim 1980–2000 using Eq. 5 with pressure (a), LH (b), DSE (c) and MSE (d) as vertical coordinates respectively. Units in Sverdrup (\(1\hbox { Sv} = 10^9\hbox { kg}\hbox { s}^{-1}\)). Positive values (solid stream lines) indicate clockwise circulation. Black dots indicate maximum and minimum values of the stream functions. Lines in ad show zonal mean pressure (a), LH (b), DSE (c), and MSE (d) in the lowermost model level respectively. e Shows meridional fluxes of LH (blue), DSE (red) and MSE (black) calculated from Eq. 6

Figure 1a shows the classical meridional overturning stream function, \({\varPsi }(y,p)\), for ERA-Interim reanalysis. It represents the total mass flux above the pressure level \(p\) at latitude \(y\) (cf. Peixoto and Oort (1992)) and is defined as

$$\begin{aligned} {\varPsi }(y,p) = \frac{1}{t_1 - t_0} \int _{t_0}^{t_1} \oint _x \int _0^p g^{-1} v(x,y,p^{\prime},t) \, dp^{\prime} \, dx \, dt. \end{aligned}$$
(4)

In practice the meridional velocities are sorted according to their pressure and then integrated from 0 to \(p\) to give \({\varPsi }(y,p)\). The result (Fig. 1a) is a three-cell structure in each hemisphere with the Hadley, Ferrel and Polar cells. The Polar cells have an amplitude of \(\sim 5\hbox { Sv}\) and do consequently not appear clearly in Fig. 1a.

Following Döös and Nilsson (2011) the definition of a meridional overturning stream function with a generalised vertical coordinate is

$$\begin{aligned} {\varPsi }(y,\chi ) = \frac{1}{t_1 - t_0} \int _{t_0}^{t_1} \oint _x \int _0^{p_s} \mu [\chi - \chi ^{\prime}(x,y,p^{\prime},t)] g^{-1} v(x,y,p^{\prime},t) ~dp^{\prime} ~dx ~dt, \end{aligned}$$
(5)

where \(\chi \) is any tracer variable (e.g. temperature or specific humidity). \(\mu [x]\) is a Heaviside function where \(\mu = 1\) for \(x \ge 0\) and \(\mu = 0\) otherwise. With this definition \({\varPsi }(y,\chi )\) only includes mass fluxes with \(\chi ^{\prime}(x,y,p^{\prime},t) \le \chi \). Note that \(\chi \) need not be monotonic with height. This allows for calculating meridional overturning stream functions using LH, DSE, and MSE as vertical coordinates (Fig. 1b–d) as also done by e.g. Czaja and Marshall (2006), Pauluis et al. (2010) and Döös and Nilsson (2011). Meridional fluxes of LH, DSE and MSE (Fig. 1e) are calculated from the meridional overturning stream functions as

$$\begin{aligned} F_\chi (y) = - \int _{\chi _{\text {min}}(y)}^{\chi _{\text {max}}(y)} \frac{\partial {\varPsi }(y,\chi )}{\partial \chi } \chi \, d \chi = \int _{\chi _{\text {min}}(y)}^{\chi _{\text {max}}(y)} ( {\varPsi }(y,\chi ) - {\varPsi }_0) \, d \chi , \end{aligned}$$
(6)

where \(\chi \) is LH, DSE, or MSE. The limits \(\chi _{\text {min}}(y)\) and \(\chi _{\text {max}}(y)\) are the values of \(\chi \) where \({\varPsi }(y,\chi ) = {\varPsi }_0\). It is assumed that \({\varPsi }_0 = 0\hbox { Sv}\).

The Hadley cell in the tropics comprises an equatorward branch in the lower troposphere and a poleward branch in the upper troposphere. The equatorward branch has high LH and low DSE and the poleward branch has low LH and high DSE (Townsend and Johnson 1985; Held and Schneider 1999; Pauluis et al. 2010). Mass fluxes in the midlatitudes are mainly associated with transient eddies where the equatorward and poleward mass transports occur at similar pressure but the latter generally have higher LH and DSE than the former (McIntosh and McDougall 1996; Döös and Nilsson 2011; Laliberté et al. 2012; Kjellsson and Döös 2012). Thus the meridional overturning circulation depends on the choice of vertical coordinate. In LH coordinates (\({\varPsi }(y,l)\), Fig. 1b) the result is two cells in each hemisphere transporting LH from the subtropics to the tropics and midlatitudes (Fig. 1e). Moreover, the midlatitude cells in \({\varPsi }(y,l)\) are stronger than the Ferrel cells in pressure coordinates. In DSE coordinates (\({\varPsi }(y,s)\), Fig. 1c) the circulation comprises only one cell in each hemisphere with two maxima each as the tropical and midlatitude cells merge. The transport of DSE is from the tropics and poleward at almost all latitudes (Fig. 1e). As a final result, in MSE coordinates (\({\varPsi }(y,h)\), Fig. 1d) the result is two cells centered in the midlatitudes where the poleward MSE fluxes peak (Fig. 1e).

The hydrothermal stream function was defined by Kjellsson et al. (2014) using a method originally developed for the ocean by Zika et al. (2012) and Döös et al. (2012). It resides in LH–DSE space and therefore has no spatial coordinate. It represents mass fluxes across a DSE surface at or below a certain LH value, or mass fluxes across a LH surface at or below a certain DSE value. The hydrothermal stream function is the global integral of all atmospheric motions where LH or DSE change, irrespectively of direction. The hydrothermal stream function is thus able to include both zonal-mean overturning circulations as in Fig. 1a–d as well as zonally asymmetric features such as the Walker circulation. This is somewhat similar to a recent study where Pauluis and Mrowiec (2013) defined an “isentropic” stream function in \(z\)\(\theta _e\) space (similar to \(z\)–MSE space).

The tendency of DSE is defined as \(ds / dt = \partial s / \partial t + \varvec{v} \cdot \nabla s = \dot{s}\). In LH–DSE coordinates \(\dot{s}\) can be written as

$$\begin{aligned} \dot{S}(l,s) = \frac{1}{t_1 - t_0} \int _{t_0}^{t_1} \int _{\varOmega } \delta [s - s^{\prime}(x,y,p,t)] \delta [l - l^{\prime}(x,y,p,t)] \dot{s^{\prime}}(x,y,p,t) \, dM \, dt, \end{aligned}$$
(7)

where \({\varOmega }\) is the whole atmosphere, \(dM\) is a mass element, and primes denote variables in \(x,y,p\) space. Throughout this paper, \(\dot{l}\) and \(\dot{s}\) are LH and DSE tendencies in \((x,y,p)\) space while \(\dot{L}\) and \(\dot{S}\) are tendencies in \((l,s)\) space. Moreover, \(\delta [x]\) is the Dirac function for which \(\int _{-\infty }^{\infty } \delta [x] dx = 1\) and \(\delta [x] = 0\) for \(x \ne 0\).

The unit of \(\delta [s-s^{\prime}]\) and \(\delta [l-l^{\prime}]\) is \(\text {kg}\hbox { kJ}^{-1}\) and \(\dot{s^{\prime}}\) has \(\text {kJ}\hbox { kg}^{-1}\hbox { s}^{-1}\), thus \(\dot{S}\) has units \(\text {kg}\hbox { s}^{-1}~(\text {kJ}\hbox { kg})^{-1}\). The LH tendency can be integrated similarly and denoted \(\dot{L}\). The hydrothermal stream function is defined as

$$\begin{aligned} \frac{\partial {\varPsi }(l,s)}{\partial l} = \dot{S}, \qquad \frac{\partial {\varPsi }(l,s)}{\partial s} = - \dot{L} \end{aligned}$$
(8)

where \({\varPsi }(l,s)\) has units of \(\text {kg}\hbox { s}^{-1}\). The two definitions above yield identical results if the tendencies are non-divergent, i.e. \(\partial _l \dot{L} + \partial _s \dot{S} = 0\). This only holds if the system conserves mass, LH and DSE perfectly which no CMIP5 model does. Furthermore, over a 20-year period decadal variabilities or global warming can give trends in LH and DSE. However, in this study this introduces only small errors and \(\partial _l \dot{L} + \partial _s \dot{S} \approx 0\) is found to be a fair approximation. The error in the hydrothermal stream function, i.e. mass fluxes outside closed stream lines, is found to be <5 % of the maximum amplitude. Mass continuity in LH–DSE space is discussed further in the Appendix. In this study the stream function \({\varPsi }(l,s)\) is calculated from \(\dot{S}\), but the results obtained from \(\dot{L}\) have been verified to be almost identical. The total tendencies \(\dot{S}\) and \(\dot{L}\) comprise both a local time derivative and an advective term. This definition differs slightly from Kjellsson et al. (2014) where an approximation of \(\dot{s}\) was done as \(\dot{s} = ds / dt = \partial s / \partial t + \varvec{v} \cdot \nabla s \approx \varvec{v} \cdot \nabla s\). This was justified since the contribution from \(\partial s / \partial t\) to \({\varPsi }(l,s)\) was found to be very small. A similar approximation is less valid in the ocean (Groeskamp et al. 2014). The full tendency as given in Eq. 7 is used throughout this paper. Hence,

$$\begin{aligned} {\varPsi }(l,s) = \frac{1}{t_1 - t_0} \int _{t_0}^{t_1} \int _{\varOmega } \mu [ l -l^{\prime} ] \delta [s - s^{\prime}] g^{-1} \dot{s^{\prime}} \, dM \, dt. \end{aligned}$$
(9)

4 Results

Fig. 2
figure 2

Hydrothermal stream function, \({\varPsi }(l,s)\), calculated from Eq. 9 for historical simulations (1980–2000) in the CMIP5 models and ERA-interim reanalysis shown in Table 1. In each figure, a thick dashed line shows the tropical mean (\(15^{\circ }\hbox {S}\)\(15^{\circ }\hbox {N}\)) profile and a thinner dashed line shows MSE for saturated air at \(z = 0 \hbox { m}\) and \(p = 1{,}013\hbox { hPa}\). The direction of the mean flow is indicated by the arrows in the bottom right figure

The hydrothermal stream function is calculated for the years 1980–2000 in the historical simulations and ERA-Interim reanalysis in Table 1 (Fig. 2). The unit of the stream function is Sverdrup (Sv) where \(1\hbox { Sv} = 10^9 \hbox { kg}\hbox { s}^{-1}\). For all models the result is a single anti-clockwise cycle as indicated by the arrows in the bottom right panel in Fig. 2. The result is very similar to Kjellsson et al. (2014) who used ERA-Interim and EC-Earth but longer periods. For all models three distinct branches are discernible along the \(-30\hbox { Sv}\) stream line. One branch extends from \((l \approx 40\hbox { kJ}\hbox { kg}^{-1}, s \approx 300\hbox { kJ}\hbox { kg}^{-1})\) to \((l = 0\hbox { kJ}\hbox { kg}^{-1}, s \approx 340\hbox { kJ}\hbox { kg}^{-1})\) with stream lines approximately parallel to the tropical-mean (\(15^{\circ }\hbox {S}\)\(15^{\circ }\hbox {N}\)) moist adiabat where MSE is almost conserved. This indicates conversion of LH to DSE by condensation of water vapour associated with cloud formation and precipitation. DSE tendencies are positive to the right of the mean tropical profile and negative to the left. This suggests that moist ascent is concentrated where MSE is higher than the tropical mean. From \((l = 0\hbox { kJ}\hbox { kg}^{-1}, s \approx 340\hbox { kJ}\hbox { kg}^{-1})\) mass is transported toward lower DSE along \(l = 0\hbox { kJ}\hbox { kg}^{-1}\), indicating radiative cooling of dry air. Some air descends in the tropics where DSE is seldom lower than \(300\hbox { kJ}\hbox { kg}^{-1}\), and some air is transported all the way to \((l = 0\hbox { kJ}\hbox { kg}^{-1}, s \approx 260\hbox { kJ}\hbox { kg}^{-1})\) where the low DSE indicates cold air near the polar regions. As the air cools and descends it is moistened either near the surface or by mixing with moisture transported upward by e.g. shallow convection. This increases LH. For lower tropospheric air, surface heating increases DSE which increases \(q_s\) (Eq. 1) so that the flow follows a Clausius–Clapeyron relation. The Clausius–Clapeyron lines in Fig. 2 are calculated from Eq. 1 using \(e_s(273\hbox { K}) = 611\hbox { Pa}\) (cf. Wallace and Hobbs (2006)) and \(T = c_p/s\) (i.e. \(z=0\hbox { m}\) and p = 1,013 hPa). This is identical to calculations by Kjellsson et al. (2014). The Clausius–Clapeyron lines thus represent saturated surface air. Testing different relative humidities shows that the actual stream lines follow a line consistent with ~80 % relative humidity. A global-mean relative humidity of ~80 % was found by e.g. Held and Soden (2006) and is not projected to change significantly with global warming. Following Kjellsson et al. (2014), the three branches are henceforth referred to as “precipitating” branch, “radiative cooling” branch and “moistening” branch after the main processes they represent.

Some differences in amplitude and shape of the hydrothermal stream function can be noted between the models (Fig. 2 and Table 1). Some models have a distinct “dip” in the precipitating branch (e.g. CCSM4, NorESM1-M, MIROC5) while some have nearly straight stream lines (e.g. GFDL-CM3, BCC-CSM1.1, IPSL-CM5B-LR). The “dip” can be explained by the MSE minimum in the tropical mid-troposphere which is due to moist air in convecting plumes detraining and mixing with the drier environment (Pauluis and Mrowiec 2013). Hence, differences in the precipitation branch between models could partly be owing to differences in parameterisations of convection and entrainment. E.g. CCSM4 and NorESM1-M display similarities in the precipitating branch and both models use the Community Atmospheric Model (CAM) but in somewhat different versions. There is also some similarity between EC-Earth and ERA-Interim which both use the IFS atmospheric model from ECMWF.

Fig. 3
figure 3

Same as Fig. 2 but for simulations of the late twenty-first century (2080–2100) following the RCP8.5 emission scenario

Fig. 4
figure 4

Tropical-mean (\(15^{\circ } \text {S}\)\(15^{\circ } \text {N}\)) vertical profiles of LH, DSE and MSE. Black lines show historical simulations and dashed lines show RCP8.5 simulations

The hydrothermal stream function is also calculated for simulations of the late twenty-first century (2080-2100) in the RCP8.5 emission scenario. The overall result is that the stream function weakens (Table 1) and widens compared to 1980–2000 (Figs. 2, 3). Increased surface air temperatures increases the surface DSE while also increasing the surface LH following the Clausius–Clapeyron relationship. The moistening branch does not move closer or further away from the Clausius–Clapeyron line but instead widens along it. This suggests that relative humidity does not change by much with global warming. The increased DSE and LH results in increased MSE almost uniformly throughout the precipitation branch which widens the hydrothermal stream function without altering its shape. Profiles of LH, DSE and MSE show that the largest increase in LH is in the lower troposphere while the largest DSE increase is in the upper troposphere (Fig. 4). This results in an MSE increase that is almost uniform throughout the tropical troposphere.

Changes in the global atmospheric circulation with global warming are assessed using mainly three measures calculated from each climate model. The width, dLH, of the hydrothermal stream function is defined as the span in LH between the cooling branch and the outermost “tip” where the moistening and precipitation branches meet. The strength of the circulation, \(\psi \), is defined as the peak amplitude of the hydrothermal stream function. Lastly, the surface air temperature, \(T_s\), is calculated by taking the temperature and pressure in the lowest model level (\(k = KM\)) and assuming constant potential temperature in the layer, i.e. \(\theta _{\text {KM}} = \theta _s\),

$$\begin{aligned} \theta _{\text {KM}} = T_{\text {KM}} \left( \frac{p_0}{p_{\text {KM}}} \right) ^{R/c_p}, \quad \theta _{s} = T_{s} \left( \frac{p_{0}}{p_s} \right) ^{R/c_p} \Rightarrow T_s = T_k \left( \frac{p_{s}}{p_{\text {KM}}} \right) ^{R/c_p} . \end{aligned}$$
(10)

Surface air temperature is calculated from \(T\) on model levels to ensure that it has the same temporal and spatial resolution as the data and to be consistent for all models and reanalysis. It should be noted that the lowermost model level mid-point is very close to the surface, typically a few meters, so \(T_s \approx T_{\text {KM}}\).

The three metrics are calculated for each model for both periods 1980-2000 and 2080-2100, and the change between the two periods are denoted \({\varDelta } \text {dLH},\, {\varDelta } \psi \), and \({\varDelta } T_s\) respectively. Linear regression of the fractional changes \({\varDelta } \text {dLH} / \text {dLH}\) and \({\varDelta } \psi / \psi \) to \({\varDelta } T_s\) shows that the hydrothermal stream function widens by \(k_{{\varDelta } \text {dLH}} = 7.1\,\%\,\hbox { K}^{-1}\) and weakens by \(k_{{\varDelta } \psi } = -5.1\,\%\,\text {K}^{-1}\) (Fig. 5).

Interestingly, the weakening is offset so that \({\varDelta } T_s \sim 2\hbox { K}\) results in almost no weakening but a clear widening. The results by Held and Soden (2006) showed that the increase in precipitation is offset in a similar manner. It should be kept in mind that the hydrothermal stream function resides in LH–DSE coordinates and a widening thus implies a shift of LH and DSE in the precipitation branch towards higher values. The widening does not imply a widening in geometrical \((x,y,p)\) coordinates, although a widening of the Hadley cells in \((y,p)\) space has been found (Lu et al. 2008). Furthermore the weakening of the hydrothermal stream function is a combination of changes in both zonal and meridional overturning circulations. Hence, it may include a strengthening of the zonal-mean meridional overturning in midlatitudes as reported by Laliberté and Pauluis (2010) and Wu and Pauluis (2013). Changes in the zonal-mean circulations in various vertical coordinates are discussed later in this paper.

Fig. 5
figure 5

Change in amplitude and width of the hydrothermal stream function from 1980–2000 to 2080–2100 shown for each model. Dashed line is a linear fit with slope \(k_{\varPsi } = -5.1\,\%\,\text {K}^{-1}\) in (a) and slope \(k_{{\varDelta } LH} = 7.1\,\%\,\text {K}^{-1}\) in (b)

As stated above LH mostly increases in the lower troposphere and DSE mostly increases in the upper troposphere such that MSE increases nearly uniformly with height (Fig. 4). The fractional increase in tropical MSE per degree of warming is denoted as \({\varDelta } h / h\), where \({\varDelta } h\) is the difference between the historical (1980-2000) and RCP8.5 (2080-2100) simulations. The fractional change in tropical near-surface MSE can be predicted from the fractional changes in DSE and LH

$$\begin{aligned} \frac{{\varDelta } h}{h}&= \frac{{\varDelta } s}{h} + \frac{{\varDelta } l}{h} \approx \left( \frac{c_p}{h} + \frac{0.07 L_v q}{h} \right) {\varDelta } T_s, \end{aligned}$$

where \(7\,\%\,\text {K}^{-1}\) is the increase in specific humidity from Eq. 2. For typical values \(q = 16\hbox { kg}\hbox { kg}^{-1}\) and \(h = 340\hbox { kJ}\hbox { kg}^{-1}\) the increase in MSE is \({\varDelta } h / h = 1.1\,\%\,\text {K}^{-1}\). Regressing tropical-mean \({\varDelta } l/l,\, {\varDelta } s/s\) and \({\varDelta } h/h\) as a function of \({\varDelta } T_s\) shows that MSE increases by \(1.02\,\%\,\text {K}^{-1}\) in the CMIP5 models, which is close to the predicted value. In the lower troposphere \(0.30\,\% \,\text {K}^{-1}\) is due to DSE (i.e. temperature) increase and \(0.72\,\% \,\text {K}^{-1}\) comes from increased LH (i.e. specific humidity). In the upper troposphere there is almost no contribution from LH increase.

Fig. 6
figure 6

Difference in mean DSE tendency in the cooling branch, \({\varDelta } \dot{S}\), and LH flux, \({\varDelta } F_{\text {LH}}\), both at \(s = 315 \hbox { kJ}\hbox { kg}\), as a function of global warming \({\varDelta } T_s\). Dashed lines show linear fits, \(k_F = 0.63 \hbox { PW}\hbox { K}^{-1}\) and \(k_{\dot{S}} = 0.48\hbox { PW}\hbox { K}^{-1}\), where only the former is statistically significant at the \(95\,\%\) level

Previous studies (Held and Soden 2006; Stephens and Hu 2010) have studied how precipitation changes with global warming and found that the fractional increase is \(2\,\% \,\text {K}^{-1}\). If the change in precipitation is estimated as the product of mass flux in the precipitation branch and LH in lower tropospheric air (similar to Held and Soden (2006)) the fractional increase of precipitation is \(-5.1\,\% \,\text {K}^{-1} + 7.1\,\% \, \text {K}^{-1} = 2\,\% ~ \text {K}^{-1}\). This consistency suggests that the estimate is reasonable and that the hydrothermal stream function provides a good measure of the strength of the global atmospheric circulation.

Increased precipitation implies increased LH flux across DSE surfaces in the precipitation branch. LH flux is calculated as

$$\begin{aligned} F_{\text {LH}}(s) = \int _{l_{\text {min}}(s)}^{l_{\text {max}}(s)} \frac{\partial {\varPsi }(l,s)}{\partial l} l ~dl = - \int _{l_{\text {min}}(s)}^{l_{\text {max}}(s)} ({\varPsi }(l,s) - {\varPsi }_0) ~dl \end{aligned}$$
(11)

where \(l_{\text {min}}(s)\) and \(l_{\text {max}}(s)\) are the minimum and maximum LH values for which \({\varPsi }(l,s) = {\varPsi }_0\). The last step in Eq. 11 is obtained by integration by parts. To assert that only closed stream lines are included, \({\varPsi }_0 = -50\hbox { Sv}\). Calculating the change of the mean LH flux across DSE surfaces between 315 and 325\(\hbox { kJ}\hbox { kg}^{-1}\) and regressing onto \({\varDelta } T_s\) shows a statistically significant (\(p < 0.05\)) increase with global warming (Fig. 6). Since the hydrothermal stream function has closed stream lines it follows that the increased DSE resulting from increased LH flux in the precipitation branch must be lost by increased radiative cooling in the cooling branch. In fact, several studies (Held and Soden 2006; Stephens and Hu 2010; Bony et al. 2013) suggest that changes in precipitation are constrained by changes in radiative cooling. In the hydrothermal stream function, radiative cooling is negative DSE tendency for very dry air. The energy that is lost through this process can be calculated by integrating the negative DSE tendencies, i.e.

$$\begin{aligned} \dot{S}^{-} = \int _s^\infty \int _{l_{\text {min}}(s)}^{l_{\text {max}}(s)} \frac{\partial {\varPsi }(l,s)}{\partial l} \mu \left[ 0 -\frac{\partial {\varPsi }(l,s)}{\partial l} \right] ~dl ~ds, \end{aligned}$$
(12)

where the Heaviside function, \(\mu \), only selects the DSE tendency where it is negative, i.e. there is diabatic cooling. As in Eq. 11 we set \(l_{\text {min}}(s)\) and \(l_{\text {max}}(s)\) as boundaries and only count \({\varPsi }(l,s)\) inside the \({\varPsi }_0 = -50\hbox { Sv}\) stream line. Regressing the change \({\varDelta } \dot{S}^{-}\) onto \({\varDelta } T_s\) does not result in a statistically significant trend (\(p \sim 0.07\)). However, changes in the cooling branch and precipitation branch (Fig. 6) result in linear trends of similar magnitude but opposite signs indicating some balance between them. Note that the changes in radiative cooling are offset to be higher than the changes in LH flux. This can be due to various reasons such as LH fluxes being underestimated as moist convection is not resolved by the data, or because DSE can increase in the precipitation branch from other factors such as radiative absorption or cloud feedbacks.

Fig. 7
figure 7

Fractional change in poleward mass transport by the meridional overturning in MSE coordinates, \({\varPsi }(y,h)\), as a function of global surface warming, \({\varDelta } T_s\). Left figure shows Southern Hemisphere (SH) and right figure shows Northern Hemisphere (NH). Both figures include a fitted dashed line with slopes \(k_{\text {SH}} = -2.0\,\%\,\hbox { K}^{-1}\) and \(k_{\text {NH}} = -5.8\,\%\,\hbox { K}^{-1}\) for SH and NH respectively

Fig. 8
figure 8

Changes in meridional fluxes of LH, DSE, and MSE as a function of latitude for each model

The hydrothermal stream function captures both the zonal and meridional overturning cells, e.g. Hadley and Walker cells as well as the midlatitude eddies. By comparing changes in the meridional overturning stream functions, \({\varPsi }(y,\chi )\), to changes in the hydrothermal stream function, it is possible to estimate to what extent the meridional component contributes to the weakening of the global atmospheric circulation. The two hemisphere-wide overturning cells in \(y\)–MSE coordinates weaken with global warming in all models studied here (Fig. 7). The Northern Hemisphere (NH) cell weakens by \(-5.8\,\%\,\hbox { K}^{-1}\) and the Southern Hemisphere (SH) cell by \(-2.0\,\%\,\hbox { K}^{-1}\). This corresponds to changes of \({\varDelta } {\varPsi }_{\text {NH}}(y,h) = -6.4\hbox { Sv}\hbox { K}^{-1}\) and \({\varDelta } {\varPsi }_{\text {SH}}(y,h) =3.4\hbox { Sv}\hbox { K}^{-1}\) for the NH and SH cells respectively. Changes in the meridional overturning in LH (\({\varPsi }(y,l)\)) and DSE \(({\varPsi }(y,s))\) coordinates are not statistically significant but do indicate a weakening of the Hadley cells in both cases, consistent with results by Lu et al. (2008) and Wu and Pauluis (2013). It should be pointed out that Laliberté and Pauluis (2010) and Wu and Pauluis (2013) find a strengthening of the meridional overturning stream function in moist isentropic coordinates in NH midlatitudes in boreal winter and SH midlatitudes in austral winter. However, the results in Fig. 7 reflect changes in amplitude of the annual-mean meridional overturning cells so a strengthening in the midlatitudes in winter may not show. Any discrepancies between Laliberté and Pauluis (2010), Wu and Pauluis (2013) and the present study can thus be explained by the different methodologies. Comparing the weakening of the meridional overturning (Fig. 7) to the weakening of the hydrothermal stream function which is \(22.8\hbox { Sv}\hbox { K}^{-1}\) (Fig. 5) the results show that only a small part of the weakening can be explained by changes in the zonal-mean circulation. This implies that most of the weakening occurs in zonal asymmetric features such as zonal overturning circulations as well as local meridional overturning circulations at different longitudes. This could be because the meridional overturning circulation is constrained by e.g. the equator-to-pole temperature gradient. Results by Vecchi et al. (2006), Vecchi and Soden (2007), and Tokinaga et al. (2012) also suggest that most of the weakening happens in the zonal overturning circulations (e.g. Walker circulation) and not the meridional overturning circulations.

The meridional fluxes of LH, DSE and MSE (Eq. 6), and their changes from 1980-2000 to 2080-2100, are calculated from the meridional overturning stream functions in LH, DSE and MSE coordinates respectively (Figs. 1, 8). LH fluxes increase at almost all latitudes, consistent with an intensified hydrological cycle as predicted by Held and Soden (2006). Furthermore, there is a distinct increase in LH flux by the SH Hadley cell by up to \(1.5\) PW in some models. The increases in LH fluxes are partly counteracted by decreases in DSE fluxes. In the tropics, increases in equatorward LH fluxes are of similar magnitude as increases in poleward fluxes of DSE which results in almost no change in MSE fluxes. In the midlatitudes, the increases in poleward LH fluxes are generally somewhat larger than the decreases in poleward DSE fluxes giving a slight increase in MSE fluxes by <0.5 PW. Czaja and Marshall (2006) and Stone (1978) suggested that the combined poleward heat transport by the atmosphere and ocean are set by the solar constant, the radius of Earth and the planetary albedo, and might thus not change by much with global warming. This study presents increases in poleward MSE flux in the atmosphere while other studies suggest a decrease in poleward heat flux by the oceans associated with a slowdown of the Atlantic Meridional Overturning Circulation Weaver et al. (2012).

5 Summary and discussion

This study diagnoses long-term changes in the global atmospheric circulation between simulations of late twentieth and twenty-first century climates from the CMIP5 archive. Using a stream function in LH–DSE coordinates that captures both Hadley and Walker cells as well as midlatitude eddies shows that the global atmospheric circulation weakens with global warming. Associated with the weakening is a widening of the stream function in LH–DSE space due to an increase of tropical MSE that is almost uniform with height. The widening and weakening results in increased LH fluxes across DSE surfaces leading to more LH being converted into DSE through precipitation. The hydrothermal stream function has closed stream lines, thus any increase in LH or DSE of air masses must be balanced by a subsequent decrease or vice versa. It was here found that LH fluxes across the DSE surfaces \(s \approx 315\hbox { kJ}\hbox { kg}^{-1}\) are in approximate balance with radiative cooling. All of the above results are robust across the set of models included in this study.

The method of using a mass stream function in thermodynamical coordinates differs from the methods used by e.g. Held and Soden (2006), Vecchi et al. (2006) and Bony et al. (2013) where vertical velocity, \(\omega \), or a sea-level pressure gradient was used as a metric for circulation strength. Both SLP gradients and \(\omega \) include both adiabatic and diabatic effects so that the flows they represent may not be associated with an overturning circulation and/or energy transport. Furthermore, this study has based all calculations on 6-h model outputs, while many other studies (Vecchi and Soden 2007; Tokinaga et al. 2012; L’Heureux et al. 2013) have used output of a lower frequency, e.g. monthly. Data of a lower frequency than daily underestimates the fluxes by midlatitude eddies thus underestimating the strength of the global atmospheric circulation. Vecchi and Soden (2007) suggested that the tropical circulation weakens by \(5-10\,\%\,\hbox { K}^{-1}\) with global warming which is in line with the \(5.1\,\%\,\hbox { K}^{-1}\) found here. The differences between the results of this study and those obtained by e.g. Vecchi and Soden (2007) could be due to some of the reasons above and more.

This study has assessed centennial changes by studying the difference between 20-year averages of 1980–2000 and 2080–2100 for a number of climate models. 20-year averages are not long enough to eliminate the effects of multi-decadal variations in ENSO or Atlantic Multidecadal Oscillations (AMO). Hence, some of the differences between models and the periods 1980–2000 and 2080–2100 could be due to natural, internal climate variations in the models. However, since internal climate variations have different periods and magnitudes as well as timings in each model, the use of several models should minimise such effects. The results of this study are thus likely due to the external forcing in the RCP8.5 scenario and not due to internal climate variations in the models.

It is of interest to speculate how the results in this study would change if global climate models were cloud-resolving and thus did not need convective parameterisation. The hydrothermal stream functions would likely look somewhat different especially in the “precipitation” branch and would most likely have larger amplitudes as the up- and downdrafts in convective plumes would be captured. The global atmospheric circulation weakens with global warming because the fractional increase in precipitation is less than that of lower tropospheric water vapour (Stephens and Hu 2010), which follows from radiative balance arguments and does not depend on model resolution. However, as argued by several previous studies (Held and Soden 2006; Stephens and Hu 2010), cloud feedbacks and aerosol effects can impact the radiative balance. With global climate models of higher resolution the global atmospheric circulation would thus still weaken, although the exact fractional decrease in mass flux could be slightly different.

Increased LH fluxes across DSE iso-surfaces above \(s \sim 310\hbox { kJ}\hbox { kg}^{-1}\) suggests increased precipitation globally. This implies increased evaporation and, combined with the increased poleward LH fluxes, implies a strengthening of the global hydrological cycle in the models studied. A strengthened hydrological cycle, i.e. increased evaporation minus precipitation, EP, patterns (“rich get richer”) and poleward moisture transports is consistent with other studies (Held and Soden 2006; O’Gorman and Schneider 2008; Durack et al. 2012). LH tendencies in the atmosphere are linked to salt tendencies in the ocean since evaporating water from the ocean to the atmosphere increases atmospheric LH and oceanic salinity. Combining the results of this study with a study of the changes in global ocean circulation in thermohaline coordinates (Zika et al. 2012; Döös et al. 2012) could make up a framework for studying the thermodynamical coupling between the atmosphere and oceans. E.g., the increased DSE and LH and thereby increased precipitation found in this study would imply that the oceanic circulation should span higher values of sea surface temperature and surface salinity as well as increased freshwater input in some regions.

Although the overall shape of the hydrothermal stream function was similar in all models and reanalysis data, there were some variations in amplitude and size. Additional calculations (not shown) with ERA-Interim data at different horizontal resolutions as well as different temporal resolutions showed that the hydrothermal stream function weakens with coarser horizontal and temporal resolution. However, in this study there is no statistically significant correlation between the amplitude of the hydrothermal stream function and native horizontal resolution for the models given in Table 1. The inter-model differences between the hydrothermal stream functions are partly due to differences in the parameterisations of convection, e.g. entrainment. EC-Earth and ERA-Interim show similarities in the precipitation branch and both share the same atmospheric model. The same goes for CCSM4 and NorESM1-M.

The CMIP5 archive includes two IPSL models, IPSL-CM5A and IPSL-CM5B, which use different versions of the atmospheric component but are otherwise identical. This study has used IPSL-CM5B but the hydrothermal stream function was also calculated for IPSL-CM5A and the stream function was 50 Sv weaker in the latter, indicating a strong dependence on model specifics. Inter-model differences of the global atmospheric circulation could be better assessed by studying simulations in the Atmospheric Model Intercomparison Project (AMIP) where different atmospheric models are run with the same prescribed orography, sea surface temperatures, and surface fluxes from 1979-2008. This could also give more insight into the mechanisms behind the recent intensification of the Walker cell (Tokinaga et al. 2012; L’Heureux et al. 2013) and the global warming “hiatus” (Kosaka and Xie 2013; England et al. 2014).