J. Meteor. Res.  2018, Vol. 32 Issue (3): 337-350 PDF
http://dx.doi.org/10.1007/s13351-018-7136-4
The Chinese Meteorological Society
0

#### Article Information

QI, Yajie, Zhongwei YAN, Cheng QIAN, et al., 2018.
Near-Term Projections of Global and Regional Land Mean Temperature Changes Considering Both the Secular Trend and Multidecadal Variability. 2018.
J. Meteor. Res., 32(3): 337-350
http://dx.doi.org/10.1007/s13351-018-7136-4

### Article History

in final form January 24, 2018
Near-Term Projections of Global and Regional Land Mean Temperature Changes Considering Both the Secular Trend and Multidecadal Variability
Yajie QI1,2, Zhongwei YAN1,3, Cheng QIAN1,3, Ying SUN4,5
1. CAS Key Laboratory of Regional Climate–Environment for Temperate East Asia, Institute of Atmospheric Physics, Chinese Academy of Sciences (CAS), Beijing 100029;
2. Institute of Urban Meteorology, China Meteorological Administration, Beijing 100089;
3. University of Chinese Academy of Sciences, Beijing 100049;
4. Laboratory for Climate Studies, National Climate Center, China Meteorological Administration, Beijing 100081;
5. Joint Center for Global Change Studies, Beijing Normal University, Beijing 100875
ABSTRACT: Near-term climate projections are needed by policymakers; however, these projections are difficult because internally generated climate variations need to be considered. In this study, temperature change scenarios in the near-term period 2017–35 are projected at global and regional scales based on a refined multi-model ensemble approach that considers both the secular trend (ST) and multidecadal variability (MDV) in the Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations. The ST and MDV components are adaptively extracted from each model simulation by using the ensemble empirical mode decomposition (EEMD) filter, reconstructed via the Bayesian model averaging (BMA) method for the historical period 1901–2005, and validated for 2006–16. In the simulations of the " medium” representative concentration pathways scenario during 2017–35, the MDV-modulated temperature change projected via the refined approach displays an increase of 0.44°C (90% uncertainty range from 0.30 to 0.58°C) for global land, 0.48°C (90% uncertainty range from 0.29 to 0.67°C) for the Northern Hemispheric land (NL), and 0.29°C (90% uncertainty range from 0.23 to 0.35°C) for the Southern Hemispheric land (SL). These increases are smaller than those projected by the conventional arithmetic mean approach. The MDV enhances the ST in 13 of 21 regions across the world. The largest MDV-modulated warming effect (46%) exists in central America. In contrast, the MDV counteracts the ST in NL, SL, and eight other regions, with the largest cooling effect (220%) in Alaska.
Key words: near-term projection     multidecadal variability     multi-model ensemble method     ensemble empirical mode decomposition (EEMD)     Bayesian model averaging (BMA)
1 Introduction

The world is facing increasingly apparent impacts of climate change on ecosystems, human lives, and economic development, so policymakers worldwide need timely and quantitative information on projections of regional climatic changes in the next 10–30 years, as this near-term timescale is important for national medium- and long-term development planning. At this timescale, the magnitudes of internally generated climate variations are comparable to or stronger than the response of the climate system to external forcing (Meehl et al., 2009; Hawkins et al., 2016). Thus, internally generated climate variations should also be considered in the near-term climate projections.

Previous climate change projections usually focused on long-term secular trends (STs) and were based on multiple climate model simulation results that used the approach of arithmetic multi-model ensemble mean (MME), which is the most common technique to cope with multi-model results and was also adopted in the Fifth Assessment Report of the International Panel on Climate Change (IPCC, 2013). However, the MME approach has been thought to be unsuitable for decadal predictions (Meehl et al., 2009; Fu et al., 2011) because internal decadal variability is filtered out when using this method due to the different phases of variations in different models, leaving only an estimate of the forced response signal as commonly adopted in detection and attribution studies (e.g., Bindoff et al., 2013). Thus, internally generated climate variations are not considered in the conventional MME approach.

