Shakun’s cauldron


Comment on Shakun et al. Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation (2012, Nature).


Samuel Schweikert*

Draft: August 31, 2016



Shakun et al. (2012) (1) found that the variations of the ‘global temperature’ lagged atmospheric CO2 partial pressure by 460 yrs ± 340 yrs during the last deglaciation – with the Southern Hemisphere temperature generally leading CO2 (by 620 ± 660 yrs) and Northern hemisphere temperature lagging CO2 (by 720 ± 330 yrs) – and repeatedly suggested (a) that this result could implicate CO2 as a major driver of deglaciations. We show that: a) the composition of the ‘global temperature’ index seemingly affects a much too large weighting factor to the Northern hemisphere (0.61) at the expense of the Southern hemisphere (0.39), which basis is not discussed, while it also could be considered surprising if most of the ocean heat content is in the Southern hemisphere; b) the reported lag for the ‘global stack’ is incompatible with the time fitting of the main inflexion points in the curves (‘function breakfit’); c) also, while the authors claim that an analysis based on function breakfit confirm their finding, it actually yields insignificant lags. In any event, the conceptual basis for the author’s conclusion regarding the possible role of CO2 is highly questionable because the very definition of a ‘mean lag’ regarding both time and space, i.e. the use of a ‘global temperature’ index, roots the approach to a linear model.

* Mechanical engineer, France.



Shakun et al. (2012) [hereafter S12] constitutes an attempt to relate the CO2 atmospheric content, during the last stadial-interstadial transition, to some variable designed to represent the warming of the entire globe surface. This study also addresses the different timing of the warming with the latitude, in the line of previous suggestions both that the Northern Hemisphere [NH] high latitudes must have warmed significantly later than the Southern Hemisphere [SH] during the last deglaciation, and that CO2, despite apparently not being the trigger of deglaciations, could have played an important role in them – e.g. Caillon et al. (2003) (2). Notwithstanding a certain (and certainly elegant) revival of the Milankovich’s thesis, since the rate of variation of the global NH ice volume, dv/dt, was considered instead of the instantaneous ice volume itself, v(t) – Edvardsson et al. (2002) (3), Roe (2006) (4) – it is also widely recognized that the warming of the Antarctic, more generally that of the SH, and of the tropical oceans, led the warming of the NH high latitudes, but the reasons for that remain largely unclear. We also note that, in the meantime, Schultz & Zeebe (2006) (5) have proposed a possible explanation to both the enigmatic lack of systematic correlation between stadial-interstadial alternations and the so-called ‘Milankovich parameter’ (insolation at 65°N in June) and to this SH lead (for Termination I).

To date, analyses of lagged correlations between CO2 and temperatures have generally been limited to Antarctica ice cores or, more rarely, to Greenland yet with Antarctic ice cores data for CO2, for different practical reasons mainly linked to accumulation rates, air trapping delay, water flow in ice cores but also, possibly, in situ decomposition of calcium carbonate (Smith et al. (1997) (6)). A few other studies focused on other places, yet in a single place or in a limited number of places or region(s), offering part conclusions or hints regarding spatiotemporal patterns, together with generally suggestive ones regarding the relation between CO2 and oceanic circulation (e.g. Shackleton (2000) (7)). One decisive common point is the limited number of locations (and phenomena) considered at the same time. This general issue is of particular concern when one wants to address the role of CO2 (either a main agent, or a more humble contributor, or essentially a consequence of the warming).

Because of energy transfers involved in lots of mechanisms – some of which highly nonlinear (some probably even unknown) -, one should not expect to find a strong correlation between atmospheric CO2 level and any local temperature. As an attempt to circumvent this problem (b), the authors have opted for what we consider a radical yet desperate way: building up a temperature index from linear combinations of a limited number of proxies, then defining a single, scalar variable representing some sort of ‘mean lead or lag’ between said index and atmospheric CO2. However the mere idea of defining a single lag already implies averaging along both time and space, i.e. to rely on such a ‘global temperature’ index but also to root one’s approach to a linear model.

In particular, a temporary stagnation, and even phases of marked decrease of the sea surface temperature [STT] in vast regions of the NH ocean, can be induced by the deglaciation process itself, because of the melting of the ice sheets and, mainly, due to surges of huge volumes of fresh water near in the Artic when icy walls collapse (see chapter ‘Global warming’, ‘global temperature’ indexes and deglaciation). In the meantime, CO2 outgasing of the SH ocean is expected after a few centuries of warming. However, because of the ‘CO2-temperature lag approach’, it is implicitly considered that the ulterior warming of the NH ocean would then have greatly needed a previous CO2 rise – not simply that latent heat or many other forms of energy are not measurable with a thermometer. In other words, the ‘CO2-temperature lag approach’, applied at a global scale, eliminates the null hypothesis in the first place.

This sort of ‘sophistry’ also brings us beyond the much more general tendency to consider that correlation means causation, because we often entirely ignore that a given variable C induces changes in both observed variables A and B whereas – even according to the authors’ thesis (see-saw via the AMOC) – it should be no surprise that, if the SH warms before the NH, some heat is afterwards transferred between them, with a given lag, while the SH oceans emits CO2, which process involves another lag (c). This S12 reasoning is not solely linked to the fact that a model cannot reproduce the observations without involving a dominating ‘green house effect’, it is beforehand a direct an unavoidable implication of the very question of the lag correlation as formulated. Yet in turn the quantitative question of the lag (see chapter ‘Lag analysis’) could reveal important because, for example, a short lag in the NH, together with a significant lead in the SH – and also between 30-60°N (see figure 2 in S12) – could probably be largely attributed to the initial cooling of the NH due to the deglaciation process itself.

Moreover, in S12, CO2 was also suggested to be an important driver of the deglaciation for Termination I because of two mistakes, one affecting the calculation of the ‘global temperature’ (see chapter ‘Global temperature stack’ and ‘hemispheric temperature stacks’), the other being an inappropriate method for the evaluation of the lag (see chapter Lag analysis) – this case illustrates the fact that ‘blind focus’ on correlation cannot lead to a decisive evidence, a fortiori when one is confronted with strongly nonlinear phenomena. Other problems affects the quality and/or representativeness of the proxies.


‘Global warming’, ‘global temperature’ indexes and deglaciation

