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 nonlinear 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, interannual, multiyear, decadal, multidecadal and centennial modes, riding a top overall warming trends. Our calculated overall rates of warming differ significantly from the estimates of the Intergovernmental Program on Climate Change [1,2]. We find the oceanic warming rate to be less than twothirds 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 econometricsbased 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 1year phase lag of global temperatures to fossil fuel burning both on land and in the sea.
Attribution of anthropogenic influence in climate change; Climate change; Climate warming; Fossil fuel burning; global surface atmospheric temperature anomalies; Global warming; Patterns of climate variability; Sea and ocean surface temperature anomalies
The Earth system absorbs the Sun’s incoming shortwave radiation and reradiates, stores or exchanges it at different rates via natural processes. For a planetary condition of thermal equilibrium, the amount of total outgoing, longwave radiation would equal the total amount of incoming, shortwave 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 peerreviewed 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 [28], 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. [9] and Ekwurzel et al. [5], have employed global energybalance coupled climatecarboncycle models and have attempted to assess the change in atmospheric CO_{2} and CH_{4} concentrations, radioactive 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 is 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 wellknown and documented surfacetemperatureanomalydata. 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 inkind. 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.
Figure 1: The monthly averaged time series of GSTA (top), GLSTA (middle) and GOSTA (bottom) January 01 1850July 31 2018 from the Climatic Research Unit and the UK Meteorology Office, Hadley Centre:
http://www.cru.uea.ac.uk/cru/info/.
The climate system consists of nonlinear and nonperiodic 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 quasistationary. “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. [10], employed the hilbert transform [11,12], 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 modemixing problem, Wu and Huang [13], 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?.
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 [14], 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 [15], the definition of a trend is stated to be “the general drift, tendency or bent of a set of data”. Chatfield [16], defines “trend” as “a longterm change in the mean”. Thus, difficulties arise in determining what “long term” means. For example, a 50year cycle would appear to be a trend in a 20year data set; but it would be shown as a 50year 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. [17]. In that paper, the EMD decomposition included the nonoscillation “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. [17], there was not an applicable definition of what constituted the overall “trend” of any continuous time series of data, especially of a nonlinear 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 [14]. With this definition of a trend, we will employ the EEMD method [13], 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.
Following the works of Wu et al., [17] and Pietrafesa and Gallagher [18], an updated analysis of global surface temperatures though 2018 is presented below. In figure 1 we presented the monthly averaged time series of global surface temperature anomaly data from 1850 through July 2018 for the landbased atmosphere (GLSTA), the Global Ocean Surface (GOSTA) and the combined (GSTA). The three time series are continuously compiled by the Climatic Research Unit (CRU), United Kingdom (UK) Meteorology Office, Hadley Centre:
http://www.cru.uea.ac.uk/cru/info/. At the beginning of the CRU, UK time series, global records of surface and ocean temperatures display a smattering of increases and decreases reflecting the relative dearth of spatial and temporal data of those early years. However, the relative amplitudes of the temperature swings on land and at the ocean surface are insync and in range, suggesting a level of coherency around the planet both in the atmosphere over the continents and across the entire planetary global ocean.
Figure 2 shows eleven IMF modes in the GSTA (also true of the GLSTA and GOSTA, so not shown) with the 11
^{th} being the overall trends of the landbased atmosphere, the ocean surface and the combination of the two. The IMF modes of the GOSTA and GLSTA are essentially the same for the ocean surface and for the atmosphere over land, save for differing amplitudes. The IMF modes have their respective physical basis. Mode 1 is several monthly variability. Mode 2 is seasonal variability. Mode 3 is the motion of the Earth on its axis of rotation. Mode 4 is the annual signal. Mode 5 is interannual variability. Mode 6 is the dominant signal of the El Nino Southern Oscillation (ENSO). Mode 7 is the quasiSolar Cycle (1012 years), centered about 11 years. Mode 8 reflects both the North Atlantic Oscillation (NAO) and of the 22year Solar Cycle. Mode 9, 6070 years is the Atlantic Meridional Overturning Circulation Belt (MOC), described by Cunningham et al., [19]. Mode 10, 105110 years represents the global thermohaline circulation conveyor belt [20,21]. Mode 11 is the gravest mode or overall trend of the 167year time series of total data (Table 1).
Figure 2: The EEMD Mode decomposition of the 1850  2017 time series of the GSTA (Shown in the Top Panel in figure 1. There are 11 IMF modes presented in the increasing period from top to bottom. The 11th mode represents the overall record length trend.
IMF

1

2

3

4

5

6

7

8

9

10

11

GLSTA GOSTAGSTA

2 M

3M

6M

1Y

23Y

57Y

1012 Y

2022Y

6070Y

105110Y

167Y Trend

Table 1: Intrinsic Mode Functions (IMFs) shown in figure 2, in Months (M) or Years (Y).
We note in figure 2, that for the GSTA, the IMF modes 2, 3 and 10 have amplitudes of 0.4
^{o}C, modes 4 and 6 have amplitudes of 0.3°C, modes 5, 7 and 9 have amplitudes of 0.2°C and mode 8 amplitude of 1°C. Mode 11 ranges between +/ 0.5°C. Thus, all the IMF modes contribute to the total time series in temperature amplitudes in a nominally equitable manner across the range of 0.1 to 0.4°C. This finding suggests that the Planet’s internal, natural modes of variability contribute all contribute significantly to monthly averaged surface temperatures and thus cannot be ignored. Thus, periods of relative warming and or cooling occur naturally by the positive and negative disposition of the natural variability, all riding atop overall atmospheric and oceanic trends and presented in the table below.
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.
Figure 3: (left panel) The overall trends of the GSTA (blue), GOSTA (red) and GLSTA (yellow) time series. The GOSTA raw data are shown in figure 2 and represents the GLSTA and the GSTA (both not shown). (Right panel) The time rates of change or 1st derivatives, of the trends of the GLSTA, GOSTA and GSTA time series.
In figure 3 (left panel), the trend of the GLSTA shows a record length total rise of 1.22°C and the GOSTA displays a total rise of 0.67°C. The GSTA or combined ocean surface and land surface rises 0.88°C. The ocean surface temperature GOSTA has risen at a much slower rate than the atmosphere over land, the GLSTA. This is an excellent example of the power of the EEMD IMF decomposition. The full time series of the GOSTA displays an overall beginning to end of the series increase of 0.75°C, the GLSTA shows a rise of 2.91°C and the GSTA indicates an overall increase of 1.25°C. If connecting lines from start to end of the three time series had been drawn, the overall trends would have been greatly overestimated, demonstrating the failing of traditional methods of computing a trend, which are to create regression lines. Thus, the IPCC rates of rise of GSTA (IPCC, 2007) overestimated the actual rises in temperatures over land by 0.72°C and by a combined airsea overestimate of 0.37°C. The IPCC straightline, regression slope estimates were 0.05°C/decade over the full GSTA temperature record, 0.07°C/decade over the prior 100 years, 0.13°C/decade over the previous 50 years and 0.18°C/decade over the latter 25 years of the record. It is of note that the GSTA temperature time series is the same time series, from the identical source, <
http://www.cru.uea.ac.uk/cru/info/> that we used in our analyses.
In figure 3 (right panel), we present the time rates of our trends (Figure 3, left panel). We find that the rates of warming from 1850 through 2017 are 7.13(102) °C/decade for the GLSTA, 4.27 (102) °C/decade for the GOSTA, and 5.29 (102) °C/decade for the combined GSTA. However, the 1st derivatives of the trend curves are all positive and curving upward indicating that the rate of warming is accelerating. More recently, from 1955 to the present the overall warming trend of the GSTA (Figure 3) has been about 8.83(102) °C/decade. Clearly, the rate of warming of air over land has globally been 167% greater than that of the surface of the global ocean. The temperature anomaly trends all range from negative values at the onset of each of the time series in 1850 to highly positive values at the end of 2017. Thus, employing straight regression lines is not the mathematically or physically plausible manner in which trends of nonlinear time series can be estimated. Only an adaptive method such as EEMD [13] can do so, as explained above [22]. 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.
In figure 4 we present the fossil fuel burning time, Carbon Emissions (CE) series from 17512014, provided via <
http://cdiac.ornl.gov/trends/emis/tre_glob.html>. As shown in figure 4, from 1751 until 1850, global CO
_{2} production was negligible, having ranged from 3MMTs in 1751 to 54MMTs 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 950MMTs. 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 US 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 causeandeffect relationship or attribution.
Figure 4: Global CO_{2} emissions from fossilfuel burning, cement manufacture and gas flaring: 17512014. confirming sources also include: T. Boden and R. Andres, The Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee 378316290, USA; and G. Marland, The Research Institute for Environment, Energy and Economics, Appalachian State University, Boone, North Carolina 286082131, USA.
We next employ the Granger Causality Test (GCT) [23] 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 “Grangercause” Y if it can be shown, through a series of ttests and Ftests on lagged values of X (and with lagged values of Y also included), that those X values provide statistically significant information about future values of Y. The idea is that if series {xt} helps to cause series {yt}, then past values of {xt} should improve predictions of series {yt}. This type of causality is established by first modeling {yt} in terms of the past values of {yt} through an Autoregressive (AR) process, then adding past values of series {xt} to create a second model. The original procedure proposed by Granger requires assuming stationarity of series {yt}. If the second model is statistically better than the first, then one has established causality in the granger sense. Testing for granger causality in global temperatures has been considered by Attanasio et al., Pasini et al. , Attanasio et al. , Triacca et al. [2427], 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 nonstationarity in the temperature series. Our approach differs in that we consider only anomaly series {yt} and one emission series {xt}, and we find a stationary AR model for the anomalies, which we justify our modeling by use of a time series goodnessoffit test [28]. 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 crosscorrelations 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 {xt} using the previous year GOSTA {yt1}, results in a pvalue 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 ({xt}) and the ocean surface temperature time series GOSTA ({yt}) 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 [29] 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 [30]. We use the goodness of fit test [28] to verify that the AR model adequately models the autocorrelation in the GOSTA series; the pvalue of 0.9464 indicates that we have found an adequately fitting stationary model for GOSTA model given in Equation (1):
????????=0.00925+0.7676????????10.2575????????2+0.4264????????3+???????? (1)
The above model explains the dynamics of the observed series solely from pastobserved 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):
????????=0.054+0.6605????????10.2772????????2+ 0.3376????????3+0.0007????????1+???????? (2)
We can test for statistical significance of each fitted coefficient using the asymptotic normality of the Gaussian maximum likelihood estimators [30]. In particular, we conclude with a pvalue of 0.00009 that the coefficient of xt1 is nonzero. 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 20th, i.e. in years 1900, 1910, 1920, 1930, 1940 and 1950, resulted in pvalues, 0.00005, 0.0002, 0.00002, 0.0002, 0.003, 0.009, for the coefficient of xt1, respectively. In each case we conclude the model using the previous year CE is statistically better than the AR (3) model alone; the carbon emission series significantly improves predictions for GOSTA over the model based solely on past GOSTA values. Regardless of starting time, in the 20^{th} century, one finds statistical evidence of a causal relationship between CE and GOSTA.
The coefficient for xt1 in model (2) is as follows. For each additional onemillion metric tons of Carbon Emissions (CE), we estimate an increase in Global Ocean Temperature (GOSTA) of 0.0007 °C. The average increase in carbon emissions per year for years 1950 through 2014 is about 130 metric tons per year. Our model estimates an increase of 0.007+ °C for each 10 metric ton increase in carbon emissions. The highly statistically significant coefficients of xt1 for land based atmospheric temperatures and total global surface temperatures are 0.0117 and 0.0087, respectively.
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.
Figure 5: Global surface temperature anomalies from 19492014 with granger causality predicted (X) versus actual temperatures (O) for 2015 and upper and lower 95% confidence limits. Top left panel is for the GOSTA. Top right panel is for the GLSTA. Bottom panel is for the GSTA.
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.0007°C/ Million Metric Tons of CE. GLSTA is increasing at the rate of 0.00117°C/Million Metric Tons of CE. Thus, the GSTA is increasing at the rate of 0.00087°C/Million Metric Tons of CE. Presently we are on a trajectory to reach 104 Million Metric Tons/year of CE by years 2020 to 2022. While these numbers are less foreboding than the straightline estimates of both the IPCC and NIPCC, they are nonetheless very noteworthy.
Since the latter part of the 19
^{th} 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 nonperiodic 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 interannual to multiyear to decadal to multidecadal to centennial. We find that both the ocean surface and the air over land display nonlinear trends depicting initially multidecadal periods of cooling from the mid to late 19
^{th} century and then persistent warming throughout the 20
^{th} century and into the 21
^{st} 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 surfacetemperatureanomalytime series in both the global ocean, the air above land and the combined time series emerge from our reduction of the nonstationary, nonlinear 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 oneyear phase lag between global carbon loading via fossil fuel burning and planetary surface temperature rise, different from Ricke and Caldeira [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 [3241].
In 2007, the IPCC was awarded a Nobel Prize for its visualcorrelativecomparison of fossil fuel burning and temperature rise. However, a visual correlation does not prove “causality”. In 2003 CWJ 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 [23] and this study definitively links the overall trends in planetary surface temperature rise, to the overall upward trend in fossil fuel burning.