Multidecadal variability (MDV), which covers seve-ral decades of time in the literature such as 65–70 yr in Schlesinger and Ramankutty (1994), 40–130 yr in Delworth and Mann (2000), and 65–80 yr in Enfield et al. (2001), is believed to at least partly reflect the internal variability of the climate system (Knight et al., 2005; Semenov et al., 2010; DelSole et al., 2011; Wu et al., 2011; Tung and Zhou, 2013; Zhang et al., 2013). MDV in the instrumental record of the land surface air temperature (TAS) can be traced back to the pre-industrial era (Tung and Zhou, 2013) and is likely an internal variability associated with the Atlantic Multidecadal Oscillation (AMO), which is possibly driven by the variability of oceanic thermohaline circulation (Delworth and Mann, 2000; Knight et al., 2005; Sutton and Hodson, 2005; Tung and Zhou, 2013). The Interdecadal Pacific Oscillation (IPO) or Pacific Decadal Oscillation (PDO) has also been suggested to modulate the global mean and regional temperatures on decadal to multidecadal timescales, particularly in the recent “global warming hiatus” period when the IPO switched from a warm phase to a cold phase corresponding to a La Niña-like cooling pattern over the eastern Pacific (e.g., Dai., 2013; Kosaka and Xie, 2016). In addition, the MDV may also be partly modulated by external forcing, notably from anthropogenic aerosol forcing on the climate in North Atlantic (Booth et al., 2012) or on the global climate (Wilcox et al., 2013), although this is still a topic of debate (Zhang et al., 2013). For example, Booth et al. (2012) used the HadGEM2-ES climate model and claimed that the bulk of the MDV in the Atlantic Ocean was forced by anthropogenic factors, such as aerosol forcing; however, Zhang et al. (2013) argued that there were major discrepancies between the simulations from the model used by Booth et al. (2012) and the observations.

The sea surface temperature (SST) changes in multiple oceans have been suggested to be tightly connected and driven by complicated trans-basin mechanisms (Luo et al., 2012; McGregor et al., 2014; Chikamoto et al., 2015; Li et al., 2016; Yao et al., 2017). For example, Luo et al. (2012) showed that the recent rapid warming in the Indian Ocean played an important role in modulating the climate changes in the Pacific Ocean. Li et al. (2016) suggested that the Atlantic Ocean plays a key role in initiating the tropical-wide teleconnection associated with the La Niña-like cooling pattern. As anthropogenic radiative forcing has affected SST changes, this makes the mechanisms of MDV rather complicated. Nevertheless, MDV exists not only in the global mean temperature but also, more profoundly, in regional temperatures (e.g., Schlesinger and Ramankutty, 1994; Zhang et al., 2007; Wu et al., 2011; Tung and Zhou, 2013; Qian, 2016) and reconstructed data (e.g., Zheng et al., 2015; Ge et al., 2017). For instance, Wu et al. (2011) suggested that the rapid warming in the global mean temperature series during the 1980s and 1990s was a result of the concurrence of the upward swing with MDV and a secular warming trend, and the former explained one third of the observed warming. MDV was also found to considerably modulate Northern Hemispheric (Zhang et al., 2007) and regional warming trends (e.g., Qian and Zhou, 2014; Gao et al., 2015; Qian, 2016; Qi et al., 2017). Therefore, changes in the MDV phase in the coming decades may considerably modulate the projected warming trend and should be considered in the near-term projections. In general, if both the ST and MDV are properly considered, we can obtain plausible near-term climate projections.

Statistical models have been proposed that consider MDV in near-term projections via a direct analysis of observations and then empirically extend the projection over the next several decades (Fu et al., 2011; Wei et al., 2015; Wu and Qian, 2015). In terms of numerical model simulations, an initialized decadal climate prediction was developed to simulate the regional climate for 1–30 yr in the future (Meehl and Teng, 2012; Taylor et al., 2012). It was recognized that the effect of initialization was the better phase of internal variability, and the skill of the annual or seasonal mean temperature predictions for the first few years was mainly contributed from the observation-based initialization (Kirtman et al., 2013; Newman, 2013; Wu et al., 2015; Xin et al., 2018). However, initialization is still in its infancy (Kirtman et al., 2013). Other than the uncertainty of individual model simulations, the performance of climate models varies among different timescales of the climate system and across regions. For instance, one model can simulate the ST relatively well; however, its performance in simulating MDV may not necessarily be good (Fu et al., 2011; Qi et al., 2017). A proper multi-model ensemble mean approach is needed to provide near-term future climate change scenarios and uncertainty estimates for practical uses that consider both MDV and STs based on the currently available modeling results. Qi et al. (2017) proposed an “alternative multi-model ensemble mean” (AMME) method for projecting near-term temperature changes. The method maximizes the strengths and minimizes the weaknesses of indivi-dual climate models by considering the capacity of climate models to simulate specific timescale components. It was inspired by Chandler (2013), who proved the validity of such a method from a mathematical point of view. This novel method was found to be skillful when applied to an example in eastern China. However, the method of determining the weights for the multi-model ensemble mean in that study was simple. Some weighting methods have been proposed in the literature, stemming from the belief that the capacities of various climate models are different; thus, these models should not be trusted equally, and some better models should be gi-ven more weight during the process of model combination (e.g., Giorgi and Mearns, 2002; Tebaldi et al., 2004; Schmittner et al., 2005; Chen et al., 2011).