The notion of a ‘global temperature’ is very often used in the climate research literature, despite the fact that any such thing is inevitably an index, based on numerous conventions – temperature being intensive ; lack of physical basis for the choice of an averaging method (even when shifting from microscopic to mesoscopic scale); strong non-linearities associated to phase changes; turbulence; chaotic internal variability; etc (8).

Notably in this field, one should add to the list more prosaic issues, among which the extremely limited number of samples (proxies), generally large uncertainties on both time and amplitude, and various factors affecting the proxies representativeness for such an exercise (d). In this case, we note that the CO2 curve exhibits some long holes, notably at the inflexion around 14,5 kyrs BP, where one cannot reasonably evaluate a lag of several centuries only. As for the temperature proxies, where some of them show large holes, this apparently is not accounted for in terms of uncertainties because the holes are filled in the first place, by interpolation every 100 yrs (‘Methods/stacking’ chapter in S12). See the grey curves in Figure S1 in the SM. We also note that two of them only start within the range 20-10 kyrs BP.

But although very important, the latter category of issues should not invite us to forget the former. Even greater care should logically be taken when one compares two climatic states as remote one from the other as modern climate and the peak of the previous glaciation. To this respect, it is hard to give credit to S12 argument, purportedly supporting their (area-weighted) averaging out of 80 proxies, that variations of the ‘global temperature’ can be successfully reproduced for the period 1880-2007 A. D. (chapter 4.1 in S12 SM) – another, perhaps less important issue, here, is that the proxies are assimilated to instrumental records at locations and altitudes necessarily distant from the proxies sites (e.g. at the nearest meteorological stations), and some sea or lakes surface temperatures to near surface air temperature.

Following S12 linear approach, it is implicitly assumed that, if one observed variable A leads another, B, then B is a consequence of factor A, and the amplitude of its effect is also implicitly proportional to the delay. We first note that, in certain circumstances and for a certain range of lags, the concept of warming is then somehow posed as a partly contradicting thing to that of deglaciation. Indeed, the huge amounts of calories involved in turning ice into water are not available for immediate production of a water or ice temperature increase. Besides, at least in the general case also, the input of huge volumes of cold fresh water in the NH oceans are not expected to make their temperature rise. This picture based on thermodynamic aspects alone gives a very poor presentation of the complex deglaciation process. In particular, very strong outbursts of freshwater are widely known to have been emitted from proglacial lakes, at times, through melting ice barriers, notably from Lake Aggasiz, with surges of up to ∼ 104 km3 in probably a few years (e.g. Teller et al. (2002) (9)), which have been modelled (Manabe & Stouffer (1997) (10)) to exhibit, in certain cases, a very strong power of reducing the thermohaline circulation and probably even stopping the formation of North Atlantic Deep Water for centuries, resulting in a strong decrease of the STT. But inverse effects have also been found in the case of a southern freshwater path, to the Mexico gulf: Clark et al. (2001) (11) – we note that S12 didn’t address such distinction when it comes to modelling (e). We also note that this approach is only apparently coherent with that proposed by Lindzen (1993b) (12) – where the characteristic ocean delay τ is proportional to the gain g, with g = 1 / (1 – f) where f is the total feedback factor – precisely because the latter, focusing on short or medium-term in modern context, ignores phenomena such as the melting and collapsing of ice sheets.

It also immediately appears that most regions are not represented at all, e.g. large parts of the SH ocean around 60°S, Russia (Figure 1 in S12). Among other problems with the proxies used by S12: the ‘global stack’ and the other temperature indexes involve averages between see surface temperatures and near air temperature over land; the proxies for temperatures in Greenland and in the Antarctic are implicitly accounted for as local ones when it comes to area-weighting (meaning a small weight) while they partly depend on phenomena which take place over oceans at lower (and much lower) latitudes. Also, at a significantly lower latitude than Greenland, S12 involves four pollen proxies in Alaska, which effectively yield comparatively late temperature increases (see S17 in S12 SM), but one could expect vegetation re-colonization to largely depend on the presence of an ice sheet. Let us also remark that some locations could be subjected to very large temperature variations due to specific local conditions. For instance: a relatively large number of proxies are close to New Zealand glaciers (f); the CH69-09 series, which effectively shows a singular curve, was taken at a place where oceanic currents, SST and even sea ice proximity could be highly variable (see figures S1 and S3 in the SM).


‘Global temperature stack’ and ‘hemispheric temperature stacks’

S12 built linear combinations of proxy temperatures records, notably: a ‘Global temperature stack’ [GS], a ‘Northern hemisphere stack’ [NHS] and a ‘Southern Hemisphere stack’ [SHS]. The area-weighting method is not addressed – S12 simply indicate (‘Methods/Stacking’ chapter) that the data is projected onto a 5° x 5° grid, then combined as area-weighted averages. However this should also be an important aspect, especially when the data is very poorly and unequally distributed. A few tests of robustness are presented in the SM but this cannot be a compensation for this lack of reporting. The SM (chapter 3) also points to an analytical code to be published in what apparently became Marcott et al. (2013) (13) but it is not clear if it concerns the area-weighting; however no further indication is provided, in this ulterior source or in its supplement, regarding this area-weighting method.

Nevertheless, we found that the GS can be simply re-obtained by a linear combination of the NHS and SHS with respective weights of 0.612 and 0.388 (constant). That is, with a much larger factor for NHS than for SHS. Even if one neglects the general issue of the ‘global temperature’ index (briefly addressed above), this seems a surprising composition: one could expect a satisfying area-weighting method to yield GS = (NHS + SHS) / 2 instead. Also, while temperatures don’t physically add up, one could try and reason on the base of heat content (even if this form of energy is obviously not generally conserved as such) and then consider that, because roughly 60% of the global ocean is in the SH, a more acceptable index for the ‘global temperature’ would correspond to GS = 0.4 NHS + 0.6 SHS (even if only a minority of the proxies were marine, contrarily to S12).

Figure 1 shows how GS evolves with those weighting factors. Superimposed are the C02 from Monnin et al. (2001) (14) and Lemieux-Dudon (2010) (15) data as well as a translation of S12 GS by 500 years (S12 460 yrs nominal lag is aimed, here, but the sampling of the ‘temperature stacks’ in S12 is every 100 yrs). We indicate the determination coefficients (R²), for the period 20-10 kyrs BP, between each of the different ‘global stacks’ and the C02 curve. We note that, while R² is 0.948 for the indicated GS, and climbs up to 0.965 when a 500 yr lag is accounted for, its value (without lag) also reaches 0.967 for the 40%-60% option and is as high as 0.960 for the 50%-50% option. See also S3 in the SM.

