Characterisation of Dansgaard–Oeschger events in palaeoclimate time series using the matrix profile method

Palaeoclimate time series, reflecting the state of Earth's climate in the distant past, occasionally display very large and rapid shifts showing abrupt climate variability. The identification and characterisation of these abrupt transitions in palaeoclimate records is of particular interest as this allows for understanding of millennial climate variability and the identification of potential tipping points in the context of current climate change. Methods that are able to characterise these events in an objective and automatic way, in a single time series, or across two proxy records are therefore of particular interest. In our study the matrix profile approach is used to describe Dansgaard–Oeschger (DO) events, abrupt warmings detected in the Greenland ice core, and Northern Hemisphere marine and continental records. The results indicate that canonical events DO-19 and DO-20, occurring at around 72 and 76 ka, are the most similar events over the past 110 000 years. These transitions are characterised by matching transitions corresponding to events DO-1, DO-8, and DO-12. They are abrupt, resulting in a rapid shift to warmer conditions, followed by a gradual return to cold conditions. The joint analysis of the δ¹⁸O and Ca²⁺ time series indicates that the transition corresponding to the DO-19 event is the most similar event across the two time series.

1 Introduction

Palaeoclimate time series reflect Earth's climate in the distant past based mainly on proxies derived from records such as sediments, ice cores, or speleothems. One of the most extensively studied palaeoclimate time series is that of oxygen isotope ratios (δ¹⁸O) retrieved from Greenland ice cores in the context of the North Greenland Ice Core Project (NGRIP) (Andersen et al., 2004; Rasmussen et al., 2014). The concentration of δ¹⁸O, as measured in the ice, serves as an indirect proxy of the air temperature at the ice core location.

The ice cores retrieved from the Greenland ice sheet revealed the occurrence of rapid warming events, which occurred over a few decades. Such abrupt transitions are designated Dansgaard–Oeschger (DO) events, during which climate conditions alternated between fully glacial (so-called stadial) and relatively mild (interstadial) conditions (Dansgaard et al., 1982; Johnsen et al., 1992; Rasmussen et al., 2014). These abrupt transitions, on average, are approximately 12 °C, with a range of 5 to 16 °C (Kindler et al., 2014). They exhibit a distinctive sawtooth shape, with rapid, decadal, increases in δ¹⁸O from Greenland stadial (GS) to Greenland interstadial (GI) conditions, followed by slow relaxations back to GS conditions on timescales of centuries or millennia. However, it should be noted that not all DO events have the same shape or duration (Rousseau et al., 2017; Lohmann, 2019).

DO events are particularly observable in the Greenland ice core records, but similar transitions were identified in diverse palaeoclimate records (e.g. Rousseau et al., 2017; Corrick et al., 2020; Bagniewski et al., 2023; Held et al., 2024), and thus DO events have been used for “wiggle matching” of records with no accurate dating information (e.g. Peterson et al., 2000; Lauterbach et al., 2011). The objective identification of DO events, other than by visual inspection, is then of critical importance.

The original identification of DO events was conducted by visual inspection of the NGRIP δ¹⁸O time series. These “canonical DOs” (Rousseau et al., 2022) and further transitions were also identified visually by Rasmussen et al. (2014), who resorted to not only the δ¹⁸O record, but also the concurrent Ca²⁺ proxy record. An algorithm was developed by Lohmann and Ditlevsen (2019) for the characterisation of DO events based on their sawtooth shape. However, the approach does not identify the events themselves. Rather, it relies on the previous visual identification by Rasmussen et al. (2014). A method based on the Kolmogorov–Smirnov (KS) test was developed by Bagniewski et al. (2021), which allows the identification of the canonical DO events in addition to other transitions that were previously unidentified by visual inspection. Although the method allows for the detection of individual jumps in the records towards either cold–warm or dry–wet conditions, it does not provide any information on their magnitude. Recurrence plots (e.g. Marwan et al., 2007) and measures of recurrence quantification analysis such as the recurrence rate (Marwan et al., 2013) allow on the other hand for the identification of the dominant changes in a record's characteristic timescale (Bagniewski et al., 2021; Rousseau et al., 2023a, b).

