1 Introduction

The carbon dioxide in the Earth’s atmosphere has increased monotonically over the past 100 years (IPCC 2013) and is considered the main cause of global warming for that period. Although this increase in atmospheric carbon dioxide has been associated with increasing temperatures, the global mean surface air temperature (SAT) has leveled off since 1998, a phenomenon known as a warming hiatus (IPCC 2013; Kosaka and Xie 2013; England et al. 2014; Meehl et al. 2014; Trenberth et al. 2014; Dai et al. 2015; Steinman et al. 2015; Guan et al. 2015) (Fig. 1a). Furthermore, there was another significant warming hiatus during the period of 1940–1975 (Fig. 1a; Beamish et al. 1998; Klyashtorin 1998; Klyashtorin and Lyubushin 2007; Van Loon et al. 2007; Wyatt and Curry 2014), which was a longer hiatus period than the recent warming hiatus. The warming hiatus has attracted great attention worldwide because it appears to contradict the theory of human-induced global warming (IPCC 2013). As the previous projection by the scientific community, the global warming will go on in the first two decades of twenty-first century (IPCC 2007). However, the unexpected warming hiatus declared a failure of the projection and a doubt to the credibility of climate science. Hitherto, many explanations on the warming hiatus have been presented by the scientific community, which increased the knowledge of hiatus process.

Fig. 1
figure 1

a The SAT anomaly time series relative to the 1961–1990 global annual mean for land and ocean (black), land (red), ocean (blue), and SST (green). The marks at the line indicate the recent warming hiatus period. b Same as a but for the NH cold season; only the SAT anomalies over land (red) and ocean (blue) are shown, and the black dashed lines are the linear trend for 2001–2013. Note that the boreal warm season runs from May to September, while the cold season runs from November to the following March

Full size image

Previous studies have suggested that the warming hiatus was mainly induced by internal climate variability (ICV), especially the oceanic ICV. Steinman et al. (2015) concluded that the hiatus has been induced primarily by a strong negative-trend in Pacific multidecadal variability, with only a small contribution from Atlantic variability. However, the Pacific multidecadal variability defined by Steinman et al. (2015) did not follow the classic ICV modes that are generally used to represent the Pacific variability (Kravtsov et al. 2015). The suggested classic Pacific variability related with the recent warming hiatus mainly were the La-Niña-like cooling induced by accelerated trade winds (Kosaka and Xie 2013; England et al. 2014), and a negative phase of the Pacific Decadal Oscillation (PDO), or (more generally) the Interdecadal Pacific Oscillation (IPO) (Dai et al. 2015; Meehl et al. 2014; Trenberth et al. 2014). The Atlantic variability represented by the Atlantic Multidecadal Oscillation (AMO) plays a significant role in the variability of the Northern Hemisphere (NH) or global mean temperature (Kravtsov and Spannagle 2008; Wyatt et al. 2012; Tung and Zhou 2013; Wyatt and Curry 2014). As suggested by Wyatt et al. (2012) and Wyatt and Curry (2014), the AMO signal culminates in an opposite signal in the NH temperature about 30 years later via a sequence of atmospheric and lagged oceanic teleconnections. These studies indicate that the previous warming hiatus (about 1940–1975) was significantly related with the AMO signal. The atmospheric ICV over the Atlantic is generally represented by the North Atlantic Oscillation (NAO) (Huang et al. 1998; Higuchi et al. 1999; Wyatt et al. 2012), which is suggested to have contributed to the recent warming hiatus (Guan et al. 2015). In addition, from an energy balance perspective, the hiatus may have been caused by heat transported to deeper layers of the Atlantic and Southern oceans (Chen and Tung 2014).

Exactly, most of these studies have mainly focused on individual oceanic ICV modes and their influence on global cooling and do not fully describe how these ICV modes have influenced SAT over NH land. The recent cooling over NH land (Cohen et al. 2012) was reversed from the previous enhanced warming, which was probably related to the circulation change over extratropical NH (Huang et al. 2012; Wallace et al. 2012; He et al. 2014; Huang et al. 2016a; Guan et al. 2015; Huang et al. 2016b). As suggested by Wallace et al. (1995) and He et al. (2014), the land-sea thermal contrast can excite feedbacks related with circulation and induces the transition between the pattern of “cold ocean and warm land” (COWL) and “warm ocean and cold land” (WOCL). In addition, some studies suggest the recent cold winters in Eurasian and North America were influenced by the accelerated Arctic warming and sea ice decline (Frolov et al. 2009; Outten and Esau 2012; Zakharov 2013; Cohen et al. 2014; Mori et al. 2014; Screen and Simmonds 2014; Xie et al. 2015). Response to these thermal forcing variation, the circulation over extratropical NH represented by the atmospheric pressure field has also changed, which was suggested to contribute much to the previous enhanced wintertime warming and recent warming hiatus (Wallace et al. 2012; Guan et al. 2015).

This body of literature greatly promoted the understanding of ocean’s effect on the process of hiatus, but the hiatus is also a significant phenomenon over land, especially in mid-to high-latitude of North Hemisphere. So how the cooling over continental NH developed is a key question to understand the recent warming hiatus. In this study, we explore the dynamics of the warming hiatus from a hemispheric perspective and examine the combined effects of ICV modes and Arctic climate on SAT over NH land. The recorded sea surface temperatures (SST) suggest a bias toward colder temperatures in recent decades, and which has induced the corresponding SAT over the oceans to also become colder (Karl et al. 2015). However, the bias in SAT over land was small, especially when averaged over a large scale (Karl et al. 2015; Jones 2016). Our results show that approximately one-half of the recent hiatus in global mean SAT was accounted by cooling over land. We attempt to address how the ICV modes and Arctic amplification induced the hiatus and examine how the IPO or the PDO, the AMO, and other decadal ICV modes influence the hiatus, especially the cooling over the continental NH.

The paper is arranged as follows. The details of the datasets and the methodology used are given in Sect. 2. In Sect. 3, a theoretical model used in this study is introduced. The global pattern of the warming hiatus and the relative contribution to hiatus are presented in Sect. 4. Section 5 discussed the decadal modulated oscillation. Section 6 examined the circulation change responding to thermal forcing during the warming hiatus. Section 7 presents circulation abrupt shift under thermal forcing change. Discussion and conclusion are presented in Sect. 8.

2 Data and methods

2.1 Temperature