One basic difficulty comes from the fact that, while the paper indicate that S12 source for the CO2 record was Monnin et al. (2001) – with age correction following Lemieux-Dudon (2010) – every figures in S12 paper or SM apparently show Monnin et al. (2004) (16) data instead. Besides, contrarily to the 80 temperature proxy records, the CO2 record effectively used is not included in the S12 SM (Excel file). As the 2004 version of Monnin et al. essentially added data at around 12,5 kyrs BP and later, including one point in the curve inflexion around 12,5 kyrs BP, and as we deliberately neglect the late deglaciation period (see chapter ‘End of the warming’), while Monnin et al. (2001) version was already rather rich up to around 9 kyrs BP, we have chosen to use the 2001 version – only said additional point in the 2004 version, near 12,5 kyrs BP, was considered in this paper (see Figure 1).

Also, while the dispersions were treated by S12 in the context of Monte Carlo tests, it is difficult to address uncertainties in a manner directly comparable to the reported lags, because figures 2a, 3 and 4 in S12 affect vertical error bars only on the temperature ‘stacks’ and time error bars only on the CO2 curve. We postulate that the time error bars somehow account for the pairs of compared curves; as for the vertical ones, they seem to concern the temperatures ‘stacks’ only – it is not clear (SM, chapter 3) what was done in the case of the CO2 curve. In the following, we will consider the ± 340 yrs range reported by S12 as the relevant incertainty (for all GSs considered).

Figure 1b in S12 seems to indicate that the southern ocean is strongly under-represented – apart from proxies concentrated around 40°-45°S, the only ‘compensation’ for them is near the equator, that is, where the temperature variations are expected to be relatively small, and the area-weighting factor high. But Figure 2b should already inform the reader about a possible mistake: clearly, the ‘global’ domain is far from centred between the ‘SH’ and the ‘NH’ domains, and 460 is far from the mean of -620 and 720.

We also observe, though this is perhaps a coincidence, that the number of proxies for the NH and the SH are respectively 50 and 30, i.e. 62.5% and 37.5% of the set, which suggests that this possible problem could be linked, directly or indirectly, to the unequal distribution of the samples between the two hemispheres.

We finally note that an equally NH-oriented weighting is to be found in the modelled temperatures. From the curves we could reproduce, we found NH and SH respective weights of ∼0.62 and ∼0.38 in the composition of the modelled ‘global temperature’, consistent with the values inferred from the ‘stacks’. This suggests that, like for the ‘stacks’, the ‘global’ and ‘hemispheric’ temperatures are not the combination of a very large number of equally distributed local results covering all the map but, rather, that of a sampling with only 80 poorly and unequally distributed locations.


Lag analysis

Besides, when looking at trends shifts around 14,5 kyr BP and around 12,5 kyr BP (S12 C02 curve timing) on Figure 1, we note that the 500-yr translated GS show trends shifts at far distance in time from those of the C02 curve, while each of the ‘temperature stacks’ with no lag seem almost synchronized with CO2. This observation already indicates that the S12 method for evaluation of the lag doesn’t mean to fit any particular pair of inflexion points or any other mid term feature. However S12 make no comment about this (and, more generally, no comment about the choice and the implications of their method). As long as the correlation coefficient is maximized, it seems that they consider it to reflect a physically sounded lag.


S12 figure1 - correl 20-10 kyrs BP

Figure 1 – Atmospheric CO2 concentration (yellow dots), S12 ‘global temperature stack’ as such (blue) and with a 500-yr translation (red);  alternative ‘global stacks’ with respectively 50%-50% (green triangles) and 40%-60% (purple squares) weighting factors for NHS and SHS. Determination factor R² between each ‘stack’ and the CO2 curve, for the period 20-10 kyrs BP.  (g)


Pedro et al. (2012a) (17), a contemporary study which involved a partly similar correlation approach, made several remarks about some of the relevant methodological aspects (18). First, they justified their choice of a lag evaluation based on a maximization of the Pearson’s coefficient (for the curves and also for their derivatives), by the fact that they ‘consider this approach to be more robust than comparing dates of events (e.g. exclusively the start of the deglaciation) because its estimation is based on a much greater number of data points (see, e.g. Mudelsee, 2001)’, but it seems that, in their case, this is not found to be at the expense of ‘misidentification of the pre- and post-transition levels’, i.e. that their maximization of the correlation coefficient is apparently not incompatible with time fitting of the main inflexion points, contrarily to the result underlined by us above. Pedro et al. (2012a) also commented that ‘[a] caveat associated with [their] lag determination method is that it implicitly assumes that the maximum of the lagged correlation does in reality provide a valid estimate of the lag. Simple linear models of this form are widely used in previous studies of temperature and CO2 phasing (e.g. Ahn et al., 2004; Siegenthaler et al., 2005; Shakun et al., 2012); however, the result may somewhat change systematically if other models are used.’

We note that, in the case of S12 (and also for Pedro et al (2012a), except for the derivatives), one deals with nearly monotonous curves, contrarily to e.g. Caillon et al. (2003). In the latter case, lag-correlation – i.e. maximizing rxy – and ‘function breakfit’ – time fitting of the main inflexion points – would yield similar results, contrarily to S12 case. In addition to the previous observations, and still pointing to Figure 1, we note a that a ‘blind’ lag correlation, here, is particularly sensitive to the parting of the curves between ∼18-15 kyrs BP, where the CO2 curve first rapidly climbs before showing a negative inflexion, which pattern doesn’t appear in the GS curve, whereas it appears in the Antarctic stack (figure 2a in S12) and more generally in the SHS – this is where the discrepancies between the two hemispheric stacks are the highest (see figure 4b in S12), and we also note that, in this same period, another discrepancy linked to the choice of the radiocarbon date scale between IntCal04 and IntCal09 could account for a 110 yrs change of S12 evaluation of the lag (see chapter 5.5, notably figure s23, in the S12 SM).

[For pre-reviewers and other readers, here: especially if you’re interested in witchcraft or in mathematics…, surely you would not miss the final quiz! See Figure X4. In short, it seems that maximizing rxy simply yields absurd lag results when one deals with quasi-monotonous curves…]