Given that DO events can be considered a recurring pattern in a palaeoclimate time series, algorithmic methods for the extraction of similar patterns from a time series can be applied to characterise these particular patterns. In this study the matrix profile approach (Yeh et al., 2016) is employed to describe DO patterns in the Greenland ice core records. The methodology and data are described in Sects. 2 and 3, respectively. The results of the matrix profile approach applied to palaeoclimate time series are presented in Sect. 4. Methodological and interpretative constraints are discussed in Sect. 5, and concluding remarks are given in Sect. 6.

2 Method

The matrix profile approach is described in detail in Yeh et al. (2016, 2018) and the references therein. This overview provides only a brief summary of the methodology. Readers are referred to the original references for further details.

For a real-value time series $\mathit{T}=\mathit{\{}{t}_{\mathrm{1}},\mathrm{\dots },t_n\mathit{\}}$ of length n, a sub-sequence T_i,m of length m is a continuous subset of the values of T: ${\mathit{T}}_{i,m}=t_i,t_{i+\mathrm{1}},\mathrm{\dots },t_{i+m-\mathrm{1}},\mathrm{1}\le i\le n-m+\mathrm{1}$. The matrix profile is an ordered vector of the Euclidean distances between the most similar (shortest distance) sub-sequences, considering all possible sub-sequences of T obtained by sliding a window of length m across T. The distance is measured using the Euclidean distance between z-normalised sub-sequences with a mean of 0 and a standard deviation of 1 (Agrawal et al., 1993). The matrix profile is defined as the distance between a sub-sequence and its most similar sub-sequence, regardless of its location. That location is stored in the profile index. The profile index is a companion vector that stores the position of the nearest neighbour of each sub-sequence. The matrix profile and the matrix profile index are two data structures that annotate a time series with the distance and location, respectively, of the nearest neighbours of all its sub-sequences (in itself or in another time series).

The brute-force obvious algorithm to compute the matrix profile requires the computation of a large number of Euclidean distances equal to ${\sum }_{i=\mathrm{1}}^{n-m+\mathrm{1}}i$. While this complexity 𝒪((n−m)²) is not a problem in the case of the typical short lengths of palaeoclimate time series, it can easily become prohibitive for even modestly sized time series. Scalable and fast algorithms have been developed to tackle this issue, reducing the complexity of matrix profile calculations to 𝒪(n) (Yeh et al., 2018). The scalability and computational efficiency of the matrix profile are not what makes it an attractive tool in itself, but the variety of useful analytical tasks that can be performed based on the matrix profile is (Zhu et al., 2020).

The matrix profile enables the identification of motifs in the time series, corresponding to sub-sequences that are highly similar to each other. The (tying) lowest points in the matrix profile correspond to the locations of the optimal time series motif pair, i.e. the pair of sub-sequences that are most similar (Mueen et al., 2009): $\forall a,b,i,j\in [\mathrm{1},\mathrm{2},\mathrm{\dots },n-m+\mathrm{1}]$ and dist$({\mathit{T}}_{a,m},{\mathit{T}}_{b,m})\le$ dist$({\mathit{T}}_{i,m},{\mathit{T}}_{j,m})$, with $i\ne j,a\ne b,|i-j|\ge w,|a-b|\ge w,w>\mathrm{0}$, and dist the z-normalised Euclidean distance between the sub-sequences. The parameter w imposes a constraint on the relative positions of the sub-sequences in the motif, ensuring a gap between the sub-sequences of at least w values. This exclusion zone is intended to exclude trivial motifs as described in Lin et al. (2002). The exclusion zone thus defined enables the avoidance of trivial matches of a sub-sequence with itself or of largely overlapping sub-sequences. For a top motif of the most similar sub-sequences, which are separated by D_1,min, neighbouring similar sub-sequences can be obtained by considering the sub-sequences within $D_{\mathrm{1},\min }\times R$ distances, where R is a small integer value (typically ≤2). The definition of a first-order motif can be generalised to that of a kth-order subsequent motif (where k is an integer number) by considering the following nearest pair of sub-sequences after excluding all sub-sequences belonging to prior motifs.

For a single time series T, the matrix profile stores the Euclidean distance from a sub-sequence T_i,m to its most similar sub-sequence. The concept can be extended to two time series T_A and T_B (which may have different lengths), with the join matrix profile containing the Euclidean distances between all sub-sequences in T_A and their most similar sub-sequence in T_B. The join motifs reflect the most similar patterns between the two time series (Yeh et al., 2017, 2018).

