Global surface temperature variability and trends and attribution to carbon emissions

We investigate the patterns of monthly time series of global ocean surface temperature and global air temperatures over the land surface, and combination thereof, from 1850 to 2018. By employing an empirical mathematical methodology, we decompose the non-linear global temperature time series and confirm patterns of frequency and amplitude modulated, discreet internal modes of variability within the two basic time series and the third, combined, time series. We find periods of warming and cooling on both the surfaces of the ocean and atmosphere over land with prominent seasonal, annual, inter-annual, multi-year, decadal, multi-decadal, and centennial modes, riding atop overall warming trends. Our calculated overall rates of warming differ significantly from the estimates of the Intergovernmental Program on Climate Change (Bernstein et al ., 2008; Stocker et al ., 2013). We find the oceanic warming rate to be less than two-thirds of surface air over land, making the ocean a regulator or comparative heat capacitor and the dominant player in determining global surface temperatures. The ocean has an enormous heat capacity relative to that of the atmosphere. Next, by employing an econometrics-based statistical formula, we establish a causal relationship between fossil fuel burning and global surface temperatures, which definitively links the overall trends in planetary surface temperature rise, to the overall upward trend in fossil fuel burning. Our study also found that there is a 1-year phase lag of global temperatures to fossil fuel burning both on land and in the sea.


Introduction
The Earth system absorbs the Sun's incoming short-wave radiation and re-radiates, stores or exchanges it at different rates via natural processes. For a planetary condition of thermal equilibrium, the amount of total outgoing, long-wave radiation would equal the total amount of incoming, short-wave radiation. If outgoing radiation does not equal incoming radiation, then the planet's global body temperature will vary both spatially and temporally; so, "climate variability". Many peer-reviewed literature studies have definitively shown that the climate has been warming over the past century and they collectively attribute the warming climate to fossil fuel burning as recent studies have found links between fossil fuel burning and climate warming  [1,2,11,17,26,34,35] , focusing on what corporations were responsible, thus the case for attribution, for having produced increases in greenhouse gases. Thus, the relationships were visually correlated with causality assumed. Numerical model studies, such as Millar et al (2016) [27] and Ekwurzel et al (2017) [11] , have employed global energy-balance coupled climate-carbon-cycle models and have attempted to assess the change in atmospheric CO 2 and CH 4 concentrations, radiative forcing, and Global Mean Surface Temperatures to the emissions traced to 90 major global carbon producers; again seeking attribution for the warming by percentage of greenhouse gas production by industry. Alternatively, our study is based entirely upon the global temperature time series and the fossil fuel burning time series and thus are data based. The public media reports that 97% of polled scientists accept this as a form of virtual attribution, however visual correlation per se, does not accomplish such a direct relationship. Attribution, by definition, links cause to effect. In the study reported on below, we will revisit the temperatures of the global surface of the Earth (upper panel in Figure 1), the atmosphere above land (middle panel in Figure 1), and the surface of the global ocean (lower panel in Figure 1) from 1850 through the first half of 2018. We will employ the well-known and documented surface-temperature-anomaly-data. While much attention has focused on the global surface atmospheric temperature record, as greenhouse gases have built up in the atmosphere, far less attention has been paid to the global ocean surface temperature record in-kind. We will investigate the variability of land and ocean surface temperatures and determine if we can shed new insights into what is revealed regarding the present state of the Earth's climate system. The climate system consists of non-linear and non-periodic processes, so we will utilize an empirical, mathematical, data adaptive technique to decompose the data. Our mathematical decomposition methodology produces internal, intrinsic modes of variability buried within the temperature time series that can be replicated with sine functions and thus can be considered quasi-stationary. "Nonlinear" (NL) is defined as a time series of systems governed by equations more complex than the linear form, = + , where 'a' and 'b' are constants. We will account for the patterns revealed by the internal modes of variability, relate them to naturally occurring physical phenomena and then compute the overall trends which we find does not have a natural causal basis. Consequently, we will employ a statistical hypothesis test for determining whether one time series, such as fossil fuel burning, can be useful in forecasting another, specifically global surface temperatures, and thus to predict climate warming, past, present and future; thus causality.  Huang et al. (1998) [19] employed the Hilbert Transform (Gabor, 1946; Van der Pol, 1946) [14,37] , in the development of the Empirical Mode Decomposition (EMD) methodology. In the application of the HT and the EMD, higher to lower frequency modes are removed from the time series in both a time and amplitude varying decomposition, until what remains is what was originally referred to as the "residue" or "gravest mode". The gravest mode or residue could go up or down or down and up or up or down. Thus, the residue or gravest mode was not in and of itself an oscillation. The three global surface temperature anomalies time series presented in Figure 1, display very a very strong sense of NL temporal variations. It is of note that a modemixing problem existed in the EMD decomposition in which successive IMFs were discovered to occasionally mix or contaminate each other. To address the mode-mixing problem,   [19] added white noise to the various time series, ensembles were created and the mean IMFs were found to stay within the natural dyadic filter windows, preserving their dyadic properties, leading to stable decompositions of frequency and amplitude modulated internal modes of variability in the record length data. However, what of a "trend" of a time series?