In this study, the weighting method in Qi et al. (2017) is improved by applying the Bayesian model averaging (BMA) approach, which has been found to be able to produce more reliable predictions than other multi-mo-del ensemble techniques (Raftery et al., 2005; Yang et al., 2012). By using this ensemble method (referred to as AMME-BMA) and the Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations, the present study aims to provide a plausible scenario of near-term projections of the changes in TAS during 2017–35 for global, hemispheric, and 21 Giorgi–Francisco land regions across the world (Giorgi and Francisco, 2000) while considering both the ST and the MDV. Although the simulations from CMIP5 are not designed to tie to the observations since their initial conditions are random, the usage of the AMME-BMA approach can be interpreted as a process to mimic “initialization” through selecting the simulations that match the observations. The details of the data and methods used are given in Section 2. Section 3 presents the feasibility of the new method and related results, followed by discussion in Section 4 and a summary in Section 5.

2 Data and methods 2.1 Basic data and pre-processing for specific regions

The observed temperature data for the period 1901–2016 are from the HadCRUT.4.5.0.0 dataset (Morice et al., 2012). The simulated TAS data are from the 34 CMIP5 models listed in Table 1. Most of the models end in 2005, so for the historical period, we use the model data from 1901 to 2005. More details about the simulations can be found at http://www-pcmdi.llnl.gov/. For the future period, we focus on the TAS projection under the “medium” Representative Concentration Pathways (RCP) scenario, i.e., RCP4.5, which corresponds to a radiative forcing of 4.5 W m–2 by 2100. The number of runs differs for different models; thus, to avoid the chance of a model with more runs receiving more weight, the different runs of a model are usually first averaged to a represent single model simulation while making full use of as many members as possible from each model in the ensemble mean procedure. However, in this work, we cannot use the simple average of different runs to represent a single model simulation since the MDV signal could be weakened by using multi-member ensembles with different phases from historical simulations that started from random initial conditions. Therefore, for this study, we consider one ensemble member (“r1i1p1”) from each model. The TAS series of the historical run up to 2005 is concatenated with the series corresponding to the RCP4.5 experiment for 2006–35, as done in other studies (e.g., van Oldenborgh et al., 2013).

Table 1 The CMIP5 models with historical runs and the RCP4.5 emission scenarios included in this study
 Number Model name Organization Resolution 1 ACCESS1.0 CSRIO/BOM, Australia 192 × 145 2 ACCESS1.3 CSRIO/BOM, Australia 192 × 145 3 BCC-CSM1.1 BCC, Beijing, China 128 × 64 4 BCC-CSM1.1(m) BCC, Beijing, China 360 × 160 5 BNU-ESM BNU, Beijing, China 128 × 64 6 CANESM2 CCMA, Victoria, Canada 128 × 64 7 CCSM4 NCAR, Boulder, USA 288 × 192 8 CESM1-BGC NCAR, Boulder, USA 288 × 192 9 CESM1-CAM5 NCAR, Boulder, USA 288 × 192 10 CMCC-CM INGV/CMCC, Italy 480 × 240 11 CMCC-CMS INGV/CMCC, Italy 192 × 96 12 CNRM-CM5 CNRM/CERFACS, France 256 × 128 13 CSIRO-Mk3.6.0 CSIRO/QCCCE, Australia 192 × 96 14 EC-EARTH Europe 320 × 160 15 FGOALS-g2 IAP, Beijing, China 128 × 60 16 GFDL-CM2.1 NOAA/GFDL, USA 144 × 90 17 GFDL-CM3 NOAA/GFDL, USA 144 × 90 18 GFDL-ESM2G NOAA/GFDL, USA 144 × 90 19 GFDL-ESM2M NOAA/GFDL, USA 144 × 90 20 GISS-E2-H NASA/GISS, USA 144 × 90 21 GISS-E2-H-CC NASA/GISS, USA 144 × 90 22 GISS-E2-R NASA/GISS, USA 144 × 90 23 GISS-E2-R-CC NASA/GISS, USA 144 × 90 24 HadCM3 MOHC, Exeter, UK 96 × 73 25 INM-CM4 INM, Moscow, Russia 180 × 120 26 IPSL-CM5A-LR IPSL, Paris, France 96 × 96 27 IPSL-CM5A-MR IPSL, Paris, France 144 × 143 28 IPSL-CM5B-LR IPSL, Paris, France 96 × 96 29 MIROC5 JAMSTEC, Japan 256 × 128 30 MIROC-ESM JAMSTEC, Japan 128 × 64 31 MIROC-ESM-CHEM JAMSTEC, Japan 128 × 64 32 MPI-ESM-LR MPIM, Hamburg, Germany 192 × 96 33 MPI-ESM-MR MPIM, Hamburg, Germany 192 × 96 34 MRI-CGM3 MRI, Tsukuba, Japan 320 × 160

Regional-scale climate information is needed for policy-makers to create targeted engagement and adaptation efforts. Here, the 21 regions defined by Giorgi and Francisco (2000) are used in the analysis (Fig. 1) to provide an overview of the temperature projections from the current models at a regional scale. These Giorgi– Francisco regions have also been used in many other studies, so the present results facilitate comparisons. Since the TAS changes at global or hemispheric scales are known indicators of global climate change, the cases for global land (GL), Northern Hemispheric land (NL), and Southern Hemispheric land (SL) are also analyzed.

 Figure 1 A diagram showing the 21 regions analyzed in this study. The regions are presented by colored boxes according to continent: Australia (red), South America (green), North America (blue), Africa (magenta), Europe (purple), and Asia (black). Giorgi and Francisco (2000) provided the definitions of the regions with latitude and longi-tude information.