The Python code stumpy (Law, 2019) is employed to perform the matrix profile analysis. The R statistical software (version 4.1.3, R Core Team, 2022) is employed for the visualisation of the matrix profile results.

3 Data

The analysis presented here applies to the 20-year-resolution time series of δ¹⁸O and calcium ion (Ca²⁺) concentrations from the NGRIP ice core on the GICC05modelext timescale (Rasmussen et al., 2014; Seierstad et al., 2014). The time series encompasses the interval from 107.6 to 10.2 ka b2k (before AD 2000), with 4869 data points in total. While the oxygen isotope ratio is a proxy for air temperature, calcium is a proxy of atmospheric dust (Fuhrer et al., 1999), reflecting changes in dust sources and transport pathways (Fischer et al., 2007). The time series are presented in Fig. 1, with the Ca²⁺ record represented on a reverse logarithmic scale, in line with the approach outlined by Rasmussen et al. (2014). A small number of missing values in the Ca²⁺ time series, <1%, were linearly interpolated.

Figure 1 Time series of δ¹⁸O and Ca²⁺ concentrations on the GICC05modelext timescale. Note that Ca²⁺ is on the inverse logarithmic scale. The vertical solid and dashed lines represent, respectively, the start and end times of canonical DO events (Dansgaard et al., 1993; Rasmussen et al., 2014).

4 Results

The results of the matrix profile analysis of the palaeoclimate time series are presented initially in terms of the characterisation of DO events from a single time series, the δ¹⁸O Greenland record, in Sect. 4.1, and the Ca²⁺ record, in Sect. 4.2. Subsequently, the joint analysis of DO events from multiple time series is documented in Sect. 4.3 using both the δ¹⁸O and Ca²⁺ proxy records.

4.1 Matrix profile analysis of the δ¹⁸O record

The matrix profile of the δ¹⁸O time series is computed using the stomp (Scalable Time series Ordered-search Matrix Profile) algorithm (Zhu et al., 2016) with a sub-sequence length (window size m) of 2500 years (125 data points). The exclusion zone $w/m$, defined as a proportion of the sub-sequence length, is set to 1, which is equal to the size of the sub-sequence length, in order to exclude trivial matches (similarity of a sub-sequence to another sub-sequence with data values in common). For each ith sub-sequence consisting of the time series values from i to $i+m-\mathrm{1}$ with $i=\mathrm{1},\mathrm{2},\mathrm{\dots }n-m+\mathrm{1}$, its z-normalised Euclidean distance to every other jth sub-sequence is computed. The matrix profile is defined as the smallest value from that set of distances, while the profile index records the location of the most similar (non-trivial) sub-sequence, or in other words the value of j. Table 1 provides an illustrative example of the matrix profile and the corresponding profile index computed for the δ¹⁸O time series. The first column contains the data point corresponding to the start of the sub-sequence, i.e. the 3044th value of the time series (and the corresponding time) and so forth. The middle column contains the matrix profile, i.e. the minimum value of all the distances between sub-sequence 3044 and all other sub-sequences in the time series. The third column gives the profile index, which for the first row indicates that the most similar sub-sequence to sub-sequence 3044 is the sub-sequence starting at the 3250th value of the time series – all other sub-sequences in the time series are separated from sub-sequence 3044 by a value higher than 2.37.

Table 1 Snippet of the matrix profile and profile index of the δ¹⁸O time series. Values corresponding to the global minimum of the matrix profile are represented in bold.

Figure 2 depicts the complete matrix profile for the δ¹⁸O time series, which represents the distance of each sub-sequence to its most similar (non-trivial) sub-sequence (same as the middle column in Table 1). The global minimum value of the matrix profile is 2.34, occurring at 71.14 ka (indicated by the vertical dashed line in Fig. 2 and the bold numbers in Table 1). Of all the sub-sequences in the δ¹⁸O time series, the one starting at 71.14 ka exhibits the shortest distance to its most similar sub-sequence, i.e. the one starting at 75.26 ka (Table 1). It should be noted that this sub-sequence does not necessarily correspond to the next minimum value of the matrix profile. Rather, it is the sub-sequence that is most similar to the one extracted from the global minimum and not necessarily the sub-sequence with the second-lowest matrix profile distance.