S12 also provided an evaluation of the lag based on ‘function breakfit’ (chapter 6.2 in the SM). In the main paper, they rather surprisingly conclude that this independent method confirms their main analysis, although it neither yields comparable results (nor does it confirms the ‘mean results’ of Figure S25b in S12 SM) nor asks the same questions. As we have shown, their lagged correlation analysis is precisely not based on fitting inflexion points; it simply admits any shifting as long as it yields a higher mean value of the maximum correlation coefficient through a very wide window. Besides, S12 write in the SM (chapter 6.2) that ‘[t]he results [of the second method] indicate that CO2 was either synchronous within error or led temperature at the starts of the deglaciation, Bølling, and Younger Dryas, and Holocene (Figure S26; Table S2)’ but this is another ‘blind’ conclusion (and handling of averages): the calculated lag (CO2 each time leading) is respectively 200 ± 60 yrs, 190 ± 120 yrs, 200 ± 140 yrs and 490 ± 90 yrs, that is, each time within the range of uncertainties, except at 10 kyrs BP (for the latter time point, see chapter ‘End of the warming’). (h)

Yet an other method was considered by S12, which we address in the following, though no conclusion was drawn from it. See chapter 6.1 and Figure s25 in S12 SM. This supplementary lagged correlation analysis consists in: 1) setting the window size successively at discrete values of 3, 4, 5, 6 and 7 kyrs; 2) selecting the maximum value of the Pearson’s coefficient between the results obtained with every given lag every 100 yrs between -1000 yrs and +1000 yrs. Practically, this method is finally involved in a Monte Carlo analysis for quantification of what is presumably to be associated to time uncertainties around the ‘mean lag’. We first note that the choice of a [-1 kyr / + 1 kyr] lead/lag range was not discussed. Neither was the choice of discrete values (every 1 kyrs) for the windows sizes. However, figure 3a1 indicates that this double convention leads to very different results at the onset of CO2, at ∼17,5 kyrs BP. While one could ‘see’ a ∼1 400 yrs temperature lead for the mean value along window sizes between 1,4-3,4 kyrs: 1) the mean result for this same range of window sizes would automatically switch to an insignificant yet positive value, i.e. a CO2 lead, of 300 yrs; 2) using discrete values of 3, 4, 5, 6 and 7 kyrs for the window size instead yields an insignificant 120 yrs temperature lead. Figure S25b in S12 SM show different results around this time point, but as it is meant to be illustrating the method, it is not clear if it supposed to stand for the final version of the data set, methods and results. In any event, the extreme sensitivity of the results at this time point clearly appears in both figures. Another oddity in this additional S12 method is that apparently results beyond the [-1 kyr / + 1 kyr] are counted, in the average, as being + or -1 kyr exactly (figure S25b in S12 SM).

As we show in Figure 2 and in Figures S2 and S4 in the SM, it appears that, in many cases, very high correlation coefficients spread along a very large range of lags. Besides, particularly in the case of large windows, selecting only the maximum value of the coefficient is not a straightforward choice, because it implies the hypothesis of a nearly constant, direct, linear relation between two physically related variables, which variables could be considered almost independent of other variables. However, especially in the case of small window sizes, but even for very large windows, it gets clear that there is nothing close to a constant time response (at least, from what can be inferred from a lagged correlation analysis with only the max(rxy) considered) – see figure S25b in S12 SM, and figure 2 in this paper. In figure s25b in S12 SM, we note very large variations of the results: 1) with time, sometimes with negative lags (temperature leads) for older points and higher, positive lags, in the last phase of the deglaciation and later; 2) with the size of the window, notably sharp inversions when it is fixed at 3 or 4 kyrs.

It is also important to note that, in the present case, contrarily to most studies involving correlation analysis, rxy is generally very high in a very large range of tested lags. It so happens that, for large windows like in S12, accepting every lag result for, say rxy > 0,8 (which would be considered a strong correlation in other geological studies) would yield any lag in most or all of the [-1 kyr / + 1 kyr] band considered by S12. The range of accepted lags would only shrink when the window size gets smaller, but then more and more regions also appear where no lag is found. Moreover, one finds that the maximum rxy regions largely shift in the lead/lag scale (and obviously also in the time scale, i.e. location of the centre of the window) when the window size changes (see Figure S25b in S12 SM). In turn, blind averages based on the maximum values of rxy here and there picks lags values even if rxy is locally very low – possibly even negative, depending on the filter – (see Figure S25b in S12 SM) and affects to such lags the same weight that to lags associated to very large values of rxy, which is another caveat of such a method – also reminding that no correlation analysis should ever been done ‘blindly’.


S12 lag VS window size - 3pts - corr

Figure 2 – Pearson’s coefficient calculated around: a) 17,5 kyrs BP, b) 14,5 kyrs and c) 12,5 kyrs BP VS lag VS window size, 1) as S12, 2) with GS = (NHS + SHS) / 2 and 3) with GS = 0,4 NHS + 0,6 SHS. The dotted lines correspond to the maximum value of the Pearson’s coefficient, max(rxy). The full black line in figure 3a1 correspond to the max(rxy) in the range [-1 kyr / + 1kyr] of the lag, illustrating the effect of the arbitrary limitation of the tested lag in S12. The lag result, based on the positions of max(rxy), is indicated in the boxes (in yrs): in figures 1), the mean lag for window sizes of 3, 4, 5, 6 and 7 kyrs (S12 convention – marked by the black horizontal lines); in figures 1), 2) and 3), the lags obtained (respectively for the 50-50% and the 40%-60% compositions of the ‘global stack’), for the larger window size into which no more than one inflexion point enters, corresponding to the position of the white circles – the thin white line indicates the limit of the domain into which this condition is fulfilled (see also Figure S4 in the SM).


S12 themselves point out (chapter 6.3 in the SM) that ‘detrending records with ramp structures can increase the resolving power of lag correlations to determine lead-lag relationships’. Figure s27d, while not entirely clear (the scale for lags is not defined) show more concentrated lag results (though lower maximum values of the correlation coefficient) for with detrended (test) curves (and different maxima positions). And figure s28 in the SM indicate smaller lags (or leads) for the detrended curve, notably 380 yrs instead of 460 yrs for the nominal lag obtained for GS. However, for some reason, S12 finally based their publish results on the undetrended approach. They also judge (chapter ‘Methods/Robustness of results’) that ‘[t]he lead–lag relation between CO2 concentration and the global temperature stack is not significantly changed by detrending the time series to remove the deglacial ramp in each quantity’.

