This study examines at the process level the climate difference between 2002–13 and 1984–95 in ERA-Interim. A linearized radiative transfer model is used to calculate the temperature change such that its thermal radiative cooling would balance the energy flux perturbation associated with the change of an individual process, without regard to what causes the change of the process in the first place. The global mean error of the offline radiative transfer model calculations is 0.09 K, which corresponds to the upper limit of the uncertainties from a single term in the decomposition analysis.
The process-based decomposition indicates that the direct effect of the increase of CO2 (0.15 K) is the largest contributor to the global warming between the two periods (about 0.27 K). The second and third largest contributors are the cloud feedback (0.14 K) and the combined effect of the oceanic heat storage and evaporation terms (0.11 K), respectively. The largest warming associated with the oceanic heat storage term is found in the tropical Pacific and Indian Oceans, with relatively weaker warming over the tropical Atlantic Ocean. The increase in atmospheric moisture adds another 0.1 K to the global surface warming, but the enhancement in tropical convections acts to reduce the surface warming by 0.17 K. The ice-albedo and atmospheric dynamical feedbacks are the two leading factors responsible for the Arctic polar warming amplification (PWA). The increase of atmospheric water vapor over the Arctic region also contributes substantially to the Arctic PWA pattern.
Observational evidence indicates that the global surface temperature (GST) has exhibited signs of a muted or stalled warming since the beginning of the twenty-first century (Easterling and Wehner 2009; Knight et al. 2009; Solomon et al. 2010; Liebmann et al. 2010; Cowtan and Way 2014). Several mechanisms have been proposed to account for the muted global warming, including external cooling forcing due to solar radiation (Lean and Rind 2009; Hansen et al. 2011; Kaufmann et al. 2011) and internal climate variability (Kosaka and Xie 2013; England et al. 2014; Maher et al. 2014; Watanabe et al. 2014). It has been suggested that the primary source of the climate variability responsible for the slowdown of warming is the increasing deep-ocean heat uptake during this period (Hansen et al. 2011; Meehl et al. 2011; Kuhlbrodt and Gregory 2012; Balmaseda et al. 2013; Guemas et al. 2013; Chen and Tung 2014; Drijfhout et al. 2014; Kosaka and Xie 2013, 2015). There are several recent studies indicating that the volcanic and aerosol forcing may contribute to the recent slowdown of warming (e.g., Solomon et al. 2011; Santer et al. 2014; Haywood et al. 2014; Schurer et al. 2015). Most recently, the work of Karl et al. (2015) and Hausfather et al. (2017) seems to suggest that the global warming rate in the last two decades could be significantly underestimated owing to biases in ocean measurements.
The recent studies of Wu et al. (2007), Semenov et al. (2010), DelSole et al. (2011), and Wu et al. (2011) present statistical evidence indicating that the warming phase of a multidecadal oscillation associated with the thermohaline is the main mechanism responsible for the unprecedented warming pace occurring in the last two decades of the twentieth century. Therefore, both the fast warming pace in the last two decades of the twentieth century and the slowdown of warming in the first decade of the twenty-first century can be interpreted by the combined effect of external forcing (Trenberth and Fasullo 2009) and internal climate variability (Fyfe et al. 2013; Maher et al. 2014; Schmidt et al. 2014). Trenberth and Fasullo (2013) further pointed out that the apparent hiatus could be attributed to an increasing trend in heat storage into deep oceans in the past decade. The main objective of this paper is to perform a process-based decomposition of the surface temperature difference between two decades: one at the beginning decade and the other at the end of the fast warming since the 1980s. Specifically, based on the energy balance principle, we attempt to quantify how much of the warming between the two periods is directly due to the increase of CO2 alone and how much of the additional warming is associated with oceanic heat storage and other climate feedback processes, such as the surface albedo, water vapor, clouds, surface evaporation, and sensible heat flux feedbacks, as well as changes in atmospheric circulation. Because the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim) does not include information of natural and/or anthropogenic aerosols, our analysis presented below does not include information of the aerosol effect.
The remaining part of this paper is organized as follows. Section 2 describes the data and presents the key features of the difference in the climate mean states of the two periods. Section 3 describes the method for the process-based decomposition. The geographical distributions of process-based decompositions are given in section 4. Section 5 presents the results of process-based decompositions for the global mean of the surface temperature change and its global pattern. Section 6 discusses regional decompositions. The conclusions are given in section 7. The appendix discusses the main sources of the error term in our process-based decomposition analysis, corrections, and implications.
2. Data and the decadal mean surface temperature anomalies
All the data used in this study are obtained from the ERA-Interim (Uppala et al. 2008; Dee et al. 2011), as in Deng et al. (2012), Park et al. (2012, 2013), and Deng et al. (2013). The variables considered by us include temperature, specific humidity, ozone mixing ratio, cloud amount, cloud liquid/ice water content, downward solar energy flux at the top of the atmosphere (TOA), and surface albedo, surface sensible, and latent heat fluxes. The time series of the annual mean CO2 concentration from 1984 to 2013 is downloaded from the Earth System Research Laboratory website (http://www.esrl.noaa.gov/gmd/ccgg/trends/).
Figure 1 shows the annual mean and 5-yr running mean of global mean surface temperature anomalies (STAs) from 1979 to 2014, indicating a strong positive trend from the mid-1980s to the late 1990s but only a weaker positive trend in the first decade of the twenty-first century. The results vividly support the consensus in the literature that there are two distinct warming periods from the 1980s to the present (Knight et al. 2009; Liebmann et al. 2010). Based on Fig. 1, we will use the slowdown-warming decade (i.e., 2002–13) as the end decade of the fast warming and 1984–95 as the beginning decade. Therefore the difference between the two periods roughly represents the net warming during the fast warming period from the beginning of the 1990s (corresponding to the middle of 1984–95) to the middle of the 2000s (close to the middle of 2002–13). We calculate the temporal means of each variable for the periods of 1984–95 and 2002–13, referred to as the mean states of 1984–95 and 2002–13, respectively. The symbol Δ denotes the difference in the mean states between 2002–13 and 1984–95.
Figures 2a and 2b are the decadal mean surface temperature anomalies for the periods of 1984–95 and 2002–13, respectively (the anomalies are defined with respect to the 30-yr mean between 1981 and 2010). It is obvious that the decade of 2002–13 is still much warmer than 1984–95 over most places around the globe, despite the warming rate during the period 2002–13 being weaker. The difference between 2002–13 and 1984–95 (Fig. 2c) exhibits two salient patterns. The first one is the pronounced polar warming amplification (PWA) pattern over the Arctic and a weak PWA over the Antarctic. The PWA is one of the most recognized and well-defined global warming patterns in the trends derived from station observations and all reanalysis products (Polyakov 2002; Johannessen et al. 2004), as well as in all climate model simulations radiatively forced by anthropogenic greenhouse gases (e.g., Holland and Bitz 2003; Gillett et al. 2008; Bindoff et al. 2013; Bracegirdle and Stephenson 2012, 2013; Flato et al. 2013). The second is the relatively stronger warming over the tropical portion of the Indian, Pacific, and Atlantic Oceans than over midlatitude oceans (Hartmann et al. 2013). However, there are three regions that had colder surface temperature anomalies in 2002–13 than in 1984–95. As reported in IPCC Fifth Assessment Report (AR5; Bindoff et al. 2013), the first is the eastern tropical Pacific, the second over the Southern Ocean along the Antarctic Circumpolar Current, and the third over the west coastlines of North America and South America (Fig. 2a vs Fig. 2b; Fig. 2c). To the west of the cooling over the eastern tropical Pacific is a warming pattern over the western tropical Pacific, referred to as a La Niña–like pattern (Meehl et al. 2011; Sohn et al. 2013). The cooling trend over the Southern Ocean (Fan et al. 2014; Armour et al. 2016; Kostov et al. 2016) is surrounded by a weak warming trend in the Antarctic (Comiso 2000; Bertler 2004; Monaghan et al. 2008; Nicolas and Bromwich 2014) and southern subtropics, defined as the warm–cold–warm (WCW) pattern.
3. The process-based decomposition method
The majority of existing methods for studying climate feedback and sensitivity are built on the premise that the temperature change is the response to radiative energy exchange with outer space, and therefore they consider primarily the radiation perturbations at the TOA. Accordingly, the climate forcing and feedback agents are mechanisms that directly affect the radiative budget at the TOA, namely, the incoming solar energy flux at the TOA, changes in the atmospheric composition (e.g., CO2), water vapor, cloud, surface albedo, and lapse rate of atmospheric temperature. The popular TOA-based climate feedback analysis methods include the partial radiative perturbation (PRP) method (Wetherald and Manabe 1988), the cloud forcing analysis method (Cess et al. 1990), and the radiative kernel method (Soden and Held 2006; Soden et al. 2008). Readers may refer to Bony et al. (2006) for a comprehensive review on these methods and Flato et al. (2013) for their applications in assessing climate sensitivity and climate feedbacks in CMIP5 climate simulations. Boer and Yu (2003) extend the PRP method to study spatial patterns of climate sensitivity. However, none of these TOA-based methods calculates the partial temperature changes due to individual feedback processes. In addition, the TOA-based feedback analysis framework cannot explicitly take into consideration the internal nonradiative processes (e.g., oceanic heat storage, convective and large-scale atmospheric energy transport, and surface turbulent fluxes such as evaporation and sensible heat fluxes). Bates (2007) further pointed out that the usage of an equation similar to the amplification formula of electronics in the PRP method is only figurative since the climate system does not follow the physical principles invoked for the electronics application.
The online feedback suppression method is designed to calculate partial temperature change due to a specific feedback process. It requires running a climate model twice: one with all processes on (original climate system) and the other with one specific process turned off (virtual climate system). The difference between such a pair of simulations corresponds to the partial temperate change due to the process under consideration (Hall and Manabe 1999; Schneider et al. 1999). However, the difference between two different systems inferred from the online feedback suppression method does not reflect exactly the effect of the suppressed feedback in the original climate system because the difference also includes the compensating effects from the other feedbacks (Cai and Lu 2009). As a result, the partial temperature change due to a specific feedback inferred by the online feedback suppression method also includes the difference in the other feedbacks between the original and virtual climate system.
To quantitatively separate the additive contributions to the global mean warming and their spatial pattern shown in Fig. 2c from the external forcing and various climate feedback processes (both radiative and nonradiative feedback processes), we have adopted the coupled climate feedback response analysis method (CFRAM) used in Deng et al. (2012), Sejas et al. (2014), and Hu et al. (2016). According to Lu and Cai (2009), given a vertical profile of nonradiative or non-temperature-induced radiative heating perturbations Δ(X) associated with a process denoted X, in a horizontal grid point, the CFRAM method allows us to calculate the vertical profile of its associated temperature changes Δ(X) by requiring the corresponding vertical profile of linearized thermal radiative cooling perturbations (∂/∂)Δ to be in radiative equilibrium balance with Δ(X):
where (∂/∂) is the Planck feedback matrix whose jth column represents the vertical profile of (linearized) thermal radiative cooling perturbation due to 1-K warming in the jth layer alone. Note that, as in Lu and Cai (2009), all the vertical profiles are denoted in the form of a vector with the bottom component corresponding to the surface layer and the remaining components corresponding to the atmospheric layers from the top layer to the layer next to the surface layer.
Using the data derived from ERA-Interim, we can evaluate each of the vertical profiles of non-temperature-induced radiative heating perturbations between the two periods listed on the right-hand side (RHS) of
where and denote the vertical profiles of the solar energy absorbed and the net longwave radiation emitted (or cooling rate) in each layer, respectively. In (2), the superscripts CO2, WV, CLD, AL, and O3 after the symbol Δ denote the partial radiative heating perturbations due to changes, respectively, in the concentration of CO2 alone, in atmospheric water vapor alone, in atmospheric cloud properties alone (which includes cloud area and ice and water cloud concentrations), in the surface albedo alone, and in ozone alone. These vertical profiles of partial radiative heating perturbations are calculated as the difference between two calculations of a radiative transfer model: one with inputs taken from the mean state of 1984–95 and the other with identical inputs except the field denoted by the superscript, which is taken from the mean state of 2002–13. As in Deng et al. (2012), all of the radiative heating calculations are made with the Fu–Liou radiative transfer model (Fu and Liou 1992, 1993). We note that the annual mean incoming solar radiative flux at the TOA used in ERA-Interim does not vary in time. As a result, it has no contribution to (2).
Besides these non-temperature-induced partial radiative heating perturbations, we have also calculated partial radiative heating perturbations due to changes in the vertical profile of temperature, which is denoted as Δ(Temp). We have verified the following:
where ∂/∂ is evaluated using the inputs taken from the mean state of 1984–95 and Δ is the vertical profile of the difference in atmospheric and surface temperatures between 2002–13 and 1984–95. Equation (3) implies that partial radiative heating perturbations due to changes in the vertical profile of temperature can be calculated under the linear approximation. The total radiative heating perturbation Δtotal( − ) can also be calculated in a similar fashion—namely, as the difference of two calculations of the Fu–Liou radiative transfer model: one with all input fields taken from the mean state of 1984–95 and the other with all input fields from the mean state of 2002–13. We have also verified the following (not shown here):
This is equivalent to saying that the changes in radiative heating rates between the two mean states of the two periods are small enough that they can be partitioned as the sum of individual partial radiative heating perturbations due to changes in one variable–parameter alone. This is the basis of our decomposition at the process level.
Besides radiative heating perturbations, changes in atmospheric and surface temperatures are also due to nonradiative heating perturbations and heat storage terms. ERA-Interim only has very limited nonradiative heating fields, such as surface latent and sensible heat fluxes. As in Deng et al. (2012), we estimate the sum of nonradiative heating perturbations and storage terms based on the energy balance equation:
Note that the surface component of Δnon_rad represents the sum of changes in surface latent and sensible heat fluxes, land surface heat storage (plus some heat loss or gain due to runoff and snow/ice melting or freezing) when the grid point is land. Over a grid point in oceans, it represents the sum of changes in surface latent and sensible heat fluxes, the net convergence of energy transport by oceanic motions into the entire oceanic column, and the oceanic heat storage term. For easy reference, we simply denote the surface component of Δnon_rad as ΔQs and the remaining components of Δnon_rad for the atmospheric layers are all zero (referred to as surface processes hereafter). Similarly, we denote the atmospheric components of Δnon_rad as Δ(ATM) whose surface component is zero and the atmospheric components equal to their counterparts of Δnon_rad Because the atmospheric heat storage is very small when taking the decadal mean, Δ(ATM) represents atmospheric dynamic processes, including latent heat perturbations due to changes in convections and large-scale atmospheric motions and changes in dry static energy flux convergence in each atmospheric layer by vertical and horizontal atmospheric motions as well as due to changes in sensible heat fluxes that enter the atmosphere from the surface below.
To have an accurate estimate of Δnon_rad, we need to derive the term Δtotal( − ) directly from the ERA-Interim dataset instead of our offline calculations. However, only the surface component of Δtotal( − ), which is equal to the change in the net downward radiative fluxes at the surface, is available from the ERA-Interim dataset. Therefore, to best utilize the available information from ERA-Interim for estimating the net nonradiative heating fields, we define
where and are the net solar and longwave radiative fluxes at the surface derived from the ERA-Interim dataset whereas Ss and Rs are the net solar and longwave radiative fluxes at the surface derived from our offline calculations. It follows that represents the errors in estimating the net nonradiative heating field at the surface in the offline radiative heating calculations. Maps of and the associated errors in our estimating partial temperature changes will be shown in the appendix. We will also present evidence suggesting that a large portion of the offline errors comes from our estimate of partial radiative energy flux perturbations associated with cloud changes.
Using the surface latent heat and sensible heat fluxes, QLH and QSH, provided by ERA-Interim, we can further decompose changes in net nonradiative heating at the surface layer based on the surface energy balance equation:
Note that the sign convention of QLH and QSH is that negative values indicate upward energy fluxes, which cause cooling at the surface. The term ΔQOCH represents changes in land surface heat storage when the surface is land. Over oceans, it represents the sum of the net convergence of energy transport by oceanic motions into the entire oceanic column and the storage term at a given grid point. Since the term ΔQOCH is small over land compared to that over oceans, we simply refer to this term as the oceanic dynamic and heat storage term.
Replacing Δ(X) in (1) with each of the five terms on the RHS of (2), we obtain the vertical profile of partial temperature changes Δ(X), where X denotes the superscript of the corresponding term. Replacing Δ(X) with Δ(ATM), whose surface component is zero and atmospheric components equal their counterparts of Δ(non_rad), we obtain Δ(ATM) by solving (1). We obtain Δ(OCH) by letting the surface component of Δ(X) be ΔQOCH and setting its atmospheric components to be zero. Similarly we obtain Δ(LHFC), Δ(SHFC), and Δ(ERR) by letting the surface component of Δ(X) be ΔQLH, ΔQSH, and , respectively, with zero values in the atmospheric components. Through the procedures outlined above, we obtained a total of 10 vertical profiles of partial temperature changes 5 for radiative processes, 4 for nonradiative processes, and 1 associated with errors in our offline radiative heating calculations. The surface components of these 10 vertical profiles of Δ(X) correspond to the partial surface temperature changes associated with the process X [denoted as hereafter].1 Note that is purely associated with the change in the downward thermal radiation resulting from changes of air temperatures in response to Δ(ATM) since its surface component is zero (Fig. 3i).
Here, we wish to add that like any other offline feedback analysis methods, our feedback analysis is a postprocessing diagnostic that cannot predict changes in other fields. The majority of other offline methods also need the information of the total temperature change as the input field for the feedback analysis. However, the CFRAM method allows us to explicitly calculate partial temperature changes associated with individual processes without requiring the information of the total temperature change. The sum of these partial temperature changes can be directly compared to the total temperature change in observations. It is seen in Fig. 2d that the spatial pattern of and their numerical values are highly consistent with the actual temperature changes (Fig. 2c) derived directly from ERA-Interim. The small differences between and are due to errors introduced in linearization of the radiative transfer model, namely, approximating the RHS of (4) with the terms on the left-hand side. Furthermore, numerical values of − are very close to zero and much smaller than individual . Also, the map of − does not have a coherent spatial structure and explains little spatial variation of (not shown). This validates the linearization of the radiative transfer model as far as the mean difference between the two periods is concerned.
Because ERA-Interim does not include information of natural and/or anthropogenic aerosols, our method does not allow us to make a direct estimate of the aerosol effect, as we have done for other variables. According to Cai and Kalnay (2005), a reanalysis made with a frozen model is still capable of capturing the temperature trend as well as its variations due to an anthropogenic forcing at its full strength (at least 95% level) after a few hundred analysis cycles, even though the forcing itself is never directly assimilated in the analysis. Therefore, the temperature fields derived from ERA-Interim can still capture the temporal cooling due to volcanic forcing and the cooling trend due to anthropogenic aerosol forcing. In our analysis, we obtain all of these partial temperature changes using non-temperature fields (e.g., atmospheric water vapor) derived from ERA-Interim except for three terms: the atmospheric dynamics, oceanic heat storage/dynamics, and offline calculation error terms. Because (i) these three terms use ERA-Interim’s temperature information and (ii) the sum of all partial temperature changes is convergent to the (total) temperature change derived from the reanalysis, we believe that the signals from the volcanic and anthropogenic aerosol forcings are blended in these three terms.
4. Direct forcings, air temperature feedbacks, and surface temperature changes
Displayed in Fig. 3 are the energy flux perturbations at the surface due to the direct effect of individual processes.2 As discussed in Cai and Tung (2012) and Sejas et al. (2016), the enhancement of the downward longwave radiative flux at the surface due to a spatially uniform increase CO2 tends to be stronger over the regions where atmospheric water vapor and/or clouds are scarce, such as polar regions and high elevation areas (Fig. 3a). The reduction of the solar energy flux reaching to the surface (Fig. 3b) is due to the upward trend of the stratosphere ozone (Fig. 4b). Because of a stronger increase of the stratosphere ozone in the tropics as well as more available incoming solar energy, the reduction in the downward solar energy flux at the surface is stronger in the tropics. The melting of sea ice over the Arctic Ocean, as indicated by a substantial reduction of surface albedo (Fig. 4b), results in more solar energy absorbed there (Fig. 3f). The increase of sea ice coverage over most regions along the Antarctic continental shelf leads to a reduction of solar energy absorption due to the increase of surface albedo. The direct effect of an increase (decrease) in atmospheric water vapor is the strengthening (reduction) of the downward longwave radiative fluxes at the surface (Fig. 3g vs Fig. 4c). The relationship of the direct effect of cloud changes strongly depends on the spatial distribution of the incoming solar energy flux. In general, the increase (decrease) of clouds in high latitudes where the incoming solar energy flux is much weaker leads to a stronger (weaker) downward longwave radiative flux at the surface (Fig. 3h vs Fig. 4d). Because of the stronger incoming solar energy flux outside high latitudes, cloud-albedo-induced changes of downward shortwave fluxes at the surface are stronger than their longwave effects. As a result, in the tropics and midlatitudes, the net of the clouds’ direct effect is positive (negative) over the regions where clouds decrease (increase). The remaining four panels in Fig. 3 (i.e., Figs. 3c,d,e,i) are the direct effect of the nonradiative processes and are obtained directly from ERA-Interim. Note that the direct effect of the atmospheric dynamics term (Fig. 3i) is zero everywhere because the energy transport by the atmospheric motions by itself, by definition, has no direct contribution to energy fluxes at the surface.
Following Sejas and Cai (2016), the indirect effect of a process X, which is defined as the downward longwave radiative flux perturbations at the surface induced by the air temperature changes associated with X, can be evaluated using
where (∂Rs/∂Tj, j = 1, 2, …, M) is the surface version of the temperature kernel and [, j = 1, 2, . . , M] is the vertical profile of partial air temperature changes associated with Δ(X) in (1). Akin to the TOA version of the temperature kernel put forward in Soden and Held (2006), the jth element of the surface version of the temperature kernel ∂Rs/∂Tj represents the downward longwave radiative flux perturbation at the surface caused by a temperature change of 1 K at the jth atmospheric layer. As discussed in Sejas and Cai (2016), [, j = 1, 2, . . , M] is determined with (1) by taking into consideration not only the vertical profile of energy flux convergence perturbations Δ(X) in the atmosphere and its value at the surface but also the associated surface temperature change . The association of [, j = 1, 2, . . , M] with is through the temperature feedback loop via the thermal radiative coupling between air temperatures and the surface temperature, as elicited in Sejas and Cai (2016). Therefore, corresponds to the surface radiative flux change due to [, j = 1, 2, . . , M], or “the air temperature feedback.”3 It is the upward longwave radiative flux perturbation emitted from the surface due to (Fig. 5) that is in balance with the sum of (Fig. 6) and its counterpart shown in Fig. 3.
To explain the linkages among individual direct forcings at the surface (Fig. 3), their air temperature feedbacks (Fig. 6), and surface temperature changes (Fig. 5) via the temperature feedback loop, we consider two special scenarios first and then the general situation [readers may consult Sejas and Cai (2016) for more detailed discussions on the importance of the temperature feedback loop]. The first special case is that most of the direct heating perturbation (without temperature feedbacks) is at the surface. Examples of the first special scenario include the change of the solar energy absorbed by the atmosphere and surface due to the change in the surface albedo Δ(AL), ocean heat storage and dynamics term Δ(OCH), and surface sensible Δ(SHFL) and latent heat Δ(LHFL) fluxes Δ(SHFL) and Δ(LHFL). In this case, nonzero values of [, j = 1, 2, . . , M] are mostly due to the temperature feedback loop. As a result, the spatial patterns of the direct forcings at the surface, the associated air temperature feedback, and partial surface temperature changes (Figs. 3c, 5c, and 6c for X = OCH; Figs. 3d, 5d, and 6d for LHFL; Figs. 3e, 5e, and 6e for SHFL; and Figs. 3f, 5f, and 6f for AL) are nearly identical with the perfect positive map correction. It is of importance to point out that not only does the air temperature feedback in the first special case have the same sign as the direct forcing but also the magnitude of the air temperature feedback is much stronger than the direct forcing. Therefore, the air temperature feedback can amplify the original direct forcing at the surface by more than doubling. The exact amount of amplification by the air temperature feedback depends on the amount of the longwave absorbers in the atmosphere as well as the vertical profile of temperatures, explaining why the ratio of the air temperature feedback to the original direct forcing varies spatially.
The second special scenario is solely related to the atmospheric dynamics term, in which the energy flux convergence perturbations are caused by changes in atmospheric convections and large-scale atmospheric advective processes and its direct forcing at the surface is zero by definition. Obviously, changes in the surface temperature in this case are entirely through the thermal radiative coupling between the atmosphere and the surface or the air temperature feedback. The comparison of Fig. 4d and Fig. 6i indicates that in the tropics, negative (positive) values of the dynamics-induced air temperature feedback tend to be over the regions where cloud liquid/ice water content increases (decreases). An increase (decrease) of clouds is indicative of the intensification (weakening) of tropical convections, which acts to transport more (less) energy upward, causing warming (cooling) anomalies in the upper atmosphere at the expense of cooling (warming) in the lower atmosphere. As illustrated in Lu and Cai (2009) and Cai and Tung (2012), the dynamics-induced cooling (warming) in the lower troposphere results in a reduction (enhancement) of the downward thermal energy fluxes at the surface, giving rise to negative (positive) values of and cooling (warming) at the surface (Fig. 5i). However, part of the negative values of in the tropics, which includes the effect of changes in tropical convections and poleward energy transport, can also be indicative of the strengthening of the poleward energy transport. By explicitly isolating changes in due to changes in atmospheric convections from those due to changes in atmospheric poleward energy transport, Lu and Cai (2009) and Cai and Tung (2012) showed that large-scale dynamics-induced has positive values in high latitudes but negative values in low latitudes. Therefore, the pattern of the dominance of negative values in tropics but large positive values over the Arctic in Fig. 6i indicates a strengthening of the poleward energy transport in the Northern Hemisphere. It is the enhanced downward thermal radiative flux [or positive values of ] due to the strengthened poleward energy transport that is responsible for the dynamical polar warming amplification, as elicited first by Cai (2005, 2006) and Cai and Lu (2007) with a simple four-box radiative–transportive climate model and Lu and Cai (2010) in a dry GCM model.
For the general situation in which both atmospheric and surface components of a direct forcing Δ(X) are not zero, such as , , Δ(WV), and Δ(CLD), changes in air temperatures consists of two parts: the direct response to the atmospheric components of Δ(X) and the indirect response to the surface component via the temperature feedback loop (or the thermal–radiative coupling between the atmosphere and the surface). However, the main contributors to are air temperature changes in the lower troposphere because ∂Rs/∂Tj is very small above the lower troposphere. Therefore, the spatial pattern of can be distinctly different from the surface component of Δ(X) only when values of Δ(X) in the lower troposphere are very large. Because the surface components of all , Δ(WV), and Δ(CLD) are much larger than their atmospheric components, (Fig. 6b), ΔIR(WV) (Fig. 6g), and ΔIR(CLD) (Fig. 6h) exhibit similar spatial patterns as their counterparts shown in Fig. 3. As a result of the overlapping of the absorption bands of CO2 and water vapor (and liquid/ice water), the elevation of the maximum heating perturbation due to an increase of CO2 is highest (about 850 hPa) in the tropics where both atmospheric water vapor and clouds in the mean state are abundant (Cai and Tung 2012; Taylor et al. 2013). The elevation of the maximum heating perturbation then decreases with latitude as the atmospheric water vapor becomes less abundant. This explains why (Fig. 4a) has its maximum value in the midlatitudes, large values in the tropics, and minimum values in the polar regions, which is distinctly different from its counterpart shown in Fig. 3. Again as a result of the amplification by the temperature feedback loop, all , , ΔIR(WV), and ΔIR(CLD) have larger amplitude than their counterparts shown in Fig. 3. This explains why the spatial patterns of the associated surface temperature changes (Fig. 5) are more similar to ΔIR(X) (Fig. 6) than (Fig. 3).
5. Contributions to global mean and pattern
Figure 7a shows the global mean values of these individual , their sum, and (the hatched portion of the last bar on the right). The global mean of the sum of individual is 0.28 K, very close to 〈〉 = 0.27 K (〈⋅〉 denotes the global mean). As discussed in section 3, the difference of 0.01 K between 〈〉 and corresponds to the global mean error due to linearization of the radiative transfer calculation, which is 4% of 〈〉.
It is seen that the partial surface temperature change due to the increase of CO2 from 1984–95 to 2002–13 alone (i.e., without any feedbacks) explains nearly half (0.15 K) of the observed global mean warming between the two periods. The largest term for the global mean of the surface temperature change turns out to be from the oceanic heat storage term (note that the global mean of oceanic dynamic term is zero), contributing as much as 1.29 K to 〈〉. However, the enhancement of upward surface turbulent heat fluxes, mostly associated with the strengthening of evaporation over oceans, has a cooling contribution of −1.18 K to 〈〉, which cancels out a large fraction of 〈〉. The net contribution from the oceanic heat storage and ocean surface processes, denoted as 〈〉, is 0.11 K. Changes in clouds contribute 0.14 K to 〈〉. The water vapor feedback adds an additional 0.09 K to the global mean warming. The increase of O3 contributes a cooling of −0.13 K to the global mean surface temperature change. The atmospheric dynamic processes also contribute negatively (−0.17 K) to 〈〉. The offline radiative transfer model calculation error 〈〉 = 0.09 K, which is still quite large even after we have removed the portions of the offline radiative transfer model calculation errors that are correlated with the cloud-induced radiative energy flux perturbations at the surface (see the appendix for more details).
The contribution analysis based on the global–regional mean warming cannot capture the strength of since the values of are not uniformly positive over the globe or a given region. Next, we wish to quantify the relative contributions of these individual shown in Fig. 5 to the global–regional spatial pattern and amplitude of shown in Fig. 2c. Following Cai and Tung (2012) and Deng et al. (2012), we calculate the pattern-amplitude projection to over the area A from each according to the following:
where ϕ and λ are latitude and longitude, respectively, and a is the mean radius of the earth. For easy reference, we refer to the term
as the amplitude of the spatial pattern of over the area A. One can easily verify that if were uniformly positive or uniformly negative over A, the results of the spatial mean contribution analysis would be the same as the pattern–amplitude analysis, as far as the sign and the relative magnitude are concerned.
For the global pattern contribution analysis, the area A is the entire global surface . The amplitude of the global pattern of is 0.63 K, more than twice as large as the global mean of . This indicates that the time mean difference of the surface temperature between the two periods has a large spatial variation, as shown in Fig. 2c. The largest positive contribution to the global pattern of comes from the surface albedo feedback , whose global pattern projection is as large as 0.19 K. The water vapor and cloud feedbacks contribute positively to the global pattern of by 0.11 and 0.15 K, respectively. According to Fig. 7b, also contributes positively to the global pattern of [ = 0.07 K], although its pattern contribution is relatively smaller than its global mean contribution. The net contribution of the oceanic heat storage and ocean surface processes to the global pattern of is about 0.05 K. The pattern projection contributions from and are negative, similar to their contributions to the global mean of . The value of is small because (Fig. 5b) is relatively small everywhere compared to .
6. Contributions to regional patterns
Our regional contribution analysis is focused on the four key differences between the two periods discussed in section 2: (i) relatively large warming over tropical oceans, (ii) a La Niña–like pattern over the tropical Pacific basin, (iii) the pronounced PWA pattern in the northern extratropics, and (vi) the WCW pattern in the southern extratropics.
The amplitude of the spatial pattern of over the tropics (30°S–30°N) is 0.39 K (Fig. 8a), which is only slightly smaller than the domain-averaged warming because is positive nearly everywhere in the tropics except over the eastern tropical Pacific. As indicated in Fig. 5c, large positive values of are found mainly over tropical oceans. It is the leading term that greatly contributes to the tropical spatial pattern of . The large positive values of are collocated with strong cooling associated with the enhancement of surface latent heat fluxes (Fig. 5d) as well as surface sensible heat fluxes (Fig. 5e). Therefore, surface turbulent heat fluxes are the main processes that carry the excessive heat away from tropical oceans entering the atmosphere. According to Fig. 8a, the net effect of ocean surface processes is the largest contributor (about 0.22 K) to the warming pattern in the tropics. Accompanied with the tropical warming are increases in both clouds and moisture in the tropical atmosphere. In addition, (Fig. 5i) also acts to offset the excessive tropical ocean warming, contributing negatively to the tropical warming pattern by as much as −0.32 K.
Some recent studies suggest that the La Niña–like trend of SST in the tropical eastern Pacific can be regarded as a combined effect of the internal climate variability and global warming (Zhang et al. 2011; Dong and Zhou 2014). The OCH process (Fig. 5c) shows a weak cooling effect in the equatorial central Pacific, which seems to be consistent with the strengthening of the upper-ocean meridional overturning circulation proposed by Zhang et al. (2011). Figure 4d shows that clouds have an upward trend in the central and eastern tropical Pacific from 1984–95 to 2002–13. The net effect of the increase in clouds acts to reduce the warming in the central and eastern tropical Pacific (Fig. 5h), thereby contributing positively to the La Niña–like trend of [ = 0.22 K]. The water vapor feedback (Fig. 5g) is the other key process responsible for the La Niña–like pattern of over the tropical Pacific [ = 0.14 K]. The cooling anomalies in the tropical central and eastern Pacific from water vapor feedback are due to the drier troposphere above the eastern half of the tropical Pacific (Fig. 4c). The pattern of (Fig. A2e) is highly negatively correlated with over the tropical Pacific.
b. Northern extratropics
The dominant feature of over the northern extratropics is the PWA pattern. The contribution analysis to the spatial pattern over the northern extratropics is made by applying (9) to the area north of 30°N. The amplitude of the spatial pattern of over the northern extratropics is 1.01 K (Fig. 8b).
The four three positive contributors are , , , and , accounting for 0.33, 0.20, 0.16, and 0.12 K, respectively, of the spatial pattern–amplitude of over the northern extratropics. Their large positive contributions are due to their large projections to the Arctic PWA pattern. In accordance with previous studies (Polyakov 2002; Holland and Bitz 2003; Johannessen et al. 2004; Screen and Simmonds 2010), the main contribution to the Arctic PWA is from caused by the reduction of surface albedo via sea ice melting. The strengthening of the atmospheric dynamical energy transport has also been regarded as a major contributor to the Arctic PWA (Moritz et al. 2002; Graversen et al. 2008). The water vapor feedback also adds additional warming in the Arctic region (Fig. 5g). These three processes are the main factors leading to the PWA in the Arctic. The spatial pattern of also exhibits a very clear PWA pattern (Fig. 5a), as in observations. The net effect of the surface processes has a strong negative contribution (about −0.17 K) to the spatial pattern of over the northern extratropics.
c. Southern extratropics
As mentioned in section 2, the spatial pattern of over the southern extratropics is characterized by the WCW pattern, namely a pronounced cooling trend over the Southern Ocean but a strong warming trend in the Antarctic and southern subtropics. The pattern of contributes positively to the Antarctic PWA signal. Because the WCW pattern is more pronounced than the Antarctic PWA signal, contributes little to the pattern of over the southern extratropics. As summarized in Fig. 8c, the two leading factors showing pronounced positive contributions to the WCW pattern (0.51 K) over the southern extratropics (30°–90°S) are and . The former contributes 0.28 K and the latter 0.27 K. The positive contribution of the albedo feedback mainly comes from both the Antarctic, where contributes strongly to the Antarctic PWA, and over the southern ocean next to Antarctic continent, where the growth of sea ice leads to < 0. The large contribution of to over the southern extratropics implies that it is the most dominant factor responsible for the WCW pattern. As in other regions and over the globe, has a large contribution to the pattern of over the southern extratropics [ = 0.18 K]. Unlike that over the Arctic, has a negative contribution to the Antarctic PWA pattern (Fig. 6i).
7. Discussion and conclusions
This study decomposes the climate difference between the periods of 2002–13 and 1984–95 into partial temperature differences due to external forcing and various internal climate processes by applying the climate feedback–response analysis method (CFRAM) to ERA-Interim. Despite that the warming rate during 2002–13 is much lower, the global mean surface temperature of the decade of 2002–13 is 0.27 K warmer than 1984–95. The surface temperature difference between the two periods is characterized by four distinct regional patterns: (i) a pronounced polar warming amplification (PWA) pattern over the Arctic and a relatively weaker PWA over the Antarctic, (ii) the relatively stronger warming over the tropical portion of the Indian, western Pacific, and Atlantic Oceans than over the midlatitude oceans, (iii) a La Niña–like pattern over the tropical Pacific, and (iv) a large-amplitude cooling trend over the Southern Ocean surrounded by a warming trend in the Antarctic and southern subtropics (referred to as the WCW pattern). All these features are consistent with other observational studies and IPCC AR5 climate simulations (Hartmann et al. 2013; Bindoff et al. 2013; Flato et al. 2013).
Our process-based decomposition indicates that the direct effect of the increase of CO2 (0.15 K) is the largest contributor to the global warming between the two periods (about 0.27 K). Changes in clouds contribute about 0.14 K to the global mean warming. Although the energy released by the oceanic heat storage term is partially canceled out by the enhancement of evaporation and sensible heat flux from the surface, the combined effect of these three surface processes still contributes 0.11 K to the global mean warming. The increase in atmospheric moisture adds another 0.1 K to the global surface warming, but the enhancement in tropical convections acts to reduce the surface warming by 0.17 K.
The oceanic heat storage term is responsible for most of the warming in the tropical oceans as well as the cooling over the southern oceans. More atmospheric water vapor and more deep clouds above the western tropical Pacific act to further strengthen the warming there, whereas the reduction of atmospheric water vapor and increase of clouds above the eastern tropical Pacific lead to a cooling trend there. The combined effects of water vapor and cloud feedbacks contribute positively to the La Niña–like pattern.
The ice-albedo feedback is the leading factor responsible for the Arctic PWA pattern, the warming along the coastline of the West Antarctic and Antarctic Peninsula, and the cooling along the coastline of the remaining portion of the Antarctic continent. The atmospheric dynamical feedback is the second largest contributor to the Arctic PWA pattern. The increase of atmospheric water vapor over the Arctic region also contributes substantially to the Arctic PWA pattern.
The period between 1984–95 and 2002–13 exhibits the fastest warming in the past three decades. Our process-based decomposition presented in this study suggests that the ocean heat storage term contributes most of the warming between 2002–13 and 1984–95. The results support the statistical analysis of the unprecedented warming pace occurring in the last two decades of the twentieth century to the ocean heat storage term reported in Wu et al. (2007), Semenov et al. (2010), DelSole et al. (2011), and Wu et al. (2011). Here we further conjecture that the heat release from the vast ocean surface, particularly from the tropical oceans, may have already reached a plateau in the first decade of the twenty-first century (i.e., 2002–13). As a result, the rate of heat release from oceanic storage is likely smaller or even turned negative in the last decade. The conjecture based on the results presented in this paper is consistent with the findings by Kosaka and Xie (2013, 2015), Chen and Tung (2014), and Drijfhout et al. (2014), which suggest that the increasing deep-ocean heat uptake is the predominant mechanism for the slowdown of the warming trend in the last decade.
Finally, we wish to comment on potential uncertainties in our decomposition analysis (see the appendix for more details). Besides having no explicit consideration of aerosol forcings, there are two main sources of errors in our decomposition analysis. One is associated with offline radiative transfer model calculations, which use the time mean fields as the inputs of the radiative transfer model instead of the instantaneous fields. In the appendix, we have explicitly identified the offline calculation errors and attributed them to be mainly associated with the use of longtime mean clouds in the offline radiative transfer calculations instead of the longtime mean radiative fluxes calculated from instantaneous clouds, as reported in Wetherald and Manabe (1988), Taylor and Ghan (1992), Kato et al. (2011), Song et al. (2014a,b), and Sejas et al. (2014). The spatial pattern of the offline calculation errors is negatively correlated with the partial radiative flux perturbations associated with clouds, indicating that the offline calculation overestimates the cloud feedback. In the appendix, we have applied the regression analysis to remove the portion of the offline calculation errors that are correlated with the partial radiative flux perturbations associated with clouds. However, the residual offline calculation errors are still relatively pronounced whose global mean is as large as 0.09 K. The second main source of errors is the uncertainty in some atmospheric components, variables, or surface fluxes provided by ERA-Interim. This type of error has little impact on our decomposition analysis as far as the convergence of individual partial temperature changes to the observed temperature change is concerned. However, this type of error could cause an overestimate of the effect of one individual process at the expense of an underestimate of other processes.
We are grateful for the editor (Mingfang Ting) and four anonymous reviewers’ insightful comments and constructive suggestions that led to significant improvements in the presentation. HXM, LYN, and YS are supported by the National Key Research Program of China (2014CB953900), the National Natural Science Foundation of China (41375081), China Special Fund for Meteorological Research in the Public Interest (GYHY201406018), and the Special Funds of Guangdong Province of China (YCJ2013-196). MC is supported by grants from the National Science Foundation (AGS-1354834) and NASA Interdisciplinary Studies Program (NNH12ZDA001N-IDS). YD is supported by the National Science Foundation (AGS-1354402 and AGS-1445956). Calculations for this study were supported by the High-Performance Grid Computing Platform of Sun Yat-sen University and the China National Supercomputer Center in Guangzhou. The datasets used in this study are all freely available on the official websites (http://apps.ecmwf.int/datasets/data/interim-full-mnth/).
Quantifying Errors of Offline Radiative Transfer Calculations
We here wish to provide additional comments on sources of without regard to errors in the ERA-Interim dataset (i.e., taking ERA-Interim as the ground truth). We have checked that there is little difference in the surface upward longwave radiative flux between the offline calculation and ERA-Interim. This is because the surface upward longwave radiative flux more or less solely depends on surface temperature. Therefore, is associated with the errors in the downward longwave radiative fluxes from the atmosphere and the net downward solar radiative fluxes at the surface. The factors that contribute to the errors in the offline calculations, or , are (i) using the decadal mean fields as the inputs of the radiative transfer model in our offline calculations, instead of taking the decadal mean of the radiative heating rates calculated using instantaneous fields, (ii) differences in radiative transfer models, and (iii) excluding changes in other processes, such as CH4 and aerosols, in our offline calculations. Therefore, represents the sum of the errors in all associated with radiative processes [i.e., the terms on the RHS of (2)] and .
Following Sejas et al. (2014), we can identify the dominant terms that are responsible for by comparing the contributions to the difference between the net downward radiative fluxes at the surface from the shortwave and longwave fluxes. Specifically, we evaluate
where the downward arrow indicates downward radiative fluxes, the upward arrow indicates upward radiative fluxes, and the remaining notations are the same as before. The approximation in (b) of (A1) reflects the fact that the upward longwave radiative flux emitted by the surface in the offline calculation is nearly identical to that directly derived from ERA-Interim , as mentioned above. The sum of (a) and (b) measures the errors in the net downward radiative fluxes at the surface by our offline calculations [i.e., the term ].
Displayed in Figs. A1a,c,e are, respectively, the terms (a), (b), and [(a) + (b)] defined in (A1). It is seen that (i) the pattern of (a) is highly (or almost perfectly) correlated with (b) negatively and (ii) the range of numerical values of (a) is larger than that of (b). The collocation of the opposite sign of (a) and (b) with the magnitude of (a) greater than (b) is the trademark of radiative forcing of clouds at the surface, namely that an increase in clouds would reduce solar energy flux (due to more reflection) but increase longwave flux (due to cloud greenhouse effect) reaching to the surface and vice versa with shortwave effect greater than longwave.
To confirm that the source of offline errors is mainly due to using the decadal mean cloud fields as the inputs of the radiative transfer model in our offline calculations (instead of taking the decadal mean of the radiative heating rates calculated using instantaneous cloud fields), we also plot the partial radiative fluxes at the surface due to the decadal mean difference in the cloud fields obtained from our offline calculations. It is seen that the main centers of large values of both shortwave (Fig. A1b) and longwave (Fig. A1d) fluxes as well as the net flux (Fig. A1f), representing the key regions where the decadal changes in cloud fields are pronounced, also appear in Fig. A1a,c,e with the opposite polarity. Therefore, it is mainly the error in cloud fields in our offline calculations that is responsible for. The spatial pattern correlation between Figs. A1a and A1b is as high as −0.82 and that between Figs. A1c and A1d is −0.77. The large negative correlations suggest that the use of time mean clouds fields tend to overestimate both shortwave and longwave component of the cloud feedback, as reported in Sejas et al. (2014). The large negative correlation (−0.55) between Figs. A1e and A1f indicate that the use of time mean clouds fields in the offline radiative model calculation also overestimates the net effect of the cloud feedbacks. The dominant contribution to the offline radiative flux calculations by using mean cloud fields instead of instantaneous cloud fields has been confirmed in Song et al. (2014a,b), who explicitly compared the difference in radiative fluxes between offline calculations using mean cloud fields and the otherwise same calculations but using instantaneous clouds.
To address the issue of overestimating the cloud feedback, we project the shortwave (Fig. A1a) and longwave (Fig. A1c) portions of the offline error fields to their counterparts (Fig. A1b and Fig. A1d) of the partial radiative fluxes at the surface due to changes in clouds derived from our offline calculations according to the following:
where fA1a(x, y), fA1b(x, y), fA1c(x, y), and fA1d(x, y) are the fields shown in Figs. A1a, A1b, A1c, and A1d, respectively. The differences between the original error terms and their counterparts of the projected patterns [i.e., fA1a(x, y) − correction_SW in Fig. A2a for shortwave and fA1c(x, y) − correction_LW in Fig. A2c for longwave] correspond to the new error terms that do not have global pattern correlations with the shortwave and longwave cloud feedbacks. The projected patterns are shown in Fig. A2b (correction_SW) and Fig. A2d (correction_LW) and each of them has a perfect correlation with its counterpart in Fig. A1, by design. We refer to the partial temperature change calculated from the sum of Figs. A2a and A2c using (1) as the partial temperature changes due to offline errors (Fig. A2e). The partial temperature calculated from the sum of Figs. A2b and A2d corresponds to a correction term (Fig. A2f) to the partial temperature change calculated from the original ΔCLD( − ) using (1). In other words, plotted in Fig. 5h, is the sum of Fig. A2f and the partial temperature change calculated from the original ΔCLD( − ).
We note that still have a signal of cloud fields even after the portions that are correlated with the surface radiative energy flux perturbations associated with cloud changes have been removed. There are at least two reasons for it. First, radiative energy flux perturbations induced by cloud changes also exist in the atmosphere. We could not consider it because of the lack of total radiative heating data for the atmosphere in ERA-Interim. Second, only the global patterns that are correlated with the surface radiative energy flux perturbations have been removed from the original offline error term. At regional scales, the (residual) offline errors are still correlated with clouds (e.g., Fig. A2a vs Fig. A2b and Fig. A2c vs Fig. A2d).
It should also be pointed out that may be caused by using mean fields of other variables or parameters (such as water vapor, temperature, and albedo) and the absence of other greenhouse gases and aerosols in our offline radiative transfer model calculations. We here also wish to add that errors in some atmospheric variables or surface turbulent latent/sensible heat fluxes provided by ERA-Interim may affect our decomposition in the same fashion. Unlike introduced in offline calculations, the errors in input fields do not affect since we use the residual to estimate nonradiative heating terms. For example, the errors in estimating ΔQLH by ERA-Interim would be assigned to ΔQOCH with the opposite sign by the residual method, resulting in no net errors in . We have also verified both that the absolute value of ΔQOCH is greater than that of ΔQLH (as well as the sum of ΔQLH and ΔQSH) and that ΔQOCH is negatively correlated with ΔQLH almost everywhere over oceans. This indicates that the energy release from oceans (ΔQOCH > 0) gets into the atmosphere via enhanced evaporation ΔQLH and enhanced thermal radiation (ocean surface warming) and vice versa.
It should be stressed here that the phrase “ associated with X” does not necessarily imply a causal relation. In essence, (1) is just a definition of , which requires the thermal radiative cooling of to balance with (in a linear sense) the energy flux perturbation associated with X, without regard to what causes the change in X in the first place.
The direct effect at the surface of a process is referred to as the corresponding surface energy flux perturbation without considering its associated changes in air and surface temperatures.
The air temperature feedback is akin to the lapse rate feedback, which is defined as the upward thermal radiative energy flux perturbation at the TOA due to the difference between the air temperature and surface temperature changes and commonly used in the PRP method. Here, the air temperature feedback is defined as the downward thermal radiative energy flux perturbation at the surface due to changes in air temperatures, which is one of the terms that directly affect the surface energy balance.