Figure 2 Matrix profile of the NGRIP δ¹⁸O time series using a window size of 2500 years (125 data points). The dashed vertical line indicates the minimum value of the matrix profile corresponding to the top motif: same age model as in Fig. 1.

The top motif pair, representing the two most similar sub-sequences in the time series starting at 71.14 and 75.26 ka, respectively, is displayed in Fig. 3. This motif pair corresponds to the canonical DO events DO-19 and DO-20, which are most similar to each other. While the sub-sequence starting at 75.26 ka is most similar to the top sub-sequence corresponding to the matrix profile minimum, at 71.14 ka, other sub-sequences, while not the most similar ones, are still not far off. All the sub-sequences for which the distance to the top sub-sequence at 71.14 ka is within R times the distance from the top motif pair can be considered neighbouring motifs. Table 2 presents the motifs obtained for the δ¹⁸O time series using values of R equal to 2 and 3. Setting a value of R=2 enables the extraction of three neighbour motifs to the top motif, while setting R=3 allows for the additional extraction of five further neighbours. It should be noted that the distance represented in Table 2 is the distance of the motif sub-sequence to the top motif and not its matrix profile distance. An alternative approach to the use of a fixed value of the parameter R is to set a maximum distance between a neighbour motif and the main motif based on the variability of the matrix profile. In this study, neighbouring motifs are considered to be separated by less than twice the global standard deviation of the matrix profile (equal to $\mathrm{2}\times \mathrm{2.38}=\mathrm{4.76}$), which in this case corresponds to the neighbours that would be extracted by considering R=2 (Fig. 3). This aspect of selecting neighbour sequences to the main motif is further discussed in Sect. 5. The last column of Table 2 indicates the canonical DO event that coincides with each motif. The seventh motif (12.98–10.48 ka) does not correspond to a DO event, but it does include the Younger Dryas cooling event (GS-1) and the transition to the Holocene.

Figure 4 depicts the normalised sub-sequences corresponding to the motifs extracted from the δ¹⁸O time series. The top motif represents the patterns corresponding to the canonical DO-19 and DO-20. These aforementioned sub-sequences and their neighbour motifs represent abrupt transitions to warmer conditions, preceded by an approximately stable stadial level and followed by a slow decrease towards stadial conditions.

Figure 3 Top motifs for the NGRIP δ¹⁸O time series and neighbour motifs. The horizontal top numbers indicate the canonical DO events: same age model as in Fig. 1.

Figure 4 Normalised motifs for the NGRIP δ¹⁸O time series. The index is the order of the δ¹⁸O values in each sub-sequence of 125 values (2500 years). The top motif pair corresponds to DO-19 and DO-20, and the neighbour motif corresponds to DO-12, DO-8, and DO-1.

Table 2 Motifs extracted from the matrix profile of the δ¹⁸O time series using a window size of 2500 years. The columns display the value of the radius parameter R for neighbouring sub-sequences, the type of motif, the end and start times of the motif, the distance to the top motif, and the corresponding canonical DO event.

4.2 Matrix profile analysis of the Ca²⁺ record

A comparable analysis is conducted for the Ca²⁺ record employing the same methodology as that employed for the δ¹⁸O time series. Figure 5 depicts the matrix profile for the Ca²⁺ record. While the overall pattern of the Ca²⁺ series is similar to that of the matrix profile of the δ¹⁸O record in Fig. 2, the peaks are sharper due to the lower noise level of the Ca²⁺ series. Moreover, the minimum value of the matrix profile for the Ca²⁺ series, indicated by the vertical dashed line in Fig. 5, occurs at a different time (at 37.14 ka) than for δ¹⁸O. Consequently, the sub-sequence with the shortest distance to its most similar sub-sequence in the Ca²⁺ record starts at 37.14 ka, corresponding to DO-8. From the profile index, the most similar sub-sequence, which constitutes the top motif pair for the Ca²⁺ time series, starts at 45.78 ka, corresponding to DO-12. The aforementioned top motif and its neighbour motifs are summarised in Table 3. The criterion of selecting as neighbour motifs the sub-sequences separated from the main motif by less than twice the matrix profile's standard deviation ($\mathrm{2}\times \mathrm{2.50}=\mathrm{5}$) yields all neighbour motifs for a radius of R=2 and the first neighbour motif for a value of R=3. This is illustrated in Fig. 6.