Background and Methodology
In discussing the trend of any time series, we first need to consider the definitions and the methodologies of computing a trend, as this is especially important for temperature time series. As such, no conventional simple averaging process can be utilized to reflect what information is buried in the multiple temperature time series. This underscores the importance of clearly defining what constitutes a trend. Granger (1966) [15] presented an insightful definition of a trend as "a trend in mean comprises all frequency components whose wavelength exceeds the length of the observed time series". In the mathematical treatise of definitions of James and James (1976) [24] , the definition of a trend is stated to be "the general drift, tendency or bent of a set of data". Chatfield (1977) [8] defines "trend" as "a long-term change in the mean". Thus, difficulties arise in determining what "long term" means. For example, a 50-year cycle would appear to be a trend in a 20-year data set; but it would be shown as a 50-year cycle when the length of the data was longer than 50 years. In speaking of a "Chatfield trend" we must take into account the number and span of observations available and then make a subjective assessment as to what constitutes "long term". For NL, NS datasets, none of these definitions is mathematically applicable leading to employment of the definition based on the work developed in Wu et al. (2007) [19] . In that paper, the EMD decomposition included the non-oscillation "residue", after the higher to lower frequency modes, the Intrinsic Mode Functions (IMFs), in the global surface temperature anomalies (the GSTA) were removed in the time varying decomposition. The authors then called that residue the trend. Thus, given our discussion above, prior to the revelation of Wu et al. (2007) [19] there was not an applicable definition of what constituted the overall "trend" of any continuous time series of data, especially of a non-linear NL time series. As such, a well mathematically defined methodology of how to compute a trend of a time series had not existed in the literature. This work presented a logical definition of a trend as an intrinsically determined function within the temporal span of the data, in which there can be at most one extremum. Being intrinsic, the method to derive a trend was adaptive; that is, it suited the span of the data. Thus, the definition of a trend in this publication presumed that a time scale exists that is logically and mathematically related to the span of the data, and was in keeping with the definition suggested previously by Granger (1966) [15] . With this definition of a trend, we will employ the EEMD method   [19] to decompose the time series of global ocean surface temperature and air temperature above land surface, from 1850 through July 2018, and identify all of their intrinsic modes and trends.

Global ocean and land based atmospheric temperature anomaly time series
Following the works of Wu et al. (2007) [19] , and Pietrafesa and Gallagher (2017) [12] , an updated analysis of global surface temperatures though 2018 is presented below. In Figure 1 we   [10] . Mode 10, 105-110 years represents the Global Thermohaline Circulation Conveyor Belt (Rahmstorf, 2003;Sabra et al., 2016) [30,32] . Mode 11 is the gravest mode or overall trend of the 167-year time series of total data (Table 1). We note in Figure 2, that for the GSTA, the IMF modes 2, 3 and   In Figure 3, the overall trends of the GOSTA, GLSTA and GSTA are presented relative to each other, from their anomaly starting points January 01, 1850 to their end values December 31, 2017. In Figure 3  Only an adaptive method such as EEMD   [19] can do so, as explained above   [19] . We also note from Figure 3 that from 1850 to 1895 the air temperature trend over land actually decreased slightly, so there was a slight global cooling over land globally that persisted for about 45 years. Then in 1895, the global air temperature over land reached its prior value in 1850 and then continued to rise through 2017. The scenario was essentially flat for the global ocean with warming commencing about 1880. The combined atmosphere and ocean time series follows closely to that of the ocean and moved from slight cooling to warming about 1880. The air over land was relatively cooler than that of the surface of the ocean until about 1915 when the trend lines cross and air temperatures on land rose at a more rapid rate than did those of the ocean. The air over land and ocean surface trend curves have diverged significantly up through 2017.

