A novel statistical decomposition of the historical change in global mean surface temperature

The station surface air temperature series have been quality controlled and homogenized (Xu et al 2018, Li et al 2021), and then converted to a 5° × 5° latitude and longitude grid dataset using the climate anomaly method (CAM) (Jones et al 1997), based on the 1961–1990 period. Then this China global land surface air temperature dataset (C-LSAT) was merged with NOAA's ERSSTv5 (Huang et al 2017), to develop a new CMST dataset and new GMST series, with 5%–95% uncertainty ranges evaluated (Yun et al 2019, Li et al 2020). The CMST dataset/series covers the monthly in situ data from January 1880 to December 2019. Comparative studies (Li et al 2020, 2021) showed that CMST has similar trends of global surface temperature and have similar uncertainties compared to other existing datasets. While the significant trends during the so called 'hiatus' period are identified by CMST, HadCRUTem4-Hybrid (infilled by the polar surface air temperature) (Cowtan and Way 2014) and ERA-5 (Simons et al 2017), etc.

Anthropogenic forcing data includes carbon dioxide (CO₂), other greenhouse gases, atmospheric aerosols, snow albedo changes and land use (Hansen et al 2007, Myhre et al 2013, Shindell et al 2013), as well as the natural forcing data including volcanic eruptions and solar irradiation (Hegerl et al 2007, Bindoff et al 2013, Myhre et al 2013, Shindell et al 2013). Data for all these (best evaluation) were downloaded from Pangaea database and KNMI's Climate Explorer website (https://doi.pangaea.de/10.1594/PANGAEA.919662) (CMST), (C-LSAT2.0) and http://climexp.knmi.nl/selectindex.cgi?Id=someone@somewhere (radiative forcing), all downloaded in Jan 2020). The time range of these forcing data are from 1850 to 2017 as shown in figure S1 (supplement information (SI) (available online at stacks.iop.org/ERL/16/054057/mmedia)).

Previous IPCC assessment reports stated that climate change signals can be decomposed into externally forced responses and internal climate (unforced) variability. Anthropogenic activities are becoming much clearer as the main contributors to surface air temperature warming, natural external forcing (such as volcanic eruptions, changes in solar output, represented by corresponding radiative forcing) also contribute to the climate change (Hegerl et al 2007, Bindoff et al 2013), and these contributions have also generally used a linear relationship expressed from at least one General Circulation Model through OPF method (Haywood et al 1997, Ramaswamy and Chen 1997). Therefore, a multiple linear regression (MLR) model is first used in this study, trying to separate the contribution of climate changes due to external forcing and internal variability.

Taking the GMST as the dependent variable, two variables (total anthropogenic forcing, natural forcing) are selected as the independent variables in the MLR model. We established the following MLR models of GMST change for the period 1880–2017:

$$\begin{equation}y = {a_0} + {a_1}{x_1} + {a_2}{x_2} + \varepsilon \end{equation} \tag{ 1 }$$

where, ${x_1}$, ${x_2}$ are the total natural and anthropogenic forcing, respectively, $\varepsilon$ is the residual series of the MLR model. ${a_0},{a_1}{\text{ and }}{a_2}$ are the regression coefficients evaluated by the ordinary least square (OLS) method. It should be noted that the recent OFP uses the total least square (TLS) method, which can be interpreted as the OLS with noise elimination (Bindoff et al 2013). However, in terms of linear regression in this paper, the OLS method will provide the best linear unbiased estimator (BLUE) of the regression coefficient according to the Gauss–Markov Theorem (Theil 1971), while the TLS needs more assumptions.

Normally, if equation (1) is strictly a BLUE, its physical meaning is very clear, the MLR represents the part of the dependent variable (GMST) affected (caused) by the independent variables (the total natural and anthropogenic forcings), and $\varepsilon$ also represents internal variability in annual GMST. Since the residual analysis shows that there is no collinearity between natural forcing and artificial forcing, we only need to investigate if the residual has obvious negative or positive autocorrelation, which leads to the linear model no longer being BLUE, and the statistical significance of the variables is meaningless, that is, the simulation (fitting) of the model is invalid. In this case, equation (1) will need to be combined with ARIMA models (including the autoregressive (AR), moving average (MA), ARMA and ARIMA models, see Li (2020), and the supplemental information (SI)).