Figure 5 Matrix profile of the NGRIP Ca²⁺ record: same conventions as in Fig. 2.

Figure 6 Top motifs for the NGRIP Ca²⁺ record: same conventions as in Fig. 3.

Table 3 Motifs extracted from the matrix profile of the Ca²⁺ record: same conventions as in Table 2.

Figure 7 displays the normalised motifs for the Ca²⁺ time series. These motifs correspond to the same canonical DOs as the motifs of the δ¹⁸O time series, with the addition of DO-21. However, the ordering of the motifs differs. In the case of the Ca²⁺ record, considered on the reverse logarithmic scale, the recurring patterns, constituted by the top motif and its neighbour motifs, represent abrupt decreases in the atmospheric dust concentration. These decreases are preceded by slightly decreasing stadial normalised values and are followed by an approximately stable interstadial level.

Figure 7 Normalised motifs for the NGRIP Ca²⁺ record. The conventions are the same as in Fig. 4, but the top motif pair corresponds to DO-8 and DO-12 and its neighbour to DO-19, DO-1, DO-21, and DO-20. The original data are on the reverse logarithmic scale.

4.3 Join matrix profile analysis of the δ¹⁸ and Ca²⁺ records

Figure 8 depicts the join matrix profile of the δ¹⁸O and Ca²⁺ time series, which is defined as the distance between each sub-sequence in the δ¹⁸O record and its most similar sub-sequence in the Ca²⁺ record. In this join case, the necessity for an exclusion zone to avoid trivial matches is negated by the fact that the sub-sequences originate from different time series. The join matrix profile (Fig. 8) differs from the previously defined individual profiles but exhibits similarities to the matrix profile for the single δ¹⁸O time series (Fig. 2). This is not unexpected given that the temporal variability of the Ca²⁺ series closely follows that of the δ¹⁸O series. For illustrative purposes, a snippet of the join matrix profile and the profile index is presented in Table 4.

Figure 8 Join matrix profile for the δ¹⁸O and Ca²⁺ time series using a window size of 2500 years (125 data points). The dashed vertical line indicates the minimum value of the matrix profile corresponding to the top motif (sub-sequence of the δ¹⁸O time series most similar to a sub-sequence in the Ca²⁺ time series). Same age model than in Fig. 1.

Table 4 Snippet of the join matrix profile and the profile index of the δ¹⁸O and Ca²⁺ time series. Values corresponding to the global minimum of the matrix profile are represented in bold.

The minimum value of the join matrix profile indicates the location of the top motif, i.e. the sub-sequence in the δ¹⁸O time series which is most similar (distance-wise) to a sub-sequence in the Ca²⁺ time series (whatever it is). The global minimum of the join matrix profile occurs at 71.28 ka, which is in close proximity to the previously defined minimum for δ¹⁸O at 71.14 ka. Table 4 indicates that the sub-sequence in the δ¹⁸O time series starting at 71.28 ka exhibits the greatest similarity to the sub-sequence in the Ca²⁺ record also starting at 71.28 ka. The most analogous sub-sequences across the δ¹⁸O and Ca²⁺ records are presented in Fig. 9, still with Ca²⁺ displayed on the reverse logarithmic scale. In this particular case, the location (time) of the top motif is the same for the two records, and thus normalised motifs are represented on the temporal scale rather than plotted as a function of the index value from m=1 to m=125, as in Figs. 4 and 7. This motif is identical to the previously identified motif in the individual analysis of the δ¹⁸O time series (Sect. 4.1). This transition, occurring approximately 70 000 years ago (the canonical DO-19), exhibits the most analogous patterns of warming and cooling across δ¹⁸O and less dusty or dustier Ca²⁺ records. Additionally, it bears the closest resemblance in shape to the transition in the δ¹⁸O record immediately after, approximately 73 000 years ago.

Figure 9 Normalised top motif across the NGRIP δ¹⁸O and Ca²⁺ time series: same age model as in Fig. 1.

5 Discussion