The TAS anomalies relative to the 1961–90 mean climatology are calculated for both simulations and observations to diminish the regional biases of individual models. The anomalies from the simulations have various resolutions and are interpolated bilinearly to a 5° × 5° grid to facilitate comparisons with the grid observations. Bilinear interpolation is also used in the IPCC (2013, page 1314) for surface air temperature. Grids without observations in time and space are masked for the model simulations. The area-weighted averages of the annual TAS anomaly series over the land regions are calculated for subsequent analyses.

2.2 Methods

The steps to obtain the AMME-BMA projection are detailed in the Appendix. In brief, the ensemble empiri-cal mode decomposition (EEMD) filter (Huang and Wu, 2008; Wu and Huang, 2009) was used to obtain the ST and MDV components from the observations as well as from the CMIP5 simulations. EEMD is an adaptive method for the analysis of nonlinear and nonstationary data that can be decomposed into a set of quasi-periodic components possessing different timescales from short to long periods and a residual called the nonlinear trend (Wu and Huang, 2009; Qian et al., 2011). The method is different from the Fourier spectra analysis that results in fixed time periods and frequencies. The number of oscillatory components is determined as the integer part of ( ${\rm{lo}}{{\rm{g}}_2}N - 1$ ), where N is the length of the target data. The ST that is extracted adaptively from the data probably reflects the underlying physics better than the linear trend and distinguishes the trend from the fluctuations superimposed upon it (Wu et al., 2011; Qian, 2016). After obtaining the ST and MDV components, we construct the AMME-BMA simulation for the MDV, ST, and the sum of MDV and ST (ST + MDV) (see details in the Appendix) and estimate the future TAS change and its range of uncertainty.

We construct the AMME-BMA models in different land regions for the historical period of 1901–2005 and apply the models to the period 2006–16 for validation by comparing the corresponding parts of the observations. The temperature changes are then projected for 2017–35. We also compare the results with those obtained by the commonly used MME method to mimic a comparison between our results and those from the most recent IPCC report (IPCC, 2013). The MME method gives each climate model equal weight. When using the MME method, how well or how poorly each model performs is overlooked (Flato et al., 2013). In contrast, BMA converts the model ensembles into calibrated and sharp predictive probability density functions (PDFs) (e.g., Raftery et al., 2005; Greene et al., 2006; Furrer et al., 2007; Fraley et al., 2010; Hu et al., 2016). Via BMA, we can obtain a weighted average of the PDFs centered on the individual bias-corrected models. The absolute bias (BI) is used to measure the performance of the AMME-BMA and MME models. The BI is the average of the absolute values of the differences between the observations and the model results over the target period and is defined as:

 ${\rm{BI}} = \frac{{\sum\nolimits_{{{t}} = 1}^{T} {\left| {f_{{t}}^{\rm obs} - f_{{t}}^{\rm model}} \right|} }}{T}.$ (1)

Here, $f_{{t}}^{\rm obs}$ and $f_{{t}}^{\rm model}$ are the observed and model-simulated values, respectively. A small BI indicates that the model performs relatively well.

To measure the importance of MDV relative to the stepwise low-frequency variation in the near-term period (2017–35), we define the ratio (R) as the warming increment induced by MDV divided by that induced by ST + MDV. A larger R indicates a greater dominance of MDV in the long-term climate evolutions.

3 Results 3.1 Evaluation and projection of the ST and MDV for the globe and two hemispheres

The ST, MDV, and ST + MDV time series of the observed annual TAS anomaly over the GL area are shown in Fig. 2. The ST shows a warming tendency from 1901 to 2005, with a warming increment of 0.99°C. This warming is slightly greater than the increment calculated from the linear trend (0.85°C) of the GL series over the same period. MDV enhances the ST from 1910 to the 1930s and from 1980 to 2005, whereas it counteracts the ST from the 1940s to the 1970s. The (ST + MDV) well captures the long-term temperature variations over the past century (Fig. 2), which corresponds to the two accelerations over the GL in the early and late 20th century and one slowdown during the mid-20th century (England et al., 2014; Kosaka and Xie, 2016; Yao et al., 2017). According to Yao et al. (2017), the SST changes in an individual ocean basin have distinct impacts on the GL warming rates between the acceleration and slowdown periods: for the two acceleration periods, the SST warming across all the oceans acted jointly to produce extensive warming over the GL and explained large portions of the GL warming rates and the global mean temperature; for the slowdown period, the compensation effect among the different oceans played an important role, especially for the strong sea surface cooling that occurred in the tropical Pacific (Kosaka and Xie, 2013, 2016; England et al., 2014; Yao et al., 2017).

 Figure 2 The time series of the global mean land surface air temperature (TAS) anomalies (black line, relative to the 1961–90 base period) and the ST (green line) and MDV (blue line) components, together with their combination (ST + MDV), during the period 1901–2005.