Generally, the residual of a successful MLR model should have the following characteristics: (a) the residual obeys the normal distribution with zero mean; (b) the residual is a normal distribution with equal variances; (c) the residual sequence is independently distributed. That is, the residuals cannot be linearly fitted with new independent variables, leaving only periodic and random components. According to the MLR model, the residual represents the random error term of the model, which is the error between the fitted value and the actual value excluding the constant term in the independent variable interpretation space. In the sense of climate change, when the fitted external forced 'signal' is subtracted from the GMST anomaly series, the remaining series (i.e. the residual) can be considered to be mainly composed of the internal variability component ('noise' of climate change). We can know from the section 3 that the residual of MLR ('noise' or internal climate variability) still explains about 8.4% of the total variance of the GMST series. Due to the relatively more important contribution of internal climate variability has been claimed in the shorter term GMST variations in recent studies (Dai et al 2015, Frankcombe et al 2015, Folland et al 2018, Li et al 2020), we further decomposed the periodic signals in the residual series by an ensemble empirical mode decomposition (EEMD) method.

EMD model is a method of decomposing complex signals into a finite number of intrinsic mode functions (IMFs) through empirical mode decomposition:

$$\begin{equation}x\left( t \right) = \mathop \sum \limits_{i = 1}^n im{f_i}\left( t \right) + {r_n}\left( t \right)\end{equation} \tag{ 2 }$$

where, imf_i (t) is the ith IMF obtained by EMD; r_n (t) is the signal residual component after n IMFs are decomposed and screened.

Each IMF component contains local characteristic signals of different time scales of the original signal. Compared with Fourier transform, wavelet decomposition and its basis function are derived from the data itself. Therefore, this method is intuitive, direct, posterior, and adaptive. EEMD is a noise-aided data analysis method based on the shortcomings of EMD method in modal aliasing (Huang and Wu 2008). The principle is that when a signal is added with a uniformly distributed white noise background, signal regions of different scales are automatically mapped to the appropriate scale related to the background white noise. Since the noise is different in each independent test, when the overall average of enough tests is used, the noise will be eliminated, and the only lasting and stable part is the signal itself (Huang et al 1998, Huang and Wu 2008, Wu et al 2011).

PLSR is an approach that is relatively new multivariate statistical data analysis tool. Because it can solve some problems that conventional ordinary multiple regression methods cannot, it has attracted the attention of researchers. PLSR provides a method for regression modeling of multiple dependent variables to multiple independent variables especially when there is a high degree of correlation among the variables tool, the analysis model is more reliable and the overall conclusion is most robust. It can also perform regression analysis, data structure simplification (PCA), and correlation analysis between two groups of variables (canonical correlation analysis) comprehensively. It is suitable for establishing regression models when there is a high degree of multiple correlation among the variables, which improves the deficiencies of traditional multiple regression analysis, making the analysis conclusions more reliable and robust. It is also suitable for regression modeling when the sample size is less than the number of variables (see the detailed introduction in text S2 in SI). In the past four decades, PLSR has been developed rapidly in both theory and practice, and its application field has also rapidly expanded to more social and natural science fields (Wang 1999, Li et al 2008). A PLSR model has also been used here to investigate the individual external forcings responses to the global and regional land precipitation variations like Li (2020).

Taking the GMST as the dependent variable, two other variables (the total anthropogenic forcing and the natural forcing) are selected as the independent variables to construct the MLR model. We established the following combined model of MLR and ARIMA models of GMST variations for the periods of 1880–2017:

$$\begin{align}y & = - 0.438 + 0.048*{x_1} + 0.410*{x_2} \nonumber\\ & \quad + 1.014*AR\left( 4 \right) + \varepsilon \end{align} \tag{ 3 }$$

The analysis of the MLR model shows that neither of the above two variables (natural forcing ${x_1}$; anthropogenic forcing ${x_2}$) can be eliminated (both are important). The adjusted square of the multiple correlation coefficient (total variance explained by the model) of the combined model is 91.6%; the analysis of variance shows that at the 5% level, the dependent variable and the two independent variables have a clear linear relationship, which proves the combined 'MLR and ARIMA' model is statistically meaningful (note that only the AR model is needed here); and both of the coefficients of the two independent variables (the anthropogenic and external natural forcing) have significant meaning. The variance inflation factors are 1.000, which also indicates that there is no co-linear relationship between the two independent variables; the Durbin–Watson (DW) test of the residuals also shows there is no first order autocorrelation in equation (3) (DW = 2.026, also see the table S1 of SI). From the above analysis and statistics, the combined model (3) are successful (BLUE) of GMST with the two independent variables, which also shows that the contribution of these two independent variables to the GMST variation is quite clear.