While the matrix profile is an algorithmic approach that enables the extraction of recurring patterns in a time series, its results are dependent on parameters that are prescribed empirically, often through a trial-and-error process, and they are dependent on the specific application and purpose of the analysis. These constraints are discussed in Sect. 5.1, while Sect. 5.2 discusses the extraction of the most similar pattern across two time series. In the case of short time series, such as the ones analysed here, the utility of the matrix profile approach is more obvious when applied to two time series. However, finding motifs in a single long and high-frequency time series is a very common need in data-mining contexts, for which the matrix profile of a single time series is considered to be an extremely useful tool (Zhu et al., 2020).

5.1 Matrix profile parameters

The matrix profile is dependent on a single parameter, the sub-sequence length m. Consequently, it is a method that can be readily applied to any time series, as it only requires the appropriate specification of a window size (see Sect. 4.1). However, the matrix profile results are highly contingent upon the value selected for m. There are no universal formulas or rigorous criteria for selecting the window size, as this depends on the objective of the analysis and the type of motif being investigated. The sub-sequence length is typically determined by considering the length of the patterns of interest in the time series being analysed. In the absence of prior knowledge regarding the length of the motif of interest in a time series, an extension of the matrix profile, which computes nearest-neighbour information for all sub-sequences of all lengths, can be considered (Madrid et al., 2019).

In our study, a window size of 2500 years was selected as adequate for the extraction of motifs with durations typical of DO events. The matrix profile obtained with this window size is compared in Fig. 10 with the matrix profile for window sizes of 3000 and 3500 years. The matrix profile values increase with increasing values of m, yet the results remain remarkably consistent across this range of window sizes, exhibiting a similar overall pattern and the lowest distance at around 70 ka. The highest values of the matrix profile (at around 90, 66, and 26 ka) reflect the flatter parts of the δ¹⁸O time series, which lack discernible features and thus comprise more disparate sub-sequences. The motifs obtained for the different window sizes are presented in Fig. 11, which provides further evidence of the robustness of the recurring patterns extracted in the considered window range.

Figure 10 Matrix profile of the NGRIP δ¹⁸O time series for different window sizes. The dashed vertical line indicates the minimum value of the matrix profile: same age model as in Fig. 1.

Figure 11 Motifs extracted from the matrix profile of the NGRIP δ¹⁸O time series for different window sizes: same age model as in Fig. 1.

The top motif is derived directly from the lowest values of the matrix profile and therefore depends only on the specified sub-sequence length. However, the extraction of neighbouring sequences to the top motif and of other motifs depends not only on the window size, but also on the tolerance (distance-wise) with which a sub-sequence is considered to match a pattern given by the value of parameter R. Once more, there are no specific criteria for specifying R, other than it should be a small integer value. The larger the value of R, the greater the distance between sub-sequences considered to match the motif and therefore the lower the similarity between potential matches. The use of smaller values of R reduces the potential for considering as matching sub-sequences that are dissimilar. This is demonstrated in Fig. 12, which depicts the normalised motifs for the NGRIP δ¹⁸O time series using R=2 and the supplementary motifs obtained using R=3. To facilitate visualisation, the first two neighbour motifs for R=3 are plotted separately from the subsequent neighbour motifs (see Table 2). The discrepancy between the top motif pattern and the neighbour motifs assigned increases with the higher-order motifs indicating a value of R=3 is not optimal in this case. Notwithstanding, the neighbour motifs 6 and 7, derived from the use of R=3, remain in close proximity to the primary motif. An alternative approach to enhance the flexibility of the extraction of neighbour motifs is to consider a maximum distance calculated not from a radius R but rather from a metric computed from the overall matrix profile variability. In this study we adopt this approach, whereby neighbour motifs are selected as the patterns that are separated from the top motif by less than twice the standard deviation of the matrix profile. For the δ¹⁸O time series, the motifs obtained using this criterion are identical to those obtained by setting R=2. However, this approach yields the same neighbour motifs for the Ca²⁺ time series as for R=2, in addition to the first neighbour motif corresponding to R=3.