The GISTEMP dataset from the National Aeronautics and Space Administration Goddard Institute for Space Studies (GISS) has a resolution of 2° × 2° and covers the period from 1880 to the present day (Hansen et al. 2010). The Extended Reconstructed Sea Surface Temperature (ERSST) dataset is provided by National Oceanic and Atmospheric Administration (NOAA)/OAR/ESRL PSD, from their Web site at http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html. These data are version 4, and begins in January 1854 continuing to the present (Huang et al. 2015). The Nino3.4, PDO, and AMO indexes were calculated using the ERSST data. Furthermore, the Merged Land–Ocean Surface Temperature Analysis (MLOST) (Smith et al. 2008) and HadCRUT4 (Morice et al. 2012) datasets were compared with GISTEMP dataset to assess the uncertainty of SAT data (Fig. S1). The SAT time series for low- and mid- latitudes were very consistent among the three datasets, especially for the period from about 1960 onward (Fig. S1 a, b). The discrepancies among the three datasets were larger for high-latitudes, because the MLOST and HadCRUT4 datasets have more missing data over high-latitudes than GISTEMP (Fig. S1c). Therefore, only the results of the GISTEMP dataset were presented in the paper. The SST data from ERSST v3b (Smith et al. 2008) and COBE (Ishii et al. 2005) datasets were compared with ERSST v4 datasets, and showed a colder bias for both the ERSST v3b and COBE datasets during the recent warming hiatus (Fig. S2). The ERSST v4 datasets have been recently corrected and updated (Karl et al. 2015), so only the results from the ERSST v4 dataset were presented in the paper. For more details about these SAT and SST datasets, see the Supplementary Material.

2.2 Oscillation indexes

The Nino3.4 index was the SST averaged over [5°S–5°N, 170°W–120°W] as defined by Trenberth (1997). The PDO index was first empirical orthogonal function (EOF) pattern and associated time series of monthly SST over field [20°N–70°N, 110°E–100°W] with subtracted global mean SST series (Mantua et al. 1997). The AMO index was SST anomalies averaged over field [0°–70°N, 0°–80°W] with the global mean SST was removed (Trenberth and Shea 2006). Furthermore, the oscillation indexes calculated from other SST datasets and definitions were compared with ours to access the uncertainty of the indexes (Fig. S3). The Nino 3.4 indexes calculated from both the HadISST (Rayner et al. 2003) and COBE datasets were very consistent with the results from the ERSST v4 dataset (Fig. S3a). Two methods—removing the global-warming fingerprint represented by the global mean SST or the linear trends series—used in calculations of PDO and AMO indexes were compared (Fig. S3 b, c). The PDO indexes calculated by both methods from the ERSST v4, HadISST combined with OISST (Reynolds et al. 2007), and COBE datasets were also very consistent (Fig. S3b). The AMO indexes calculated by both methods from the ERSST v4, Kaplan SST (Kaplan et al. 1998), and COBE datasets had some differences at the interannual time scale (Fig. S3c). However, they were very consistent at decadal and multi-decadal time scales (Fig. S3c), and therefore had little influence on our results. In addition, the AMO index calculated by a method that removed the spatio-temporal signal of the global-warming fingerprint was also similar when calculated by the method that removed the linear trend series (Kravtsov and Spannagle 2008). Therefore, we only presented the results from the ERSST v4 dataset and used the global mean SST removal method to calculate the PDO and AMO indexes. More details about these oscillation indexes and SST datasets are provided in the Supplementary Material.

2.3 Reanalysis products

The wind speed, geopotential height data were from NCEP/NCAR Reanalysis I (Kalnay et al. 1996), NCEP-DOE Reanalysis II (Kanamitsu et al. 2002), and ERA-interim (Dee et al. 2011) reanalysis product sets. All these reanalysis products have a resolution of 2.5° × 2.5°. But the ERA-interim has many different resolution choices. The NCEP/NCAR Reanalysis I was from 1948 to present, and the other two products were from 1979 to present. Because the reanalysis products are model-based, the modelled outputs may differ from real observations, especially on decadal-plus time scales (Wyatt and Peters 2012; Kravtsov et al. 2014). To access some of the uncertainties in the results based on these reanalysis products, all the results associated with these three products were presented by the ensemble mean of them. Furthermore, the wind speed and geopotential height data from the NOAA–CIRES Twentieth Century Reanalysis V2c (Compo et al. 2011) and the ERA-20C (Hersbach et al. 2015) reanalysis product sets were used to examine the atmospheric dynamics during the previous warming hiatus. The NOAA–CIRES Twentieth Century Reanalysis V2c products have a resolution of 2° × 2°, and are available from 1851 to 2012 (Compo et al. 2011). The ERA-20C products has a resolution of 2° × 2°, but has many different resolution choices, and is available for the period 1900–2010 (Hersbach et al. 2015). All the results associated with these two reanalysis products were also presented by the ensemble mean of them.

2.4 EEMD decomposition

The ensemble empirical mode decomposition (EEMD) method was used in our study. EEMD is an adaptive one-dimensional data analysis method; it is temporally local and therefore can reflect the nonlinear and nonstationary nature of climate data. Climate variability can be split into different oscillatory components with intrinsic timescales, including interannual, decadal and multidecadal extents. The steps followed in the EEMD method are taken from Ji et al. (2014) and are as follows:

  1. (1)

    Add a white noise series with an amplitude 0.2 times that of the standard deviation of the raw data to the raw data series x(t);

  2. (2.1)

    Set \(x_{1} (t) = x(t)\) and find the maxima and minima of x1(t), and obtain the upper envelope e u (t) and lower envelope e l (t) using cubic splines to connect the maxima and minima, respectively;

  3. (2.2)

    Find the local mean \(m(t) = [e_{u} (t) + e_{l} (t)]/2\), and then determine whether m(t) is close to zero (equivalent to the symmetric of the upper and lower envelopes with respect to the zero line) at any location based on the given criterion;

  4. (2.3)

    If yes, stop the sifting process; otherwise, set \(x_{1} (t) = x(t) - m(t)\) and repeat steps 2.1 to 2.2;

  5. (2.4)

    In this manner, we obtain the first intrinsic mode function (IMF), and by subtracting it from x(t), we obtain a remainder. If the remainder still contains oscillatory components, we again repeat steps 2.1 to 2.2 but with the new x1(t) as the remainder.

So, each time series is decomposed into different IMFs, which can be expressed as:

$$x(t) = \sum\limits_{j = 1}^{n} {C_{j} (t) + R_{n} (t)}$$

where c j (t) represents the jth IMF, which is an amplitude-frequency-modulated oscillatory component, and R n (t) is the residual of data x(t), which is either monotonic or contains only one extreme.

  1. (3)

    Repeat steps 1 and 2 again and again but with different white noise series added each time and obtain the (ensemble) means for corresponding IMFs of the decompositions as the final result.