Figure 3 shows the AMME-BMA and MME simulations for the MDV, ST, and ST + MDV over the GL. The observed data generally lie in the 5th–95th percentile range of the CMIP5 simulations, implying that the current climate models have the capacity to roughly capture the observed changes. The AMME-BMA results track the observed variations more closely than the MME-based results during both the historical and validation periods (Fig. 3). The BIs of the MDV, ST, and ST + MDV over the GL from the AMME-BMA results are smaller than those from the MME results during the training period, by 83.76%, 98.25%, and 77.37%, respectively (Fig. 4). For the validation period of 2006–16, the MDV, ST, and ST + MDV over the GL from the AMME-BMA simulation also have much smaller BIs than the corresponding results from the MME simulation, by 90.24%, 97.42%, and 43.55%, respectively. Thus, we have more confidence in the AMME-BMA projection of the near-term TAS changes.

 Figure 3 The (a) MDV, (b) ST, and (c) ST + MDV in the global mean land surface air temperature (TAS) anomaly series from 1901 to 2005, compared between the observations (solid black line) and the multi-model ensemble results based on the AMME-BMA approach (solid red line) and the MME approach (solid blue line). The dotted lines are the corresponding components during the validation period 2006–16 and the future period 2017–35. The shaded part denotes the 5th–95th percentile range of the model simulations of MDV, ST, and ST + MDV.

For 2017–35, the MDV over the GL from the AMME-BMA simulation switches to a cooling tendency after reaching its peak around the year 2025 (Fig. 3a), which will thereby partly offset the ST during 2025–35, although the value of the MDV is very small. The AMME-BMA projection of the MDV-modulated GL temperature increase from 2017 to 2035 is 0.44°C; the 90% uncertainty range of which is 0.30°C–0.58°C. In contrast, the MME-based projection displays an increase of 0.57°C (Fig. 3c), which is 0.13°C higher than that of the AMME-BMA projection.

 Figure 4 (a) Box plots of the BI for MDV (1901–2005) in the GL, NL, and SL regions. Each box represents the interquartile range (IQR) and contains 50% of the results of 34 CMIP5 models for a region: the upper edge of the box represents the 75th percentile [upper quartile (UQ)], while the lower edge is the 25th percentile [lower quartile (LQ)]. The horizontal black line within the box is the median. The “+” signs represent outliers [either > (UQ + 1.5 × IQR) or < (LQ – 1.5 × IQR)]. The vertical dashed lines indicate the range of the non-outliers. The blue squares represent the BIs of the MME-based results, and the red circles represent those of the AMME-BMA results for the historical period. (b) and (c) are the same as (a), but for ST and ST + MDV, respectively.
 Figure 5 (a, c, e) are the same as in Fig. 3, but for the Northern Hemispheric land; (b, d, f) are also the same as in Fig. 3, but for the Southern Hemispheric land.

For the NL and SL, the AMME-BMA simulations also result in much smaller BIs than the results of the individual models and the MME in terms of MDV, ST, and ST + MDV (Figs. 4a, b, c) during the historical period. According to the results of the AMME-BMA projections, MDV considerably counteracts the ST in the NL and SL from 2017 to 2035 because it switches to a cooling tendency after the peak around the year 2021/2022 (Fig. 5a). The percentage of the MDV cooling effect is 8.5% for the NL and 5.9% for the SL during 2017–35. The MDV-modulated temperature change in the NL displays a warming increment of approximately 0.48°C (90% uncertainty range from 0.29 to 0.67°C), and that of the SL is approximately 0.29°C (90% uncertainty range from 0.23 to 0.35°C) for the period 2017–35. For comparison, the MME approach produces larger warming increments in both the NL (0.66°C) and the SL (0.47°C) (Fig. 5).

3.2 Evaluation and projection of the ST and MDV for 21 regions

In the observations, strong MDV signals during 1901–2005 exist mostly in the Northern Hemisphere, specifically in Alaska (ALA), Greenland (GRL), Central North America (CNA), Eastern North America (ENA), Northern Europe (NEU), and North Asia (NAS) (Fig. 6a). The observed STs differ among the 21 regions (Fig. 6b): the ST-induced warming increment from 1901 to 2005 reaches up to 2.56°C for Alaska and only 0.45°C for Central America (CAM).

 Figure 6 (a) The MDV in the Giorgi–Francisco regions for 1901–2005 and 2006–35 under RCP4.5. The red curves indicate the AMME-BMA results and the blue curves indicate the MME results. The thin solid black lines show the observed MDV based on the HadCRU4 data. The gray shading shows the 5th–95th percentile range of the multi-model simulations. (b) and (c) are the same as (a), but for ST and ST + MDV, respectively. The axes are colored according to the continents, which are defined in the caption of Fig. 1.

