Oceanic forcing of the global warming slowdown in multi‐model simulations

Concurrent with the slowdown of global warming during 2002–2013, the wintertime land surface air temperatures over Eurasia, North America, Africa, Australia, South America, and Greenland experienced notable cooling trends. The oceanic effects on the continental cooling trends are here investigated using two sets of uncoupled experiments with six different climate models. Daily and annually varying sea ice is prescribed for both sets of experiments, while daily and annually varying SST is used in the first set (EXP1) and daily and annually repeating climatological mean SST in the second set (EXP2). All six models capture the slowdown of global‐mean land surface air temperature during 2002–2013 winters in EXP1 only. The slowdown concurs with a negative phase of the Pacific Decadal Oscillation (PDO), indicating that PDO plays an important role in modulating the global warming signal. Not all ensemble members capture the cooling trends over the continents, suggesting additional contribution from internal atmospheric variability.


| INTRODUCTION
The slowdown in the global-mean surface air temperature (GMST) trend during the beginning of the 21st century, concurring in time with a strong increase in atmospheric greenhouse gas concentrations, led to an intense public and scientific interest in the cause for the so-called global warming hiatus (Easterling & Wehner, 2009;Kosaka & Xie, 2016). Although under debate, it has been suggested that the central and eastern tropical Pacific surface cooling played an important part in the hiatus period (1998Kosaka & Xie, 2013;Trenberth, Fasullo, Branstator, & Phillips, 2014). In boreal winter, the warming hiatus manifested as the land surface cooling and cold extremes over Northern Hemisphere mid-latitudes (Cohen, Furtado, Barlow, Alexeev, & Cherry, 2012;Johnson, Xie, Yu, & Li, 2018).
The prominent continental cooling trends in boreal winter have aroused great research attention. The relative importance of tropical Pacific SST and radiative forcing (Kosaka & Xie, 2013), the Arctic sea ice and internal atmospheric variability (Li, Stevens, et al., 2015a), and the tropical Pacific SST and internal atmospheric variability (Deser et al., 2017), on temperature change over Eurasia, have all been examined.
In this study, the effects of SST and sea ice on the regional cooling trends on all continents (black frames in Figure 1c,d) are investigated using two sets of multimodel simulations. Six uncoupled climate models with observation-based SST and sea ice as the surface boundary are utilized. The reason why we choose uncoupled models is that the coupled models' simulations usually overestimate warming during the beginning of the 21st century (Kosaka & Xie, 2013;Medhaug, Stolpe, Fischer, & Knutti, 2017). The novelty of this study is the robust evidence from multi-model simulations which supports the oceanic effects on the global warming slowdown. The rest of the paper is organized as follows. Data and methods are described in section 2. The observed slowdown of continental warming is illustrated in section 3. Discussion on the simulated land surface air temperature trend in the multi-model simulations are presented in section 4. Finally, a summary is given in section 5.
Two coordinated experiments  are performed to investigate the effects of global SST and sea ice on the climate change. Five independent uncoupled (atmosphere-only) climate models, composed of 20-member ensembles each, are adopted: CAM4  (Reynolds et al., 2007). The first set of experiment (EXP1) is forced by both daily and annually varying SST and sea ice, while the second set (EXP2) is forced by daily and annually varying sea ice but daily and annually repeating climatological mean SST (Screen, Simmonds, Deser, & Tomas, 2013). For evaluating the robustness of the results , the sixth model (AFES; 1.5 × 1.5 with 56 vertical levels up to 0.09 hPa; Ohfuchi et al., 2004) has the same experiment set up as the other five models, but was prescribed with monthly mean SST and sea ice from (Hurrell, Hack, Shea, Caron, & Rosinski, 2008). Thirty ensembles are performed for AFES of both experiments. No systematic deviations are found between results derived from AFES and the other five models .
The extended winter 1982 refers to November and December in 1982 and January, February, and March in 1983. The starting year of hiatus in multi-model simulations (2002) is one year later than the observed (2001) (Figure 1a,b vs. Figure 2a -f,h-m,o-t), the choice of the hiatus period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) is therefore contrived to achieve consistency between model simulations and observations. Only land surface air temperatures are considered in this study since we focus on temperature trends over land areas and that the modelled temperatures over the oceans in the uncoupled models are nearly identical to observations. The global-mean land surface air temperature (GMLST) is calculated by area-weighted averaging. It should be noted that there is a clear pattern similar to the negative phase of the Pacific Decadal Oscillation (PDO|−; Wang & Miao, 2018) but no apparent signal of the positive phase of the Atlantic Multidecadal Oscillation (Li, Orsolini, Wang, Gao, & He, 2018a) in the surface air temperature (SAT; Figure 1d). In order to study the regional characteristics of the continental cooling trends during 2002-2013, we choose the six sub-regions (Eurasia: 50 -135 E, 20 -65 N; North America: 90 -165 W, 35 -70 N; Africa: 10 -40 E, 25 S-5 N; Australia: 120 -155 E, 30 -10 S; South America: 40 -80 W, 20 S-10 N; and Greenland: 20 -60 W, 60 -75 N) as indicated schematically by the black frames in Figure 1c,d.