In the EEMD calculation as outlined in Ji et al. (2014), the noise added to data has an amplitude that is 0.2 times the standard deviation of the raw data, and the ensemble number is 400. The number of IMFs is 6. A MATLAB EMD/EEMD package with the above stoppage criteria and end treatment is available for download at http://rcada.ncu.edu.tw/research1.htm (Wang et al. 2014).

3 Theoretical model

The model used in this study is governed by the nondimensional quasi-geostrophic barotropic vorticity equation for a sphere (Källén 1981; Charney and DeVore 1979):

$$\frac{\partial }{\partial t}\nabla^{2} \psi = J(\nabla^{2} \psi + h, \, \psi ) - 2\frac{\partial \psi }{\partial \lambda } + k\nabla^{2} (\psi^{*} - \psi )$$
(1)

where \(\psi\) is the nondimensional streamfunction, \(\psi^{*}\) is the nondimensional streamfunction forcing, t is nondimensional time, \(\nabla^{2}\) is the nondimensional horizontal Laplacian operator, J is the nondimensional Jacobian operator, \(\lambda\) is the longitude, \(\mu\) is the sine of the latitude \((\mu = \sin (\varphi ))\), h is the nondimensional topographic height parameter and k is the nondimensional dissipation rate.

A purely zonal component and two wave components are taken into account in the low-order model for simplicity; thus, we select three spherical harmonics

$$F_{A} = P_{2}^{0} (\mu ),\quad F_{K} = \cos (2\lambda )P_{3}^{2} (\mu ),\quad \, F_{L} = \sin (2\lambda )P_{3}^{2} (\mu )$$
(2)

and truncate the expansions for \(\psi\), \(\psi^{*}\) and \(h\) as follows:

$$\psi = \psi_{A} F_{A} + \psi_{K} F_{K} + \psi_{L} F_{L}$$
(3)
$$\psi^{*} = \psi_{A}^{*} F_{A} + \psi_{K}^{*} F_{K}$$
(4)
$$h = h_{0} F_{K}$$
(5)

The standardized associated Legendre function is defined as