In the model simulations, the observed TAS varia-tions fall within the 5th–95th percentile range of the CMIP5 simulations for most regions (Fig. 6). However, the MDV, ST, and ST + MDV results in the 21 regions from the CMIP5 models have generally larger BIs than those for the globe and the two hemispheres, as previous studies have shown that the biases of the model simulations generally increase with decreasing spatial scale (Masson and Knutti, 2011; Räisänen and Ylhäisi, 2011; Zhang and Yan, 2014). The simulated TAS changes (MDV, ST, and ST + MDV) based on the AMME-BMA method are generally more consistent with the observations than those based on the MME method for both the historical and validation periods, mainly due to the better phases in most regions (Fig. 6). Compared to the results of the MME method, the BIs from the AMME-BMA results decline by more than 27% for MDV, 69% for ST, and 15% for ST + MDV during 1901–2005 among the 21 land regions. For the validation period of 2006–16, the AMME-BMA results also have smaller BIs than those of the MME for a majority of the regions. According to the assessments so far, the AMME-BMA simulations capture the major components of the observations in the 21 regions. Thus, we have more confidence in applying the AMME-BMA approach to the near-term projections of regional TAS changes.

 Figure 7 Box plots of warming increments (ΔT) due to (a) ST and (b) ST + MDV from 2017 to 2035. The blue lines represent the MME-based results and red lines represent the AMME-BMA results.

For the near-term projections, the AMME-BMA results of MDV, ST, and ST + MDV in terms of the temperature anomalies relative to the average of these results for the 30-yr reference period of 1961–90 are also shown in Fig. 6. If only the ST is considered, the regional warming increments from 2017 to 2035 range from 0.14°C [in Central America (CAM)] to 0.86°C in [Central Asia (CAS)]. For comparison, the ST-induced warming increments from the MME projections are larger for most regions, especially Central America (CAM), Central North America (CNA), Eastern North America (ENA), and Alaska (ALA), by up to 0.24°C (Figs. 6b, 7). If both the MDV and ST are considered, the regional warming increments projected by the AMME-BMA from 2017 to 2035 range from 0.18°C (90% uncertainty range from –0.35 to 0.72°C in Alaska) to 0.90°C (90% uncertainty range from 0.66 to 1.15°C in North Asia). Central North America (CNA), Greenland (GRL), Northern Europe (NEU), Sahara (SAH), Central Asia (CAS), Tibet (TIB), and North Asia (NAS) exhibit larger warming increments (above 0.5°C) than elsewhere (Figs. 6c, 7).

 Figure 8 The ratio (R) of the warming increment induced by the MDV to that induced by the sum of MDV and ST for the GL, NL, SL, and the 21 Giorgi–Francisco regions during the near-term period (2017–35).

The contribution of MDV is prominent for the near-term period in Central America (CAM), Central North America (CNA), Eastern North America (ENA), and Alaska (ALA), and the R values are larger than 30% during 2017–35 (Fig. 8). The four regions are all adjacent to the North Atlantic Ocean, where the AMO signal originates. The 50–88-yr oscillation contributed to more than half of the variance of the detrended temperature record in that region and was suggested to be responsible for the 65–70-yr temperature oscillations in many Northern Hemisphere continents (Schlesinger and Ramankutty, 1994). Then, how will this important MDV modulate the secular warming trend in the near-term period in different regions?

According to the AMME-BMA results, MDV enhan-ces the ST in 13 of the land regions for 2017–35, especially in Central North America (CNA), Central America (CAM), and Eastern North America (ENA), as noted above (Fig. 8). The warming enhancement contributed by MDV occurs in Central North America (CNA) (0.18°C or 35%), Central America (CAM) (0.12°C or 46%), and Eastern North America (ENA) (0.14°C or 42%), owing to its upward swing since the 1980s (Figs. 6, 7, 8). The total warming increment induced by the ST + MDV during 2017–35 is 0.52°C (90% uncertainty range from 0.42 to 0.62°C) for Central North America (CNA), 0.26°C (90% uncertainty range from 0.098 to 0.43°C) for Central America (CAM), and 0.33°C (90% uncertainty range from 0.083 to 0.58°C) for Eastern North America (ENA), while the TAS increments induced by the ST are only 0.34, 0.14, and 0.19°C, respectively (Fig. 7).

In contrast, MDV counteracts the ST in 8 regions, including Alaska (ALA), Australia (AUS), Western North America (WNA), the Mediterranean Basin (MED), Southern Africa (SAF), East Asia (EAS), Central Asia (CAS), and Tibet (TIB), for the near-term period, owing to a downward swing of the MDV (Figs. 6, 7, 8). The largest cooling effect (–0.40°C or 220%) of the MDV occurs in Alaska (ALA), with a turning point from a warming tendency to a cooling tendency in the MDV occurring in approximately 2017. The warming increment from ST + MDV in this region is only 0.18°C (90% uncertainty range from –0.35 to 0.72°C) during 2017–35, which is the weakest among all the regions, and that of ST is 0.58°C (Figs. 6, 7, 8). In addition, the second and third largest percentages of the MDV cooling effect exist in East Asia (19%) and Tibet (18%) during 2017–35, with the turning point from a warming tendency to a cooling tendency of the MDV occurring in approximately 2023/24 and 2018/19, respectively, in these two regions (Figs. 6, 7, 8).