Following Ahn et al. (2004) (19) and Pedro et al (2012a), we also apply a lagged-correlation analysis to the derivatives – although only the nominal values are considered (like for the previous results). The same 3 different couples of weighting factors are considered for the composition of the GS. First we use the raw values for CO2 and for GS, then we apply on each curve, before calculation of the derivatives, running mean filters of 300 yrs, 500 yrs, 700 yrs, 900 yrs and 1 100 yrs. No lag is found, whatever the composition of GS, for the unfiltered curves and for a 300 yr filter. An increasing lag – with GS lagging CO2 – is obtained for an increasing degree of smoothing, yet with seemingly insignificant lags, appart from the original GS (as S12) combined with large filters of 900 or 1 100 yrs (as mentionned above, we consider the ± 340 yrs range reported by S12 as the relevant incertainty range, for all GSs considered). See Figure 3 and also Figure S6 in the SM.


Figure 3 – Lagged correlation based on the derivatives of CO2 and GS curves, for the 20-10 kyrs BP period. Number of occurences for which max (rxy) correspond to a given lag, for different compositions of the ‘Global stak’ (GS), and for different filters (smoothing applied on both CO2 and GS curves). Negative values of the lag mean that CO2 leads GS.


End of the warming

While the main method for the lag analysis in S12 consists in a direct lag-correlation between 20 and 10 kyrs BP (then neglecting events later than 10 kyrs BP), the second one, which results are purportedly consistent with that of the main method, actually yields insignificant lags, except at 10 kyrs BP. For this latter result, events as late as 8 kyrs BP were involved. See chapter 6.2 and Figure S26 in S12 SM.

However, not only the ‘global temperature’ (S12 GS) was already reaching its maximum around 10 kyrs BP, even if the deglaciation itself was not over yet (judging by the ∼28% remaining NH ice-sheet area), but the correlation coefficient had sharply dropped long ago, reaching zero near ∼11.5 kyrs BP and being even negative since ∼10.5 kyrs BP. See figures 3a, b, and d in S12, and S26 in S12 SM. The only thing that this independent method confirms, then, is that accounting for a CO2 warming effect even once there is no warming left to occur, and blindly mobilizing ulterior variations in a positive correlation research, in the corner of very large windows, while they are themselves persistent negative correlations, has the effect of increasing the apparent lag in an earlier phase.

We also observe that the averaging along time, in S12 supplementary method (chapter 6.1 in S12 SM), accounts for a long period posterior to 10 kyrs BP, after the global warming was over. Here also, events as late as around 8 kyrs BP have been considered, (judging from Figure S25(b) in the SM). However, as figure also S25(b) shows, it is clear that this ‘late’ component strongly draws up the S12 ‘composite lag’ in the direction of a longer temperature lag behind CO2.

Pedro et al. (2012a) note that ‘[t]here is a suggestion in both the Byrd and Siple Dome records of a larger lag at the onset of the Holocene’ but decided that ‘this period is not included in the time window of [their] lag assessment as [they] assume it is likely influenced by processes acting on the carbon cycle outside the Southern Ocean.’ In any event, one must raise the question why the generally increasing curves ceased to increase – instead of going on with what could have been a runaway process, in the case of a two-ways relation. At the very least, simply asking this question should be enough to prohibit the inclusion of period after ∼10 kyrs BP in the considered time range, because it clearly shows that the underlying hypothesis (constant, linear relationship between two independant variables) no longer stand, at least for this late period: this would prove that the observed variables (sometimes, at least) are themselves dominated by others (some of which probably unknown).

However, it also appears that, in the S12 model, while the influence of the atmospheric CO2 concentration on the ‘global temperature’ is expected to be by far dominant (figure 3e in S12), the direct influence of ‘global SST’ on the CO2 level seems to be very small and largely overwhelmed by other factors. In such a context, i.e. that of an almost one-way relation, S12 could have found themselves authorized to include a large post-deglaciation period. But one is confronted to a circular reasoning here, if the quantification of the CO2 effect is not brought from an independent, strongly based, analysis. Furthermore, once a post-deglaciation period is considered, one could also wonder why not to include an even more recent phase, from around 7 kyrs BP to modern times, during which the correlation is also negative with, this time, an increasing CO2 curve – EPICA Dome C, Monnin et al. (2004) – and a decreasing ‘global temperature’ (e.g. GISP2). Regarding this point, we also note that the long persistence of a relatively high CO2 level after the temperature peak is coherent with what was found by Solé et al. (2007) (20) for every type-B events in their classification of DO-events.


Underlying mechanisms

The use of such a simple model (linear regression), not to mention ‘blind analysis’, is certainly even more problematic if one considers that, based on entirely independent studies of physical mechanisms (Henry’s law, among others), the causal relation between the two variables should be (could be) considered to be both-ways (each way involving a given range of time responses). Obviously we don’t intend to propose a better model when dealing with such spatiotemporal scales and complex phenomenon, let alone when the relation could be both-ways; we simply note that S12 raise absolutely no question about the choice of their own model, while publishing claims about leads or lags for ‘global variables’ at such a scale.

Yet if ones now tries and addresses the causes (even with almost no clue), it seems rather surprising that, while the southern hemisphere oceans, representing about 60% of the oceans of the planet, are reported to having being warming since between around two millennia and a few centuries before the CO2 increased (while the range of uncertainty allows for the case of CO2 leading the SHS, this is not the case when local or regional temperatures are considered), it can be claimed that this CO2 increase accounts for a dominating part of the ‘global’ warming. At least, reconciling both claims obviously implies a two-way link. Simply observing a correlation is not enough; physical mechanisms must be proposed together with a (non blind) correlation, which should possibly explain not only the observed correlation but also the observed variations of the one variable which is considered as being mainly the cause and not the effect of the other. For example, if the greenhouse effect of the CO2 and its feedbacks is suggest to be the dominant factor, one still has to explain, in particular, where the CO2 emissions came from, in such quantities at such times. Then, either this explanation must show that the warming of the oceans is a minor cause of those emissions, or the correlation is no longer to be analysed under the hypothesis of a mostly one-way relation. Besides, if the outgasing of warming oceans cannot be ruled out as a strong agent in the process, one would once again have to consider the fact that the southern hemisphere oceans represent about 60% of the oceans of the planet. In such a context, one wonders what can be the meaning of using a correlation to an index affecting a 61%-39% factors distribution for the NH and the SH. Even under the additional hypothesis that another agent (e.g. strongly varying populations of marine species or terrestrial vegetation) strongly compensate said process anywhere, this NH-SH distribution remains simply because it is an obvious implication of the correlation as it is expressed.