The identification of higher-order motifs, other than the top (global) motif corresponding to the matrix profile minimum, is more challenging due to the dependence on the constraints that must be set in terms of the maximum distance for which sub-sequences are taken as matching. As illustrated in the preceding section, varying the radius parameter yields both matching and dissimilar patterns. Therefore, we have adopted a more flexible criterion than the radius parameter R in order to ensure appropriately constrained results and interpretable motifs. While the matrix profile calculation is only dependent on the sub-sequence length and the top motif is objectively obtained from the global minimum value of the matrix profile, subsequent motifs or neighbours depend on the adopted distance constraints. To circumvent this limitation, this study focuses on the extraction of the top motif from the palaeoclimate records. In general, the extraction of motifs has to take into account the purpose of the analysis and is very much problem-dependent in terms of what the patterns of interest are in a time series (which may differ, even for the same time series, depending on the goal of the analysis). In the case of very large datasets for which information on the patterns of interest is not available, the statistical significance of motifs can be assessed to evaluate which additional motifs are significant, other than the top motif (Castro and Azevedo, 2011).

The comparison of sub-sequence similarity is performed here based on the Euclidean distance metric. An alternative would be to consider instead dynamic time warping (DTW) as a more robust distance measure for time series (Keogh and Ratanamahatana, 2005). However, empirical comparisons showed that the Euclidean distance is competitive with or superior to more complex measures (Ding et al., 2008). Therefore, the Euclidean distance between z-normalised sub-sequences is employed as the similarity measure in this study.

Here we only discuss the influence on matrix profile results of methodological options in terms of parameters and metric selection. A further extension would be to assess how differences in the time series data would impact the matrix profile results, e.g. by applying surrogate time series methods such as the approach of Nakamura and Small (2005). Extending the methodology to account for uncertainty in dating and proxy values would be a further methodological extension particularly relevant for palaeoclimate time series.

Figure 12 Normalised motifs for the NGRIP δ¹⁸O time series: top motif and neighbour motifs 3 to 5 (see Table 2) corresponding to R=2 (a), the top motif and neighbour motifs 6 and 7 (b), and the top motif and neighbour motifs 8 to 10 (c): same convention as in Fig. 4.

Figure 13 Various time series of Ca²⁺ concentration. (a) Original values, (b) version shifted by 10 ka with the same size, and (c) version trimmed by 10 ka with a smaller size: same age model as in Fig. 1.

Figure 14 Join matrix profile for the δ¹⁸O time series and the Ca²⁺ series in Fig. 13 using a window size of 2500 years. (a) Original Ca²⁺ record, (b) shifted version, and (c) trimmed version: same age model as in Fig. 1.

5.2 Extraction of the dominant motif across the two time series

The most analogous pattern (top motif) across two time series does not need to occur at the same time in the two series. This fact is illustrated by considering the same δ¹⁸O and Ca²⁺ proxy records but artificially changing the Ca²⁺ time series in two distinct ways. The first case introduces an artificial shift in time of 500 data points (10 ka) by adding to the beginning of the record 500 points with the same value as the mean of the Ca²⁺ time series. The second case removes the first 500 data points of the Ca²⁺ series, thus obtaining a shorter record of a different length than the δ¹⁸O series. The time series are displayed in Fig. 13, and the join matrix profile between the δ¹⁸O series and these versions of the Ca²⁺ series is presented in Fig. 14. The outcomes are largely analogous, with only minor discrepancies in the join matrix profile for the different versions of the Ca²⁺ time series. In the case of the shifted version of the Ca²⁺ series (Fig. 14b), the difference between this version and the original is observed in the highest values of the join matrix profile, which are flatter. The matrix profile peaks correspond to parts of the series with no discernible temporal structure, typically a featureless noise level. In such cases the most similar sub-sequences are the stable level values (equal to the mean of the time series) that were introduced at the beginning of the record. In the case of the trimmed version of the Ca²⁺ series (Fig. 14c), it should be noted that the join matrix profile has the same length, despite the reduced length of the Ca²⁺ series. This is because the join matrix profile contains the distance of every sub-sequence in the δ¹⁸O series to the trimmed Ca²⁺ series (of a smaller size). Consequently, the number of sub-sequences (in the δ¹⁸O series) remains constant, yet the distances are calculated between each of these sub-sequences and a smaller number of sub-sequences in the Ca²⁺ series. The principal distinction between the join matrix profile values in this particular case is observed at the beginning of the record, as the initial sub-sequences of the δ¹⁸O series are most similar to the initial sub-sequences in the Ca²⁺ series, which no longer exist, having been matched to other sub-sequences in the Ca²⁺ time series.