4 Discussion

Near-term predictions are constructed in this paper by using the AMME-BMA approach to select models from a probability perspective. An earlier paper (Qi et al., 2017) evaluated the models in terms of MDV by its amplitude and period, not by phase, since historical runs were started from random initial conditions. If a model has an amplitude and period identical to those of the observations and only a phase mismatch, the model would still be suitable for prediction as long as a run could be found that has a phase similar to the observations. However, for practical users of the near-term climate projections who use only the model outputs, it is necessary to choose those runs that match the observed phase. In addition, only a limited number of good performance models were used in the AMME-BMA method, so the cancellation of the MDVs from different models is not as poor as that from the MME approach, which used all models, as indicated in Fig. 3.

Since low-frequency temperature variations over Alaska, especially in the winter season when Alaska is affected very little by solar radiative forcing, are closely linked to the phase changes of the PDO (e.g., Hartmann and Wendler, 2005; Bieniek et al., 2014), the peculiar weak warming in the near-term period in Alaska reported in this study might be a result of phase transition in the IPO from its warm phase to its cold phase around 2000 (Dai, 2013; Kosaka and Xie, 2013; England et al., 2014). This will deserve more regional studies.

In most of the regions analyzed in this study, the multi-model simulations roughly captured the MDV (Figs. 3, 5, 6). However, as noted by previous studies, many models tend to underestimate the inter-basin influence of the Indian and Atlantic Oceans on the Pacific Ocean (Luo et al., 2018), and some models have simulated weaker trends for the acceleration and slowdown periods than the observations (Yao et al., 2017). Although the AMME-BMA approach was used to select models and ensemble with limited models, the projections will have greater confidence if the models are further improved.

5 Summary

In this study, an improved version of the multi-model ensemble analysis of CMIP5 simulations was introduced for near-term climate projections. Compared to the MME method, which mainly results in an estimation of the response of the climate signal to external forcing, this new ensemble approach emphasizes not only the ST but also the MDV. The ST and MDV are extracted from the observations and individual model results over the histori-cal period (1901–2005) via EEMD, and the signals are reconstructed via BMA with observation-based validation. The AMME-BMA method is applied to produce near-term temperature change projections for the global lands, Northern and Southern hemispheric lands, and 21 land regions under the RCP4.5 scenario for the period 2017–35. The main results are summarized as follows.

(1) The AMME-BMA method well captures the ST and MDV components of the observations for both the training and validation periods and is able to reproduce the observed MDV feature for different regions.

(2) The results of the AMME-BMA near-term projections suggest that the MDV in the GL, NL, and SL all switch to a cooling tendency after a peak during 2017–35, although the value of the MDV is very small. The warming increments from the ST + MDV are smaller than those from the conventional MME method. For the 21 land regions, the lowest regional warming increment from ST + MDV occurs in Alaska, while the largest warming increment occurs in North Asia from 2017 to 2035. MDV will enhance the near-future secular warming in 13 land regions due to its upward swing, especia-lly in Central North America, Central America, and Eastern North America, where the impacts of MDV are prominent. Among these regions, Central America will have the greatest MDV-modulated warming effect (46%). In contrast, MDV will counteract the ST in the other 8 land regions, including Alaska, Australia, the Mediterranean Basin, Western North America, Southern Africa, East Asia, Central Asia, and Tibet, due to a downward swing of MDV. Among these regions, Alaska will exhibit the largest cooling effect from MDV due to its turning point in approximately 2017 from an upward to a downward swing, which will markedly slow down the warming trend.

Acknowledgments. The authors thank the reviewers and the Editor for their comments and suggestions to help improve the manuscript. We also thank the Met Office Hadley Center and Climate Research Unit for pro-viding the observed HadCRUT4 data used in this work and the Program for Climate Model Diagnosis and Intercomparison (PCMDI) for collecting and archiving the model data.

Appendix: Steps of the AMME-BMA Projection

A.1 Extracting the ST and MDV using the EEMD filter

White noise at an amplitude of 0.3 times the standard deviation of the target series was added to the target series with an ensemble size of 1000 to execute the EEMD tool. Five (six) oscillatory components and a nonlinear ST can be obtained from an observed series of annual TAS anomaly for the period 1901–2016 (a simulated series for the period 1901–2035). The observed MDV is obtained from the fifth oscillatory component, with two or more extrema during the historical period 1901–2005 and a mean period of approximately 40–100 yr over the different regions analyzed in this study. For the longer time series of the simulations, the combination of the fifth and sixth oscillation components is considered as the MDV.