Coherently with this manner of reasoning of S12, in S12 model the flow of fresh water, i.e. the melted ice quantity, is not a result but is treated like a free variable, that is, as an input which own causes do not need to be determined – as noted in our introduction, not only the amplitude of the flows but also their path should be decisive. As for the causes of both the temperatures variations and the CO2 variations, we know that they also remain poorly understood; we simply note that, in short, S12 model has no ‘primary driver’ itself. This lack of understanding, a fortiori touching, of the causal relation(s) between mechanisms, is obviously an important limitation when trying to analyse an observed correlation. Yet the CO2 itself intervenes merely as an input in the model. While at the very first glance (see Figure 4 in S12), the model exhibits an impressive ability to reproduce the ‘observed’ hemispheric ‘temperature stacks’, this should not come as a surprise if the dominant temperature-CO2 correlation were to be observed with a only a small lead or lag, and if CO2 is hypothesized as the main driver of the temperature. The former possibility is confirmed in this paper. The latter hypothesis appears clearly in figure 3e in S12, where the CO2 factor alone (together with its modelled feedbacks) is apparently expected to account for an overwhelming part of the global warming.

Unfortunately, S12 online material does not include the modelled temperatures. We reproduced them from screen grabs, in order to evaluate the lag between modelled and ‘observed’ temperatures. Figures S5 in the SM indicates that: the warming in S12 climate model is very late (> 1 kyr) at the onset of the deglaciation (Figure S5a); the modelled ‘global temperature’ also lags the ‘observed’ ‘global temperature’ by several centuries around 12,5 kyrs BP (Figure S5c), and lags CO2 by ∼600 yrs; whereas it leads both CO2 and the ‘observed’ GS at the time of the main fresh water surge, around 14.5 kyrs BP (see Figure S5b in the SM).

This large modelled lag at ∼17,5 kyrs and ∼17,5 kyrs is seemingly coherent with the very large sensitivity of the climate model used by S12 – the reversed tendency around ∼14,5 kyrs is apparently linked to a violent phenomena involving a huge fresh water input and a very strong oscillation of the modelled AMOC (see figures 3e, 3a and 4f in S12).


Supplementary material


 s12 - vue globale 81 séries 11

S12 CO2 VS NH series

Figure S1 – Temperatures and CO2 proxies in S12 (bottom: NH series only – if readers want to observe a given series VS CO2).


S12 - correl stacks S12 C1 et C2 entre 20-10 BP et 20-11.5 kyrs BP

Figure S2 – Lag for different ‘stacks’ and time window.

The grey area indicates the approximate range of uncertainty.


S12 spots courants 2

Figure S3 – SST and currents (left), wave height and ice proximity (sea ice in grey) near ‘CH69-09’ place, on March 20th, 2016.


sh12 lags correl coeff S12 2-5 kyrs 2

Figure S4 – Lag VS time and window size (S12 global stack).

The full black curves indicate the domains where the lag can be evaluated.


  S12 lag VS window size - 3pts - corr - model - annotations

Figure S5 – Lag VS time and window size for the modelled temperatures in S12. Pearson’s coefficient calculated around: a) 17,5 kyrs BP, b) 14,5 kyrs and c) 12,5 kyrs BP. 1) as S12, 2) with GS = (NHS + SHS) / 2 and 3) with GS = 0,4 NHS + 0,6 SHS. Full black lines: max(rxy) for the modelled temperatures; dotted lines: max(rxy) for the ‘measured’ temperatures.


derivées 2

Figure S6 – Lagged-correlation based on the derivatives. Pearson’s coefficients for different degrees of smoothing on the curves before calculation of the gradient. GSO is GS as in S12; GS1 = 0.5 NHS + 0.5 SHS; GS2 = 0.4 NHS + 0.6 SHS. See also Figure 3.


A few other figures for pre-reviewers and for my readers here


S12 for figure 1

Figure X1 – As Figure 1 but with curves superimposed.

Y axis: normalized according to std. dev. between 20-10 kyrs BP


S12 for figure 1b

Figure X2 – As Figure X1 but in parametric view (don’t miss Figure X4).


sh12 lags view 2

Figure X3 – Evaluation of the uncertainties regarding time, for CO2 only. Top: from (the so stricking) Figure S25(b) in S12 SM; uncertainties by me (see bottom). Bottom: visual evaluation of said uncertainties (+ ‘global stack’).


s12 test ABC 2

Figure X4 – Final quizzzz!

Maybe if you follow that, you’d agree it could be added in the paper (SM)…

In sum: take 3 ‘quasi-monotonous curves’, then try and put them in Shakun’s cauldron… Witch one lags witch other?…



The following notes are for pre-reviewers and other readers here. They are destined to be removed from the final draft – yet some of these notes could be kept if pre-reviewers suggest it.

