Changyu LI , Jianping HUANG* , Yongli HE , Dongdong LI , and Lei DING Key Laboratory for Semi-Arid Climate Change of the Ministry of Education, College of Atmospheric Sciences,Lanzhou University, Lanzhou, Gansu 730000, China

ABSTRACT

Key words: warming slowdown, ocean heat content, box model, heat transport

1. Introduction

As a result of industrialization, the carbon dioxide in the Earth's atmosphere has increased continuously over the past 100 years, which is considered as the main reason behind global warming (IPCC, 2013). However, the global mean atmospheric temperature (especially in the surface layer) leveled off over the first decade of the 21st century, in contrast to the rapid warming during the late 20th century (Fig.1). This phenomenon has been referred to as an atmospheric warming slowdown or hiatus (e.g., Easterling and Wehner, 2009; Foster and Rahmstorf, 2011; Meehl et al.,2011; England et al., 2014; Guan et al., 2015; Steinman et al., 2015; Huang et al., 2017, He et al., 2018), and has attracted great attention worldwide owing to its ostensible contradiction of the human-induced global warming theory (Otter? et al., 2010).

Fig. 1. Time series of the global mean temperature anomaly relative to the climatology of 1961—90 (a) and trends of global mean temperature for 1983—98 (RWP, red) and 1998—2013 (WSP,blue) (b). The black solid and dashed lines represent the surface temperature derived from GISSTEMP and HadCRUT4,respectively. The colored lines represent the atmospheric temperature derived from ERA-Interim for the near surface(brown), lower atmosphere (green, from the surface to 500 hPa), upper atmosphere (blue, from 500 hPa to 100 hPa), and the whole atmosphere (red, from the surface to 100 hPa).

Improving comprehension of the atmospheric warming slowdown can provide a better projection of future climate(Yao et al., 2017), and the energy balance of the climate system may help us further understand the mechanism behind the phenomenon (Trenberth et al., 2014a). According to observations, greenhouse gases (GHGs) have continued to increase in the 21st century (WMO, 2014), which means that the radiative forcing induced by GHGs is increasing and the energy continued to be trapped in Earth's climate system during the atmospheric warming slowdown. Thus, there must have been a redistribution of heat from the atmosphere to elsewhere in the climate system to fit the slowdown of the atmospheric warming. This points toward the vital role played by the ocean in this slowdown owing to its huge heat capacity,and, more importantly, there is a continuous “global warming” if we focus on the whole climate system rather than just the atmosphere (Cheng et al., 2018). Indeed, many studies have suggested that the oceanic decadal modes, which are tightly linked to the changes in ocean heat content(OHC), may be the main cause of the atmospheric warming slowdown. Kosaka and Xie (2013) concluded that the current warming slowdown is tied specifically to La-Ni?a-like decadal cooling in the Pacific, induced by accelerated trade winds. Some researchers suggest that the phenomenon occurs during the negative phase of the Pacific Decadal Oscillation, or the Interdecadal Pacific Oscillation, which can broadly explain the Pacific OHC changes in the slowdown period (Trenberth and Fasullo, 2014b; Dai et al., 2015). Besides, the Atlantic Multidecadal Oscillation has been proposed as a predictor or key factor for the mean surface temperature over the Northern Hemisphere or globe (Wyatt et al.,2012; Chylek et al., 2014). The atmospheric teleconnection over the Atlantic is generally represented by the North Atlantic Oscillation, which has been suggested to have contributed to the recent warming slowdown (Guan et al., 2015).Huang et al. (2017) suggested that these internal climate variability modes can excite a decadal modulated oscillation,which can explain up to 95% of the temperature variance in the boreal cold season over the continental Northern Hemisphere during the warming slowdown. Furthermore, Liu and Zhou (2017) explored the atmospheric footprint of the warming slowdown, and found it to be coupled with changes in the air—sea latent heat flux. Chen and Tung (2014) also proposed that the warming slowdown may have been caused by oceanic heat transport.

2. Data and methods

Although the warming slowdown has been partially explained well by internal climate variability, its physical mechanism still needs to be further studied. In this study, the characteristics of atmospheric warming slowdown are explored from the global energy balance perspective. We investigate the energy redistribution between the atmosphere and ocean in different latitudes and depths based on analysis of observations, reanalysis products, and the simulations of a coupled atmosphere—ocean box model. We attempt to provide a better understanding of the oceanic energy transport and reveal the fundamental factors that influence this process.