Fossil Fuel Burning, Carbon Emissions and CO2
In Figure 4 we present the fossil fuel burning time, Carbon Emissions (CE) series from 1751-2014, provided via <http://cdiac.ornl.gov/trends/emis/tre_glob.html>. As shown in Figure 4, from 1751 until 1850, global CO2 production was negligible, having ranged from 3 MMTs in 1751 to 54 MMTs in 1850. In the latter half of the 19 th century, carbon emissions increased considerably, with a slight decrease during World War II, followed by a dramatic upsurge around 1949 up to 2014, reaching a value of 950 MMTs. The CE was essentially flat from 1751 through the mid latter half of the 19 th century when the change of carbon burning began to be spiky. Reasons for that are the spiky advances in the industrialization of different parts of Europe and the emergence of industry in the U.S. One could proceed here with cross correlations between the CE and GLSTA, GOSTA and/or the GSTA. However, although the temperature curves all display strong visual correlation with the fossil fuel burning curve, visual correlation does not create a cause-andeffect relationship or attribution.

Granger causality relating carbon burning to global surface temperature anomalies
We next employ the Granger Causality Test (GCT) (Granger, 1980) to obtain evidence for a causal relationship between the carbon burning time series and yearly temperature anomaly time series. GCT statistically tests for determining whether one time series is useful in forecasting another. Generally regressions reflect mere correlations, but Granger argued that by measuring the ability to predict the future values of a time series using prior values of another time series, causality in economics could be tested. A time series X is said to " Granger- [3] , Attanasio et al. (2013) [4] , and , amongst others. The above papers consider more time series than those considered here so they can for example separate out anthropogenic forcing and they also allow for non-stationarity in the temperature series. Our approach differs in that we consider only anomaly series {y t } and one emission series {x t }, and we find a stationary AR model for the anomalies, which we justify our modeling by use of a time series goodness-of-fit test (see Fisher and Gallagher, 2012) [12] . One advantage of our methodology is that our final fitted model can be used to predict the impact of current (future) carbon emissions on temperature. We will approach this in a systematic process. First, we conducted a series of cross-correlations between surface temperature anomalies and fossil fuel burning. We started with the Earth's surface temperatures (GLSTA, GOSTA and GSTA) leading the fossil fuel burning curve (CE). We found, by way of example, that the lagged value of GOSTA provides no statistically significant improvement in predicting the CE. Therefore, one cannot establish that global surface temperature changes (such as the GOSTA) "causes" the CE change outcome in the same manner as suggested above and discussed next below. For example, a fourth order autoregressive model, which fits the data starting in 1950 to {x t } using the previous year GOSTA {y t-1 }, results in a p-value of only 0.65 of the significance of the coefficient for GOSTA. We would definitely not conclude that GOSTA is causing outcomes for CE. Now we move onto demonstrating the effect or lack thereof that CE has on the GOSTA and by association on the GLSTA and GSTA (both not shown).
We first consider establishing a causal relationship between the carbon burning time series ({x t }) and the ocean surface temperature time series GOSTA ({y t }) in the latter half of the 20 th century. Here we model on the monthly scale by taking the average temperature anomaly for each year as y and the average monthly carbon emission created by taking the yearly value divided by 12. Our first model relates GOSTA to past values of GOSTA. Using Akiake's Information Criterion (Akiake, 1978) we select the order of the auto regression to be three, meaning that each year's GOSTA is related to the values from the last 3 years. The model is fitted using Gaussian maximum likelihood (Brockwell and Davis, 1991) [7] . We use the goodness of fit test (Fisher and Gallagher, 2012) [12] to verify that the AR model adequately models the autocorrelation in the GOSTA series; the p-value of 0.9464 indicates that we have found an adequately fitting stationary model for GOSTA model given in Equation (1): (1) = 0.00925 + 0.7676 −1 − 0.2575 −2 + 0.4264 −3 + The above model explains the dynamics of the observed series solely from past-observed values. We note that we could improve the predictive ability the above model by adding a linear trend component and would do so if our goal was to estimate the rate of change of temperature. However, to establish the (Granger) causal relationship we instead add the previous year CE value to the model. The resulting model is given by Equation (2): (2) = −0.054 + 0.6605 −1 − 0.2772 −2 + 0.3376 −3 + 0.0007 −1 + We can test for statistical significance of each fitted coefficient using the asymptotic normality of the Gaussian maximum likelihood estimators (Brockwell and Davis, 1991) [7] . In particular, we conclude with a p-value of 0.00009 that the coefficient of x t-1 is non-zero. In other words, the predictive model for GOSTA is significantly statistically better if the previous year CE is included. Model (2) explains the changes in temperature through a combination of autocorrelation (AR) and the impact of CE. We have thus established causality in the Granger sense. In this analysis, we selected the data beginning in 1950, since the warming trend in the latter half of the 20 th century is well established (cf. Figure 3). However, a similar analysis could be conducted starting at any point in the past. For example, the same analyses beginning each decade in the first half of the 20 th , i.e., in years 1900, 1910, 1920, 1930, 1940, and 1950, resulted  We next employ our fitted models to predict the global sea surface anomaly (GOSTA), the land atmosphere (GLSTA) and combined global surface temperature anomalies (GSTA) for 2015. Figure 5 shows the observed anomalies for years 1950 through 2014. The prognostic values are remarkably accurate. In the plots, we mark the predicted value from our fitted model, the upper and lower 95% prediction limits and the observed value for Year 2015. We find that the model, which uses both past values of the anomaly series and the previous year's carbon emission, provides excellent predictions for the observed anomalies of the ocean surface, atmosphere on land and the combination therein for 2015. This provides further empirical evidence of Granger Causality relating previous year carbon emissions to sea surface, atmosphere on land and the combination temperatures on a global scale. We utilized the 2014 and 2015 data because 2014 is the last year for which we were able to obtain Fossil Fuel Burning data from the website provider at the time that this study was conducted.
The results presented above provide a pathway and portend to future increases in global surface temperatures given anticipated increases in fossil fuel burning, GOSTA, GLSTA and GSTA as a function of CE. GOSTA is increasing at the rate of 0.