Shakun et al. (2012) is under paywall: However:

  • Sylvestre Huet (see note (a), one of our most-beloved French warmists, graduated in litterature, sorta ‘journalist’ in a journal created by Sartre and owned (at that time) by a famous banker, apparently felt free to upload some of his foes’ copy. So… Here you are:
  • The supplementary material and Excel File for S12 is online on Nature website.

(a) This suggestion was even interpreted as a demonstration in the editorial of the same Nature release, and subsequently in various articles in grand public journals worldwide – e.g. S. Huet in Libération, France (this journalist had already drawn a similar conclusion from… Caillon et al. 2003):

(b)  See Shakun’s PhD dissertation, ref. given in footnote h.

(c)  The latter remark was made by N. Shaviv online:

(d)  Regarding this ‘secondary’ kind of issues, some important remarks have been made by Don J. Easterbrook online: and

(e)  Which is even more surprising, knowing that Peter U. Clark is 2nd author here, and was director of Shakun’s PhD (see ref. in footnote h).

(f)   See also the remarks made by Don J. Easterbrook online – ref. indicated in footnote d.

(g)  Mudelsee’s method (Break function regression. A tool for quantifying trend changes in climate time series. Eur. Phys. J. Special Topics 174, 49–63 (2009)) should probably be used for ages determination, before submission.

(h)  It its perhaps also worth noting that, in what seems to be an earlier version of the S12 paper – Chapter 4.4 in Shakun, J. D. & Clark, P. U., The Role of CO2 During the Last Deglaciation, in preparation for submission to Science, in J. D. Shakun, PhD dissertation submitted to the Oregon State University (presented May 21, 2010) -, J. D. Shakun and P. U. Clark wrote: ‘We use lag correlations to address the phasing of CO2 and global temperature and find CO2 exhibited a small but insignificant lead of 45 ± 130 years from 19-11 ka, increasing to 835 ± 240 years when considering the entire time interval of ice sheet melting from 19-6.5 ka Figure 4.3)’. Apparently the set of 77 temperature proxies was roughly similar, but the temperature records were ‘combined as unweighted averages’ (also the IntCal09 and MARINE09 calibrations were used for radiocarbon dates instead of IntCal04). However, according to chapter 5.2 in the S12 SM, this difference should probably have changed the variance of the ‘global stack’ (by 7-8%, but not after ∼11 ka) but not its structure. This nevertheless confirms that accounting for the ‘late period’ (see chapter ‘End of the warming’) leads to a very strong increase of such a ‘mean lag’. It is also interesting that, despite of an even less significant lag reported in this draft (before 11 kyrs BP), the two authors of this published draft concluded that the result ‘implicate CO2 as a major driver of global deglaciation and highlight its importance in controlling global temperature’, which seems a strikingly surprising conclusion, not only because observed correlation is obviously no evidence for causation, and not only because this result is that of an insignificant lag, which the two authors themselves explicitly acknowledged, but also because, even under the hypothesis of an actually positive lag, a low value of the time response doesn’t seem strongly compatible with a high value of so-called ‘climate sensitivity’, and because, in this case also, the moderate lag of ∼300 yrs in the NH (together with a ∼700 yrs lead in the SH) could probably be largely attributed to the initial cooling of the NH due to the deglaciation process itself (see next chapter). This is not clear, however, how many proxies were actually considered at this stage, because figure 4.1 in this draft shows only 68 proxies. What is clear, though, is that the main authors should have been aware that their ‘mean lag’ result is probably very sensitive to certain parameters in the lagged correlation computation, especially to the window size but probably also to the inclusion of certain proxy records in the composition of the data set. However, it does not seem that this important lack of robustness in their own standards was properly addressed.



This work was not funded. The author wishes to thank in advance every people interested in an helpful review of this draft paper, before possible submission to Nature or another journal.



(1)   Shakun, J. D., Clark, P. U., He, F., Marcott, S. A., Mix, A. C., Liu, Z., Otto-Bliesner, B., Schmittner, A. & Bard, E. Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation. Nature 484, 49-55 (2012).

(2)  Caillon, N., Severinghaus, J. P., Jouzel, J., Barnola, J. M., Kang, J. & Lipenkov, V. Y. Timing of Atmospheric CO2 and Antarctic Temperature Changes Across Termination III. Science 299, 1728-1731 (2003).

(3)  Edvardsson, S., Karlson, K. G & Engholm, M. Accurate spin axes and solar system dynamics: Climatic variations for the Earth and Mars. Astronomy & Astrophysics 384, 689-701 (2002).

(4)   Roe, G. In defence of Milankovitch. Geophysical Research Letters, 33 L24703 (2006).

For my readers. Edvadsson et al. (2002) (top); Roe (2006) (bottom). As Lindzen put it, ‘Roe had it (…) never seen such a nice correlation in this field’ (his words as I remeber them):

Milankovitch Edvarsson Roe 2

(5)   Schultz, K. G. & Zeebe, R. E., Pleistocene glacial terminations triggered by synchronous changes in Southern and Northern Hemisphere insolation: The insolation canon hypothesis. Earth and Planetary Science Letters 249, 326–336 (2006).

(6)  Smith, H.J., Wahlen, M. Mastroianni, Taylor, K. C,  The CO2 concentration of air trapped in GISP2 ice from the last glacial maximum-holocene transition. Geophysical Research Letters, 24,1-4. (1997).

(7)   Shackleton, N., The 100,000-Year Ice-Age Cycle Identified and Found to Lag Temperature, Carbon Dioxide, and Orbital Eccentricity, Science 289, 1897 (2000).

(8)  Essex, C., McKitrick, R. & Andresen, B. Does a global temperature exist? Journal of Non-Equilibrium Thermodynamics (2006).

(9)  Teller, J. T., Leverington, D. W., and Mann, J. D., Freshwater outbursts to the oceans from glacial Lake Agassiz and their role in climate change during the last deglaciation, Quaternary Science Reviews 21, 879–887 (2002).

(10)   Manabe & Stouffer (1997)    (ref. tbi)

(11)   Clark, P. U., Marshall, S. J., Clarke,, G. K. C., Hostetler, S. W., Licciardi, J. M. and Teller, J. T., Freshwater Forcing of Abrupt
Climate Change During the Last Glaciation, Science 293, 283-287 (2001).

(12)  Richard S. Lindzen, Constraining Possibilities Versus Signal Detection, Natural Climate Variability on Decade-to-Century Time Scales, U. S. National Research Council, 1995.

(13)  Marcott, S. A., Shakun, J. D., Clark, P. U. & Mix, A. C., A Reconstruction of Regional and Global Temperature for the Past 11,300 Years, Science 339, 1198 (2013).

(14)  Monnin et al (2001)    (ref. tbi)

(15)  Lemieux-Dudon (2010)    (ref. tbi)

(16)  Monnin et al (2004)    (ref. tbi)

(17)  Pedro, J. B., Rasmussen, S. O. and van Ommen, T. D., Tightened constraints on the time-lag between Antarctic temperature and CO2 during the last deglaciation, Climate of the Past 8, 1213–1221 (2012).

(18) The review process of both papers largely covered each other. The arguments in Pedro et al. (2012a) which we recall here correspond to objections we already had found during our own analysis. There are cited here mainly in order to illustrate the fact that, contrarily to S12, some authors in the field openly raised such questions about the method for lag evaluation.