The paper is arranged as follows. Section 2 describes the details of the datasets and the methodology. Section 3 introduces the coupled atmosphere—ocean box model used in this study. The characteristics of the atmospheric temperature anomaly and OHC anomaly related to the warming slowdown period (WSP) and rapid warming period (RWP) are presented in section 4. Section 5 examines the heat transport between the atmosphere and ocean at different depths and latitudes during the WSP and RWP in the box model. Discussion and conclusions are presented in section 6.

2.1. Data sources

The analysis of OHC was based on the Institute of Atmospheric Physics (IAP) gridded ocean temperature dataset and the ECMWF's Ocean Reanalysis System 4 (ORAS4)dataset. The IAP dataset features global coverage of the oceans at 41 vertical levels from the surface to 2000 m, and monthly temporal resolution from 1940 to present (Cheng et al., 2016, 2017). The ORAS4 dataset, which assimilates a variety of observational information in an ocean model (Balmaseda et al., 2013), extends to a depth of 5000 m, but only 0—2000 m is used here to compare with the IAP dataset. The analysis of global surface temperature was based on the GISTEMP and HadCRUT4 datasets. The GISTEMP dataset, which consists of global surface temperature and covers the period from 1880 to the present day (Hansen et al.,2010), was used to calculate the linear trend of the surface temperature, and HadCRUT4 (Morice et al., 2012) was used for comparison with GISTEMP. The analysis of air temperature for the whole atmosphere was based on ERA-Interim(Dee et al., 2011).

2.2. Ensemble Empirical Mode Decomposition method

The Ensemble Empirical Mode Decomposition(EEMD) method was used in our study. It is an adaptive 1D data analysis method that can reflect the nonlinear and nonstationary nature of climate data. Therefore, by using this method, the climate variability can be split into different oscillatory components with intrinsic time scales, including interannual, decadal, and multidecadal extents. In this article, the steps of EEMD were taken from Ji et al. (2014) and Huang et al. (2017), which are described in detail in the Supplementary Material.

2.3. Division of the ocean

In this study, the ocean was divided into five regions:the Pacific Ocean, Atlantic Ocean, Indian Ocean, Southern Ocean, and others. The Southern Ocean was separated from other basins by 35°S, and the boundary between the Indian Ocean and the Pacific Ocean was 120°E. The ranges of each ocean basin are presented in Fig. S1 in Electronic Supplementary Material.

3. Coupled atmosphere-ocean box model

A coupled atmosphere—ocean box model was used in this study to investigate the characteristics of energy transport in the climate system during the RWP and WSP. It takes regional differences by latitude into consideration and divides the atmosphere and ocean into a small number of isothermal reservoirs. The ocean model is based on Schneider and Thompson (1981) and further developed in Danny Harvey and Schneider (1985). Detailed derivations can be found in these publications.

In the model, the atmosphere has only one layer and is divided into two boxes; the ocean has three layers—namely,the upper layer, intermediate layer, and bottom layer—and each layer is also divided into two boxes (Fig. 2). There is heat transport among these boxes of the atmosphere—ocean system. The boxes at low latitude and those at high latitude are linked at 35°N, where the zonal-mean net radiative forcing is close to zero and the poleward atmospheric heat transport is near its peak (Nakamura et al., 1994). The zonalmean net radiative forcing is positive (negative) over the south (north) of 35°N. The equations for this system are as follows:

Here, TA, TU, TIand TBrepresent the temperature of the atmosphere, the upper layer ocean, the intermediate layer ocean,and the bottom layer ocean, respectively; CA, CU, CIand CBare the corresponding heat capacities, which are equal to the depth of each ocean layer multiplied by the heat capacity of a unit of water; A is the net incoming radiation and B is the climate feedback parameter, where A-BTArepresents the net radiative forcing at the top of atmosphere; F is the human-induced radiative forcing; R represents the energy exchange occurring at the surface of the ocean, the perturbation of which is of vital importance in our study and will be discussed in detail below; DUis the depth of the upper layer ocean and Cwis the volumetric heat capacity of water; τ is the turnover time of upper-layer water with intermediate layer water;is the thermohaline mass flux and σ is the area of the ocean; χ is the efficiency of atmospheric heat transport. The subscript “1” represents the low latitudes and the subscript “2” represents the high latitudes. All parameters used in this study are listed in Table 1. Their values are based on previous box-model studies (Schneider and Thompson, 1981; Danny Harvey and Schneider, 1985;Marotzke and Stone, 1995; Krasovskiy and Stone, 1998;Yang et al., 2016).

The equilibrium states of the box model can be obtained by setting all of the time derivations to zero (dTi/dt =0 for each box). For the atmosphere, we obtain:

Note that we do not take the human-induced radiative forcing F into consideration now. According to the values listed in Table 1, we can determine the equilibrium states of the model. The temperature of the atmosphere is 28.2°C in low latitudes (TA1) and 1.2°C in high latitudes (TA2), with a global-averaged temperatureof 14.7 °C. The meridional temperature contrast is equal to 27.0 °C, which is consistent with the result of Yang et al. (2015). However,it should be emphasized that the temperature of some ocean boxes becomes the same when the model reaches the equilibrium state (TU2=TI2=TB1=TB2). Thus, the model is often used to study the response (the perturbation) to external forcing rather than the final state of the climate system.

Now let us proceed to consider the response (or the perturbation) of the system to the human-induced greenhouse effect. We assume that the concentration of CO2 increases linearly with time (Kim et al., 1992), and its radiative forcing can be expressed as

where r = 0.03 W m-2yr-1, which is equivalent to the greenhouse effect caused by a doubling of the CO2concentration in 140 years; H(t) represents the impulse function; t represents the time. The additional heat into the climate system as a result of the greenhouse effect can be calculated using Eq.11 and 12.

Fig. 2. Schematic plot of the coupled atmosphere—ocean box model. H0 is the net energy gains at the top of the atmosphere; R represents the energy exchange occurring at the surface of the ocean; Fo1—Fo6 illustrate qualitatively the heat transports among different ocean boxes; Fa is the meridional atmosphere energy transport.

Table 1. Parameters used in the box model.

There is another key point when we calculate the perturbation of energy transport after taking CO2forcing into consideration, and that is how to deal with the heat uptake in the air—sea interfacePrevious studies have proposed many parameterizatiΔons of the energy transport at the Earth's sea surface (Warren and Schneider, 1979; Danny Harvey and Schneider, 1985; Morantine and Watts, 1990), but it seems quite difficult to reproduce the warming slowdown with these schemes. Besides, although over 90% of this positive energy imbalance is manifested in increased OHC (von Schuckmann et al., 2016), there is still not enough clarity about how the heat is taken into the ocean, due to large uncertainty in observational air—sea heat flux (Josey, 2011). In this case, we handle it on the basis of energy conservation rather than applying these schemes. The purpose of the box model is to investigate the characteristics of energy transport in the climate system during the WSP, so we directly fix the change rate of atmospheric temperature, which is derived from ERA-Interim during the WSP (shown in Fig.3c), in the box model. The changes of corresponding atmosphere heat content can thus be calculated. Due to the energy conservation of the model atmosphere—ocean system,we can now get the OHC change under the warming slowdown scenario by subtracting the atmosphere's heat content change from the total increased energy of the model atmosphere—ocean system. Similarly, the characteristics of heat distribution in the system under the rapid warming scenario can also be investigated with this method.

Assume any variable can be expressed as A+ΔA, where A is the equilibrium state and ΔA is the perturbation. Considering Eqs. (9) and (10), the functions of perturbation in the model system (only the equations of the atmosphere part are shown here) can be written as

where the terms in Eqs. (13) and (14), from left to right, represent the change of atmospheric heat content, the climate feedback, the meridional atmospheric heat transport, the ocean heat uptake, and the radiative forcing, respectively.The ocean heat uptake or the energy transported into the ocean can be easily determined according to the given change in atmospheric temperature.

4. Change in atmospheric temperature and OHC

The characteristics of both the global atmospheric temperature anomaly and OHC anomaly during the WSP are evidently different from those during the RWP.

For the atmosphere, both observational datasets (GISSTEMP and HadCRUT4) show that the warming trends of the global mean surface temperature anomaly during the WSP are obviously smaller than those during the RWP (Fig.1b). The characteristics of the atmospheric temperature anomaly derived from ERA-Interim are consistent with the two observational datasets, though their amplitudes have some differences. It is apparent that the whole atmosphere (from the surface up to 100 hPa) experiences a warming slowdown in the 21st century, compared with the rapid warming in the 1980s and 1990s (Fig. 1; details of atmospheric temperature trend distribution in different pressure levels shown in Fig. S2).Comparing the trends of atmospheric temperature during the two periods, it is evident that the atmosphere accommodates less heat in the WSP than the RWP: while the atmospheric heat content in total increases about 0.62 × 1021J in the period of the warming slowdown, this increase is much less than that in the RWP (approximately 1.37 × 1021J), which reveals the different features of heat transport between the two periods.