A.2 Constructing the AMME-BMA projection model

We construct the AMME-BMA model for MDV, ST, and ST + MDV during the training period 1901–2005, validate the model in 2006–16, and make the ensemble projection of TAS changes under the RCP4.5 scenario during 2017–35 by applying the weights from the training period. The BMA weights reflect the contributions of the individual models to the overall predictive skill over the training period. The forecast PDF by BMA is:

 ${\text{p}}\left( {{{y|}}{{{f}}_{\text{1}}},...{f_k}} \right) = \sum\nolimits_{k = 1}^K {{w_k}{g_k}\left( {y|{f_k}} \right)} ,$ (A1)

where ${w_k}$ is the weight of the kth model, and ${g_k}\left( {y|{f_k}} \right)$ is the predictive PDF based on ensemble member ${f_k}$ . The conditional PDF ( ${g_k}\left( {y|{f_k}} \right)$ ) can be approximated by a normal distribution centered at a linear function of the predictor, ${a_k} + {b_k}{f_k}$ . Therefore, the BMA predictive mean is the conditional expectation of y, given the forecast, and can be computed as:

 ${\rm{E}}\left[ {{{y|}}{{{f}}_{\rm{1}}},...{f_k}} \right] = \sum\nolimits_{k = 1}^K {{w_k}\left( {{{\rm a}_k} + {{\rm b}_k}{f_k}} \right)} ,$ (A2)

where ${{{\rm a}}_k}$ and ${{{\rm b}}_k}$ can be calculated by the regression between ${{{f}}_k}$ and y in the training period, and the BMA weight ${{{w}}_k}$ is based on the relative performance of model k in correctly replicating the observed historical changes (MDV or ST) in the surface temperature in the training period 1901–2005. The ${{{w}}_k}$ values are probabilities; thus, they are nonnegative and are assumed to sum to 1. Following the recommendation of Raftery et al. (2005), we used the expectation–maximization (EM) algorithm to obtain the weight ${{{w}}_k}$ for each model. After obtaining ${{{\rm a}}_k}$ , ${{{\rm b}}_k}$ , and ${{{w}}_k}$ , the deterministic BMA model can be determined. More details of the BMA process can be found in Raftery et al. (2005).

A.3 Calculating the warming increment and its uncer- tainty range

To represent the TAS change from 2017 to 2035, we calculate the warming increment induced by the ST, MDV, and the combination. The ST-induced warming increment is defined as:

 ${\rm {Increment}}{_{\rm {ST}}} = {\rm {ST}}\left( {\rm {last}} \right) - {\rm {ST}}\left( {\rm {first}} \right),$ (A3)

where IncrementST indicates the accumulation of long-term warming (Ji et al., 2014), ST(last) is the ST value for the year 2035, and ST(first) is the value for the year 2017 in the ST series. The MDV-induced and MDV-modulated warming increments (IncrementMDV and IncrementST+MDV, respectively) are defined similar to Eq. (4), but for the MDV and ST + MDV series, respectively.

The uncertainty range of the MDV-modulated warming increment via the AMME-BMA model is estimated as follows. First, the 90% prediction interval of ST + MDV (ST, MDV) through the BMA predictive PDF (Raftery et al., 2005; Yang et al., 2012) is obtained. Second, the 95th and 5th percentiles from the first step are used to represent the upper and lower uncertainty limits of the ST + MDV (ST, MDV) in each region and at each time:

 ${\rm{{Upper}_{ST + MDV}}} = {\rm {\left( {ST + MDV} \right)}}^{95{\rm {th}}},$ (A4)
 ${\rm{Lowe}}{{\rm{r}}_{\rm {ST + MDV}}} = {\left( {{\rm {ST}} + {\rm MDV}} \right)^{5{\rm {th}}}}.$ (A5)

Third, the upper and lower uncertainty limits of the MDV-modulated warming increment in the near-term projection are calculated as:

 \begin{align} & {\rm{Incremen}}{{\rm{t}}_{{\rm{ST + MDV\underline {\,\,\,} upper}}}} \hfill \\ & \,\,\, {\rm{ = Uppe}}{{\rm{r}}_{{\rm{ST + MDV}}}}\left( {{\rm{last}}} \right) - {\rm{Lowe}}{{\rm{r}}_{{\rm{ST + MDV}}}}\left( {{\rm{first}}} \right), \end{align} (A6)
 \begin{align}& {\rm{Incremen}}{{\rm{t}}_{{\rm{ST + MDV\underline {\,\,\,} lower}}}} \\&\,\, {\rm{ = Lowe}}{{\rm{r}}_{{\rm{ST + MDV}}}}\left( {{\rm{last}}} \right) - {\rm{ Uppe}}{{\rm{r}}_{{\rm{ST +MDV}}}}\left( {{\rm{first}}} \right). \end{align} (A7)
References