(19) Ahn, J., Wahlen, M., Deck, B. L., Brook, E. J., Mayewski, P. A., Taylor, K. C., and White, J. W. C. A record of atmospheric CO2 during the last 40,000 years from the Siple Dome, Antarctica ice core, J. Geophys. Res., 109, D13305, doi:10.1029/2003JD004415 (2004).

(20)  Solé, J., Turiel., A. & Llebot, J. E. Classification of Dansgaard–Oeschger climatic cycles by the application of similitude signal processing. Physics Letters A 366, 184–189 (2007).


3 Réponses à “Shakun’s cauldron

  • Don’t miss what Nature wrote in its editorial:

    […] as Jeremy Shakun at Harvard University in Cambridge, Massachusetts, and his colleagues show on page 49, carbon dioxide does drive atmospheric warming. Uncontroversial stuff, perhaps, yet the link continues to be questioned by people who would play down the risks of human greenhouse-gas emissions. The queries re-emerged in 2006, when former US vice-president Al Gore showed a graph of historical carbon dioxide levels and temperature in his movie, An Inconvenient Truth, and was accused of glossing over the relationship between the two. So let there be no confusion now: the new study confirms that, as Earth emerged from the last ice age some 19,000 to 10,000 years ago, rising global temperatures were preceded by increased global carbon dioxide in the atmosphere — a result that emphasizes the role of carbon dioxide in driving global change in the present day. This relationship is a foundation stone of climate science and of policies to regulate greenhouse-gas emissions, and it is solid.
    Quelle surprise! the climate sceptics may shout — scientists find proxy data and use a computer model to get the answer they wanted, to seal the conspiracy. Then let the second paper this week show that modern science does anything but offer self-serving results to support existing ideas. […]

  • Stabilizing the draft for reviewers

    At this stage (September 7, 2016), I’ll try and keep the draft unchanged, for a few months maybe. One reason is obviously to provide a stable version for reviewers (in particular, references and notes must keep their numbers), but another one is that, while it’s an endless stuff…, I would also certainly be tempted to remove or strongly alter some developments whenever I feel there are of secondary importance, etc. Let the reader judge that. Of course, every modification, removal or addition can be suggested. Don’t be afraid or ashamed: the paper we’re talking about is mainly bs anyway (somebody at WWWT even said it’s unfalsifiable science). In any event, my suggestions for modifications, like yours, will only be made in the comments; at a certain stage, I’ll post a new stable version, in a new page.


    Supplementary material for Shakun et al. (2012) (text and Excel file):
    ‘Temperature Stacks’ and proxies in the Excel file. CO2 missing (as indicated). My Excel sheets can be provided on requests, but I have many versions so try and be precise about what you need.

    Shakun’s PhD compilation of papers and drafts, mentionned in note (h), can be found here:;sequence=1
    See chapter 4 (also have a look at Figure 3.2 in chapter 3).

    Shakun’s list of publications here:

    Nature editorial mentionned in note (a) and quoted my previous comment:

    See also the FAQs here :

    Some precision is made in there, regarding the freshwater path(s) in their model (see my remark associated to note (e)) – emphasis mine:

    Was the climate model tuned to match the proxy data?
    Mostly not. The climate model was driven by four factors: (1) variations in the earth’s orbit around the sun (known precisely), (2) the retreat of ice sheets (known reasonably well from geologic records), (3) rising greenhouse gases (known from ice-core records), and (4) freshwater poured into the ocean associated with the melting of the ice sheets. We know approximately how much freshwater went into the ocean during the deglaciation from sea-level records, but it is less certain where exactly the freshwater went in through time (e.g., Antarctica vs. North America; Hudson R. vs. Mississippi R.). The location is important because it controls if/how the AMOC responds to the freshwater addition. Several reasonable freshwater scenarios were tried, and the one that yielded North Atlantic climate oscillations in best agreement with the Greenland ice core records was used. It is worth reiterating that while freshwater-AMOC changes affect how heat is redistributed around the planet, they have relatively little impact on the global mean temperature. ”

    Also (emphasis mine):
    What does the global temperature lag behind CO2 mean?
    The climate does not equilibrate instantly with rising CO2, because there are a couple big sources of thermal inertia – the ice sheets and ocean, which take a long time to melt away and warm up. Therefore, it would make sense that global temperature lags behind CO2 during the last deglaciation. This is the same reason that even if CO2 concentrations were held constant today, the earth would still continue warming – the so-called “warming in the pipeline”.
    As I put it in the draft, the NH lag would then probably make sense anyway…, even if CO2 had no effect at all. But see how Shakun turns its own ‘explanation’! I especially like this one…

  • Missing references (14), (15) and (16).

    For those, of course I was essentially interested in the data (I’m not trying and judge their methods).

    (14) Monnin et al (2001)

    Ref. Monnin, E., Indermühle, A., Dällenbach, A., Flückiger, J., Stauffer, B., Stocker, T. F., Raynaud, D. and Barnola, J. M., Atmospheric CO2 concentrations over the last glacial termination, Science, 291, 112-114 (2001).


    (15) Lemieux-Dudon (2010)

    My mistake, it’s Lemieux-Dudon et al. (2010).

    Ref. Lemieux-Dudon, B., E. Blayo, J.-R. Petit, C. Waelbroeck, A. Svensson, C. Ritz, J.-M. Barnola, B.M. Narcisi, and F. Parrenin. 2010. Consistent dating for Antarctic and Greenland ice cores. Quaternary Science Reviews, Vol. 29, Issues 1-2, pp. 8-20 (2010).


    (16) Monnin et al (2004)

    Ref. Monnin, E., Steig, E. J., Siegenthaler, U., Kawamura, K., Schwander, J., Stauffer, B., Stocker, T. F., Morse, D. L., Barnola, J.-M., Bellier, B., Raynaud, D. and Fischer, H., Evidence for substantial accumulation rate variability in Antarctica during the Holocene, through synchronization of CO2 in the Taylor Dome, Dome C and DML ice cores. Earth and Planetary Science Letters, 224, 45-54 (2004).


    NB: don’t miss the Readme’s (and references) associated to the data.

    Note that, as late as in 2015, Shakun was still uncapable (probably he still is today) to correctly cite Monnin et al. (2004) when using its data; he still indicates Monnin et al. (2001) instead. See figure 2d (light-blue curve = Monnin et al. (2004) data) and reference 63 in:

Répondre à Sam Annuler la réponse

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *