Anomaly Detection in Paleoclimate Records using Permutation Entropy
Joshua Garland, Tyler R. Jones, Michael Neuder, Valerie Morris, James W. C. White, Elizabeth Bradley
AArticle
ANOMALY DETECTION IN PALEOCLIMATERECORDS USING PERMUTATION ENTROPY
Joshua Garland ∗ , Tyler R. Jones , Michael Neuder , Valerie Morris , James W. C.White and Elizabeth Bradley Santa Fe Institute University of Colorado at Boulder Institute for Arctic and Alpine Research University of Colorado at Boulder Department of Computer Science ‡ These authors contributed equally to this work. * Correspondence: [email protected]: December 3, 2018;
Abstract:
Permutation entropy techniques can be useful in identifying anomalies in paleoclimate datarecords, including noise, outliers, and post-processing issues. We demonstrate this using weightedand unweighted permutation entropy of water-isotope records in a deep polar ice core. In one regionof these isotope records, our previous calculations [1] revealed an abrupt change in the complexityof the traces: specifically, in the amount of new information that appeared at every time step. Weconjectured that this effect was due to noise introduced by an older laboratory instrument. In thispaper, we validate that conjecture by re-analyzing a section of the ice core using a more-advancedversion of the laboratory instrument. The anomalous noise levels are absent from the permutationentropy traces of the new data. In other sections of the core, we show that permutation entropytechniques can be used to identify anomalies in the raw data that are not associated with climatic orglaciological processes, but rather effects occurring during field work, laboratory analysis, or datapost-processing. These examples make it clear that permutation entropy is a useful forensic tool foridentifying sections of data that require targeted re-analysis—and can even be useful in guiding thatanalysis.
Keywords: paleoclimate; permutation entropy; ice core; anomaly detection.
1. Introduction
Paleoclimate records like ice and sediment cores provide us with long and detailed accounts ofEarth’s ancient climate system. Collection of these data sets can be very expensive and the extractionof proxy data from them is often time consuming, as well as susceptible to both human and machineerror. Ensuring the accuracy of these data is as challenging as it is important. Most of these cores areonly fully measured one time and most are unique in the time period and region that they “observe,”making comparisons and statistical tests impossible. Moreover, these records may be subject to manydifferent effects—known, unknown, and conjectured—between deposition and collection. Thesechallenges make it very difficult to understand how much information is actually present in theserecords, how to extract it in a meaningful way, and how best to use it while not overusing it.An important first step in that direction would be to identify where the information in a proxyrecord appears to be missing, disturbed, or otherwise unusual. This knowledge could be used to flagsegments of these data sets that warrant further study: either to isolate and repair any problems or,more excitingly, to identify hidden climate signals. Anomaly detection is a particularly thorny problemin paleorecords, though. Collecting an ice core from polar regions can cost tens of millions of $US.Given the need for broad geographic sampling in a resource-constrained environment, replicate coresfrom nearby areas have been rare. This has restricted anomaly detection methods in paleoscience tothe most simple approaches: e.g., discarding observations that lie beyond five standard deviationsfrom the mean. Laboratory technology is another issue. Until recently, for instance, water isotopes a r X i v : . [ phy s i c s . d a t a - a n ] N ov of 15 in ice cores could only be measured at multi-centimeter resolution, a spacing that would lump yearsor even decades worth of climate information into each data point. The absence of ground truth is afinal challenge, particularly since the paleoclimate evidence in the core can be obfuscated by naturalprocesses: material in a sediment core may be swept away by an ocean current, for instance, and icecan be deformed by the flow of the ice sheet. These kinds of effects can not only destroy data, but alsodeform it in ways that can create spurious signals that appear to carry scientific meaning but are infact meaningless to the climate record.Thanks to new projects and advances in laboratory techniques, the resolution challenge is quicklybecoming a thing of the past [2]. In addition to the recently measured West Antarctic Ice Sheet (WAIS)Divide ice core (WDC), many additional high-resolution records are becoming available, such as theSouth Pole Ice Core (SPC) [3] and the Eastern Greenland Ice Core Project (EGRIP) [4]. Replicate data ison the horizon as well, which may solve some of the statistical challenges. The SPC project, for instance,will involve dual analysis (i.e. two replicate sticks of ice) from three separate subsections of the ice core.However, replicate analyses of deep ice cores beyond a few hundred meters of ice will not occur anytime soon, let alone multiple cores from a single location . Thus, a rigorous statistics-based treatmentof this problem is still a distant prospect for deep ice cores. In the meantime, information-theoreticmethods—which can work effectively with a single time-series data set—can be very useful. In previouswork, we showed that estimates of the Shannon entropy rate can extract new scientific knowledgefrom individual ice-core records [1,7]. There were hints in those results that this family of techniquescould be more generally useful in anomaly detection. This paper is a deeper exploration of that matter.To this end, we use permutation entropy techniques to study the water-isotope records in theWAIS Divide ice Core (WDC). The resolution of these data, which were measured using state-of-the-artlaboratory technology [2,8], is an order of magnitude higher than traditionally measured deep ice corepaleoclimate records(e.g.,[9]), and also representative of new datasets that are becoming available. Weidentify abrupt changes in the complexity of these isotope records using sliding-window calculationsof the permutation entropy [10] (PE) and a weighted variant of that method known as weightedpermutation entropy (WPE) that is intended to balance noise levels and the scale of trends in thedata [11]. Via close examination of both the data and the laboratory records, we map the majorityof these abrupt changes to regions of missing data and instrument error. Guided by that mapping,we re-measured and re-analyzed one of these segments of the core, where the PE and WPE resultssuggested an increased noise level and the laboratory records indicated that the processing had beenperformed by an older version of the analysis pipeline. The PE and the WPE of this newly measureddata—produced using state-of-the-art equipment—are much lower, and consistent with the valuesin neighboring regions of the core. This not only validates our conjecture that permutation entropytechniques can be used to identify anomalies, but also suggests a general approach for improvingpaleoclimate data sets in a targeted, cost-effective way.Permutation entropy is a complexity measure: it reports the amount of new informationthat appears, on the average, at each point in a sequence. In the context of a time series, thistranslates to a measure of how information propagates forward temporally. This has implications forpredictability [12,13] among other things; indeed, these measures have been shown to converge to theKolmogorov-Sinai Entropy under suitable conditions [14,15], as described in Section 2.2. The goal hereis not to measure the complexity of paleoclimate data in any formal way, however. Our intent, andthe underlying conjecture behind our approach, are more practical: abrupt changes in a quantity likepermutation entropy suggest that something changed, either in the system or the data. The followingsection describes the data and methods that we use to demonstrate the efficacy of that approach. There are some replicate data available from shallow snow pits [5] or closely drilled shallow ice cores [6], but these do notextend beyond the past few hundred years, at most. of 15
Figure 1.
Water-isotope records for δ O (top panel) and δ D (bottom panel) taken from the WestAntarctic Ice Sheet (WAIS) Divide ice core (WDC).
2. Materials and Methods
As a proof of concept for the claim that information theory can be useful in detecting anomaliesin paleorecords, we focused on water-isotope data from the West Antarctic Ice Sheet (WAIS) Divideice core (termed WDC hereafter) [8]. These data, which consist of laser absorption spectroscopymeasurements of the isotopic composition of the ice, are proxies for local temperature and regionalatmospheric circulation during the past 67,000 years. The specific variables that we consider are theratios of the heavy and light isotopes of hydrogen ( H / H ) and oxygen ( O/ O). Time-series tracesof these ratios, which we will identify as δ D and δ O in the rest of this paper, are shown in Figure 1.The details of the experiment and data processing are as follows. The water-isotope data wererecorded at a rate of 1.18 Hz (0.85s intervals). Ice samples were moved through a cavity ring-downspectroscopy continuous flow analysis (CRDS-CFA) system [2] at a rate of 2.5 cm/min, yieldingmillimeter resolution. Registration of the laser that tracks the depth along the core can be a challenge,particularly in the brittle zone, a section of ice from approximately 577 − ≈ − (cid:104) by the International Atomic Energy Agency (IAEA): δ = ( R sample / R VSMOW − ) where R is the isotopic ratio. The δ O and δ D symbols refer to fractional deviations from VSMOW,normally expressed in parts per thousand (per mille or (cid:104) ).Converting the δ D and δ O data from the depth scale to an age scale is a major undertaking,as well as a significant source of uncertainty. Constructing the age-depth model for the WAISDivide core required dozens of person-months of effort. In the upper ≈ of 15 by combining annual dating of high-resolution ( < eradicate the annual variations below a certain depth. Similar issuesarise in age models for sediment cores, where marine organisms mix the material from different timeperiods. These, and other effects that deform the timelines (e.g., brittle zone deformities introducedinto the depth scale during the registration process) are critically important considerations when oneis doing any kind of time-series analysis on paleoclimate data sets.The final stage in the data-processing pipeline addresses the uneven temporal spacing of thesamples. Because of compression, 0.5 cm of ice—the width of a sample—represents roughly 1/40th ofa year of accumulation near the top of the WDC and roughly 1.4 years near the bottom . One couldcertainly calculate PE or WPE on such a sequence, but the timeline of the results would be deformedbecause of the temporally irregular spacing of the data points. Since we are specifically interested inthe time scales of the changes in complexity, we needed to perform our calculations on data that areevenly spaced in time, so we interpolated the δ D and δ O data to an even 1/20th year spacing using acombination of downsampling and linear splines. While this is standard practice in the paleoclimatefield, it is not without problems if one plans to use information-theoretic methods on the results, asdescribed in more detail in Section 2.2. Indeed, irregular temporal spacing is a fundamental challengefor any kind of time-series analysis method, and the effects of interpolation methods used to regularizethe temporal spacing of these data sets must be considered very carefully if one is doing sophisticatednonlinear statistics on the results. The following section describes how this plays out in the context ofWPE calculations on ice cores; several recent papers offer more-general treatments of this matter forother kinds of paleoclimate data [21–25].
Permutation entropy (PE)—a measure of how much new information appears at each timestep, on the average, across the computation window—was originally proposed as a “naturalcomplexity measure for a time series” [10]. Since its conception, PE has been shown to convergeto the Kolmogorov-Sinai Entropy [14,15], or the metric Entropy [26,27] under a variety of conditions,and depending on features (e.g., stationary, ergodic) of the underlying generating process. It has alsobeen shown to correlate with intrinsic predictability in a variety of time series [12,13]. For the purposesof this paper, we simply view permutation entropy as a measure of the temporal complexity of a timeseries and we treat abrupt changes in that complexity as a signal of possible anomalies, naturallyoccurring or data related.Permutation entropy calculations combine traditional information-theoretic probability estimateswith ordinal analysis, a process by which time-ordered elements of a time series, also known as delayvectors [ x t , x t + τ , . . . , x t +( m − ) τ ] ≡ X m , τ t , are mapped to an ordinal pattern , or “permutation,” of thesame size. If the first, second, and third points in a time series were X = [ x , x , x ] = [ − ] ,for instance, the corresponding permutation would be φ ([ − ]) =
231 since x ≤ x ≤ x .We will refer to this mapping as φ : R m → S m , which accepts an m -dimensional delay vector and Water isotope records in ice cores are affected by diffusion [20], which smooths and attenuates the signals. Firn diffusion inthe upper ≈ This spacing is also affected by changes in yearly snow accumulation, obviously. That effect actually turns out to be a majoradvantage, as we describe in [1]; together with the actions of diffusion, it creates a link between the accumulation rate andthe information content that can be used to back the accumulation record out of WPE calculations. of 15 returns the corresponding permutation π ∈ S m , the set of permutations of length m . After using φ to convert the time-series data into a series of permutations, one computes the permutation entropyusing the formula: H ( m , τ ) = − ∑ π ∈S m p ( π ) log ( π ) (1)where p ( π ) is the estimated probability of each permutation π ∈ S m occurring in the time series,calculated as follows: p ( π ) = (cid:12)(cid:12) { t | t ≤ N − ( m − ) τ , φ ( X m , τ t ) = π } (cid:12)(cid:12) N − ( m − ) τ (2)Here, | · | is set cardinality.Notice that the mapping φ , as defined above, does not distinguish between [ − ] and [ − ] ; both will be mapped to the same 231 permutation. This may be inappropriate ifthe observational noise is larger than the larger-amplitude trends in the data, or if there is meaningfulinformation in the amplitude of the signal [11]. The now-standard way to address these concerns , weighted permutation entropy (WPE) [11], effectively emphasizes permutations that are involved in“large” features and de-emphasizes those whose amplitudes are small relative to the features of thetime series. This is accomplished by weighting each permutation by its variance: w ( X m , τ t ) = m m ∑ j = (cid:16) x t +( j − ) τ − X m , τ t (cid:17) (3)where X m , τ t is the arithmetic mean of the values in X m , τ t . The weighted probability of a permutation isdefined as: p w ( π ) = ∑ t ≤ N − ( m − ) τ w ( X m , τ t ) · δ ( φ ( X m , τ t ) , π ) ∑ t ≤ N − ( m − ) τ w ( X m , τ t ) (4)where δ ( x , y ) is 1 if x = y and 0 otherwise. The final calculation of WPE is similar to Equation (1)above, but with the weighted probabilities: H w ( m , τ ) = − ∑ π ∈S m p w ( π ) log p w ( π ) . (5)Following standard practice, we normalize both PE and WPE results by dividing by log ( m ! ) , whichcauses them to range from 0 (low complexity) to 1 (high complexity) [12].The temporal resolution of a permutation entropy analysis is dictated by the span of the dataacross which the sum is taken. Since we are interested in the changing patterns in the δ D and δ Otraces described at the end of Section 2.1—not overall results that are aggregated across the entiretraces—the calculations in this paper employ Equations (1) and (5) in sliding windows across thosetraces. Choosing the size W for those sliding windows is not trivial. The minimum window width iseffectively dictated by the value of m because the number of data points that one needs in order togather representative statistics on the permutations grows with the length of those permutations. If,in expectation, one wants 100 chances for each of the m ! possible permutations to be discovered ina given window, there must be at least 100 m ! data points in the calculation. But even this minimumwindow size can be problematic in practice, leading to high variances in the probability estimatesfrom Equations (2) and (4). For our purposes, this is a real issue: if W is chosen too low, the highvariances can eradicate interesting features in the PE and WPE curves; too-large W values will reducethe resolution of the analysis, thereby washing out short-time-scale anomalies. To work around this, itis often helpful to increase W until the results stabilize, and that is the approach taken here. Since thiswill depend on the underlying statistics of the signal, careful PE and WPE calculations require some of 15 case-by-case hand-tuning and/or testing. (This is true of many other data-analysis methods, of course,though that is not widely appreciated in many scientific fields.)Unfortunately, the literature offers little rigorous advice regarding how to choose m ; the generalrecommendation is 3 ≤ m ≤
6, without any formal justification. The convergence proofs mentionedin the first paragraph of this section require that m → ∞ ; that would require an infinitely long timeseries (viz., W ≥ m !) and is obviously unreasonable for real-world data. In practice, this choiceis a balance between detail and data length: short permutations cannot capture the richness of thedynamics, but long permutations require long calculation windows. If m =
2, for instance, then theonly dynamics that are captured are “up” and “down.” In a time series with even moderate complexity,the probabilities of these two events are usually roughly similar, so both H ( τ ) and H w ( τ ) saturatenear 1 and neither PE nor WPE is informative. The richer variety of permutations offered by longerword lengths can more-accurately capture the complexity of the time series, but—as described in theprevious paragraph—that will drive up the window size and lower the temporal resolution of theanalysis. In the face of this, a useful practical strategy is to vary m and observe the effects on thefeatures in the PE and WPE curves. For real-world time series, those features often stabilize at very low m (hence the loose recommendation cited in the first sentence of this paragraph). In the case of the datastudied here, the features of the PE and WPE curves stabilized at m =
3. To be conservative, we used m = m and W values around these specific choices,confirming that the features in the PE and WPE curves remained the same. There is, of course, noformal guarantee that this tuning procedure will converge to good choices for this key free parameter,but similar “persistence” approaches are used in many different data-analysis approaches (e.g., delaysand dimensions for delay-coordinate embedding [28] or length scales for topological data analysis[29]).The third and final free parameter in PE and WPE calculations is the delay, τ . In the literature, itis customary to fix τ = original δ D and δ O measurements were spacedevenly in depth but unevenly in time. Because timeline deformation will obfuscate the results ofpermutation entropy calculations, one must transform the measured data to an even timeline, asdescribed in the last paragraph of Section 2.1, before invoking Equation (1) or (5). Linear interpolation,the standard practice in paleoclimate data analysis, can be problematic in this context, as the repeating,predictable patterns introduced by interpolation can skew the distribution of the permutations andthereby lower the PE and WPE values . This effect will generally worsen with depth because moreinterpolation is required lower in the core, where the temporal spacing between the measuredsamples is larger. As mentioned in Section 2.1, the uncertainty of the age scale also comes intoplay here—another effect that worsens with depth [23]. For all of these reasons, one cannot compareWPE values across wide ranges of the timeline of an ice core.Permutation entropy techniques have been used successfully to detect a variety of changes intime series: e.g., epileptic seizures in EEG signals [30], bifurcations in the transient logistic map [30],voiced sounds in a noisy speech signal [10], and market inefficiencies in financial records [31–33].They have also been used to assess predictability of time series in a variety of fields [12,13] and revealvarious interesting effects in paleoclimate data [1,7,34,35]. The bulk of this work has used the weighted variant of the technique. As we will show in this paper, though, PE and WPE work together in acomplementary fashion, detecting different kinds of anomalies in paleoclimate records. One could also downsample the original measurements to a fixed temporal spacing, but that would involve discarding largeamounts of the data in the top layers of the core, and thus would greatly reduce the resolution of the resulting WPE analysis. of 15
Figure 2.
Permutation entropy of the δ O and δ D data from Figure 1, calculated in rolling 2400-pointwindows (i.e., W = m = τ =
1. Top panel: weighted andunweighted permutation entropy of δ O in cyan and blue, respectively. Bottom panel: WPE and PE of δ D in orange and red, respectively.
3. Results
Figure 2 shows the weighted and unweighted permutation entropies of the δ O and δ D data fromFigure 1. There are a number of interesting features in these four curves, many of which have scientificimplications, as discussed in [1,7]. Here, we will focus on apparent anomalies in the complexity of theseclimate signals: i.e., discontinuities or abrupt changes in the curves. The most obvious such feature isthe large jump in all four traces between ≈ δ D and δ O depend less on their previous values during this period than in surrounding regions—i.e.,that a large amount of new information is produced at each time step by the system that generatedthe data. There are two potential culprits for such an increase: data issues or a pair of radical, rapidclimate shifts at either end of this time period. Since no such shifts are known or hypothesized by theclimate-science community, we conjectured that this jump was due to processing-related issues in thedata.Recall that PE and WPE are designed to bring out different aspects of the information in thedata. In view of this, the fact that both calculations pick up on this particular feature is meaningful: itsuggests that the underlying issue is not just noise—which would be at least partially washed out bythe weighting term, causing the jumps in WPE to be smaller than those in PE. Rather, there may beother causes at work here. Going back to the laboratory records and the original papers, we found thatan older, less precise, instrument was used to analyze this section of ice, and that this region of the corereceived the poorest possible quality score [18,36]. This part of the core is from the aforementionedbrittle zone, which shatters when removed from the ice sheet. This not only makes the timeline difficultto establish, as mentioned in Section 2.1, but can also affect the δ D and δ O values, as the drillingfluid can penetrate the ice and contaminate the measurements [18]. This could easily disturb the of 15
Figure 3.
Top panels: the re-measured δ D (red) and δ O (blue) data points in the range ≈ amplitude-encoded information in the core. Unfortunately, permutation entropy techniques cannottell us what the underlying causes of this anomaly are; that requires expert analysis and laboratoryrecords. Even so, their ability to flag problems, and help experts form scientific hypotheses about theircauses, is a major advantage of these information-theoretic techniques.As a case in point, we tested our hypotheses about the grey-shaded jump in Figure 2 byre-measuring the 331 m segment of the WAIS Divide ice core that corresponds to this time period.After obtaining the archived ice at the National Science Foundation Ice Core Facility (NSF-ICF), were-measured the isotope data with state-of-the-art equipment, repeating the depth-to-age conversionand temporal regularization processes described in Section 2. Plots of these new traces appear inFigure 3, with the old values shown in black. Visual examination of these traces makes two thingsapparent: lower levels of high-frequency noise in the new signals and small phase offsets betweenthe old and new data. The lower noise levels were, of course, the point of the resampling. The minorphase disparities are a consequence of updates in the laboratory pipeline and the difficulties posedby this particular section of the core. In the years since the creation of the first generation of thissystem—which was used to measure the original data—there have been a number of updates andimprovements. As a result, there are subtle differences (on the order of a few centimeters) in the depthregistration between the re-measured and the original data. This is a particular challenge in the brittlezone, since the broken or shattered ice pieces are difficult to piece back together and can settle alongfractures, reducing the length of the ice stick by a few centimeters. This can even cause completedata loss in short segments. Re-measuring this ice using state-of-the-art technology has allowed us toimprove the data in several important ways, by reducing the overall noise levels, improving the depthregistration, and even filling in missing parts of the original record that we obtained obtained uponremeasuring the archived ice.Given the kind of issues that were present in the data, such as increased small-scale variance, itis worth considering whether simpler approaches—less computationally complex than permutation of 15 Figure 4.
A simple anomaly detection algorithm. Here σ is estimated on a rolling 2400-pointoverlapping window for δ D (top) and δ O (bottom). Neither the feature in the grey-shaded boxin Figure 2, nor the other anomalies described later in this paper, are brought out by this calculation. entropy, and with fewer free parameters—would be equally effective in flagging the grey-shaded jumpin Figure 2. Because of the inherent data challenges (unknown processes at work on the data, lack ofreplicates, laboratory issues, etc.), this community has been traditionally limited to fairly rudimentaryapproaches to anomaly detection: e.g., discarding observations that lie beyond five standard deviationsfrom the mean. Accordingly, we performed a rolling-variance calculation on the same isotope recordsusing a 2400-point window. These results, shown in Figure 4, do not bring out the ≈ δ D and δ O records of Figure 1 with there-measured values from Figure 3 and repeat the PE and WPE calculations. The results are quitestriking, as shown in Figure 5: the large square waves are completely absent from the PE and WPEtraces of the repaired data set. This experiment—the longest high-resolution re-measurement andre-analysis of an ice core that has been performed to date—not only confirms our conjecture that thisanomaly was due to instrument error, but also allowed us to improve an important section of the datausing a targeted, highly focused effort.The repaired segment of water isotope data ( ≈ − underlyingmechanism for abrupt changes. However, their ability to identify a region of the core that shouldbe revisited—because of issues that would only have been apparent with a laborious, fine-grained Figure 5.
PE and WPE of δ O (top panel) and δ D (bottom panel) using re-measured data from1037-1368m ( ≈ − traditional analysis of the data—is a major advantage. This result highlights the main point that wewish to make in this paper: PE and WPE are useful methods for identifying anomalies in time-seriesdata from paleoclimate records. This is particularly useful in contexts where the data are difficultand/or expensive to collect and analyze. In these situations, a tightly focused new study, specificallytargeted on the offending area using permutation entropy, can maximize benefit with minimum effort.Another obvious feature in Figure 5 is the downward spike around 58 kybp in all four curves.Alerted by this anomaly in WPE and PE, we returned to the laboratory records and found that 1.107m (110.1 yr) of ice was unavailable from the record in this region and so a span of ≈ δ D and δ O traces were filled in by interpolation. This series of points—a linear ramp withpositive slope—translated to a long series of “1234” permutations. This causes a drop in the PE curvesas the calculation window passes across this expanse of completely predictable values. Indeed, forcalculations with τ = W = ≈ − δ D) or vice versa (e.g., 25.6, 30.1, and 47.5kybp for δ D and δ O). This brings out the differential nature of these two techniques, and the powerof the combined analysis: each can, independently, detect types of anomalies that the other misses.In the remainder of this section, we dissect a representative subset of the anomalies in Figure 5 toillustrate this claim, beginning with the case where WPE, but not PE, suggests data issues. Manualre-examination of the data in the regions around 25.6, 30.1, and 47.5 kybp revealed that the δ D and δ Osignal in these regions had been compromised in very small segments of the record ( < m of ice). Thetop left panel of Figure 6 shows the region around 47.5 kybp. The signal is visibly different in the smallshaded region, combining sharp corners and high-frequency oscillations—unlike the comparatively Figure 6.
Top row: δ D data from 47.5-47.8 kybp before and after removal of and interpolation over(black dashed line) a range of faulty values (shaded in gray). Bottom row: WPE calculated fromthe corresponding traces. The width of the square wave in the lower left plot is the size of the WPEcalculation window (2400 points at 1/20th year per point) plus the width of the anomaly. The horizontalshift between the earliest faulty value and the rise in WPE is due to the windowed nature of the WPEcalculation. smooth signal in the regions before and after this segment. WPE is much better at picking up this kindof anomaly because there has been a drastic shift in the information that is encoded in the amplitude ofthe signal. PE, in contrast, does not flag this segment because the distribution of permutations in thisregion is similar to that of original signal. However, the variance of each delay vector is quite differentthan in the surrounding signal: an effect to which PE is blind, but that WPE is designed to emphasize.There are many potential mechanisms that could be causing this distortion in the amplitude:extreme events, for instance, such as large volcanic eruptions, or the same kinds of data anddata-processing issues that are discussed above. Returning again to the laboratory records, wediscovered that the CRDS-CFA system had malfunctioned while analyzing this region. We could,as before, confirm that this was the cause of the anomaly by re-measuring this region of the core.If the outliers disappeared as a result, that would resolve some significant data issues with only aminimal, targeted amount of effort and expense. If they did not disappear, then the permutationentropy calculations would have identified a region where the record warranted further investigationby paleoclimate experts. To date, we have not carried out this experiment, but plan to request archivedice from the NSF-ICF for re-measurement and re-analysis. Access to this limited and irreplaceable
Figure 7.
Isotope data near 38.7 kybp that produced a spike in PE but not in WPE. According to the labrecords, the graphical user interface froze during the analysis of this segment of the ice, compromisingthe results. resource is, justifiably so, highly guarded—particularly in regard to the deeper ice. In lieu of thatexperiment, we re-processed the data from the region surrounding 47.5 − δ D record byremoving the offending data points and interpolating across the interval. As shown in the bottom-rightpanel of Figure 6, this removes the small square wave from the WPE trace. Without being able tore-measure this ice, of course, we cannot narrow down the cause of this anomaly. Even so, the WPEresults are useful in that they allow us to do some targeted reprocessing of the data in order to mitigatethe effects.Spikes in PE that are not associated with equally strong spikes in WPE are indications of effectsthat are dominated by small-scale noise. Figure 7 shows an example of the isotope record in one of theseregions. In contrast to the situation in Figure 6, which is dominated by high-frequency oscillations,this anomaly takes the form of a strong low-frequency component with a small band of superimposednoise. The weighting strategy in WPE causes it to ignore these features, but PE picks up on them veryclearly . To diagnose the cause of this anomaly, we again returned to the lab records, finding that thegraphical user interface froze while analyzing this particular ice stick, and that the data-processingroutines involved in distance registration and isotope content were compromised.Finally, there is a single spike (near 16.6 kybp in the WPE of δ D) where visual inspection of theraw data by expert paleoscientists does not suggest any issues, nor are there any recorded problems inthe laboratory records. This is not part of any sections of the core that are known to be problematic,like the brittle zone, and most likely warrants further investigation.
4. Discussion
This manuscript illustrates the potential of information theory as a forensic tool for paleoclimatedata, identifying regions with anomalous amounts of complexity so that they can be investigated Indeed, the frequency shift may magnify the impact of noise on that calculation.3 of 15 further. As a proof-of-concept demonstration of this claim, we used permutation entropy techniquesto isolate issues in water-isotope data from the WAIS Divide core, then re-measured and re-analyzedone of the associated regions using state-of-the-art laboratory equipment. The results of thisexperiment—the longest segment of replicate ice-core analysis in history—verified that the issueflagged by these information-theoretic techniques in that region were caused by outdated laboratorytechniques, as we had conjectured. We also showed that permutation entropy techniques identifiedseveral other smaller anomalies throughout the record, which we showed were caused by variousminor mishaps in the data processing pipeline. Finally, unlike other anomaly-detection studies thatuse either PE or WPE, we showed that using these two techniques together, rather than in isolation,allows detection of a much richer landscape of possible anomalies.This re-analysis approach is already providing valuable insights that reach well beyond thescope of this work. Compared to the prohibitive cost and time commitment of collecting a new core,re-measuring ice from an archived core can be comparatively easy and inexpensive if the regionrequiring re-analysis can be precisely and clearly defined. This forms the basis for what we believewill be an effective general targeted re-investigation methodology that can be applied not only to thewater-isotope record from the WAIS Divide ice core, but to other high-resolution paleoclimate records.Indeed, several such records are being finalized at the time of this publication, which makes this aperfect time to start developing, applying, and evaluating sophisticated anomaly-detection methodsfor these important data sets. Excitingly, as a result of the study reported in this manuscript, PE-basedtechniques are now being added to the quality control phase of the analysis pipelin for at least one ofthese new high-resolution records.The information in paleoclimate records contains valuable clues and insights about the Earth’s pastclimate, and perhaps about its future. While these records have a lot of promise, it is crucial—especiallyin an era of rampant climate-change denial—to be careful in extracting the useful and meaningfulinformation from these records while simultaneously identifying regions that are problematic. Fora multitude of reasons, distinguishing useful information from a lack thereof can be a particularlychallenging task in paleodata. This is especially true in parts of these records that are hard to analyze,such as the brittle zone of the WDC. Our work has identified a number of intervals in this core wherethe data require a closer look, including several that contain information corresponding to the timeperiod of the dawn of human civilization.
Author Contributions: conceptualization, J.G.; methodology, J.G., T.J., E.B., M.N. and J.W.; software, J.G. andM.N..; validation, J.G., T.J., E.B., M.N. and J.W.; formal analysis, J.G., E.B. and T.J.; investigation, J.G., T.J, M.N.,E.B.and V.M.; resources, T.J, J.W. and V.M.; data curation, V.M., J.W. and T.J.;writing—original draft preparation J.G.and E.B.; writing—review and editing, J.G., E.B. and T.J.; visualization, J.G. and M.N.; supervision, E.B. and J.W..;project administration, J.G. and E.B.; funding acquisition, J.G. and T.J.
Funding:
This research was funded by US National Science Foundation (NSF) grant number 1807478. JG was alsosupported by an Omidyar Fellowship at the Santa Fe Institute. For original ice core data, this work was supportedby NSF grants 0537593, 0537661, 0537930, 0539232, 1043092, 1043167, 1043518 and 1142166. Field and logisticalactivities were managed by the WAIS Divide Science Coordination Office at the Desert Research Institute, USA,and the University of New Hampshire, USA (NSF grants 0230396, 0440817, 0944266 and 0944348). The NSFDivision of Polar Programs funded the Ice Drilling Program Office (IDPO), the Ice Drilling Design and Operations(IDDO) group, the National Science Foundation Ice Core Facility (NSF-ICF), the Antarctic Support Contractor,and the 109th New York Air National Guard.
Acknowledgments:
Water isotope measurements were performed at the Stable Isotope Lab (SIL) at the Instituteof Arctic and Alpine Research (INSTAAR), University of Colorado. We wish to thank Bruce Vaughn for his effortsin designing and implementing the laser-based, continuous flow measurement system for water isotopes, as wellas Mark Twickler, Geoff Hargreaves, and Richard Nunn for their assistance planning and executing ice samplingat NSF-ICF. We thank Bedartha Goswami for his advice. We would also like to thank the reviewers of this paperfor their valuable feedback—particularly Peter Ditlevsen.
Conflicts of Interest:
The authors declare no conflict of interest.
Abbreviations
The following abbreviations are used in this manuscript:CFA Continuous Flow AnalysisCRDS-CFA Cavity Ring-Down Spectroscopy Continuous Flow AnalysisNSF-ICF National Science Foundation Ice Core FacilityPE Permutation EntropyWPE Weighted Permutation EntropyWAIS West Antarctic Ice SheetWDC WAIS Divide Ice Core
References
1. Garland, J.; Jones, T.R.; Bradley, E.; Neuder, M.; White, J.W. Climate entropy production recorded in a deepAntarctic ice core. arXiv preprint arXiv:1806.10936 .2. Jones, T.R.; White, J.W.C.; Steig, E.J.; Vaughn, B.H.; Morris, V.; Gkinis, V.; Markle, B.R.; Schoenemann, S.W.Improved methodologies for continuous-flow analysis of stable water isotopes in ice cores.
AtmosphericMeasurement Techniques , , 617–632.3. Casey, K.; Fudge, T.; Neumann, T.; Steig, E.; Cavitte, M.; Blankenship, D. The 1500 m South Pole ice core:Recovering a 40 ka environmental record. Annals of Glaciology , , 137–146.4. eastgrip.org.5. Münch, T.; Kipfstuhl, S.; Freitag, J.; Meyer, H.; Laepple, T. Regional climate signal vs. local noise: Atwo-dimensional view of water isotopes in Antarctic firn at Kohnen Station, Dronning Maud Land. Climateof the Past , , 1565–1581.6. Jones, T.; White, J.; Popp, T. Siple Dome shallow ice cores: A study in coastal dome microclimatology. Climate of the Past , , 1253–1267.7. Garland, J.; Jones, T.; Bradley, E.; James, R.; White, J.W.C. A first step toward quantifying the climate’sinformation production over the last 68,000 years. Proceedings of the 12th International Symposium onIntelligent Data Analysis; Springer: Stockholm, Sweeden, 2016; Lecture Notes in Computer Science.8. Jones, T.R.; Roberts, W.H.G.; Steig, E.J.; Cuffey, K.M.; Markle, B.R.; White, J.W.C. Southern Hemisphereclimate variability forced by Northern Hemisphere ice-sheet topography. Nature , , 351–355.9. Johnsen, S.J.; Dahl-Jensen, D.; Gundestrup, N.; Steffensen, J.P.; Clausen, H.B.; Miller, H.; Masson-Delmotte,V.; Sveinbjörnsdottir, A.E.; White, J. Oxygen isotope and palaeotemperature records from six Greenlandice-core stations: Camp Century, Dye-3, GRIP, GISP2, Renland and NorthGRIP. Journal of QuaternaryScience: Published for the Quaternary Research Association , , 299–307.10. Bandt, C.; Pompe, B. Permutation entropy: A natural complexity measure for time series. Physical ReviewLetters , , 174102.11. Fadlallah, B.; Chen, B.; Keil, A.; Príncipe, J. Weighted-permutation entropy: A complexity measure fortime series incorporating amplitude information. Physical Review E , , 022911.12. Garland, J.; James, R.G.; Bradley, E. Model-free quantification of time-series predictability. Physical ReviewE , , 052910.13. Pennekamp, F.; Iles, A.; Garland, J.; Brennan, G.; Brose, U.; Gaedke, U.; Jacob, U.; Kratina, P.; Matthews,B.; Munch, S.; Novak, M.; Palamara, G.M.; Rall, B.; Rosenbaum, B.; Tabi, A.; Ward, C.; Williams, R.; Ye,H.; Petchey, O. The intrinsic predictability of ecological time series and its potential to guide forecasting. Preprint available at bioRxiv 350017 .14. Keller, K.; Unakafov, A.M.; Unakafova, V.A. On the relation of KS entropy and permutation entropy.
Physica D: Nonlinear Phenomena , , 1477 – 1481.15. Bandt, C.; Keller, G.; Pompe, B. Entropy of interval maps via permutations. Nonlinearity , , 1595.16. Epstein, S.; Buchsbaum, R.; Lowenstam, H.A.; Urey, H.C. Revised carbonate-water isotopic temperaturescale. Geological Society of America Bulletin , , 1315–1326.17. Mook, W.; Rozanski, K. Environmental isotopes in the hydrological cycle. IAEA Publish , .18. Sigl, M.; Fudge, T.J.; Winstrup, M.; Cole-Dai, J.; Ferris, D.; McConnell, J.R.; Taylor, K.C.; Welten, K.C.;Woodruff, T.E.; Adolphi, F.; Bisiaux, M.; Brook, E.J.; Buizert, C.; Caffee, M.W.; Dunbar, N.W.; Edwards, R.;Geng, L.; Iverson, N.; Koffman, B.; Layman, L.; Maselli, O.J.; McGwire, K.; Muscheler, R.; Nishiizumi, K.; Pasteris, D.R.; Rhodes, R.H.; Sowers, T.A. The WAIS Divide deep ice core WD2014 chronology – Part 2:Annual-layer counting (0–31 ka BP).
Climate of the Past , , 769–786.19. Buizert, C. The WAIS Divide deep ice core WD2014 chronology—Part 1: Methane synchronization (68–31ka BP) and the gas age–ice age difference. Climate of the Past , , 153–173.20. Jones, T.R.; Cuffey, K.M.; White, J.W.C.; Steig, E.J.; Buizert, C.; Markle, B.R.; McConnell, J.R.; Sigl, M. Waterisotope diffusion in the WAIS Divide ice core during the Holocene and last glacial. Journal of GeophysicalResearch: Earth Surface , , 290–309. 2016JF003938, doi:10.1002/2016JF003938.21. Sakellariou, K.; McCullough, M.; Stemler, T.; Small, M. Counting forbidden patterns in irregularly sampledtime series. II. Reliability in the presence of highly irregular sampling. Chaos: An Interdisciplinary Journal ofNonlinear Science , , 123104.22. McCullough, M.; Sakellariou, K.; Stemler, T.; Small, M. Counting forbidden patterns in irregularly sampledtime series. I. The effects of under-sampling, random depletion, and timing jitter. Chaos: An InterdisciplinaryJournal of Nonlinear Science , , 123103.23. Goswami, B.; Heitzig, J.; Rehfeld, K.; Marwan, N.; Anoop, A.; Prasad, S.; Kurths, J. Estimation ofsedimentary proxy records together with associated uncertainty. Nonlinear Processes in Geophysics , , 1093–1111.24. Boers, N.; Goswami, B.; Ghil, M. A complete representation of uncertainties in layer-counted paleoclimaticarchives. Climate of the Past , , 1169–1180.25. Eroglu, D.; McRobie, F.H.; Ozken, I.; Stemler, T.; Wyrwoll, K.H.; Breitenbach, S.F.; Marwan, N.; Kurths, J.See–saw relationship of the Holocene East Asian–Australian summer monsoon. Nature communications , , 12929.26. Amigó, J.M.; Kennel, M.B.; Kocarev, L. The permutation entropy rate equals the metric entropy rate forergodic information sources and ergodic dynamical systems. Physica D: Nonlinear Phenomena , , 77 –95. doi:https://doi.org/10.1016/j.physd.2005.07.006.27. Amigó, J.M. The equality of Kolmogorov–Sinai entropy and metric permutation entropy generalized. Physica D: Nonlinear Phenomena , , 789 – 793.28. Bradley, E.; Kantz, H. Nonlinear time-series analysis revisited. Chaos: An Interdisciplinary Journal ofNonlinear Science , , 097610.29. Robins, V. Towards computing homology from finite approximations. Topology Proceedings , , 503–532.30. Cao, Y.; Tung, W.; Gao, J.; Protopopescu, V.; Hively, L. Detecting dynamical changes in time series usingthe permutation entropy. Physical Review E , , 046217.31. Zunino, L.; Zanin, M.; Tabak, B.M.; Pérez, D.G.; Rosso, O.A. Forbidden patterns, permutation entropy andstock market inefficiency. Physica A: Statistical Mechanics and its Applications , , 2854 – 2864.32. Zunino, L.; Zanin, M.; Tabak, B.M.; Pérez, D.G.; Rosso, O.A. Complexity-entropy causality plane: A usefulapproach to quantify the stock market inefficiency. Physica A: Statistical Mechanics and its Applications , , 1891 – 1901.33. Zanin, M. Forbidden patterns in financial time series. Chaos: An Interdisciplinary Journal of Nonlinear Science , , 013119.34. Saco, P.M.; Carpi, L.C.; Figliola, A.; Serrano, E.; Rosso, O.A. Entropy analysis of the dynamics of ElNiño/Southern Oscillation during the Holocene. Physica A: Statistical Mechanics and its Applications , , 5022–5027.35. Balasis, G.; Donner, R.V.; Potirakis, S.M.; Runge, J.; Papadimitriou, C.; Daglis, I.A.; Eftaxias, K.; Kurths, J.Statistical mechanics and information-theoretic perspectives on complexity in the Earth system. Entropy , , 4844–4888.36. Souney, J.M.; Twickler, M.S.; Hargreaves, G.M.; Bencivengo, B.M.; Kippenhan, M.J.; Johnson, J.A.; Cravens,E.D.; Neff, P.D.; Nunn, R.M.; Orsi, A.J.; others. Core handling and processing for the WAIS Divide ice-coreproject. Annals of Glaciology ,55