Discussion and Conclusions
Since the latter part of the 19th century up to the present, the reported overall rise in global surface temperatures has been viewed largely as an atmospheric phenomenon. However, we show that the global ocean is an important component in determining global surface temperatures; basically, the ocean is the dog and the atmosphere is the tail. Via an empirical, mathematical methodology, we are able to decompose the nonlinear and non-periodic data sets, and reveal the patterns of buried internal modes of variability of planetary temperatures over the past 167 years, from 1850 through 2018. We find periods of both cooling and warming, both in the ocean and the atmosphere over land, with natural variability ranging from seasonal to annual to inter-annual to multi-year to decadal to multi-decadal to centennial. We find that both the ocean surface and the air over land display non-linear trends depicting initially multi-decadal periods of cooling from the mid to late 19th Century, and then persistent warming throughout the 20th Century and into the 21st Century. Our calculated overall trends of the rates of warming differ significantly from the estimate of the IPCC or the NIPCC, with the oceanic rate 60% less that in the atmosphere. It is special note here that while the overall trends of planetary surface trends indicate a persistent and increasing warming (IMF Mode 11 in Figure 2), the 10 higher frequency modes of variability can modulate the overall temperature record from seasonal to centennial time scales. Nonetheless, while the cars on the train may be oscillating with different modulated frequencies (periods) and amplitudes, the train is moving forward. Empirical relationships between billions of tons of fossil fuel burning and the overall trends of the global surface-temperatureanomaly-time series in both the global ocean, the air above land and the combined time series emerge from our reduction of the non-stationary, non-linear data. Mathematical relationships between fossil fuel burning and surface temperatures in the oceans and over land are presented. The statistical relationship curves reveal strongly that there is a one-year phase lag between global carbon loading via fossil fuel burning and planetary surface temperature rise, different from Ricke and Caldeira (2014) [31] who proposed that planetary temperatures changed about a decade after burning. In fact, robust relationships are presented between the GLSTA, GOSTA and GSTA, and their past time series, and the CE time series. Via Granger Causality, these surface temperatures are predicted very accurately from fossil fuel burning a year earlier. Thus, the conclusion we reach is that we have proven that there is "attribution" between fossil fuel burning and climate warming. In 2007, the IPCC was awarded a Nobel Prize for its visualcorrelative-comparison of fossil fuel burning and temperature rise. However, a visual correlation does not prove "causality". In 2003 C.W.J Granger was awarded a Nobel Prize for his work in econometrics theory and applications. We invoked "Granger Causality" to attribute the overall trends in global surface warming to fossil fuel burning and carbon emissions (https://en.wikipedia.org/wiki/Clive_Granger). The conclusion reached in this study of overall planetary temperature trends versus the time series of fossil fuel emissions is that if present fossil fuel burning is not curtailed there will be continued warming of the planet in the future. We find that while the global ocean has an enormous capacity to absorb and retain heat, and thus act as a planetary heat reservoir, it is releasing stored heat in increasing amounts over time. Our findings are causally correlative in the sense of Granger Causality (Granger, 1980) and this study definitively links the overall trends in planetary surface temperature rise, to the overall upward trend in fossil fuel burning.