| OBSERVED COOLING TRENDS IN MOST OF THE GLOBAL CONTINENT
Specifically, the observed regional-mean land surface air temperatures (RMLSTs) over the six sub-regions have all shown a decrease after the early 2000s ( Figure 3). Note that the RMLST is defined as the ratio between the regional area-weighted averaging and the summation of

| RESULTS FROM MULTI-MODEL SIMULATIONS
We first address whether the six models can reproduce the realistic GMLST variability. During 1982-2013, the correlations between the observed GMLST and the simulated GMLST in the ensemble-mean of CAM4, WACCM, IFS, IAP4, LMDZOR, and AFES simulations in the EXP1 range from 0.77 to 0.84 (significant at 99% confidence level; Table 3(a)). In the EXP2, the linear correlations with observations are lower than those in the EXP1 but still significant at 99% confidence level (Table 3(a)). The simulated GMLST in the EXP1 is also closer to the observed GMLST in magnitude than that in the EXP2 (Figure 2a-f vs. h-m). Although the simulations suggest the contribution of both SST and sea ice to the GMLST, it is interesting to note that global SST and sea ice variations jointly account for 59.3-70.6% of the inter-annual variability of GMLST, approximately double of that explained by sea ice alone (Table 3(b)). As shown in Figure 2g,n, EXP1 reproduces the slowdown of warming trend in the GMLST which has not been reproduced by the EXP2. The oceanic effect is further estimated through subtracting the simulations in the EXP2 from EXP1 (referred to as EXP1-minus-EXP2; Screen, Deser, & Simmonds, 2012). As shown in Figure 2o-t, EXP1-minus-EXP2 effectively depicts the time evolution of the observed GMLST, with the highest correlation coefficient of 0.80 (  (Figure 2u and Table 1(a)). The successful simulation of inter-annual variability and recent warming slowdown of the GMLST by most of the models implies the potential influence of global SST and sea ice on the global continental temperature variability, with SST being the dominant factor. Regionally, cooling trends over North America, Africa, Australia, and South America are generally reproduced by the EXP1 and EXP1-minus-EXP2 (Figures 3a-f and 4a-f,m-r). While the simulated F I G U R E 4 (a-f) Simulated global surface temperature trend patterns ( CÁdecade −1 ) in winter for 2002-2013 from ensembles' means of (a) CAM4, (b) WACCM, (c) IFS, (d) IAP4, (e) LMDZOR, and (f) AFES in EXP1, respectively. (g-l) and (m-r) same as (a-f), but in EXP2 and EXP1-minus-EXP2, respectively Eurasian winter temperature trend shows diversity among models. WACCM, IAP4, LMDZOR, and AFES ensembles show a slight cooling of the Eurasian winter climate, but CAM4 and IFS capture a slight warming trend in the high latitudes of Eurasia (Figures 3a and  4m-r). The EXP1-minus-EXP2 displays broadly similar trend patterns with the EXP1, except for the absent warming over the Barents-Kara Seas (Figure 4a-f vs. m-r). It can be seen that the spatial distribution of the forced tropical cooling and extratropical warming over the Pacific Ocean resemble the PDO|−, accompanied with weak cooling trends above the North Atlantic Ocean (Figure 4m-r). Moreover, the PDO|− signal is the most dominant feature in the trend of the daily varying global SST which is prescribed as one of the boundaries in the EXP1 ( Figure 5). The presence of the PDO|− in both the SAT and SST trend fields indicates the atmospheric response to SST variability in numerical simulations and thus suggests the potential influence of the PDO|− on the slowdown of warming trend in the GMLST, consistent with a previous numerical study (Trenberth et al., 2014). Also, the PDO|− pattern is complemented with strong easterly winds in the tropical Pacific and remarkable change of the sea level pressure in the subtropics in both reanalysis ( Figure 6a) and numerical simulations (Figure 6b-g). The observed upper tropospheric global teleconnection wave patterns from the tropics to extratropics (Figure 6h) are proposed as possible mechanisms through which the PDO|− causes regional climate anomalies during the hiatus (Trenberth et al., 2014). Here, the wave trains have also been generally replicated (Figure 6i-n), supporting the hypothesis that the PDO|− could lead to continental cooling.
It is worth noting that the model ensemble mean reduces the amplitude of internal variability relative to prescribed forcing . A recent study by Sung et al. (2019) emphasized the role of internal climate variability on the remarkable North American cooling. In addition to the fundamental effect of PDO (Sung, Kim, Baek, Lim, & Kim, 2016), we further examine the potential role of atmospheric internal dynamics through removing the ensemble mean of each model from individual members in the EXP1. It shows that six members from CAM4, six members from WACCM, three members from IFS, seven members from IAP4, six members from LMDZOR, and ten members from AFES (from which the ensemble mean has been removed) did reproduce the North American cooling trend (identified by the threshold above 0.5 CÁdecade −1 over North America), though the temperature trend over the ocean becomes considerably weak (Figure 7a-f). The accompanied intensification of the anticyclonic ridge near Alaska and cyclone to the south (Figure 8a-f), which show a basin-scale northsouth atmospheric circulation over the North Pacific, resemble the negative phase of the North Pacific Oscillation (NPO|−; Lee, Hong, & Hsu, 2015;Rogers, 1981). The overlying atmospheric circulation trend in the upper troposphere (Figure 8g-l) corresponds to the negative phase of the Pacific North America (PNA|−) pattern (Wallace & Gutzler, 1981). It suggests that the internal atmospheric variability might also be a possible driver for the North American cooling (Sung et al., 2019). Therefore, our results suggest that the internal atmospheric variability still has impacts on climate change though the contribution of external forcing is dominant.
Apparent temperature trends mainly occur in high latitudes of Northern Hemisphere in the EXP2, showing much weaker amplitude over North America relative to EXP1-minus-EXP2 (Figure 4g-l vs. m-r). Strong warming emerges over the Barents-Kara Seas and subarctic Russia, while the observed surface warming in the East Siberian-Chukchi Seas, which is suggested to be responsible for North American cold winters (Kug et al., 2015), is not clear in the EXP2 (Figure 4g-l). The feature of Eurasian temperature trend in the multi-model simulations forced by sea ice variations is also interesting. Generally, cooling in two (Figure 4g,l) and warming in four (Figure 4h-k) sets of simulations over Eurasia make it difficult to attribute the Eurasian cooling to the reduction of sea ice in this study. A recent study by Ogawa et al. (2018), which is based on the same model simulations, has focused on the impacts of Arctic sea ice loss on the Eurasian winter cooling trend. Only a small number of members in both experiments can simulate Eurasian cooling trend with the observed amplitude   fig. 3), leading to the conclusion that atmospheric internal dynamics instead of sea ice changes might play an important role in the observed Eurasian cooling. In this study, some individual members of EXP1, from which the ensemble mean has been removed, have  (Figure 7g-l). Correspondingly, an intensified surface high pressure over Northern Eurasia (Honda, Inoue, & Yamane, 2009) is found in these individual members (Figure 8m-r), consistent with the reanalysis (Figure 5a). High-latitude internal variability (Figure 8m Kosaka, Watanabe, Nakamura, and Kimoto (2019) concluded that~44% of the Eurasian cooling trend in recent decades is due to sea ice loss in the Barents-Kara Seas by using a hybrid analysis of observations and multi-model ensembles from atmospheric general circulation models. Both viewpoints can get supporting evidences from our multi-model simulations (Figures 4 and 7), which explains why it has been a controversial issue whether the Arctic sea ice loss Liu et al., 2012;Mori et al., 2014Mori et al., , 2019 or internal atmospheric variability (Kosaka & Xie, 2013;Li, Stevens, et al., 2015a;McCusker, Fyfe, & Sigmond, 2016;Ogawa et al., 2018;Sun, Perlwitz, & Hoerling, 2016) can influence Eurasian winter climate.
Quantitatively, RMLST trends simulated by the EXP1-minus-EXP2 are in broad agreement with the observations, except for those over Eurasia (Figure 3). In response to the oceanic forcing, the simulated RMLST in North America exhibits a cooling trend of −0.073, −0.089, −0.059, −0.069, −0.039, and −0.080 CÁdecade −1 during 2002-2013, approximate to the observed −0.077 CÁdecade −1 (Figure 3b and Table 1(c)). Over the F I G U R E 8 (a-f) Simulated sea level pressure (shading) and 850-hPa vector winds (vectors) and (g-l) 300-hPa streamfunction trend patterns (10 6 m 2 Ás −1 Ádecade −1 ) in winter for 2002-2013 from the ensemble members (from which the model ensemble mean has been removed) that have simulated the North American cooling trend in (a, g) CAM4, (b, h) WACCM, (c, i) IFS, (d, j) IAP4, (e, k) LMDZOR, and (f, l) AFES in EXP1, respectively. (m-x) Same as (a-l), but for the ensemble members that have simulated the Eurasian cooling trend other four sub-regions, the simulated RMLST trends are close to the observed ones (Figure 3c-f and Table 1(d)-(g)).
Note that the observed Greenland cooling is not well reproduced by two models in the EXP1-minus-EXP2 (Figures 3f and 4p-q), which is likely associated with the atmospheric teleconnection (e.g., the North Atlantic Oscillation; figure not shown) as suggested by previous studies (Hurrell, 1995;Hurrell, Kushnir, & Visbeck, 2001). In general, the high consistency between the RMLST trends in observations and numerical simulations (Figure 3), which targets potentially the response to dominate SST change, strengthens the hypothesis that the PDO|− largely explains the surface temperature cooling trends contributing to the slowdown of warming trend in GMLST.
The comparisons of atmospheric response between simulations forced with observed daily-varying SST and daily-varying sea ice and those forced with climatological SST and daily-varying sea ice have suggested the potential oceanic contribution to the slowdown of global warming. The results are consistent with previous studies (such as Kosaka & Xie, 2013). However, it might be noted that the effects of SST and sea ice have not been totally separated in this study due to the potential interaction between SST and sea ice. Further study can focus on the atmospheric response to the observed daily-varying SST and climatological sea ice, which will further improve our understanding of oceanic impacts on climate. Additionally, the simulations in this study are not suitable to investigate the role of external forcing on the continental cooling, which, however, has been revealed by many studies (Fyfe, von Salzen, Cole, Gillett, & Vernier, 2013;Huber & Knutti, 2014;Marotzke & Forster, 2015;Santer et al., 2014). For example, a modelling study of Fyfe et al. (2013) suggested that increasing stratospheric aerosol since the late 1990s has a global cooling impact of about 0.07 CÁdecade −1 . Such knowledge on the global warming slowdown will improve our understanding of climate change in the future.

| SUMMARY
Observations indicate that the land surface air temperatures over Eurasia, North America, Africa, Australia, South America, and Greenland exhibit remarkable cooling trends, concurrent with the slowdown of global warming during 2002-2013. The inter-annual variability and recent warming slowdown of GMLST are successfully reproduced by the multi-model simulations. Global SST and sea ice anomalies in combination play a major role, and that oceanic contribution dominates. The observed global teleconnection wave trains, which are suggested associated with the PDO|−, have been reproduced by the multimodel simulations. We thus hypothesize that the wintertime cooling trends in North America, Africa, Australia, South America, and Greenland during 2002-2013 might be the response to the shift of the PDO to its negative phase. Additionally, the internal atmospheric variability has also contributed to the cooling trend over Eurasia and North American. It might be noted that there is an apparent La Niña-like decadal cooling trend which may also contribute to the global warming hiatus. However, the multi-model simulations in the present study cannot identify the impact of regional SST. Further investigation is needed in this regard, and numerical simulation by multi-models is an essential way to understand the climate response to the natural forcing.