In today’s post, I’m going to re-examine (or more accurately, examine de novo) Ed Cook’s Mt Read (Tasmania) chronology, a chronology recently used in Gergis et al 2016, Esper et al 2016, as well as numerous multiproxy reconstructions over the past 20 years.
Gergis et al 2016 said that they used freshly-calculated “signal-free” RCS chronologies for tree ring sites except Mt Read (and Oroko). For these two sites, they chose older versions of the chronology, purporting to justify the use of old versions “for consistency with published results” – a criterion that they disregarded for other tree ring sites. The inconsistent practice immediately caught my attention. I therefore calculated an RCS chronology for Mt Read from measurement data archived with Esper et al 2016. Readers will probably not be astonished that the chronology disdained by Gergis et al had very elevated values in the early second millennium and late first millennium relative to the late 20th century.
I cannot help but observe that Gergis’ decision to use the older flatter chronology was almost certainly made only after peeking at results from the new Mt Read chronology, yet another example of data torture (Wagenmakers 2011, 2012) by Gergis et al. At this point, readers are probably de-sensitized to criticism of yet more data torture. In this case, it appears probable that the decision impacts the medieval period of their reconstruction where they only used two proxies, especially when combined with their arbitrary exclusion of Law Dome, which also had elevated early values.
Further curious puzzles emerged when I looked more closely at the older chronology favored by Gergis (and Esper). This chronology originated with Cook et al 2000 (Clim Dyn), which clearly stated that they had calculated an RCS chronology and even provided a succinct description of the technique (citing Briffa et al 1991, 1992) as authority. However, their reported chronology (both as illustrated in Cook et al 2000 and as archived at NOAA in 1998), though it has a very high correlation to my calculation, has negligible long-period variability. In this post, I present the case that the chronology presented by Cook as an RCS chronology was actually (and erroneously) calculated using a “traditional” standardization method that did not preserve low-frequency variance.
Although the Cook chronology has been used over and over, I seriously wonder whether any climate scientist has ever closely examined it in the past 20 years. Supporting this surmise are defects and errors in the Cook measurement dataset, which have remained unrepaired for over 20 years. Cleaning the measurement dataset to be usable was very laborious and one wonders why these defects have been allowed to persist for so long.
The “RCS” Chronology: Calculated vs Archived
To avoid burying the lede in further details, the diagram below shows the difference between my RCS chronology calculation from measurement data (left panel) and the “RCS” reconstruction used in Esper et al 2016 (as well as Gergis et al 2016 and many other studies), extended prior to 1000AD from the underlying chronology archived at NOAA in 1998. Versions of both are available back to 1500BC (and will be shown later), but only the past 1500 years is shown in the diagram below, which is intended to illustrate the differences.
Despite the difference in visual appearance, the two versions are very highly correlated (r=0.57 over nearly 3600 years). However, the RCS chronology shows a long-term decline, with 20th century values returning to values reached earlier in the millennium, and less than values of the first millennium. On the other hand, the Cook chronology (archived by Esper) has flattened values earlier in the millennium, such that late 20th century values appear somewhat anomalous.
Figure 1. Mt Read versions (converted to SD units): left – RCS chronology calculated from Esper et al 2016 archive, edited to remove defects as described in Postscript; right – Esper et al 2016 reconstruction for Mt Read, extended prior to 1000AD with NOAA arfilt reconstruction (see discussion below).
Cook’s RCS Description and Age Profile Curve
A tree ring chronology is, in its essence, an estimate of an annual growth index after allowing for juvenile growth (since ring widths decline with age.) “Traditional” standardization fit a growth curve to each tree individually. But if growth rates varied between century, such techniques transfer variability over time to variability between trees. This type of problem is well known in statistics as fixed and random effects – techniques long advocated at Climate Audit – but these techniques are unfamiliar to tree ring scientists, who describe such phenomena in tree rings in coarse and artisanal terms.
Cook was keenly aware of the lack of low frequency variability in “traditional” standardization ( a technique that he had used in his original publication of Tasmanian data in 1991 and 1992). The issue also concerned Briffa, who, in two influential publications (Briffa et al 1991, 1992), advocated the use of a single growth curve for each site (rather than individual growth curves for each tree) as a means for retaining centennial variability, but, like Cook, being unable to express the issue in formal statistics. The single-growth-curve technique was subsequently labeled as “RCS” standardization. The technique has problems if the dataset is inhomogeneous between sites (as many are). There are numerous CA posts structuring chrobnology development in terms of fixed and random effects.
Back to Mt Read: Cook et al 2000 clearly and unambiguously stated that they used a single age profile curve (“RCS”) to allow for juvenile growth in chronology development:
In an effort to preserve low-frequency climatic variance in excess of the individual segment lengths, we applied the regional curve standardization
(RCS) method of Briffa et al. (1992, 1996) to the data.
Cook described his RCS technique in straightforward terms and (of considerable assistance for subsequent analysis) included a figure showing the age profile curve.
The RCS method requires that the ring-width measurements be aligned by biological age to estimate a single mean growth curve that reflects the intrinsic trend in radial growth as a function of age… The mean series declines in a reasonably systematic and simple way as a function of age. The simplicity of this ring- width decline with increasing age indicates that the RCS method may work well here…Therefore, a simple theoretical growth curve has been fit to the mean series. This RCS curve is shown in Fig. 2A as the smooth curve superimposed on the mean series.
Cook’s Figure 2A can be closely replicated from a quality-controlled version of the Mt Read dataset, as shown below, where my replication matches Figure 2A down to fine details. This match accomplishes three things. It confirms that differences between my RCS chronology and the archived Cook version do not arise from different age profile curves. Second, it also confirms the validity of my quality control editing for the defects and errors in the archive. While the edits were motivated by other inconsistencies, the age profile curve calculated on the data prior to quality control for stupid defects does not match. Third, the data is better fitted by a Hugershoff curve, a common form in dendro analysis (y= A+ B*x^D *exp(-C*x) ), than by a negative exponential (another common form). I accordingly used the Hugershoff fit in my calculation of the RCS chronology.
Figure 2. Mt Read age profiles. Top – from Cook et al 2000; bottom – calculated from archived measurement data. Blue- fit from Hugershoff curve y= A+ B*x^D *exp(-C*x). A slight discrepancy in fit on the far tail won’t effect results since so few measurements at this age.
Reverse Engineering
If Cook’s chronology wasn’t an RCS chronology, then what was it? Its lack of centennial-scale variability strongly suggests that it employed some form of fitting to each individual core. Cook’s first articles on Tasmania (Cook et al 1991, 1992) had employed smoothing splines – a technique which made it impossible to recover centennial variability, as Cook understood very clearly.
Another standardization technique employed at the time was (supposedly) “conservative” standardization with “stiff” negative exponential curves, a style of standardization used by Cook’s colleague, Gordon Jacoby – a technique which still standardized on individual cores (rather than one curve for the entire site). As a guess, I did a Jacoby-style standardization on the Mt Read data, comparing to the Cook et al 2000 Figure 3 chronology (archived as the NOAA “original” chronology) – see figure below.
Figure 3. Mt Read versions. Left- using “traditional” standardization (negative exponentials for each core); right – NOAA “original”. The NOAA “original” version was illustrated in Cook et al 2000 Figure 3
It doesn’t match exactly, but the general similarity of appearance, combined with the dissimilarity to a freshly calculated RS chronology, convinces me that something went awry in Cook’s calculations. While he clearly intended to calculate an RCS chronology, I am convinced that the Cook et al 2000 chronology was calculated using some variation of “traditional” standardization (in which curves were separately fit to individual trees/cores), a conclusion supported by the lack of centennial scale variability.
Compare RCS Chronology to Mean Ring Width Series
The RCS chronology calculated in this post is closely related to the series of mean ring widths illustrated in Cook et al Figure 1A.
Figure 4. Top – Cook et al 2000 Figure 1A showing mean ring width; bottom – mean ring width series calculated from quality-edited measurement data.
Cook observed that his method (claimed to be RCS) had “removed some long-term growth” variation, but claimed (incorrectly in my evaluation) that it had preserved “much of the century-scale information”:
Comparing the mean series in Figs. 1A and 3A, it is clear that the RCS method has removed some long-term growth variations in the standardized chronology, but much of the century-scale information has been preserved.
This observation can be dramatically re-stated in light of the preceding analysis. The RCS calculation (this post) preserves virtually all of the long-term growth variation of the mean ring width series.
The “Original” and “Arfilt” Versions
Both Cook et al 200 and the associated NOAA archive contain two closely related versions of his chronology: the “original” version (shown in Figure 3) and the “arfilt” version (shown in Figure 7). To foreclose potential issues arising from differences between these versions, I’ve shown both variations against an RCS chronology in the figure below, this time showing values back to 1500BC (as in the NOAA archive), more than doubling the coverage shown in the first figure above.
Figure 3. MT Read chronologies from -1500BC to 1991AD: left – RCS chronology (this post); middle: “original” chronology; right: “arfilt” chronology.
Relative to the RCS chronology, both the “original” and “arfilt” chronology appear very similar to one another and differ from the RCS chronology through their lack of low-frequency variability. In this view, one can discern a general correspondence in decadal features of the smoothed (yellow) version, while noticing that the Cook versions have flattened out high early values apparent in the RCS chronology.
Cook described the calculation of his “arfilt” chronology as follows:
Simple linear regression analysis was used to transform the Lake Johnston Huon pine tree-ring chronology into estimates of November-April seasonal temperatures. Prior to regression, both the tree-ring and climate series were prewhitened as order-p autoregressive (AR) processes… The tree-rings were modelled and prewhitened as an AR(3) process, with the AR order determined by the minimum AIC procedure (Akaike, 1974). The AR coefficients and explained variance of this model are given in Table 2. Note that the majority of the persistence is contained in the first coefficient (ARl = 0.397), which basically reflects an exponentially- damped response to environmental inputs such as climate. In total, this AR(3) model explained 21.8% of the tree-ring variance. In contrast, the November-April average temperatures were modelled as an AR(1) process that explained 10.5% of the variance (Table 2).
Reading this article today, it’s impossible to see any value added through the addition of the arfilt procedure to the “original” chronology, however this was actually calculated.
Even though the arfilt series is even further from an RCS chronology than the “original” chronology, this was the version used in both Gergis et al 2016 and Esper et al 2016. (This can be proven by digital comparison of the data.) Gergis et al 2016 contained 10 values for the arfilt series (1992-2001) not included in the NOAA archive and not supported by the the present measurement data archive, which ends in 1991. If the additional values arise from the inclusion of fresh tree ring data, it’s hard to understand why this data wouldn’t have some impact on the chronology up to 1991, but these values remain unchanged – even to the third decimal place.
Gergis et al 2016
As noted in the introduction, Gergis et al re-calculated tree ring chronologies from the underlying measurement data for all sites using an RCS variation (“signal-free detrending”) recently developed at the University of East Anglia:
All tree-ring chronologies were developed from raw measurements using the signal-free detrending method, which improves the resolution of medium-frequency variance…
The method has been presented as yet another recipe, with no attempt to place it in a broader statistical context.
When the method was applied to tree ring sites in the Gergis network, it frequently resulted in much enhanced 20th century relative to previously published results, as, for example, for Celery Top East (shown below), where the published version had little centennial variability, but, as re-calculated by Gergis et al, had a very pronounced increase in the 20th century.
Figure 4. Celery Top East chronologies. Left- the (unarchived) BLT (East) chronology as illustrated in Allen et al (2001); right – the Celery Top East chronology archived in connection with Gergis et al 2016. The series name includes identifiers BLT and RFR, both sites discussed in Allen et al 2001. The difference between versions appears to arise from methodology, not data.
However, for Mt Read (and Oroko), Gergis used prior chronology versions, ostensibly for consistency with published results” – though such inconsistency had not troubled them in cases where they obtained elevated 20th century values. They described their choice of older versions as follows:
The only exceptions to this signal-free tree-ring detrending method were the New Zealand silver pine tree-ring composite (Oroko Swamp and Ahaura), which contains logging disturbance after 1957 (D’Arrigo et al. 1998; Cook et al. 2002a, 2006), and the Mount Read Huon pine chronology from Tasmania, which is a complex assemblage of material derived from living trees and subfossil material. For consistency with published results, we use the final temperature reconstructions provided by the original authors that include disturbance-corrected data for the silver pine record and regional curve standardization for the complex age structure of the wood used to develop the Mount Read temperature reconstruction (Cook et al. 2006).
It is evident that their decision to use a prior version of Mt Read was only made after examining results of the fresh chronology, which must have been similar to the results calculated above i.e. with elevated values in the first millennium and early second millennium. Such an election, only made after getting (presumably) adverse results is yet another example of data torture (Wagenmakers, 2011, 2012), for which Gergis’ proffered rationale (“consistency with published results”) is both flimsy and inconsistent with handling other series.
The version used instead in Gergis et al 2016 was, as noted above, identical to Cook’s arfilt series between 1000 and 1991. Gergis did not use values in the archive from prior to AD1000, but her dataset included 10 additional values (1992-2001). Presumably these additional values were provided by Cook, but, if they were calculated from additional data, one would have expected that earlier values would at least be somewhat impact (and not identical to three decimal places).
Esper et al 2016: Reconstruction and Measurement Data
Esper et al 2016 stated that they used the “most recent” reconstruction from each site. However, as noted above, the Tasmania reconstruction (values from 1000-1991 AD) in their archive ends in 1991 and is a shortened form of Cook’s arfilt reconstruction (archived at NOAA in 1998). Values before 1000AD were not archived for Tasmania, though earlier values were archived for many other sites.
Esper et al asserted that Cook’s Tasmania reconstruction had been produced from an RCS chronology, but this assertion does not appear to have been verified by Esper and/or other lead authors.
The Mt Read measurement data in the Esper et al 2016 archive appears to be more or less identical to measurement data archived at NOAA in 2002 (ausl024). Both contain many booby traps for the unwary – all of which ought to have been corrected long ago. Because of these defects, the archive is not readable using the R package dplR. Cook’s booby traps have several levels. Level 1: Cook’s measurement data contains in two different units. Each core segment has a trailing (final) value of either -9999 or 999; cores with a trailing value of 999 are in units that are 10 times larger. I do not recall any other measurement data archive which isn’t consistent. I don’t understand why this inconsistency isn’t tidied in the archive. Level 2: in all or nearly all measurement data archives, each core segment has a separate ID through adding a single-character suffix to the core ID. The end of each segment is marked by a trailing value of -9999 or 999. However, Cook’s archive contains multiple segments with the same ID, blowing up data assimilation. In my read program, problematic cores could be picked up by seeing if the core, after assimilation, still included a trailing value. There were about 15 such cores – in each case, I manually added a suffix to one segment of the core, thereby differentiating the segments to accommodate a yield. This took a lot of time to diagnose and patch. Level three: this was the trickiest. As a cross check, I examined the distribution of average ring widths by core segment and found that average ring widths for six cores (KDH89, KDH16B, KDH85, KDH61A, KDH64, KDH83) were about 10 times the average ring widths. Using the “native” values of the Esper archive, the age profile curve didn’t match the Cook et al figure. I concluded that the values for these cores in the Cook archive needed to be divided by 10 to make sense. I did this semi-manually and saved a clean measurement dataset, which I then applied. After this cleaning, the age profile match improved considerably. Note that these cores were all in early periods (mostly first millenium). Thus the editing of the data reduced first millennium values.
Conclusion
An RCS chronology calculated according to the stated methodology of Cook et al 2000 yields an entirely different result than that reported by Cook. In my opinion, Cook, like Gergis et al 2012, did not use the procedure described at length in the article – in Cook’s case, he did not use the RCS procedure described in the article as a method to preserve low-frequency variability. In my opinion, Cook’s chronology was most likely produced using a variation of “traditional” standardization that did not preserve low frequency variability.
Cook’s chronology has been used over and over in multiproxy studies: Mann et al 1998, Jones et al 1998, Mann and Jones 2003, IPCC AR4, Mann et al 2008; most recently, Gergis et al 2016 and Esper et al 2016. Despite its repeated use, one can only conclude that no climate scientist ever looked closely at Cook’s actual chronology, a conclusion circumstantially supported by the persistence of gross errors in the Cook measurement data, even in the Esper et al 2016 version, issued more than 20 years after the original measurements.
The actual RCS chronology for Mt Read has elevated values in the late first millennium and early second millennium. Gergis et al evidently calculated such a chronology and, in another flagrant instance of ex post cherry picking, decided to use the ancient Cook chronology, which turns out to have been erroneously calculated (like Gergis et al 2012, one might add). Use of the Mt Read RCS chronology and Law Dome series would obviously lead to substantially different results in the medieval period where Gergis only used two proxies.
Postscript
A reader pointed out that Allen et al (2014) reported an update to Mt Read. I should have taken note of this, but take some consolation in the fact that Esper et al 2016 didn’t take note of this update, despite claiming to have used the most recent series from each site and having a common co-author (Cook).
Allen et al took cores from 18 trees. Their reported chronology was not an RCS chronology, but a chronology based on fitting individual series – a technique that Cook had purported to avoid in order to preserve centennial variability. Allen et al describe the use of negative exponential curves, conceding that this technique “potentially loses some centennial time scale variability” relative to Cook et al 2000:
The updated chronology was constructed from the new samples and the previously crossdated material (Buckley et al. 1997; Cook et al. 2000) and is based on individual detrended series with a mean segment length greater than 500 years (Cook et al. 2000). All samples, including those previously obtained, were standardised… using fitted negative exponential curves or linear trends of negative or zero slope, with the signal-free method applied to minimise trend distortion and end-effect biases in the final chronology (Melvin and Briffa, 2008). The use of negative exponential and linear detrending preserves low-frequency variations due to climate consistent with the ‘segment length curse’ (Cook et al., 1995b), but potentially loses some centennial time scale variability that had been preserved in the Mt Read chronology based on regional curve standardisation (Cook et al. 2000).
My surmise is that the technique used by Allen et al is probably similar to, if not identical, to the technique that Cook et al 2000 actually used. The figure below compares Allen et al 2014 Figure 2C to the RCS chronology of this post. The two series are highly correlated, but the Allen chronology has flattened out the long-term decline of the RCS chronology.
Figure ^. Left – Allen et al 2014 Figure 2C. Right – RCS chronology (this post) overlaid onto Allen Figure 2C.
Once again, it seems implausible that Allen et al 2014 did not first attempt an RCS chronology in order to directly update the Cook et al 2000 chronology, deciding to use an individually fitted chronology only after inspecting an RCS chronology. In addition, the start of their figure (400BC) neatly and perhaps opportunistically excludes high values in the immediately preceding interval.
References:
Allen et al 2014. Continuing upward trend in Mt Read Huon pine ring widths–Temperature or divergence? Quat Sci Rev. pdf
Cook et al 2000. Warm-season temperatures since 1600 BC reconstructed from Tasmanian tree rings and their relationship to large-scale sea surface temperature anomalies. Clim. Dyn. pdf
Gergis et al 2016. Australasian Temperature Reconstructions Spanning the Last Millennium. Journal of Climate.
Esper et al 2016. Ranking of tree-ring based temperature reconstructions of the past millennium. Quat Sci pdf.