$$\begin{aligned} P_{n}^{m} (\mu ) & = \sqrt {\frac{(2n + 1)(n - m)!}{{2\left( {n + m} \right)!}}} \frac{{(1 - \mu^{2} )^{{{\raise0.5ex\hbox{$\scriptstyle m$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}}} }}{{2^{n} n!}}\frac{{d^{n + m} }}{{d\mu^{n + m} }}(\mu^{2} - 1)^{n} \\ & \quad (m = 0,1,2, \ldots ,n;\;n = 0,1,2, \ldots ) \\ \end{aligned}$$
(6)

so that

$$\int_{ - 1}^{1} {P_{n}^{m} (\mu )P_{n*}^{m*} (\mu )} d\mu = \left\{ {\begin{array}{*{20}l} 1, \,\hfill & {{\text{m}} = {\text{m* }}\& {\text{ n}} = {\text{n*}}} \hfill \\ 0,\, \hfill & {{\text{m}} \ne {\text{m* or n}} \ne {\text{n*}}} \hfill \\ \end{array} } \right.$$
(7)

The flow of the atmosphere is confined in the Northern Hemisphere in the model (i.e., \(\lambda \in [0,2\pi ], \, \mu \in [ 0 , 1 ]\)). Inserting the expansions (2)–(5) into Eq. (1), we can obtain

$$\left\{ {\begin{array}{*{20}l} {\frac{{d\psi_{A} }}{dt} = - \frac{1}{6}\alpha h_{0} \psi_{L} + k\left( {\psi_{A}^{*} - \psi_{A} } \right)} \hfill \\ {\frac{{d\psi_{K} }}{dt} = \alpha \psi_{A} \psi_{L} + k\left( {\psi_{K}^{*} - \psi_{K} } \right) + \frac{1}{3}\psi_{L} } \hfill \\ {\frac{{d\psi_{L} }}{dt} = - \alpha \psi_{A} \psi_{K} + \frac{1}{6}\alpha h_{0} \psi_{A} - \frac{1}{3}\psi_{K} - k\psi_{L} } \hfill \\ \end{array} } \right.$$
(8)

where

$$\alpha = \frac{105}{128}\sqrt {10}$$

Solving for the equilibrium solution of the system given by Eq. (8), we obtained the following equation using \(\psi_{\text{A}}\)

$$\begin{aligned} & 6k\left( {\psi_{A}^{*} - \psi_{A} } \right)\left[ {(\alpha \psi_{A} + \frac{1}{3})^{2} + k^{2} } \right] \\ & \quad + \alpha h_{0} k\left[ {\psi_{K}^{*} \left( {\alpha \psi_{A} + \frac{1}{3}} \right) - \frac{1}{6}\alpha h_{0} \psi_{A} } \right] = 0 \\ \end{aligned}$$
(9)

and

$$\psi_{L} = \frac{{6k(\psi_{A}^{*} - \psi_{A} )}}{{\alpha h_{0} }}$$
(10)
$$\psi_{K} = \frac{{\alpha \psi_{A} + \frac{1}{3}}}{k}\psi_{L} + \psi_{K}^{*}$$
(11)

The cubic Eq. (9) shows that we can, at most, have three equilibria for certain values of the forcing parameters. The stability of an equilibrium solution is determined by the characteristic values of the three-order coefficient matrix of the linear perturbation equations resulting from (8), i.e., from

$$\left| {\begin{array}{*{20}l} { - (k + \sigma )} \hfill & 0 \hfill & { - \frac{1}{6}\alpha h_{0} } \hfill \\ {\alpha \overline{\psi }_{L} } \hfill & { - (k{ + }\sigma )} \hfill & {\alpha \overline{\psi }_{A} + \frac{1}{3}} \hfill \\ { - \alpha \overline{\psi }_{K} + \frac{1}{6}\alpha h_{0} } \hfill & { - \alpha \overline{\psi }_{A} - \frac{1}{3}} \hfill & { - (k{ + }\sigma )} \hfill \\ \end{array} } \right| = 0$$
(12)

4 Global distribution of the SAT trend and its relative contribution

As shown in Fig. 1, the temperature changes have varied depending on the season, land/sea, and the hemisphere (Wu et al. 2011; Huang et al. 2012; Ji et al. 2014). The change of mean SAT over land is much faster than the ocean, especially in the previous decades before the recent warming hiatus. Although both the land and ocean show the hiatus during recent decades, it shows that the land SAT have been to a higher level after the enhanced warming during 80–90s in last century. Compared with land, the SAT and SST of ocean maintain a relatively moderate warming rate before the recent hiatus period (Fig. 1a). The SAT over NH exhibits the largest warming trend before the recent hiatus, and it also exhibits the most obvious hiatus in recent decades (Fig. 1b). The SAT of land over NH takes an obvious cooling trend in the boreal cold season during recent decades, accompanied by a weak increase of mean SAT of ocean (Fig. 1b).

The linear trends of SAT in the boreal warm and cold season for 2001–2013 are presented in Fig. 2. The apparent differences in the spatial patterns of SAT changes between the warm and cold seasons during the recent hiatus are shown in Fig. 2. There are more cooling regions and larger cooling trends in the boreal cold season than in the warm season during the recent hiatus. The cooling centers in warm season are mainly located in small regions over Eurasia continents and some regions over Pacific (Fig. 2a). Most regions of NH land still show the warming trend in warm season. However, large parts of NH land exhibit the obvious cooling trend in cold season (Fig. 2b). Besides, the cooling over the equatorial central and eastern Pacific (ECEP) (Kosaka and Xie 2013; Dai et al. 2015) is also more significant in cold season (Fig. 2b). Therefore, we confine our analysis to SAT changes in the NH during the boreal cold season.

Fig. 2
figure 2

a The linear trend in the boreal warm season SAT for 2001–2013. The blue (red) number in left (right) corner is the global mean cooling (warming) trend in °C per 13 years. b Same as a but for the boreal cold season

Full size image

To quantify the spatial features of the recent hiatus, we calculated the areal weighted relative contribution of the SAT for each grid cell to the global mean cooling or warming trend during the boreal cold season (Fig. 3). The relative contribution identified in Fig. 3 was calculated using Eq. (13). The global mean cooling SAT trend was calculated as \(\sum\nolimits_{i}^{N} {K_{i} \omega_{i} }\), and N is the number of the whole grids whose SAT trends were negative. The same procedure is used to calculate the global mean warming trend. For Table 1, the relative contribution of the regional mean cooling trend to global mean cooling is the sum of the relative contribution of all the grids in a region (Fig. 3). The relative contributions of the regional mean cooling trends to NH, land or ocean were the ratio of the regional mean relative contributions to the average relative contributions over NH, land, and ocean, respectively.

$$r_{i} = \frac{{K_{i} \omega_{i} }}{{\sum\nolimits_{i}^{N} {K_{i} \omega_{i} } }}$$
(13)

where K i is the SAT trend of grid i, \(\omega_{i} = \cos (\theta_{i} \cdot \pi /180.0)\), θ i is the latitude of the grid, and N is the number of grids whose SAT trend is negative (cooling) or positive (warming), e.g., N is the number of the whole grids whose SAT trend is negative when the SAT trend of the grid i is negative.

Fig. 3
figure 3

Areal weighted relative contribution of the SAT trend in each grid cell to the global mean warming/cooling trend during the boreal cold season for 2001–2013. The relative contribution was calculated using Eq. (13). The blue (red) number in left (right) corner is the global mean cooling (warming) trend in °C per 13 years. The black squares represent the North American (NA), the Eurasian continents (EUA), and the Equatorial Central and Eastern Pacific region (ECEP), respectively

Full size image
Table 1 The relative contribution (%) of the regional mean cooling trend to the averaged Northern Hemisphere (NH), global, land and ocean cooling trend, respectively, for the boreal cold season
Full size table

Figure 3 shows that the cooling centers of the recent warming hiatus are mainly located above the Eurasian (EUA) and North American (NA) continents as well as in the ECEP regions. Meanwhile, continued warming centers have been identified over the Arctic, in some regions around the Mediterranean, and over the North Pacific Ocean. The areal weighted global mean warming trend was 0.47 °C per 13 years, which is slightly lower than the global cooling trend of −0.50 °C per 13 years. The EUA and NA continents contributed 18 % (60 %) and 43 % (25 %) to NH (terrestrial) mean cooling, respectively (Table 1). Although the oceanic ICV played a major role in explaining the hiatus (Steinman et al. 2015), the local cooling contribution by the ocean was only 52 % (Table 1), which indicates that 48 % is explained by terrestrial cooling, a figure that should not be ignored.

5 Decadal modulated oscillation

We refer to the ICV modulated components of SAT variation on the decadal to multi-decadal scale as the decadal modulated oscillation (DMO). To examine the DMO components of SAT, we used the EEMD method, which can split the evolution of SAT into long-term trends and oscillation components (Wu et al. 2011; Ji et al. 2014). The long-term warming trend is related primarily to radiative forcing (Wu et al. 2011; Guan et al. 2015), and the oscillation component is thought to be induced mainly by ICV (Wu et al. 2011; Wallace et al. 2012; Guan et al. 2015). The DMO components of the global annual mean SAT and NH cold season mean SAT are presented in Fig. 4a, b respectively. The DMO component is the sum of IMF 3, 4, 5 from the EEMD, and the long-term trend is the IMF 6. The DMO enhances or suppresses the long-term trend on decadal to multi-decadal timescale. The DMO’s magnitude is comparable to that of the long-term trend at the decadal scale (Fig. 4). When the DMO is in an upward (warming) phase, it contributes to an accelerated warming trend, as in 1985–1998. It appears that there is a downward swing in the DMO occurring since about 2000, which has balanced or reduced the radiative forced warming and resulted in the recent global warming hiatus. During the boreal cold season, the DMO’s cooling trend is greater than the long-term warming trend, which is shown as obvious cooling in the mean SAT over the NH (Fig. 4b).

Fig. 4
figure 4

a The EEMD decomposed global annual mean SAT anomalies for the long-term trend component (red), DMO component (blue), and the trend plus DMO (red-blue). b Same as a but for NH cold season

Full size image

To examine the modulated effect of the oceanic ICV modes on the NH SAT during the cold season, the multiple linear regression (MLR) and EEMD methods were used to create the DMO index. First, the EEMD method was applied to the indices of three classic oceanic modes: the PDO, Nino3.4, AMO, and Arctic Oscillation (AO), respectively. Then, the DMO component of the NH SAT during the cold season (Fig. 4b) was regressed with standardized IMF 3, 4, 5, for each index using a stepwise MLR. The regressed DMO component in Fig. 5 is expressed in Eq. (14). The regression-based approximation of DMO component for NH SAT using PDO, Nino3.4, AMO, and Arctic Oscillation (AO) indexes can explain 88 % of its variance during the boreal cold season (Fig. 5). The three classic oceanic ICV modes—PDO, Nino3.4 and AMO—contributed 22, 18 and 56 %, respectively.

Fig. 5
figure 5

a The DMO component of the NH SAT anomalies during the cold season (black line) and its regression using the decadal variability of the PDO, Nino3.4, AMO, and AO indices (bar). b The contribution of PDO (blue), Nino3.4 (grey), and AMO (red) to the regressed DMO component in a. Note that the regression in a explained 88 % of the variance in corresponding SAT anomalies, and the PDO, Nino3.4, AMO, and AO contributed 22, 18, 56, and 4 %, respectively. As the contribution of AO was little, the line that represents AO’s contribution was not shown in b

Full size image

Moreover, a similar method was applied to the temporal coefficients in the first three empirical orthogonal function (EOF) modes of the de-trended global SST (Fig. 7) and the AO index. The results are shown in Fig. 6 and expressed in Eqs. (15, 16). To obtain the results shown in Fig. 6a, the EOF analysis was first applied to the linear de-trended global SST for the boreal cold season (Fig. 7). Then, methods similar to those used to derive Fig. 5 were applied to the temporal coefficients of the first three EOF modes and the AO index. The results in Fig. 6c were derived by employing similar methods but replaced the EEMD method with the 11-year running mean method, which is expressed using Eq. (16). The first EOF mode of the de-trended global SST closely correlates with the PDO and Nino3.4 indices, with correlation coefficients of 0.69 and 0.93, respectively (Table 2). The second and third SST’s EOF modes closely correlate with the AMO and present the correlation coefficients of 0.61 and 0.64, respectively (Table 2). The relative contribution of each classic oceanic ICV mode or SST’s EOF mode and AO to the regressed DMO component was calculated following the method described by Huang and Yi (1991).

Fig. 6
figure 6

a The DMO component for the NH cold season mean SAT anomalies (black line), and its regression using decadal variability of EOF modes for global SST, and AO indexes (bar). b The contribution of first (blue), second (orange), and third (red) EOF modes to the regressed DMO component in a. c, d The same as a, b but using the 11-year running mean instead of the EEMD method. Note that the regression in a, c explained 93 and 95 % of the variance in corresponding SAT anomalies, respectively. The EOF1, EOF2, EOF3, and AO contributed 33, 45, 15, and 7 % in b, and contributed 11, 42, 43, and 4 % in d, respectively. The lines that represent AO’s contributions were not shown in b, d

Full size image
Fig. 7
figure 7

a, b The spatial pattern and temporal coefficient of the first EOF mode for the linear de-trended global SST anomaly over the 1900–2013 period with respect to the 1961–1990 reference period for the boreal cold season, respectively. c–f The same as a, b but for the second and third EOF modes, respectively. Note that the three black rectangles in a indicate the ocean regions where the PDO, AMO and Nino3.4 indices are defined

Full size image
Table 2 The correlation coefficients between the temporal coefficients in Fig. 7b, d, f and the PDO, Nino3.4, and AMO indices
Full size table
$$\begin{aligned} DMO & = 0.008 + 0.032PDO_{3} + 0.094PDO_{5} + 0.063Nino_{4} \\ & \quad + 0.015Nino_{5} + 0.206AMO_{5} + 0.016AO_{3} - 0.042AO_{4} \\ \end{aligned}$$
(14)
$$\begin{aligned} DMO & = 0.008 + 0.018EOF1_{3} + 0.012EOF1_{4} + 0.104EOF1_{5} \\ & \quad + 0.063EOF2_{4} + 0.095EOF2_{5} + 0.031EOF3_{3} \\ & \quad + 0.05EOF3_{5} + 0.024AO_{3} - 0.037AO_{4} \\ \end{aligned}$$
(15)
$$DMO = 0.002 + 0.044EOF1 + 0.092EOF2 + 0.114EOF3 + 0.029AO$$
(16)

Figures 5 and 6 suggest that the Atlantic Ocean’s decadal variability makes the largest contribution to the DMO, whereas the North Pacific and tropical Pacific’s decadal variability have the relative weaker effect. However, the sum of North Pacific and tropical Pacific’s contribution is comparable with Atlantic’s contribution. The North Pacific’s contribution to SAT variability of NH is larger than tropical Pacific for decadal timescale. Nevertheless, the recent warming hiatus was mainly the contribution of decadal variability in the Pacific Ocean (Kosaka and Xie 2013; England et al. 2014; Meehl et al. 2014; Trenberth et al. 2014; Dai et al. 2015). It seems that the AO make little contribution to the DMO for the decadal scale as suggested by the results, which showed the relative contributions of AO were 4, 7, and 4 % in Figs. 5, 6b, d, respectively. However, there are many studies indicated that the Arctic sea ice and enhanced Arctic warming influenced the SAT variability over mid-latitude NH continents significantly in last decades (Outten and Esau 2012; Cohen et al. 2014; Mori et al. 2014; Screen and Simmonds 2014). So the Arctic’s influence on recent warming hiatus will be examined further in next section. As most parts of NH are land, how the suggested oceanic modulation influence the SAT variability of NH land will be examined in next section as well.

6 Thermal forced atmospheric circulation change

Atmospheric circulation change is the dominant factor affecting the decadal variability of terrestrial SAT over the NH during the boreal cold season (Van Loon et al. 2007; Wallace et al. 2012; Guan et al. 2015), and oceanic ICV influences that atmospheric circulation via oceanic thermal forcing (Latif and Barnett 1994; Rodwell et al. 1999). To study the thermally forced atmospheric circulation change, the surface thermal forcing fields were examined firstly. Figure 8 shows the regional mean SAT time series with 11-year running mean over different latitudinal zones (Fig. 8a), and two warming and two cooling centers in the mid-latitude (Fig. 8b), respectively. Although warming over the mid- to high NH latitudes has paused over the past decades for boreal cold season, Arctic warming accelerated during the same period (Fig. 8a, b), which is consistent with Cohen et al. (2014) and Wyatt and Curry (2014). The obvious enhanced Arctic warming during the period of warming hiatus was also suggested by the Fig. 9a, which shows the zonal mean SAT time series of NH. The zonal mean SAT over Arctic regions shows the persistent increasing trends, but the mid-latitude regions have no significant trend during recent hiatus (Fig. 9a). Despite the obvious differences in varied latitudinal zones, the SAT change also shows a difference over varied longitudinal zones (Fig. 8b). As shown in Fig. 8b, the SAT over NA and EUA show significant decreasing trends, but the western EUA (regions around the Mediterranean) and Pacific regions still show the warming trend. The meridional mean SAT time series shown in Fig. 9b also confirm these results. The specific locations of these four regions are presented in Fig. 9c. The red and orange filled regions are two warming centers during the warming hiatus, and the green and blue filled regions are two cooling centers (Fig. 9c). The regional mean SAT in Fig. 8b was calculated over the color filled regions in Fig. 9c. The meridional mean SAT in Fig. 9b also was calculated over the color filled regions, but the skyblue filled regions of NA were not used in the calculation of the meridional mean.

Fig. 8
figure 8

a The 11-year running mean SAT time series for the boreal cold season averaged across 10°–30°N (green), 30°–50°N (blue), 50°–70°N (black), and 70°–90°N (red) latitudinal zones. b The 11-year running mean SAT time series for the boreal cold season averaged over North America (NA, black), Eurasia (EUA, blue), Western EUA (green), and the Pacific (red)

Full size image
Fig. 9
figure 9

a The 11-year running zonal mean SAT time series for the boreal cold season. b The 11-year running meridional mean SAT time series for the boreal cold season. The dashed lines indicate the edges of the NA and EUA continents, and the western EUA regions. c The distribution of the regions. The colored zones represent the regions in a, b. Note that both the sky blue and blue regions represent NA in Fig. 8b, but only blue regions are used to calculate the meridional mean in b

Full size image

Figure 10 shows the meridional thermal forcing (MTF) and zonal thermal forcing (ZTF) evolution over NH for boreal cold season, and the corresponding atmospheric circulation response. The MTF represents the meridional temperature gradients between the mid- and high-latitudes. The ZTF represents the asymmetry in temperatures between the extratropical large-scale warm and cold zones in the zonal direction, especially the land-sea thermal contrast. As shown in Fig. 10a, the different SAT changes over the mid-latitudes and the Arctic induced a rapidly weakened meridional temperature gradient over the NH extratropics, which induced an accelerated weakening of the extratropical MTF during the recent warming hiatus. However, the asymmetric ZTF, which is related to the land-sea thermal contrast, has heightened during the hiatus, shifting from a decreasing trend to a significant increasing trend (Fig. 10b). This enhanced extratropical ZTF was induced by the rapidly increasing discrepancy between the two cooling centers over the EUA and NA continents and the two warming centers over the western EUA and the North Pacific (also see Figs. 3, 8, 9).

Fig. 10
figure 10

a The meridional temperature gradient represented by the difference between the regional mean SAT over the 30°–50°N and 70°–90°N latitudinal zones. b Zonal asymmetry thermal forcing represented by the SAT difference between North America (NA), eastern Eurasia (EUA) and Pacific, western EUA. The solid line with arrows in a, b is the linear trend for 2001–2013. The difference in the time series is the anomaly relative to the 1981–2010 reference period. c The difference in the mean geopotential height (GPH) over the 500-hPa layer for 2002–2013 and 1990–2001. The vector indicates the corresponding wind field difference. d The difference in the mean meridional wind speeds over the 500-hPa level for 2002–2013 and 1990–2001. Note that the minimum latitude of c, d is 20°N, and the dashed circle lines are 30°N and 60°N. Only the wind vectors over the regions in which the meridional and zonal wind speed change that calculated from three reanalysis products (ERA-interim, NCEP I and NCEP II) had the same sign are presented in c. The black dots in d indicate the regions in which the signs for the results from the three reanalysis products did not agree

Full size image

The weaker MTF followed a pattern of “warm Arctic and cold land” (WACL) (Overland et al. 2011; Mori et al. 2014) (see Fig. 3), which corresponded to an obvious positive anomaly at 500-hPa geopotential height (GPH) over the Arctic (Fig. 10c). As shown in Fig. 11a, c, e, the positive GPH anomalies were also shown in other three height levels, the 200- and 850-hPa GPH, and sea level pressure (SLP) fields, which indicate the low pressure system over Arctic–referred as polar vortex in upper troposphere was quite weak in the recent warming hiatus. In addition, the mid-latitudes over the eastern EUA, NA, the Atlantic, and the Mediterranean all exhibited a negative 500-hPa GPH anomaly (Fig. 10c). Similar situations occurred above other three height levels (Fig. 11a, c, e). As shown in Figs. 10c and 11a, c, e, the changes in the GPH or SLP fields over the Arctic and the mid-latitudes are asymmetric. These asymmetries induced the weaker polar vortex and slower westerly winds over the high-latitudes as shown in the wind fields (Figs. 10c, 11a, c, e). Besides, the weakened meridional GPH gradient and slower westerly winds will induce a slower speeds and a larger amplitude in planetary waves (Screen and Simmonds 2014), which will facilitate the invasion of extremely cold air at the mid-latitudes, leading to cooling over the EUA and NA during the recent hiatus. Furthermore, during the previous warming hiatus period (1940–1975), the Arctic experienced cooling trends, and the corresponding polar vortex was strengthened (Fig. S4). Therefore, the Eurasian extratropical regions (poleward of 45°N) showed warming trends because of the stronger MTF in high-latitudes (Fig. S4). However, the Eurasian mid- and low-latitudes (southward of 45°N) experienced cooling trends, because of the weaker MTF over these regions (Fig. S4, Table 3). The North America high-latitudes showed significant cooling trends during the previous warming hiatus period, because the stronger polar vortex expanded to high-latitudes over North America (Fig. S4b).

Fig. 11
figure 11

a, b Same as Fig. 10c, d but for the 200-hPa level. c–f Represent 850-hPa and sea level, respectively. Note that the wind speed in e, f was measured at 10 m height above the surface

Full size image
Table 3 The flow patterns in the boreal cold seasons for three periods: the previous warming hiatus (1940–1975), the previous accelerated warming (1980–2001), and the recent warming hiatus (2001–2013)
Full size table

The enhanced ZTF generally followed a “warm ocean and cold land” (WOCL) pattern (Wallace et al. 1995; He et al. 2014), but one warm center was over the western EUA regions rather than over the Atlantic Ocean (see Fig. 3). The North Pacific had one of the most obvious warming centers (see Fig. 3), which corresponded to a negative PDO phase and a positive 500-hPa GPH anomaly (Fig. 10c). As shown in Fig. 11a, c, e, the positive anomalies over North Pacific are also shown in other heights. As the Aleutian Low is defined in SLP fields, it indicates a quite weaker Aleutian Low in recent warming hiatus (Fig. 11e). The positive GPH or SLP anomalies over the North Pacific enhanced the southward flow of air over NA (Figs. 10d, 11b, d, f), which caused the cooling across NA (Fig. 3). Furthermore, the positive GPH anomalies over the North Pacific were also presented during the previous warming hiatus period, which also contributed to the NA’s cooling during that period (Fig. S4, Fig. S6 a, c, e). Previous studies (e.g. Beamish et al. 1998; Klyashtorin 1998; Klyashtorin and Lyubushin 2007; Van Loon et al. 2007; Wyatt and Curry 2014) also suggested the evident positive GPH anomalies and the weaker Aleutian Low over the North Pacific during the previous warming hiatus period. In addition, these studies examined the atmospheric circulation changes throughout a period (about 1915–1975) that contains the entire evolutions form the accelerated warming (about 1915–1940) to the warming hiatus (1940–1975). By examine the Pacific Circulation Index and Atmospheric Circulation Index that represent the large-scale winds and atmospheric-mass-transfer, and the Aleutian Low Pressure Index, they suggested the GPH and Aleutian Low over the North Pacific changed significantly during both the accelerated warming (about 1915–1940) and the previous warming hiatus period (Vangenheim 1940; Girs 1971a, b; Beamish et al. 1998; Klyashtorin 1998; Klyashtorin and Lyubushin 2007; Van Loon et al. 2007; Wyatt and Curry 2014). The cyclone-like negative 500-hPa GPH anomaly over the eastern EUA also enhanced the flow of air southward over EUA, and cold Arctic air was carried by the southward airflow and the cyclone anomaly to large parts of the eastern EUA (Fig. 10c, d), which induced obvious cooling in that areas (Fig. 3). Furthermore, the similar cyclone anomaly also occurred during the previous warming hiatus period, but its location was over the Eurasian mid- and low-latitudes (Fig. S4). In addition, the opposite SLP or GPH anomaly at 850-hPa between Greenland and the mid-latitude Atlantic indicated a negative phase NAO-like pattern (Fig. 11d, f), which also contributed to the recent warming hiatus (Guan et al. 2015). However, there was a positive NAO-like pattern during the previous warming hiatus period, which probably contributed to the warming trends in extratropical regions over the Eurasian continent during that period (Fig. S4).

7 Abrupt shift in thermally forced atmospheric circulation

A theoretical dynamic model (Källén 1981; Charney and DeVore 1979) was used to confirm the responses of circulation patterns to the MTF and ZTF changes over the extratropical Northern Hemisphere suggested above. Charney and DeVore (1979) found multiple equilibrium states for a given topographic and thermal forcing using this model. Among the multiple equilibrium states, two stable equilibrium states vary distinctly: one state is a low-index flow with a strong wave component and a weaker zonal component (Fig. 12e), whereas the other is a high-index flow with a weaker wave component and a stronger zonal component (Fig. 12c). The low-index equilibrium is a colder state with more blocking events over extratropical NH, but the high-index equilibrium is a warmer state with less blocking (Charney and DeVore 1979). In addition, there is an abrupt shift between the low-index and high-index equilibrium with the change of MTF and ZTF over the extratropical NH (Charney and DeVore 1979) (see Fig. 12a).

Fig. 12
figure 12

a The equilibrium bifurcation associated with the change both in the parameter of MTF (\(\psi_{A}^{*}\)) and the parameter of ZTF (\(\psi_{K}^{*}\)). Note that \(\psi_{A}^{*}\) is the opposite of MTF, and \(\psi_{K}^{*}\) is same-phase with ZTF. All points in the orange and skyblue surface indicate high-index and low-index equilibrium states, respectively. All the points with unstable equilibria are not shown. The black line is a possible track of equilibrium states that correspond to the change of both \(\psi_{A}^{*}\) and \(\psi_{K}^{*}\). The thin black lines over the skyblue and orange surfaces in a represent the streamfunction field corresponding to the two black points. b, c The heating field and streamfunction field corresponding to the black point on the orange surface, respectively. d, e The same as b, c but for the black point on the skyblue surface. The colored background shows the topography distribution in the model, and the warm tone indicates the “land”; the cool tone indicates the “ocean”. “C” and “W” in b, d represent cold and warm, respectively. “H” and “L” in c, e represent high and low pressure, respectively. The contour intervals are all 0.1 units. Parameter values are: k = 0.05, \(h_{0} = - 0.15\), \(\psi_{A}^{*} = - 0.42\), \(\psi_{K}^{*} = 0.2\) in b and \(\psi_{A}^{*} = - 0.32\), \(\psi_{K}^{*} = 0.3\) in d

Full size image

To evaluate the thermally forced atmospheric circulation changes that occur with the simultaneous changes of both MTF and ZTF over the extratropical NH, we first separately examine MTF and ZTF to assess each pattern’s individual impact profile. The equilibrium bifurcation associated with the changes to the parameter of MTF (\(\psi_{A}^{*}\); \(\psi_{A}^{*}\) is opposite to MTF) is shown in Fig. 13. \(\psi_{A}^{*} < 0\) indicates that thermal forcing decreases from low to high latitudes. \(\psi_{A} < 0\) indicates that the large-scale flow is a westerly wind, whereas \(\psi_{A} > 0\) indicates that the large-scale flow is an easterly wind. When the absolute value of \(\psi_{A}\) is greater than the absolute value of \(\psi_{K}\) and \(\psi_{L}\), the zonal component is stronger than the wave component. We can see that when the MTF is large (\(\psi_{A}^{*} < - 0.42\)), the large-scale flow has only one possible equilibrium state, which has a large zonal component (Fig. 13a) and small wave components (Fig. 13b, c) and is a high-index flow. When the MTF is small (\(- 0.3 < \psi_{A}^{*} < 0\)), the large-scale flow has only one possible equilibrium state, which has a small zonal component (Fig. 13a) and large wave components (Fig. 13b, c) and is a low-index flow. When the MTF falls between the two cases outlined above (\(- 0.42 < \psi_{A}^{*} < - 0.3\)), there are two stable equilibrium solutions; one is a high-index flow, and the other is a low-index flow. Which stable flow regime the atmosphere is attracted to depend on the initial state of the flow and the way in which the MTF changes over time (Shabbar et al. 2001).

Fig. 13
figure 13

The equilibrium bifurcation associated with the changes in the parameter of MTF (\(\psi_{A}^{*}\)). The ordinate gives the equilibrium values \(\psi_{A}\) in a, \(\psi_{K}\) in b and \(\psi_{L}\) in c, respectively. A numerical evaluation of the eigenvalues shows stability properties as indicated by the black (stable) and red dashed (unstable) lines. Parameter values are: k = 0.05, \(h_{0} = - 0.15\), and \(\psi_{K}^{*} = 0.2\)

Full size image

The equilibrium bifurcation associated with the changes to the parameter of ZTF (\(\psi_{K}^{*}\)) is shown in Fig. 14. \(\psi_{K}^{*} > 0\) indicates that the ocean’s warming is stronger than the land’s, whereas \(\psi_{K}^{*} < 0\) means that the ocean’s cooling is stronger than the land’s. When the ZTF with a “warm ocean and cold land” pattern (WOCL, Fig. 12d) is large (\(\psi_{K}^{*} > 0.8\)), the large-scale flow has only one possible equilibrium state, which has a small zonal component (Fig. 14a) and large wave components (Fig. 14b, c) and is a low-index flow. When the ZTF is small (\(\psi_{K}^{*} < 0.2\)), or if the ZTF exhibits a “cold ocean and warm land” pattern (COWL, \(\psi_{K}^{*} < 0\)), the large-scale flow has only one possible stable equilibrium state, which has a large zonal component (Fig. 14a) and small wave components (Fig. 14b, c) and is a high-index flow. When the ZTF falls between the two cases outlined above (\(0.2 < \psi_{K}^{*} < 0.8\)), there are also two stable equilibrium solutions.

Fig. 14
figure 14

Same as Fig. 13, except that the abscissa all give the parameter of ZTF (\(\psi_{K}^{*}\)) and have a parameter value \(\psi_{A}^{*} = - 0.42\)

Full size image

Figures 13 and 14 suggest that both a decrease in MTF and an increase in ZTF favored the occurrence of a low-index flow. The equilibrium bifurcation associated with the simultaneous change in the parameters of MTF (\(\psi_{A}^{*}\)) and ZTF (\(\psi_{K}^{*}\)) is shown in Fig. 12a.

In addition, blocking corresponding to the low-index flow occurred more frequently over the EUA and NA during the recent hiatus (Fig. 15). Figure 15a shows the decadal variability of sector blocking frequency over EUA (50°E–140°E) and NA (170°E–80°W) regions. Figure 15b compares the 5-year mean difference of local blocking frequency between the beginning of hiatus (1999–2003) and present (2009–2013). The sector and local blocking was calculated following the procedure outlined in D’Andrea et al. (1998). As shown in Fig. 15a, the blocking frequency over EUA increased by 8 % in boreal cold season since 2001, which is nearly one-third of its climatology, which is about 27 %. Such a severe increase of the blocking events will definitely induce the EUA’s cooling during recent hiatus. The blocking frequency over NA also increased by 6 %, which led the NA’s cooling too (Fig. 15a). Besides, the blocking frequency over EUA and NA decreased significantly during the enhance warming period (Huang et al. 2012, He et al. 2014) before the recent warming hiatus. The frequency of blocking events also contributed to the enhanced warming. The increase frequency of the local blocking over the longitude zones of EUA and NA also indicated the more frequent blocking events during the recent hiatus (Fig. 15b). Over all, the Arctic amplification (Cohen et al. 2014; Mori et al. 2014), land-sea thermal contrast (Wallace et al. 1995; He et al. 2014), the negative phase of PDO (England et al. 2014; Dai et al. 2015), and negative phase of NAO (Luo 2005; Luo et al. 2010; Guan et al. 2015) jointly induced the weaker MTF and stronger ZTF, and then made the frequent blocking events in recent warming hiatus (Table 3).

Fig. 15
figure 15

a The 11-year running mean sector blocking frequency for the boreal cold season over EUA (blue line) and NA (red line). The shading represents the one standard deviation range of results from the ERA-interim, NCEP I and II reanalysis products. b The local blocking frequency for the boreal cold season average from 1999 to 2003 (black line) and from 2009 to 2013 (red line). The pink and grey shading represent the one standard deviation range of results from the ERA-interim, NCEP I and II reanalysis products. The blue area represents the positive difference between the average of the results for 2009–2013 and 1999–2003 over the EUA and NA regions. Note that the sector and local blocking was calculated following the procedure outlined in D’Andrea et al. (1998)

Full size image

8 Discussion and conclusion

The recent warming hiatus was identified several years ago (Easterling and Wehner 2009), however, it attracted worldwide attention since the famous scientific community-Intergovernmental Panel on Climate Change (IPCC) highlighted it in recently IPCC report (IPCC 2013). Although the IPCC has discussed some possible causes for the hiatus, they were not enough to win the confidence of people to the fact of global warming and climate science. Many papers about the hiatus have been published since then. The hiatus has stirred debate within the community of climate science, on topics such as the influence of internal climate variability on global temperature change (Kravtsov et al. 2014; Meehl et al. 2014; Dai et al. 2015; Steinman et al. 2015; Guan et al. 2015), and projection of climate in interdecadal time scale (Meehl and Teng 2014; England et al. 2015; Roberts et al. 2015). However, some important aspects of the hiatus, such as the global dynamical figure of it, are still unclear.

Although cooling over land contributed 48 % to the global warming hiatus, the manner in which ICV modes influenced this cooling remains unclear. Our research concludes that these oceanic ICV modes and Arctic influenced the terrestrial SAT by eliciting both weaker MTF and stronger ZTF, which then induced a weaker polar vortex and westerly winds, slower but amplified planetary waves, and specific variations in meridional wind speeds. Thus, extremely cold air was able to invade the mid-latitude regions easily. In addition, there will be more frequent blocking events under the weaker MTF and stronger ZTF situation. We use the term DMO to denote the modulation effect of these oceanic ICV modes and Arctic on the SAT change at decadal to multidecadal timescales. When the DMO trended upward, the warming was accelerated by an associated stronger asymmetric MTF and weaker ZTF, and vice versa. In addition, the asymmetrical thermal forcing and terrestrial SAT change may elicit positive feedbacks, which may further induce amplified regional cooling or warming (Wallace et al. 1995; He et al. 2014). The enhanced warming in the Arctic (Cohen et al. 2014; Mori et al. 2014; Wyatt and Curry 2014) and the downtrend PDO (England et al. 2014; Meehl et al. 2014; Trenberth et al. 2014; Dai et al. 2015) and NAO (Guan et al. 2015) were the main classic ICV modes responsible for the current warming hiatus.

The robust twenty-first century warming projection presented when considered only the ensemble members from the Coupled Model Intercomparison Project Phase 5, which captured the recent hiatus (England et al. 2015). Regardless of the model’s projections, our study also suggests that warming will accelerate again when the DMO enters its upward phase. However, when this change is likely to occur, and how much time we have left to prepare for the accelerated warming expected in the near future are questions that remain unanswered. Furthermore, although the contributions of the temperature change in the warm seasons contributed little to the recent warming hiatus, the warm seasons made a large contribution during the pervious warming hiatus (Fig. S7). Compared with the recent warming hiatus, the spatial patterns of the pervious warming hiatus and its atmospheric dynamics also showed some differences in the cold season (Fig. S4, S6 a, c, e; Table 3). Therefore, more studies about the previous warming hiatus should be undertaken in the future.