Figure 15 Top motif across δ¹⁸O and the Ca²⁺ time series. (a) δ¹⁸O record, (b) original Ca²⁺ time series, (c) shifted version, and (d) trimmed version: same age model as in Fig. 1.

For the original version and the two versions of Ca²⁺ under consideration, the minimum value of the join matrix profile occurs at the same time corresponding to the top motif previously identified in the δ¹⁸O record at around 71 ka. Tables 5 and 6 present excerpts of the join matrix profile and profile index for the shifted and trimmed Ca²⁺ time series, respectively. In the former case, the minimum value of the matrix profile occurs at the same index value of the δ¹⁸O time series as previously defined. However, the most similar sub-sequence in the Ca²⁺ record is no longer coincidental, although it is correctly identified at 10 ka later. A comparable outcome is observed in the latter case. The profile index indicates that the most similar sub-sequence occurs at the index value of the time series corresponding to the correct time. Figure 15 displays the top motif across the δ¹⁸O time series and the different versions of the Ca²⁺ series. It can be seen that the most similar sub-sequence across the two records corresponds to the canonical DO-19. Furthermore, the same results are obtained with shifted and trimmed versions of the Ca²⁺ time series, demonstrating the robustness of the approach.

Table 5 Snippet of the matrix profile and profile index for the shifted Ca²⁺ time series: same convention as for Table 4.

Table 6 Snippet of the matrix profile and profile index for the trimmed Ca²⁺ time series: same convention as for Table 4.

6 Conclusions

In this study, the algorithmic matrix profile approach was employed to identify recurring patterns in the well-studied NGRIP palaeoclimate record. The matrix profile is dependent on a single parameter, the sub-sequence length. This is generally set by considering the typical duration of the patterns of interest, as there are no stringent criteria for its specification. In this analysis, a window size of 2500 years was considered. Shorter patterns exist in the time series and are of interest, but short-length sub-sequences can be harder to identify in terms of recurring patterns, and thus this study focused on patterns spanning around 2500 years. Consistent patterns were obtained for window sizes of 3000 and 3500 years, indicating that the results are robust to window sizes within this range.

The objective of the matrix profile approach was not to identify the abrupt DO transitions or to determine their precise timing. Rather, the objective was to characterise the abrupt transitions in a purely data-driven manner, based on the shape of the corresponding DO patterns. For the δ¹⁸O time series, the transitions corresponding to the canonical events DO-19 and DO-20, occurring at around 72 and 76 ka, respectively, are identified as the most similar ones. These form the most prominent top motif pair in the time series, indicating that analogous mechanisms may have been responsible for these abrupt climate transitions. Further transitions corresponding to the DO-12, DO-8, and DO-1 events are established as neighbouring transitions, with a similar shape characterised by an abrupt transition to warm conditions, preceded by approximately stable stadial conditions and followed by a slow return to cold conditions. The same transitions are identified in the Ca²⁺ time series, but their ordering differs. The transitions corresponding to events DO-8 and DO-12 are identified as the most similar ones in the Ca²⁺ record. These events are distinguished by an abrupt decrease in the terrestrial dust concentration, followed by a period of stable dust conditions.

The matrix profile method has also been employed to identify the most analogous pattern across the two different δ¹⁸O and Ca²⁺ time series, despite their disparate lengths. Given the high degree of similarity between the two records, the assumption of their simultaneous change, within the 20-year resolution of the records, serves as the basis for the definition of stratigraphic events in Rasmussen et al. (2014). The join matrix profile identifies a coincident top motif, which also corresponds to the DO-19 canonical event. When considering a shifted version of the Ca²⁺ timescale, the join matrix profile is able to identify the same sub-sequence as the most similar pattern across the two time series, though they are not coincident in time. This allows for accurate identification of the time shift that was introduced. A shorter version of the Ca²⁺ time series also demonstrates the ability of the join matrix profile to correctly identify matching patterns across the δ¹⁸O and Ca²⁺ records. This indicates the potential of the join matrix profile as an objective quantitative approach for matching palaeoclimate time series as an alternative to visual wiggle matching. The identification of similarities in events across distinct proxy records can assist the investigation of the key factors affecting the marine (δ¹⁸O) and continental (dust) hydrological cycles during the last 130 000 years. This can be tested by Earth system models.