The shift in heat transport in the climate system can also be demonstrated by the trends of zonal mean atmospheric temperature (pressure versus latitude) during the RWP and WSP, which are shown in Fig. 3. It can be seen that the vertical distributions of the zonal mean temperature trend are obviously different between the two periods (Figs. 3a and b). The low-latitude atmosphere shows an apparent warming trend from the surface to around the 150 hPa level during the RWP; however, this warming trend weakens or even becomes negative during the WSP. As for the high latitudes,significant cooling is in the upper atmosphere during the RWP, but this cooling weakens over the Arctic regions and even becomes warming over the Antarctic during the WSP.The spatial distribution of the temperature trend in each pressure level can be found in Fig. S2. Overall, in the high latitudes, the atmosphere warms even faster during the WSP than during the RWP; whereas, in the low latitudes, the warming trend during the WSP shows a sharp decline, compared with the warming trend during the RWP (Fig. 3c). These results indicate that the atmospheric temperature trend in the low latitudes contributes most to the warming slowdown,and more importantly, changes in heat transport in the climate system occur during the two periods.

Figure 4 presents the time series of OHC in the global ocean and each ocean basin, from which continuous rapid increases in OHC since the 1980s due to the growth of greenhouse gasses (GHGs) can be seen. This GHG-forced longterm warming can also be found in the reanalysis dataset ORAS4 (Fig. S3), though it suffers from discontinuity due to an upgrading of the oceanic observation system (Palmer et al., 2017). That is, the substantial anthropogenic warming(GHG-forced) signal in the OHC records should be removed before investigating OHC changes related to the slowdown (i.e., natural variability of OHC). Here, we used the EEMD method, proposed by Wu and Huang (2009) and Ji et al. (2014), to separate the decadal variability from the long-term anthropogenic warming. Based on this method,we were able to decompose the atmospheric temperature anomaly and OHC anomaly into different time scales. As shown in Fig. 5, there are two main components in both OHC and atmospheric temperature series: a decadal to multidecadal oscillation (Figs. 5a and b), and the long-term upward trend (Figs. 5c and d). The long-term trend is considered to be mainly contributed by the radiative forcing(Wu et al., 2011; von Schuckmann et al., 2016; Huang et al., 2017). The oscillation component, which represents the natural variability, enhances or suppresses the long-term trend on decadal to multi-decadal time scales. For the OHC,when the oscillation is in a downward phase (Fig. 5a),which means less heat uptake in the ocean, it contributes to an accelerated warming trend in atmospheric temperature,as in the 1990s (Fig. 5b). There is an upward swing in the oscillation of OHC in 1995—96 (Fig. 5a), and thus a subsequent downward phase of the atmospheric temperature oscillation since about 2002 (Fig. 5b), which reduces the longterm upward trend in atmospheric temperature and ultimately results in the recent surface warming slowdown (Fig. 5f).It should be noted that the response of atmospheric temperature to the OHC increase seems to have a lag in time (～5 years; Figs. 5a and b). The downswing of OHC arises in 2006—07 and, after about 5 years, the atmospheric warming slowdown comes to an end with an upswing of atmospheric temperature. Similar results about this lag can also be found in Chen and Tung (2018), when they discussed the relationship between the warming slowdown and Atlantic Overturning Circulation. The sensitivity of the EEMD to the time period selected was also examined, and the results showed the time period has little impact on the pattern of the oscillation component of OHC (Fig. S4).

Fig. 3. Trend of zonal-mean atmospheric temperature during the (a) RWP(1983—98) and (b) WSP (1998—2013) derived from ERA-Interim. Trends for the atmosphere in low and high latitudes (divided by 35°N) in the two periods are shown in (c). Shaded zones are statistically significant at the 5% level.

The characteristics of OHC differ in different ocean basins (Fig. 4), which indicates the contribution of each ocean basin to the warming slowdown may be different(Yao et al., 2017). Therefore, to remove the long-term anthropogenic warming signal, as well as quantify the contribution of each ocean basin to the warming slowdown, we applied the “Trend-Diff method”, which was proposed in Cheng et al. (2018), to each ocean basin with different depths. This method uses the trend differences to reduce the impact of the substantial warming signal, assuming the long-term warming rate is constant in the research period. Although there should be an acceleration of long-term warming (non-constant in rate) due to the accumulation of GHGs,this method shows similar results with analysis of the nonconstant trend proposed in Trenberth and Shea (2006), and is thus acceptable, as mentioned in Cheng et al. (2018). Due to the 5-year prior phase of OHC shown in the EEMD(Figs. 5a and b), the trend differences were calculated between the years 1993—2008 (period of OHC decadal variability related to warming slowdown) and 1979—92 (period of OHC decadal variability related to rapid warming), and the oceans were divided into two layers vertically: the upper layer (0—200 m) and the deeper layer (200—2000 m). The results are shown in Table 2. It can be seen that both the upper and deeper layer in the Atlantic, Pacific and Indian oceans present a positive trend difference, while the Southern Ocean makes a negative contribution. According to the calculations, the largest trend differences exist in the Atlantic, especially its deeper layer (～2.44 × 1021J yr-1; ffgt; 70% in the deeper layer), followed by the Indian and Pacific oceans, with trend differences of ～2.08 × 1021J yr-1and 1.69 × 1021J yr-1,respectively. The global trend of OHC in 1993—2008 in total increases about 5.15 × 1021J yr-1compared with the reference period 1979—92, and the deeper ocean (below 200 m) contributes about 60% of the total. The results indicate that more heat is transported to the deeper ocean, such as the deeper layer of the Atlantic.

Fig. 4. The OHC anomaly (red) averaged in the global ocean, Atlantic Ocean, Southern Ocean, Pacific Ocean, and Indian Ocean, relative to the climatology from 1961 to 1990 based on the IAP dataset. The dashed lines represent the linear trend of OHC during 1983—98 and 1998—2013, respectively.

Fig. 5. EEMD global annual mean OHC from the IAP dataset (red, from the sea surface to 2000 m depth) and atmospheric temperature from GISSTEMP (black) anomalies for the oscillation component (a, b), long-term trend component (c, d), and the trend plus oscillation (e, f). The oscillation component is the sum of IMF (intrinsic mode function) 3, 4 and 5 from the EEMD, and the long-term trend is IMF 6. We mainly focus on the decadal and the multidecadal time scales, so IMF 1 and 2 are ignored, which represents high-frequency oscillation.

For further investigation, EOF analysis of the detrended OHC was also applied to each layer. The results are shown in Fig. 6. After removing the linear trend from the OHC series, the EOF analysis for the OHC of the upper layer ocean showed that the spatial pattern of the first leading mode (EOF1), which explains about 33% of the total variance, is mainly characterized in the tropical Pacific, and has an east—west structure (Fig. 6a). The time series (PC1) associated with EOF1 changes to a negative trend from 1998 when the warming slowdown occurs, while it persists in a positive trend during the 1980s and 1990s (Fig. 6c). This pattern, which has a close relationship (correlation coefficient of ～0.91; significant at the 0.01 level) with Oceanic Ni?o Index, indicates the important role of ENSO in the modulation of OHC in the upper ocean. Consistent with the EOF analysis, the trend differences of OHC in the Pacific Ocean mainly occur within a shallow layer and in the east—west direction (Figs. S5a and b), with an opposite pattern between the RWP and WSP. As for the deeper layer ocean, EOF1 explains about 15% of the total variance, and the positive patterns of OHC changes in the deeper layer are mainly located in the Atlantic and the Indian Ocean (Fig. 6b). PC1 shows a downward trend in the 1980s—90s and a subsequent upward trend at the end of the 20th century (Fig.6d). This indicates that there is a positive OHC anomaly associated with the warming slowdown and a negative OHC anomaly associated with the rapid warming, suggesting more heat is stored in the deeper layer ocean during the WSP than during the RWP (also shown in the trend differences of the OHC anomaly in the Atlantic basin in Figs. S5c and d). Theincrease of Atlantic OHC in EOF1 was also found by Lee et al. (2011) and Chen and Tung (2014).

Table 2. OHC trend differences between two periods (1993—2008 minus 1979—92).

Fig. 6. First EOF mode of the detrended OHC in the upper layers (0—300 m) (a) and the deeper layers (300—2000 m)(b). The PC1 time series associated with the EOF mode are shown in (c) and (d) for the upper layers and deeper layers, respectively. The percentage of variance explained by each mode is indicated in the right-hand corner. The second EOF modes and their associated PC2 time series are not shown. The observational OHC is based on the IAP dataset.