Figure 1(a) shows that the above regression fits well the GMST long term trend signals due to the response of the anthropogenic and natural forcing in the 138 years from 1880 to 2017. Figure 1(b) shows the residual of the combined model, and figure 1(c) shows the residual is a normal distribution with equal variances and the model is meaningful (BLUE) according to the three criteria mentioned in section 2.3. Of the two kinds of the external factors, the anthropogenic forcing has caused most of the GMST warming during the last 138 years. Before the second half of the 20th century, greenhouse gas emissions increased slowly, so the GMST also increased slowly, and it entered a period of rapid warming since the 1970s. However, we also notice that a certain contribution of the 'noise' (or internal variability), which explains about the rest 8.4% of the total variance contribution also agrees well with the previous studies (Li et al 2020). In addition, the fitted series was obviously affected by external forcing influences from five large volcanic eruptions: 1884 (Krakatau volcano, Indonesia, eruption in 1883), 1903 (Mt Pelee, Martinique, eruption in 1902), 1964 (Agung volcano, Indonesia, eruption in 1963), 1983 (Kilauea volcano, Hawaii USA, eruption in 1983), 1992 (Pinatubo volcano, Philippines, eruption in 1991) (Zhang and Zhang 1985, Hou and Li 2000).

Zoom In Zoom Out Reset image size

Based on the EEMD model, the residual series of the above regression (figure 1(b)) is decomposed into six IMFs (the fluctuations at different timescales) and the remainder (figure 2). As can be seen from figure 2, there are obvious periodic signals in the internal variability of 2.9 years, 6.0 years, 13.8 years, 33.4 years, 69.8 years and 136.1 years. Among these, the first five IMF1–IMF5 (the variance contribution explained are 25.8%, 19.0%, 17.9%, 23.5% and 13.7%, respectively) passed the significant test at the 5% level, IMF6 explains only about 0.21% of the variance, and its periodicity was not significant. In other words, all these former five periodic changes are significantly reflected in the long-term variations of GMST. The further comparison also shows that the deposition of the multi-scale periodic signals is not sensitive to whether the external forcing effect is extracted in advance.

Zoom In Zoom Out Reset image size

Figures 2(a)–(f) also show that the internal climate variability term in the residual can be well decomposed into these several periodic signals. It can be seen that the inter-annual changes and the high and low phase of the GMST anomalies are well explained. IMF1, IMF2 explain the inter-annual variations (IAV) of the GMST, IMF3, IMF4 explain the IDV of the GMST, and IMF5 well explains the MDV warmer or colder fluctuation, and these five IMFs are exactly the eigenmodes with the largest contribution of variance, explaining a total of about 99% of the total variance of the combined model's residual, which in total accounts for about 8.4% of the total variance of GMST variations. Further, numerous previous studies indicated that the ENSO (2–7 year cycle), AMO (65–80 year cycle) and PDO/IPO (inter- decadal variation) are the dominant modulation signals of these three components (Dai et al 2015, Folland et al 2018, Wei et al 2019, Wu et al 2019). This agrees well with the multi-scales variations of internal climate (unforced) variability separated by EEMD in this paper.

Although the components of the natural and anthropogenic forcing responses have been separated through the above combined model of MLR and ARIMA, it is obviously not enough. More in-depth studies often need to clarify the contributions of individual forcing (such as CO₂, aerosol, and even ozone in tropospheric and stratospheric layers), then the above models are no longer applicable. It has always been a difficult problem to separate the impacts of various anthropogenic forcings on climate change, that is, the 'confusion' between multiple independent variables. For example, the effects of urbanization and greenhouse gases on surface temperature variations are very similar (The average temperature increases, but the daily temperature range decreases). In terms of statistical language, these two variables may have a very clear correlation (Chong and Jun 2005) (this is very common in the case of multiple variables, namely multiple correlations or collinearity between factors, see table S2 in SI).

Figure 3 shows that the first six principal components (PCs) explained nearly 90% of the variance GMST, and it will not increase much with the introduction of more PCs. This shows the PLSR fits the GMST change well during the period of 1880–2017 with the first six PCs. Figure 3 presents the standardized PLSR coefficients of the independent variables. From figure 3 (and table S3 in SI), the coefficients for CO₂, black carbon on snow (snow_BC), other GHGs, Aerosol and Natural forcing are significant at the 5% level, the land use's coefficient is marginally significant at 10% level, and the rest are insignificant. The independent variables including CO2, other GHG, land use and natural forcing, contribute to the GMST warming, while the rest (including aerosol, snow_BC, Ozone, water vapor and contrails (the last three variables' coefficients are insignificant)) contribute to cooling of the GMST. The slight difference with Myhre et al (2013) occurs in the warming (cooling) effect caused by the surface albedo change related to the land use and the black carbon on the snow surface. This shows there is still a large uncertainty about the correlation between surface albedo and global climate change currently (Schwaiger and Bird 2010, Sieber et al 2019).

Zoom In Zoom Out Reset image size