These characteristics of OHC described above reflect the different mechanisms of OHC redistribution in the ocean basins. The thermohaline circulation plays an important role in the Atlantic (Rahmstorf, 2002). Therefore, the OHC anomaly contributed by the Atlantic is mainly in the deeper layer (Table 2, Fig. 6c); whereas, heat redistribution in the upper 200 m in the Pacific is strongly influenced by ENSO (Fig. 6a). The OHC changes in the Indian Ocean are dominated by the Indonesian Throughflow, the increases of which could be responsible for the positive OHC trend differences shown in Table 2 (Lee et al., 2015; Li et al., 2017).

In summary, we found that the whole atmosphere (from the surface to 150-hPa level), particularly the low-latitude atmosphere, shows an obvious warming slowdown during the WSP (1998—2013) compared with that during the RWP(1983—98). After removing the long-term anthropogenic warming signal (forced by GHGs), we found that the global changes of OHC related to the atmospheric warming slowdown increase faster than its changes in the reference period related to rapid warming, particularly the OHC in the deeper layer (200—2000 m) of the Atlantic and Indian Ocean.From the global energy balance perspective, we suggest that the recent global atmosphere warming slowdown is closely coupled with the faster increasing of global OHC.

5. OHC anomaly simulated by the box model

To further investigate the characteristics of energy transfer in the climate system, especially the heat uptake at the air—sea interface and vertical heat transport in the ocean, under different scenarios, a coupled atmosphere—ocean box model was used in this study. Its simple process can help us understand the nature of energy transportation during the WSP.Its basic idea and equations are discussed in detail in section 3.

Firstly, we simulated the changes of energy transport in the climate system when warming slowdown occurs. The trend of atmospheric temperature in the box model was set as 0.03°C (10 yr)-1in the low latitudes and 0.08°C (10 yr)-1in the high latitudes (Table 3) to correspond with results from ERA-Interim (shown in Fig. 3c). The equations of the model are solved numerically using an Euler scheme with 4-day time steps. The model runs for 100 years, and note that all of the results shown below are the perturbation state of the model system. Figure 7 presents the OHC trend in the box model and the observations during the WSP. The result of the bottom layer ocean (2000—4000 m) in the box model is not shown because its temperature changes are relatively small and observations are limited for the ocean below 2000 m. It can be seen from the observation that the OHC in both low- and high-latitude areas shows an obvious positive trend, with an increase of about 5.8 × 1021J yr-1and 3.4 ×1021J yr-1, respectively. For the box model, the simulated patterns of OHC trends in the deeper-layer and upper-layer ocean are quite consistent with the observations (Fig. 7),though their magnitudes are relatively larger. In addition, it is notable that the OHC trend in the deeper layer is larger than its trend in the upper layer, as shown both in the observation and the simulation, indicating more heat storage in the deeper ocean.

Table 3. Atmospheric temperature trends imposed in the box model.

Fig. 7. OHC trend in the observations and simulation of the box model, from the left to the right representing the global,low-latitude and high-latitude ocean, respectively. The red part represents the upper-layer and the blue part the deeper-layer ocean. The observational OHC is based on the IAP dataset.

The warming patterns in the ocean during the WSP have a tight relationship with the features of energy transport in the climate system, which is also demonstrated in the box model. The OHC changes in the upper ocean are affected by the energy taken into the ocean at the sea surface.Figures 8c and e demonstrate the contribution of related processes to the ocean heat uptake in low latitudes and high latitudes during the warming slowdown scenario, respectively.It is shown that the differences are mainly due to changes in features of energy transport in the climate system, especially the atmospheric heat transport and climate feedback.Because the warming rate of atmospheric temperature in high latitudes is larger than that in low latitudes during the WSP, the meridional atmospheric temperature gradient would be reduced, which would lead to less atmospheric heat transport from the low to high latitudes. As a result, the atmospheric heat transport facilitates the ocean heat uptake in low latitudes, but suppresses it in high latitudes (Figs. 8c and e; green lines). Meanwhile, the higher warming rate in high latitudes means more outgoing longwave radiation (outgoing energy). As a result, the climate feedback has a stronger inhibiting effect on the ocean heat uptake in high latitudes than that in low latitudes (Figs. 8c and e; blue lines).Consequently, the heat uptake and the increased OHC of the upper ocean in high latitudes is less than those in low latitudes during the WSP (Fig. 8a, Fig. 7). As for the deeper (intermediate layer) ocean, the seawater mixing and thermohaline circulation impact the deeper-ocean heat storage. During the WSP, the surface ocean heat uptake in the low latitudes increases faster than that in the high latitudes. As a result, compared with the equilibrium states (all results shown here are perturbation states), the seawater mixing brings more heat from the upper-layer to deeper-layer ocean in low latitudes than in high latitudes (Figs. 9c and e; blue lines),due to a faster increase of temperature gradient. The thermohaline circulation wells up and bring cold water from the bottom layer in the low latitudes, while it conversely sinks in the high latitudes, thus facilitating heat storage in high latitudes and suppressing it in low latitudes (Figs. 9c and e; red line). However, this cannot compensate the differences due to the thermohaline, and eventually, under the concurrent effect of these processes, the heat storage and the increased OHC in low latitudes are larger than those in high latitudes(Fig. 9a, Fig. 7).

Fig. 8. Ocean heat uptake at the sea surface (a, b) and the contribution of each term to it (c—f) during the warming slowdown scenario (left-hand panels) and rapid warming scenario (right-hand panels). The dashed,dotted and solid lines in (a, b) represent the low-latitude, high-latitude and global-averaged oceans,respectively. The four terms—AHT (green), CF (blue), AHC (red), and RF (black)—in (c—f) represent the atmospheric heat transport, climate feedback, atmospheric heat content, and radiative forcing, respectively.The ocean heat uptake at the surface is equal to the sum of these terms.

Fig. 9. Heat storage in the deeper-layer (intermediate-layer) ocean (a, b) and the contribution of each term to it (c—f) during the warming slowdown scenario (left-hand panels) and rapid warming scenario (right-hand panels). The dashed, dotted and solid lines in (a, b) represent the low-latitude, high-latitude and globalaveraged oceans, respectively. The red and blue lines in (c—f) represent the contributions of the thermohaline circulation and mixing effect, respectively. Heat storage in the deeper-layer ocean in the model is equal to the sum of the two terms.

Furthermore, we compared the characteristics of OHC changes between the scenarios of rapid warming and warming slowdown in the box model, by fixing different atmospheric warming trends for the two scenarios. A faster atmospheric warming trend, about 0.196°C (10 yr)-1in low latitudes and 0.025°C (10 yr)-1in high latitudes (Table 3),which was also derived from ERA-Interim (shown in Fig.3c), was set in the box model to simulate the rapid warming scenario. The simulated OHC difference between the warming slowdown scenario and the rapid warming scenario is shown in Fig. 10a. It is clear that, compared with the rapid warming scenario, OHC increases faster in both the upper-layer ocean and deeper-layer ocean in the warming slowdown scenario; particularly, the increased OHC in the deeper-layer ocean is greater than that in the upper-layer ocean after several simulated years. The box model ultimately simulates a total OHC trend difference of about 6.2 × 1021J yr-1between the two scenarios, and the deeper ocean contributes over half of it. This trend difference can be attributed to faster ocean heat uptake at the sea surface and faster deeperocean heat storage during the WSP than those during the RWP (Figs. 10b and c). The warming rate of atmospheric temperature in high latitudes is much smaller than that in low latitudes during the RWP, which is opposite to the result during the WSP. Thus, the atmospheric heat transport (Figs. 8d and f; green lines) and climate feedback (Figs. 8d and f; blue lines) show a much stronger inhibiting effect on the ocean heat uptake in low latitudes than in high latitudes, which are also opposite to the results during the WSP. Consequently,the perturbation of ocean heat uptake at the sea surface in low latitudes is negative, which is far less than that in high latitudes during the RWP (Fig. 8b). It indicates that more heat is transported from the ocean to the atmosphere in low latitudes during the RWP. Thus, the mixing effect shows a strong inhibiting effect on the deeper-ocean heat storage in low latitudes (Fig. 9d; blue line) owing to the reduced temperature gradient between the upper-layer and intermediate-layer ocean—opposite to the situation during the WSP. To compensate, the thermohaline circulation in low latitudes subsequently behaves as a positive effect (Fig. 9d; red line). As for the high latitudes, due to the energy conservation, there would be faster heat transport from the atmosphere to the ocean during the RWP than that during the WSP. Thus,both the mixing effect and thermohaline circulation show a moderate facilitating effect on the deeper-ocean heat storage in high latitudes (Fig. 9f). Consequently, for the global ocean (including the low latitudes and high latitudes), the surface ocean heat uptake and the deeper-ocean heat storage during the RWP (Fig. 8b, Fig. 9b, Fig. 10) are both slower than those during the WSP.

6. Summary and discussion

The warming slowdown phenomenon has attracted worldwide attention and has been intensely debated since it was proposed (Easterling and Wehner, 2009; IPCC, 2013). Studies on the warming slowdown have been published since then,and many of them have found a tight link between the warming slowdown and the oceanic energy change (Meehl et al.,2011; Chen and Tung, 2014; Steinman et al., 2015). The importance of the thermohaline circulation has been proposed by several studies, the change of which could influence energy transport in the climate system (Chen and Tung, 2018).However, the physical mechanism of how the ocean energy changes still needs be further explored and clarified.

In this study, we show a distinction between the energy transport in the climate system of the RWP and WSP.Notably, the characteristics of the atmospheric temperature change and the OHC change, as well as the ocean heat uptake, in the low and high latitudes during the WSP are almost opposite to those during the RWP. After removing the substantial anthropogenic warming signal, it can be seen that, during the period of global atmospheric warming slowdown, the global OHC increases faster than that during the period of rapid atmospheric warming. From the global energy balance perspective, we suggest the accelerated ocean heat uptake is tightly coupled with the recent global warming slowdown. Based on the EEMD method, it can be seen that the oscillation component of OHC differs between the WSP and RWP, which indicates the different “drain rate” of energy from the atmosphere to ocean. The upward swing in the oscillation of OHC implies enhancement in ocean heat uptake, thus leading to a downward phase of the atmospheric temperature oscillation, which suppresses the long-term warming trend and finally behaves as a period of atmospheric warming slowdown. Furthermore, EOF analysis showed that the dominant variability of OHC in the upper-layer ocean occurs mainly in the Pacific, while the OHC in the deeper-layer ocean varies mostly in the Atlantic, which reveals the important influence of the thermohaline circulation.

A coupled box model was used in the final part of this study to investigate the characteristics of energy transport in the climate system and confirm the OHC changes under the two scenarios. The model reproduces the accelerated increasing of global OHC under the warming slowdown scenario,the trend of which corresponds well with the results derived from the observational IAP dataset. In particular, the simulated OHC differences presented in Fig. 10 emphasize the important role of the deeper ocean in the warming hiatus, which contributes over half of the total trend differences. Furthermore, the model also revealed that the changes of atmospheric temperature are tightly linked with the characteristics of air—sea energy transport. Comparing the warming slowdown scenario with the rapid warming scenario, it can be seen clearly that the different atmospheric warming patterns in low/high latitudes under the two scenarios are accompanied by distinct outgoing longwave radiation and atmospheric meridional heat transport, as well as other related processes, thus leading to different characteristics of air—sea heat transport and demonstrating the distinct energy transport in the climate system in the two scenarios.

Some limitations of this study should be acknowledged.The real climate system is extremely complex; neither the simple box model we used (which is oversimplified) nor complex coupled general circulation models (which possess considerable uncertainty) can directly reveal the energy transfer in the atmosphere—ocean system. Therefore, the box model was mainly used to explain some phenomena qualitatively, and not for quantitative analysis. In addition, an accurate heat transport at the sea surface on the decadal time scale remains elusive, due to the large uncertainty (～10 W m-2) in observations of air—sea heat flux (Josey, 2011). Besides, the vertical and inter-basin heat redistribution in the ocean is not well understood, and estimates of OHC from oceanic datasets may be uncertain (Wang et al., 2018). The results presented in this paper are primarily derived from the IAP dataset; however, the characteristics of OHC in the deeper layer of the Indian Ocean seem to be different if other datasets are used, and this may contribute to the uncertainty of estimates of oceanic warming in the low latitudes(Fig. S6). Therefore, this study is just preliminary and heuristic. For a further understanding of the energy transfer in the climate system, more accurate observational datasets and reliable climate models should be developed in the future.

Fig. 10. OHC difference (units: J) between the warming slowdown scenario and the rapid warming scenario(former minus the latter) in the box model (a) and comparison of the global ocean heat uptake at the surface(b) and heat storage in the deeper layer ocean (c) between the two scenarios. The dark red and light red in (a)represent the deeper-layer and upper-layer ocean; the blue line and red line in (b) and (c) represent the warming slowdown scenario and rapid warming scenario, respectively.

Acknowledgments.The authors thank Professor Ruixin HUANG of the Woods Hole Oceanographic Institution for his several helpful suggestions. This work was jointly supported by the National Science Foundation of China (Grant Nos. 41521004,41575006 and 41705047), the China 111 project (Grant No.B13045), and the Foundation of the Key Laboratory for Semi-Arid Climate Change of the Ministry of Education in Lanzhou University from the Fundamental Research Funds for the Central Universities (Grant No. lzujbky-2017-bt04).

Electronic supplementary material:Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-019-8281-0.