{'fraction_non_alphanumeric': 0.037206376853666207, 'fraction_numerical': 0.019268451992161987, 'mean_word_length': 4.695831898036514, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Processes in Geophysics

The analysis of palaeoclimate time series is usually affected by severe methodological problems, resulting primarily from non-equidistant sampling and uncertain age models. As an alternative to existing methods of time series analysis, in this paper we argue that the statistical properties of recurrence networks – a recently developed approach – are promising candidates for characterising the system's nonlinear dynamics and quantifying structural changes in its reconstructed phase space as time evolves. In a first order approximation, the results of recurrence network analysis are invariant to changes in the age model and are not directly affected by non-equidistant sampling of the data. Specifically, we investigate the behaviour of recurrence network measures for both paradigmatic model systems with non-stationary parameters and four marine records of long-term palaeoclimate variations. We show that the obtained results are qualitatively robust under changes of the relevant parameters of our method, including detrending, size of the running window used for analysis, and embedding delay. We demonstrate that recurrence network analysis is able to detect relevant regime shifts in synthetic data as well as in problematic geoscientific time series. This suggests its application as a general exploratory tool of time series analysis complementing existing methods.


Introduction
Palaeoclimate proxy data representing past variations of environmental conditions can be obtained from various types of geological archives distributed over the Earth's surface. The study of time series of such proxies, i.e. data that encode the Correspondence to: J. F. Donges (donges@pik-potsdam.de) temporal variability of physical, chemical, biological or sedimentological properties, is a major source of information fostering our understanding of the functioning of the complex Earth system in the past, present, and future. However, nonequidistant sampling, uncertain age models, multi-scale, and multi-stable state variability as well as relatively high noise levels render the study of these proxy records a challenging problem for time series analysis.
Methods used for time series analysis can be roughly classified as linear or nonlinear. On the one hand, linear methods are based on the evaluation of certain classical statistical characteristics and assume the presence of an underlying linear stochastic process with eventually some superimposed deterministic (e.g. periodic) components Davis, 1991, 2002;Hamilton, 1994). Prominent examples that are frequently used for the analysis of realworld time series, including such obtained from geological archives (Schulz and Stattegger, 1997;Schulz and Mudelsee, 2002;Mudelsee et al., 2009;Rehfeld et al., 2011), are correlation functions and power spectra. On the other hand, nonlinear methods follow a dynamical systems point of view, implicitly assuming the presence of certain types of deterministic behaviour (Abarbanel, 1996;Kantz and Schreiber, 1997;Donner and Barbosa, 2008).
The vast majority of existing linear or nonlinear methods of time series analysis relies on the quantification of patterns of temporal dependences between observations x(t) made at different times t, i.e. aims to quantify functional relationships of the form on time t. In standard linear methods of time series analysis, f is often assumed to be a linear function; in this case, the parameters of f encode linear temporal correlations. More generally, one may consider arbitrary (i.e. not explicitly specified) deterministic relationships f , which may be characterized using concepts such as mutual information (Kantz and Schreiber, 1997).
In the following, we will refer to methods of time series analysis that are based on the quantification of temporal interrelationships between observations, e.g. correlation and mutual information functions or power spectra, as correlative methods. These clearly depend on how well the observation points are specified. In particular, in case of a nonuniform sampling of the considered time series, estimates of even simple linear characteristics can often not be expressed in a straightforward analytical way. For example, if one wishes to avoid interpolation (which leads to additional uncertainties), power spectra can be estimated using harmonic regression of the data (e.g. by means of the Lomb-Scargle periodogram;Lomb, 1976;Scargle, 1982), projection methods (Foster, 1996a,b), or a variety of alternative approaches Rehfeld et al., 2011). However, in the specific case of palaeoclimate data where typically not even the exact timing of the individual observations is sufficiently well known (Telford et al., 2004), correlative methods can have strong conceptual disadvantages.
In contrast to this large class of methods (which characterise time series from a more or less rigorous statistical point of view), alternative concepts such as fractal dimensions and generalisations thereof have been first developed in different mathematical disciplines and later applied to the characterisation of the properties of certain dynamical systems (Sprott, 2003). Statistical estimates of such measures can be obtained by a variety of different approaches, most of which take into account the spatial arrangement of observations in the (possibly reconstructed) phase space. From this perspective, the mentioned methods do not directly require knowledge about the timing of observations, i.e. can be considered as non-correlative or geometric methods, since they rely on geometric attractor properties in phase space rather than on dynamical information. In the case of palaeoclimate data with uncertain age models, geometric methods may provide a considerable alternative for statistical analysis. However, as a particular disadvantage, we note that the proper estimation of fractal dimensions usually requires a considerably larger amount of data than necessary for most correlative methods (Sprott, 2003), which are typically not available in palaeoclimatology.
Some fundamental relationships between the geometric properties of attractors in phase space (e.g. Hausdorff and box dimensions) and important invariants of the associated dynamics (e.g. Lyapunov exponents) are known to exist . Note that certain measures of dimensionality include both geometric and dynamical information, i.e. all Rényi dimensions D q for q > 1 including the information dimension D 1 (Sprott, 2003). However, besides fractal dimension estimates based on attractor topology there are only very few suitable and purely geometric methods available. Recently, it has been suggested to characterise the mutual proximity relationships of all pairs of state vectors from the sampled attractor in phase space by means of complex network methods (Zhang and Small, 2006;Yang and Yang, 2008;Xu et al., 2008;Marwan et al., 2009;. Among others, the concept of recurrence networks (RNs) Donner et al., 2010a,b) has been proven particularly useful for this purpose. Since such complex network representations of time series take only spatial information into account, they can be considered as important examples of geometric methods of time series analysis. RNs provide a set of nonlinear measures characterising the complexity of dynamical systems (Donner et al., , 2011a, e.g. allowing to distinguish periodic from chaotic dynamics. While recent findings demonstrate close interrelationships between certain RN properties and fractal dimensions (Donner et al., 2011b), the graph-theoretical measures can often be estimated with high confidence from much shorter time series than fractal dimensions. This warrants their application as a tool for windowbased analysis of non-stationary data Donner et al., 2011a). In contrast to the aforementioned approaches, transition networks (Nicolis et al., 2005) and visibility graphs (Lacasa et al., 2008) are correlative methods in the sense that they depend explicitly on the temporal ordering of observations. When considering network-based methods of time series analysis, only RNs Donner et al., 2011a;Hirata et al., 2011) and visibility graphs (Elsner et al., 2009) have been used to analyse geoscientific data. So far RN analysis is the only network-based technique that has been applied to investigate palaeoclimate proxy records. In this work, we discuss the application of RNs to studies of palaeoclimate records, with a special focus on the identification of structural changes in the dynamics that are not easily found when relying on simple linear statistics. As a benchmark example, we will mainly utilise three marine records of aeolian dust flux from Northern Africa during the last 5 Myr (million years) (Trauth et al., 2009;Marwan et al., 2009;Donner et al., 2011a;Donges et al., 2011). In Sect. 2, we present a detailed description of the considered data sets, the necessary preprocessing steps, and the general idea of RNs and their quantitative analysis. Application to typical nonlinear model systems with a systematic drift of the control parameters in Sect. 3 suggests that network statistics are well suited for identifying dynamical transitions from finite time series. Finally, in Sect. 4, we describe the results of our investigations obtained for the different palaeoclimate time series and discuss their robustness with respect to the fundamental parameters of our method.  (Tiedemann et al., 1994;deMenocal, 1995deMenocal, , 2004Larrasoaña et al., 2003).
sampled at times t i , where i is an index variable and N the number of samples, are available from three sediment cores: ODP 659 (Atlantic Ocean offshore subtropical West Africa) (Tiedemann et al., 1994), ODP 721/722 (Arabian Sea) (deMenocal, 1995(deMenocal, , 2004, and ODP 967 (Eastern Mediterranean Sea)   (Fig. 1). In addition, the benthic oxygen isotope (δ 18 O) record from ODP site 659 (Tiedemann et al., 1994) will be studied as a proxy for variations in global ice volume, which can be assumed to have a considerable impact on the continental aridification via a southward displacement of climate and vegetation zones. All time series are shown in Fig. 2.

Detrending
All considered time series {x i } show a nonlinear trend of increasing amplitude and variance towards the present. This trend reflects the successive aridification of North and East Africa and the intensification of Northern hemisphere glacial cycles during the Plio-Pleistocene (Trauth et al., 2009). To prevent corruption of the results of our analysis and significance test due to this nonlinear trend, we attempt to remove it to first order by subtracting from x i the mean of a sliding window of size W D (t i ) centered at t i for all time points t i , i.e.,  (Tiedemann et al., 1994), and of terrigenous dust flux from North Africa at ODP sites (B) 659 (Tiedemann et al., 1994), (C) 721/722 (deMenocal, 1995(deMenocal, , 2004, and (D) 967 . The horizontal red bars in panel (A) indicate two consecutive recurrence windows of length W * = 410 kyr and mutual offset ∆W * = 41 kyr as used in the analysis of  where for a chosen detrending window size W D , That is, the effective detrending window size decreases towards the time series' boundaries, resulting inx 1 =x N = 0. This simple approach avoids the complication of locally or globally fitting higher-order polynomials or performing highpass filtering given irregular sampling and uncertain dating of measurements to remove the nonlinear trend. Since RN analysis as our method of choice is a non-correlative technique and its results are permutation invariant (Sec. 2.6), spurious autocorrelations which may be introduced by the sliding window detrending do not pose a serious problem here. We will show in Sec. 4.2 that our results are robust with respect to a large range of reasonable choices of W D . Except of the detrending, no further preprocessing was applied to the data. Particularly, we do not resample the time series to obtain an evenly spaced record in the time domain, since the necessary interpolation would corrupt the results of the further analysis to be performed below (see, e.g., Rehfeld et al. (2011)).  (Tiedemann et al., 1994;deMenocal, 1995deMenocal, , 2004Larrasoaña et al., 2003).

Description of the data
Marine records of terrigenous dust flux from North Africa are an important source of information on the long-term aridification of the continent during the Plio-Pleistocene (Trauth et al., 2009) where i is an index variable and N the number of samples, are available from three sediment cores: ODP 659 (Atlantic Ocean offshore subtropical West Africa) (Tiedemann et al., 1994), ODP 721/722 (Arabian Sea) (deMenocal, 1995(deMenocal, , 2004, and ODP 967 (Eastern Mediterranean Sea)   (Fig. 1). In addition, the benthic oxygen isotope (δ 18 O) record from ODP site 659 (Tiedemann et al., 1994) will be studied as a proxy for variations in global ice volume, which can be assumed to have a considerable impact on the continental aridification via a southward displacement of climate and vegetation zones. All time series are shown in Fig. 2.

Detrending
All considered time series {x i } show a nonlinear trend of increasing amplitude and variance towards the present. This trend reflects the successive aridification of North and East Africa and the intensification of Northern hemisphere glacial cycles during the Plio-Pleistocene (Trauth et al., 2009). To prevent corruption of the results of our analysis and significance test due to this nonlinear trend, we attempt to remove it to first order by subtracting from x i the mean of a sliding window of size W D (t i ) centered at t i for all time points t i , i.e.
sampled at times t i , where i is an index variable and N the number of samples, are available from three sediment cores: ODP 659 (Atlantic Ocean offshore subtropical West Africa) (Tiedemann et al., 1994), ODP 721/722 (Arabian Sea) (deMenocal, 1995(deMenocal, , 2004, and ODP 967 (Eastern Mediterranean Sea)   (Fig. 1). In addition, the benthic oxygen isotope (δ 18 O) record from ODP site 659 (Tiedemann et al., 1994) will be studied as a proxy for variations in global ice volume, which can be assumed to have a considerable impact on the continental aridification via a southward displacement of climate and vegetation zones. All time series are shown in Fig. 2.

Detrending
All considered time series {x i } show a nonlinear trend of increasing amplitude and variance towards the present. This trend reflects the successive aridification of North and East Africa and the intensification of Northern hemisphere glacial cycles during the Plio-Pleistocene (Trauth et al., 2009). To prevent corruption of the results of our analysis and significance test due to this nonlinear trend, we attempt to remove it to first order by subtracting from x i the mean of a sliding window of size W D (t i ) centered at t i for all time points t i , i.e.,  (Tiedemann et al., 1994), and of terrigenous dust flux from North Africa at ODP sites (B) 659 (Tiedemann et al., 1994), (C) 721/722 (deMenocal, 1995(deMenocal, , 2004, and (D) 967 . The horizontal red bars in panel (A) indicate two consecutive recurrence windows of length W * = 410 kyr and mutual offset ∆W * = 41 kyr as used in the analysis of  where for a chosen detrending window size W D , That is, the effective detrending window size decreases towards the time series' boundaries, resulting inx 1 =x N = 0. This simple approach avoids the complication of locally or globally fitting higher-order polynomials or performing highpass filtering given irregular sampling and uncertain dating of measurements to remove the nonlinear trend. Since RN analysis as our method of choice is a non-correlative technique and its results are permutation invariant (Sec. 2.6), spurious autocorrelations which may be introduced by the sliding window detrending do not pose a serious problem here. We will show in Sec. 4.2 that our results are robust with respect to a large range of reasonable choices of W D . Except of the detrending, no further preprocessing was applied to the data. Particularly, we do not resample the time series to obtain an evenly spaced record in the time domain, since the necessary interpolation would corrupt the results of the further analysis to be performed below (see, e.g., Rehfeld et al. (2011)).  (Tiedemann et al., 1994), and of terrigenous dust flux from North Africa at ODP sites (B) 659 (Tiedemann et al., 1994), (C) 721/722 (deMenocal, 1995(deMenocal, , 2004, and (D) 967 . The horizontal red bars in panel (A) indicate two consecutive recurrence windows of length W * = 410 kyr and mutual offset W * = 41 kyr as used in the analysis of Sect. 4 and in Figs. 9-12.
That is, the effective detrending window size decreases towards the time series' boundaries, resulting inx 1 =x N = 0. This simple approach avoids the complication of locally or globally fitting higher-order polynomials or performing highpass filtering given irregular sampling and uncertain dating of measurements to remove the nonlinear trend. Since RN analysis as our method of choice is a non-correlative technique and its results are permutation invariant (Sect. 2.6), spurious autocorrelations which may be introduced by the sliding window detrending do not pose a serious problem here. We will show in Sect. 4.2 that our results are robust with respect to a large range of reasonable choices of W D . Except of the detrending, no further preprocessing was applied to the data. Particularly, we do not resample the time series to obtain an evenly spaced record in the time domain, since the necessary interpolation would corrupt the results of the further analysis to be performed below (see, e.g. Rehfeld et al., 2011).

Embedding
Univariate time series often reflect the dynamics of a higherdimensional complex system as viewed through some observation function. In typical situations it is possible to reconstruct the phase space trajectory using time-delay embedding, i.e. considering state vectors instead of the univariate observationsx i themselves (Packard et al., 1980;Takens, 1981). Due to the finite length of the available time series, the index i is now restricted to the range i = 1,...,N − (m − 1)τ . The embedding parameters embedding dimension m and delay τ have to be appropriately determined from the available data, e.g. using approaches such as the false nearest-neighbours (Kennel et al., 1992) and average mutual information (Fraser and Swinney, 1986) methods, respectively. Although there are good reasons for applying embedding techniques, it is known that this approach also has conceptual disadvantages and may induce spurious structures in recurrence plots and corresponding misleading results of recurrence quantification analysis (RQA) (Thiel et al., 2006). In contrast, many important dynamical invariants can be estimated from non-embedded time series as well, especially using recurrence plot-based methods (Thiel et al., 2004a). From here on we will use the simplified notation y i for reconstructed state vectors and assign to them the ages t i , respectively. While the standard approaches for determining the optimum embedding parameters typically provide feasible results in the case of many applications, the situation is considerably more challenging for palaeoclimate records: on the one hand, traditional embedding methods require equally spaced observations, so that interpolation of the available data might become necessary with all corresponding conceptual disadvantages. On the other hand, in the presence of dating uncertainties, even such interpolation is hardly possible and would lead to an enormous enhancement of uncertainty in the embedded record.
Given these methodological difficulties we attempt a compromise: (i) the embedding dimension m = 3 is a tradeoff given the relatively short time series forbidding larger embedding dimensions (Eckmann et al., 1992;Kantz and Schreiber, 1997) and the underlying high-dimensional dynamics as estimated by the false nearest-neighbours criterion (Kennel et al., 1992;Marwan et al., 2009). (ii) Using a Gaussian kernel-based estimator of the autocorrelation function adapted to irregularly sampled time series (Rehfeld et al., 2011), we find that the autocorrelation of all four time series has decayed markedly after 10 kyr (Fig. 3). Hence, we choose the delay τ to cover approximately the same time scale τ * = 10 kyr for all considered records, i.e. τ = τ * / T , where T is the average sampling time (Table 1) and x denotes the integer part of x. This yields τ 1 = 2 for ODP site 659, τ 2 = 5 for site 721/722, and τ 3 = 27 for site 967. A promising technique for consistent embedding of irregularly sampled time series is based on Legendre polynomials (Gibson et al., 1992) and should be explored in future studies.

Embedding
Univariate time series often reflect the dynamics of a higherdimensional complex system as viewed through some observation function. In typical situations it is possible to reconstruct the phase space trajectory using time-delay embedding, i.e., considering state vectors instead of the univariate observationsx i themselves (Packard et al., 1980;Takens, 1981). Due to the finite length of the available time series, the index i is now restricted to the range i = 1,...,N − (m − 1)τ . The embedding parameters embedding dimension m and delay τ have to be appropriately determined from the available data, e.g., using approaches such as the false nearest-neighbours (Kennel et al., 1992) and average mutual information (Fraser and Swinney, 1986) methods, respectively. Although there are good reasons for applying embedding techniques, it is known that this approach also has conceptual disadvantages and may induce spurious structures in recurrence plots and corresponding misleading results of recurrence quantification analysis (RQA) (Thiel et al., 2006). In contrast, many important dynamical invariants can be estimated from non-embedded time series as well, especially using recurrence plot-based methods (Thiel et al., 2004a). From here on we will use the simplified notation y i for reconstructed state vectors and assign to them the ages t i , respectively. Whereas the standard approaches for determining the optimum embedding parameters typically provide feasible results in the case of many applications, the situation is considerably more challenging for palaeoclimate records: On the one hand, traditional embedding methods require equally spaced observations, so that interpolation of the available data might become necessary with all corresponding conceptual disadvantages. On the other hand, in the presence of dating uncertainties, even such interpolation is hardly possible and would lead to an enormous enhancement of uncertainty in the embedded record.
Given these methodological difficulties we attempt a compromise: (i) The embedding dimension m = 3 is a tradeoff given the relatively short time series forbidding larger embedding dimensions (Eckmann et al., 1992;Kantz and Schreiber, 1997) and the underlying high-dimensional dynamics as estimated by the false nearest-neighbours criterion (Kennel et al., 1992;Marwan et al., 2009). (ii) Using a Gaussian kernel-based estimator of the autocorrelation function adapted to irregularly sampled time series (Rehfeld et al., 2011), we find that the autocorrelation of all four time series has decayed markedly after 10 kyr (Fig. 3). Hence, we choose the delay τ to cover approximately the same time scale τ * = 10 kyr for all considered records, i.e. where ∆T is the average sampling time (Tab. 1). This yields τ 1 = 2 for ODP site 659, τ 2 = 5 for site 721/722, and τ 3 = 27 for site 967 corresponding to τ * = 10 kyr. A promising technique for consistent embedding of irregularly sampled time series is based on Legendre polynomials (Gibson et al., 1992) and should be explored in future studies.

Windowed analysis
For detecting structural changes in the dynamics encoded by the time series, we slide a window over the embedded record {y i } and perform the subsequent analysis for the data contained in each window separately. However, the records under study are quite heterogeneous with respect to their basic sampling properties (Tab. 1). The average sampling time ∆T differs widely across the records. In order to assure comparability of our results uncovered from the different time series, the most natural approach is to choose windows of a fixed size W * in units of time. However, this approach has two disadvantages: The exact timing t i of the available observations is not known as is the case for most geological proxy records, and due to the non-uniform sampling rates, different windows would contain different amounts of data. While the latter is not problematic for statistical tests against homogeneity of the distribution of values in different windows, a quantitative comparison of statistical characteristics of the associated RNs (see Sec. 2.5) is not possible. Therefore, in the following, we will proceed in a different way by prescribing both the window size W and step size ∆W for RN analysis measured in units of sampling points. In order

Windowed analysis
For detecting structural changes in the dynamics encoded by the time series, we slide a window over the embedded record {y i } and perform the subsequent analysis for the data contained in each window separately. However, the records under study are quite heterogeneous with respect to their basic sampling properties (Table 1). The average sampling time T differs widely across the records. In order to assure comparability of our results uncovered from the different time series, the most natural approach is to choose windows of a fixed size W * in units of time. However, this approach has two disadvantages: The exact timing t i of the available observations is not known as is the case for most geological proxy records, and due to the non-uniform sampling rates, different windows would contain different amounts of data. While the latter is not problematic for statistical tests against homogeneity of the distribution of values in different windows, a quantitative comparison of statistical characteristics of the associated RNs (see Sect. 2.5) is not possible. Therefore, in the following, we will proceed in a different way by prescribing both the window size W and step size W for RN analysis measured in units of sampling points. In order to derive W and W from the desired quantities in units of time, W * and W * , we divide by the average sampling time, In turn, the actual window size W * (t i ) is determined by the average sampling time in the size-W window centred around Nonlin. Processes Geophys., 18, 545-562, 2011 www.nonlin-processes-geophys.net/18/545/2011/ t i . For a particular choice of W * and the associated step size W * in units of time, the resulting values of W and W , the mean window widths, and step sizes as well as the corresponding standard deviations are given in Table 1.
The simple approach for determining the window size described above guarantees that the windows cover approximately the same time span for all records and positions within the time series. While most sampling intervals take values close to the mean, there are distinct outliers, which most likely correspond to missing data due to incomplete core recovery, hiata or disturbances of the sediment such as turbidites (Fig. 4). Nevertheless, the standard deviation of window size σ (W * ) is still small in comparison to the average window size W * (Table 1), which suggests that statistical characteristics computed for different windows can still be quantitatively compared in a reasonable way.
Formally, the data series {y µ i } within the µ-th window, µ = 1,2,..., N−W W , is given by where from here on i = 1,...,W . We use the window's midpoint's timing to attach an age to the scalar network measures g µ calculated from the data within the µ-th window.

Recurrence network analysis
Recurrence in phase space is a basic property of complex dynamical systems. Since the seminal work of Poincaré (1890), it is known that under rather general conditions, dynamical systems tend to return arbitrarily close to their previous states in the long-term limit. In the last decades, the recurrence property has attracted considerable interest, since it has been shown that essential information about the main dynamical properties is contained in the temporal pattern of mutual recurrences of a state (Thiel et al., 2004b;Robinson and Thiel, 2009). Particularly, the visual representations of recurrence plots Marwan et al., 2007) have found wide use, which are most commonly expressed by a binary recurrence matrix where for the µ-th window, ε is the recurrence threshold and (·) denotes the Heaviside function. In the following we use the supremum norm · to measure distances in the reconstructed phase space of the considered observable y (see Fig. 5 for examples). The appropriate choice of the important parameter ε is discussed below.
It turned out that recurrence plots show distinct line structures, whose length distribution can be used for defining suitable measures of complexity in terms of RQA, or for J. F. Donges et al.: Identification of dynamical transitions in marine palaeoclimate records by RN analysis 5 to derive W and ∆W from the desired quantities in units of time, W * and ∆W * , we divide by the average sampling time, where x denotes the integer part of x. In turn, the actual window size W * (t i ) is determined by the average sampling time in the size-W window centred around t i . For a particular choice of W * and the associated step size ∆W * in units of time, the resulting values of W and ∆W , the mean window widths, and step sizes as well as the corresponding standard deviations are given in Tab. 1. The simple approach for determining the window size described above guarantees that the windows cover approximately the same time span for all records and positions within the time series. While most sampling intervals take values close to the mean, there are distinct outliers, which most likely correspond to missing data due to incomplete core recovery, hiata or disturbances of the sediment such as turbidites (Fig. 4A). Nevertheless, the standard deviation of window size σ(∆W * ) is still small in comparison to the average window size ∆W * (Tab. 1), which suggests that statistical characteristics computed for different windows can still be quantitatively compared in a reasonable way.
Formally, the data series {y µ i } within the µ-th window, µ = 1,2,..., is given by where from here on i = 1,...,W . We use the window's midpoint's timing to attach an age to the scalar network measures g µ calculated from the data within the µ-th window.

Recurrence network analysis
Recurrence in phase space is a basic property of complex dynamical systems. Since the seminal work of Poincaré (1890), it is known that under rather general conditions, dynamical systems tend to return arbitrarily close to their previous states in the long-term limit. In the last decades, the recurrence property has attracted considerable interest, since it has been shown that essential information about the main dynamical properties is contained in the temporal pattern of mutual recurrences of a state (Thiel et al., 2004b;Robinson and Thiel, 2009). Particularly, the visual representations of recurrence plots Marwan et al., 2007) have found wide use, which are most commonly expressed by a binary recurrence matrix where for the µ-th window, ε is the recurrence threshold and Θ(·) denotes the Heaviside function. In the following we . The distribution for the δ 18 O record at ODP site 659 is visually almost indistinguishable from that of the corresponding dust flux record and therefore not shown. The PDFs were estimated using a Gaussian kernel with bandwidth h = σ(∆T )(N − 1) −1/5 (Tab. 1) following Scott's rule (Scott, 1982). (B,C,D) Temporal variation of the sampling times for the three dust flux records.
use the supremum norm · to measure distances in the reconstructed phase space of the considered observable y (see Fig. 5 for examples). The appropriate choice of the important parameter ε is discussed below. It turned out that recurrence plots show distinct line structures, whose length distribution can be used for defining suitable measures of complexity in terms of RQA, or for estimating dynamical invariants such as correlation dimension, 2nd-order Rényi entropy, or generalised mutual information (Marwan et al., 2007). In the context of palaeoclimate research, recurrence plots and RQA have been successfully applied for tracing dynamical changes Marwan et al., 2003) and aligning records with different agedepth models (Marwan et al., 2002). RQA is a correlative method of time series analysis, as it explicitly relies on temporal structures in the form of diagonal and vertical lines.
Recently, it has been suggested to approach recurrence  Table 1) following Scott's rule (Scott, 1982). (B, C, D) Temporal variation of the sampling times for the three dust flux records.
estimating dynamical invariants such as correlation dimension, 2nd-order Rényi entropy, or generalised mutual information (Marwan et al., 2007). In the context of palaeoclimate research, recurrence plots and RQA have been successfully applied for tracing dynamical changes Marwan et al., 2003) and aligning records with different agedepth models (Marwan et al., 2002). RQA is a correlative method of time series analysis, as it explicitly relies on temporal structures in the form of diagonal and vertical lines.
Recently, it has been suggested to approach recurrence matrices from a complex network perspective by identifying (δ ij denoting Kronecker's delta) with the adjacency matrix of a complex network associated to the underlying time series  1 . This Table 1. Basic properties of the analysed palaeoclimate time series. N is the number of samples contained in the time series, T the mean sampling interval, and σ ( T ) the standard deviation of sampling intervals. For a desired window size W * = 410 kyr and step size W * = 41 kyr (as chosen in Sect. 4.1 for RN analysis), W and W give the corresponding window and step size (in numbers of observations), W * and W * the average effective window and step size, and σ (W * ) and σ ( W * ) the associated standard deviations (in units of time).   matrices from a complex network perspective by identifying (δ ij denoting Kronecker's delta) with the adjacency matrix of a complex network associated to the underlying time series  1 . This analogy implies that each sampled state vector is assigned a vertex in the RN, where two vertices are linked if the corresponding state vectors are recurrent, i.e., mutually close, in phase space (Fig. 6). According to the conventions of Sec. 2.3, each vertex i in the µ-th window has an age t µ i = t (µ−1)∆W +i attached to it. To simplify the notation when defining network measures, we will drop the window index µ in the following.
The edge density measures which fraction of the maximum theoretically possible number W (W − 1)/2 of undirected edges is present in the RN, where the number of vertices W is determined by the chosen recurrence window size. ρ(ε) is equivalent to the recurrence rate in traditional RQA. The properties of the resulting RNs (parameterised by the single parameter ε) have been shown to trace structures in phase space corresponding to dynamically invariant objects (Donner et al., 2010a, in press) as well as changes in the dynamical behaviour of arbitrary time series Donner et al., 2011). For detecting bifurcations in time series, global-scale network characteristics of complex network theory are of main interest (Newman, 2003;da Costa et al., 2007). Here we will focus on the following four measures: (i) Transitivity T : The transitivity of an unweighted and undirected network characterises the overall probability that two randomly chosen neighbours of an also randomly chosen vertex are connected analogy implies that each sampled state vector is assigned a vertex in the RN, where two vertices are linked if the corresponding state vectors are recurrent, i.e. mutually close, in phase space (Fig. 6). According to the conventions of Sect. 2.3, each vertex i in the µ-th window has an age t µ i = t (µ−1) W +i attached to it. To simplify the notation when defining network measures, we will drop the window index µ in the following. The edge density measures which fraction of the maximum theoretically possible number W (W − 1)/2 of undirected edges is present in the RN, where the number of vertices W is determined by the chosen recurrence window size. ρ(ε) is equivalent to the recurrence rate in traditional RQA.
geoscientifically relevant applications of data analysis, such as dendrograms in agglomerative cluster analysis, or nonlinear decomposition of multivariate data using isometric feature mapping (Gámez et al., 2004).
The properties of the resulting RNs (parameterised by the single parameter ε) have been shown to trace structures in phase space corresponding to dynamically invariant objects (Donner et al., , 2011b as well as changes in the dynamical behaviour of arbitrary time series Donner et al., 2011a). For detecting bifurcations in time series, global-scale network characteristics of complex network theory are of main interest (Newman, 2003;. Here we will focus on the following four measures: i. Transitivity T : the transitivity of an unweighted and undirected network characterises the overall probability that two randomly chosen neighbours of an also randomly chosen vertex are connected (Newman, 2003). In case of RNs, T serves as a measure for the regularity of the dynamics as encoded in the RN's mesoscopic structure  ). The two-dimensional graph visualisation has been obtained with the software package GUESS using a force-directed placement algorithm (http://graphexploration.cond.org). It is important to note that in this visualisation, node positions are determined by the aforementioned algorithm and do not correspond to a projection of the node coordinates in the reconstructed three-dimensional phase space. (Newman, 2003). In case of RNs, T serves as a measure for the regularity of the dynamics as encoded in the RN's mesoscopic structure . Specifically, regular dynamics (e.g., on a periodic orbit) is typically characterised by higher values of the transitivity T than chaotic dynamics. T can furthermore be interpreted as a global measure of the underlying attractive set's effective dimensionality d (Donner et al., in press), i.e., the theoretical result is T = (3/4) d when using the supremum norm in phase space. For continuous-time systems, this implies T = 3/4 for a periodic orbit and T < 3/4 for chaotic dynamics. However, for small numbers of vertices (state vectors) W as used in this work the estimated values of T will deviate from these theoretical expectations (Donner et al., in press).
When dealing with short time series (segments) as it is the case in this work, transitivity is a more robust measure than the related global clustering coefficient C (Watts and Strogatz, 1998;Newman, 2003), since the latter gives relatively more weight to sparsely sampled regions in phase space (vertices with low degree k) .
(ii) Average path length L: The average path length is defined as the mean value of the shortest path lengths l ij between all mutually reachable pairs of vertices (i,j) (measured in terms of geodesic graph distance, i.e., the minimum number of edges that have to be traversed on any path connecting the vertices i and j) (Watts and Strogatz, 1998;Newman, 2003). A pair of vertices (i,j) is called mutually reachable if there exists at least one path connecting i and j. Since for comparable values of ε, the average distances along different types of orbits typically differ significantly, changes in L can be used as sensitive indicators of dynamical transitions .
(iii) Assortativity R: A complex network is called assortative if vertices tend to connect preferentially to vertices of a similar number of connections (degree k i = j A ij ). On the other hand, it is called disassortative if vertices of high degree prefer to link to vertices of low degree, and vice versa (Newman, 2002). This assorta- . The two-dimensional graph visualisation has been obtained with the software package GUESS using a force-directed placement algorithm (http://graphexploration.cond.org). It is important to note that in this visualisation, node positions are determined by the aforementioned algorithm and do not correspond to a projection of the node coordinates in the reconstructed three-dimensional phase space.
Specifically, regular dynamics (e.g. on a periodic orbit) is typically characterised by higher values of the transitivity T than chaotic dynamics. T can furthermore be interpreted as a global measure of the underlying attractive set's effective dimensionality d (Donner et al., 2011b), i.e. the theoretical result is T = (3/4) d when using the supremum norm in phase space. For continuoustime systems, this implies T = 3/4 for a periodic orbit and T < 3/4 for chaotic dynamics. However, for small numbers of vertices (state vectors) W as used in this work the estimated values of T will deviate from these theoretical expectations (Donner et al., 2011b).
When dealing with short time series (segments) as it is the case in this work, transitivity is a more robust measure than the related global clustering coefficient C (Watts and Strogatz, 1998;Newman, 2003), since the latter gives relatively more weight to sparsely sampled regions in phase space (vertices with low degree k) (Donner et al., , 2011a.
ii. Average path length L: the average path length is defined as the mean value of the shortest path lengths l ij between all mutually reachable pairs of vertices (i,j ) (measured in terms of geodesic graph distance, i.e. the minimum number of edges that have to be traversed on any path connecting the vertices i and j ) (Watts and Strogatz, 1998;Newman, 2003 iii. Assortativity R: a complex network is called assortative if vertices tend to connect preferentially to vertices with a similar number of connections (degree k i = j A ij ). On the other hand, it is called disassortative if vertices of high degree prefer to link to vertices of low degree, and vice versa (Newman, 2002). This assortativity property can be quantified by the Pearson correlation coefficient between the degrees k i ,k j of the vertices on both ends of all L = j >i A ij edges (i,j ), where is the mean of the average edge end-point degree (k i + k j )/2 . In the RN context, R can be considered as a measure for the local continuity of www.nonlin-processes-geophys.net/18/545/2011/ Nonlin. Processes Geophys., 18, 545-562, 2011 the phase space density of state vectors .
iv. Network diameter D: the network diameter is the maximum geodesic (shortest-path) distance between all mutually reachable pairs of vertices in the network (Newman, 2003). From this definition, there are obvious relationships with the average path length L, which are expected to lead to strong correlations between both measures .
In order to apply RNs in a sliding window analysis, a reference framework is necessary. Here, we consider a dataadaptive choice of ε that guarantees for a fixed edge density ρ of 5 %, which has been found a reasonable choice in previous studies . One should note, however, that even with this choice the characteristics of RNs can only be compared in a meaningful way if the network size W is kept fixed (see Sect. 2.1). Among the considered complex network measures, T and R are mainly affected by finite-sample problems otherwise, whereas L and D explicitly depend on ε and W (Donner et al., 2010a).

Significance test
We perform a relatively simple statistical test of whether the network characteristics in a certain time interval differ significantly from the general network characteristics expected given the phase space distribution of state vectors y i from the whole detrended and embedded record and a certain recurrence window size W . The corresponding null hypothesis is that the network measures observed for a certain window are consistent with being calculated from a random draw of W state vectors from the prescribed phase space distribution induced by the complete detrended time series. We can justly assume a thus randomised embedded time series without losing essential information, because all network measures g(·) considered here are permutation-invariant when considering a fixed subset of state vectors y 1 ,...,y W . More formally, g(y 1 ,...,y W ) = g(y π(1) ,...,y π(W ) ) for arbitrary permutations π. A similar test for RQA measures requires a more advanced method (Schinkel et al., 2009). In order to create an appropriate null-model, we use the following approach: (i) draw randomly W state vectors from the embedded time series (corresponding to the window size chosen for the original data), (ii) construct a RN from this set of state vectors, and (iii) calculate the network measures of interest. Repeating this procedure sufficiently many times, we obtain a test distribution for each of the network measures and estimate its 0.05 and 0.95 quantiles that can be interpreted as 90 % confidence bounds. The proposed significance test can be viewed as a test against stationarity of the higher-order geometrical properties of the time series that are quantified by qualitatively different RN measures.

Implementation
We implemented the above described methods using the programming language Python (van Rossum and Drake, 2006), the packages NumPy (Oliphant, 2006) and SciPy (Jones et al., 2011) as well as embedded C++ code. Complex network measures have been calculated employing the Python package igraph .

Dynamical transitions in model systems
To validate the proposed methodology for detecting transitions in time series based on RNs, we apply it to the logistic map and the Lorenz system with drifting bifurcation parameter as paradigmatic examples of discrete and continuous-time dynamical systems, respectively. While step-like changes of bifurcation parameters have already been studied for discrete ) and continuous-time dynamical systems (Zou et al., 2010;Donner et al., 2011a), here we are particularly interested in the effect of transients, which are expected to be present in real-world systems and, hence, data extracted from them. We will check whether the global network quantifiers described above are able to detect transitions in the system's dynamics induced by bifurcations due to a slowly changing control parameter. For this purpose we are specifically looking for time intervals (or equivalently values of the bifurcation parameter) where one or more of the considered network quantifiers undergo sudden changes. This requires taking into account the measures' interpretation in terms of dynamical systems theory (Sect. 2.5). Furthermore, we will study how their performance and the level of resolved detail depend on the window size W . This analysis particularly shows that the window sizes W chosen for the RN analysis of terrigenous dust flux records (Table 1) are indeed appropriate for detecting bifurcations.

Logistic map
We iterate the logistic map while varying the bifurcation parameter linearly from r 1 = 3.8 to r M = 3.9 in M = 10 000 equidistant steps setting r = 1 × 10 −5 (Fig. 7), similar to Trulla et al. (1996). We analyse the resulting time series {x i } without embedding or detrending. The transition from chaotic to 3-periodic dynamics after an interior crisis at r = 1 + danger of confusing statistical fluctuations with proper dynamical changes substantially. Time series from geological archives are typically characterised by a variety of different types of nonstationarities, including (i) changes in the longterm mean or variance of the recorded proxies, (ii) variations in the amplitudes of almost periodic variability components (e.g., such attributed to Milankovich-type variations caused by periodic changes in the Earth's orbit), or (iii) even multimodal behaviour (e.g., transitions between glacial and interglacial periods). All these three types of nonstationarities are contained in our data ( Fig. 2 and Trauth et al. (2009)). While these different phenomena can be analysed using more specific methodological approaches, we propose RN analysis as a general exploratory tool for detecting time intervals containing changes in the dominating type of dynamical behaviour. In the following, we will illustrate the robustness of this approach for the four marine records introduced in Sec. 2.1 and briefly discuss the possible climatological background of the observed dynamical changes.

Time-dependence of network properties
We consider the four marine palaeoclimate records embedded in a three-dimensional reconstructed phase space with a time delay of approximately τ * = 10 kyr, resulting in the embedding parameters described in Sec. 2.3. For an initial inspection, we use recurrence windows of size W * = 410 kyr with a mutual offset of subsequent windows ∆W * = 41 kyr. Note that the latter two parameter choices correspond to those used in previous work on the ODP site 659 dust flux record Donner et al., 2011). The selection of both parameters results from a compromise between systems Donner et al., , 2011b, T and R abruptly increase to their maximum value of 1 following this transition, whereas at the same time L and D sharply decrease to their minimum value of 1. Among all the four measures, T and R most clearly detect the termination of the period-doubling cascade following the period-3 behaviour at the accumulation point r ≈ 3.849, while T , L and D highlight the merger of the subsequently formed three chaotic bands at the interior crisis at r ≈ 3.857 (Wackerbauer et al., 1994). The latter transition is only weakly visible in R. Additionally, much fine-structure is resolved by the network measures, e.g. a narrow period-4 window at r 3.89 that is most clearly indicated by an increased transitivity T across all W . Generally, the transitions appear more and more blurred as W increases, which is due to the growing number of samples from both periodic and chaotic dynamical regimes contained in the recurrence windows when sliding over the bifurcation point. In consequence, some of the narrow periodic windows appearing for r < 3.83 and r > 3.86 are only visible for small recurrence window sizes W . As a rule of thumb, we can expect a periodic/chaotic window of width w r embedded within a chaotic/periodic background to be detectable if w r W r.
Another notable feature is that both L and D show a clear tendency to increase with growing W in the chaotic parameter ranges ( Fig. 7b and d). This is theoretically expected, since both measures are extensive, i.e. they depend explicitly and nonlinearly on the number of vertices W in the RN for a general phase space distribution of state vectors as induced by chaotic dynamics ). In contrast, L and D do not change with W in the periodic windows, most notably in the large period-3 window of the logistic map ( Fig. 7b and d). We can explain this behaviour by recalling that for discrete-time systems in a p-periodic regime, the RN reduces to a set of p fully connected components . Following the definitions in Sect. 2.5, this in turn leads to L = D = 1 in any periodic regime and independent of W .

Lorenz system
To illustrate the performance of windowed RN analysis for detecting transitions in continuous-time dynamical systems, we consider the Lorenz system with a time-dependent bifurcation parameter r = r(t), high temporal resolution of the finally produced RN measures (small W * , ∆W * ) and larger statistical confidence in the results (large W * ). The choice of W * is more critical than that of ∆W * , because the former directly influences the number of vertices W in the RNs via Eq. (6). Since a formal criterion for determining an optimal choice of W * and ∆W * is not available so far, we study the robustness of our results with respect to variations in the more critical parameter W * in Sec. 4.2.
We additionally apply local detrending by removing the long-term average taken over windows of W * D = 500 kyr, where which has not been considered in the aforementioned studies. As we will show in the following, the main features recovered by our analysis are not qualitatively changed when applying detrending. However, this step appears relevant in other kinds of statistical analyses, e.g., for estimating spectrograms or time-dependent coefficients of autoregressive processes, since the data show considerable long-term trends in both mean and variance (Fig. 2). Regarding the transitivity (Fig. 9), we find a synchronous behaviour of the two geographically distinct records at ODP sites 659 and 721/722 during the Pliocene (∼5.3-2.6 Myr BP (before present)) and Early Pleistocene (2.6-1.0 Myr BP) 2 , Transitivity T , (B) average path length L, (C) assortativity R, and (D) diameter D for varying recurrence window size W for the Lorenz system (Eq. 19) with drifting bifurcation parameter r (see text) and initial condition (x 0 ,y 0 ,z 0 ) = (10,10,10). Because we are interested in the performance of our method for scalar time series, we chose the x-component of the trajectory sampled with sampling time t = 0.05 and embed it with embedding dimension m = 3 and delay τ = 15. W was varied linearly in the interval [100,600], the recurrence window step size was fixed to W = 10 samples. We varied the recurrence threshold ε to yield a fixed edge density ρ = 0.05 . Vertical dashed lines indicate the critial values of r discussed in the main text.

Fig. 8. (A)
While the system is evolving, r increases linearly from r 0 = 160 at time t 0 = 0 to r f = 170 at time t f = 500, i.e.
For consistency with the analysis of scalar palaeoclimate time series to be performed below, we use an embedding of the x-component time series for RN analysis without prior detrending (see caption of Fig. 8 for details). The RN measures indicate two major transitions towards increasingly irregular dynamics at r ≈ 161 and r ≈ 166.5 (Fig. 8). The former possibly reflects an initial transient due to the chosen initial condition. The latter agrees well with the major shift from periodic (large T , large L and D for continuous-time systems; Donner et al., , 2011bZou et al., 2010) to chaotic (small T , small L and D) behaviour which is present in the Lorenz system's non-transient bifurcation scenario at r ≈ 166 Donner et al., 2011a). On a shorter time scale, the path-based measures L and D among others detect weaker transitions at r ≈ 163.5, r ≈ 164.5 and r ≈ 166. Note that one has to be careful when comparing these results to bifurcation studies where distinct realisations of the Lorenz system with fixed parameter r (not varying in time) are studied (e.g. Donner et al., 2011a), since transients influence the results and cannot be excluded by construction when r is continuously varied in time. However, our results are consistent with the work of Trulla et al. (1996) who observed that in transient scenarios bifurcations may appear for larger bifurcation parameters than in their nontransient equivalents. The dependence of the results on the recurrence window size W is more pronounced than that described above for the logistic map. This is likely due to the fact that transients play a larger role in continuous-time systems like the Lorenz model than in discrete-time systems.

Dynamical transitions in palaeoclimate records
Our studies in the previous section demonstrated that RN analysis can be meaningfully applied for detecting dynamical transitions in non-stationary time series from different Nonlin. Processes Geophys., 18, 545-562, 2011 www.nonlin-processes-geophys.net/18/545/2011/ model systems by applying this kind of analysis to running windows. This is a necessary, but not sufficient condition for ensuring the feasibility of RN analysis for detecting regime shifts in palaeoclimate records as well. However, the application of our simple significance test (Sect. 2.6) diminishes the danger of confusing statistical fluctuations with proper dynamical changes substantially. Time series from geological archives are typically characterised by a variety of different types of nonstationarities, including (i) changes in the long-term mean or variance of the recorded proxies, (ii) variations in the amplitudes of almost periodic variability components (e.g. such attributed to Milankovich-type variations caused by periodic changes in the Earth's orbit), or (iii) even multimodal behaviour (e.g. transitions between glacial and interglacial periods). All these three types of nonstationarities are contained in our data ( Fig. 2 and Trauth et al., 2009). While these different phenomena can be analysed using more specific methodological approaches, we propose RN analysis as a general exploratory tool for detecting time intervals containing changes in the dominating type of dynamical behaviour. In the following, we will illustrate the robustness of this approach for the four marine records introduced in Sect. 2.1 and briefly discuss the possible climatological background of the observed dynamical changes.

Time-dependence of network properties
We consider the four marine palaeoclimate records embedded in a three-dimensional reconstructed phase space with a time delay of approximately τ * = 10 kyr, resulting in the embedding parameters described in Sect. 2.3. For an initial inspection, we use recurrence windows of size W * = 410 kyr with a mutual offset of subsequent windows of W * = 41 kyr. Note that the latter two parameter choices correspond to those used in previous work on the ODP site 659 dust flux record Donner et al., 2011a). The selection of both parameters results from a compromise between high temporal resolution of the finally produced RN measures (small W * , W * ) and larger statistical confidence in the results (large W * ). The choice of W * is more critical than that of W * , because the former directly influences the number of vertices W in the RNs via Eq. (6). Since a formal criterion for determining an optimal choice of W * and W * is not available so far, we study the robustness of our results with respect to variations in the more critical parameter W * in Sect. 4.2.
We additionally apply local detrending by removing the long-term average taken over windows of W * D = 500 kyr, where which has not been considered in the aforementioned studies. As we will show in the following, the main features recovered by our analysis are not qualitatively changed when applying detrending. However, this step appears relevant in Here we used a detrending window size W * D = 500 kyr, recurrence window size W * = 410 kyr and step size ∆W * = 41 kyr, embedding dimension m = 3 and delay τ * = 10 kyr. The recurrence threshold ε was chosen adaptively to yield a fixed edge density ρ = 0.05. The grey bars represent the 5% and 95% quantiles with respect to the test distribution obtained from 10,000 realisations of our null-model for each record separately. Vertical dashed lines indicate the detected epochs of transitions discussed in the main text.
including two periods of extraordinarily large values of T at about 3.45-3.05 and 2.2-2.1 Myr BP, related to pronounced clusters of vertices shown in Figs. 6B,C. The first of these periods results from a time interval of strongly suppressed and almost constant dust flux in the Mid Pliocene (see Fig. 2), while the latter one coincides with a period of almost periodic Milankovich-type variations (Trauth et al., 2009). We note that it is known (and empirically understood) that both types of dynamics typically lead to large values of T Donner et al., 2010a, in press;Zou et al., 2010), so that this result is consistent with theoretical expectations. During the Early Pleistocene, the signatures at both sites decouple from each other, which could be the result of an enhancement of the atmospheric Walker circulation (Ravelo et al., 2004). For the last about 1.5 Myr, the variations of transitivity become more similar between ODP sites 721/722 and 967, particularly highlighting the Mid Pleistocene transition between 1.2 and 0.7 Myr BP (Fig. 6A), which corre-1.0 Myr BP is not motivated stratigraphically, but climatologically, i.e., by the onset of the Mid-Pleistocene transition around 1.0 Myr BP.  other kinds of statistical analyses, e.g. for estimating spectrograms or time-dependent coefficients of autoregressive processes, since the data show considerable long-term trends in both mean and variance (Fig. 2).
Regarding the transitivity (Fig. 9), we find a synchronous behaviour of the two geographically distinct records at ODP sites 659 and 721/722 during the Pliocene (∼5.3-2.6 Myr BP; before present) and Early Pleistocene (2.6-1.0 Myr BP) 2 , including two periods of extraordinarily large values of T at about 3.45-3.05 and 2.2-2.1 Myr BP, related to pronounced clusters of vertices shown in Figs. 6b and c. The first of these periods results from a time interval of strongly suppressed and almost constant dust flux in the Mid Pliocene (see Fig. 2), while the latter one coincides with a period of almost periodic Milankovich-type variations (Trauth et al., 2009). We note that it is known (and empirically under-f dynamical transitions in marine palaeoclimate records by RN analysis d from 9, (C) an clig win-10 kyr nd deptively ent the tained sepatransif T at unced se ped and ig. 2), t peri-). We t both (Maret al., xpect both ult of (Ravons of 1/722 trancorreically, .0 Myr sponds to a change in the dominating Milankovich-type periodicity. The results obtained for the average path length L (Fig. 10) are mostly consistent with these findings, also highlighting the Mid Pliocene, Early Pleistocene, and Mid Pleistocene as periods with changes in the long-term dust flux variability. Specifically, L tends to show significant peaks at abrupt change points between regular and more erratic climate variability, as indicated by T (see Marwan et al. (2009) for a theoretical explanation of this behaviour).
The oxygen isotope anomaly obtained from the analysis of benthic foraminifera characterises a distinctively different climatic parameter (i.e., global ice volume) than terrigenous dust flux, so that it can be expected that the variability recorded by this proxy differs from that of the dust flux. An inspection of the RN properties indeed confirms this expectation. Specifically, the transitivity T does not show any systematic maxima at all (indicating time intervals with regular long-term dynamics) (Fig. 9A), which is in clear contrast to the aeolian dust flux. The average path length L shows significant maxima around 2.9 Myr BP (possibly being related to the intensification of Northern hemisphere glaciation at around this time), between 1.8 and 1.3 Myr BP (consistent with the corresponding results for the dust flux records, suggesting a high-latitude mechanism behind the large-scale climatic changes during this time period), and after about 900 kyr BP (possibly resulting from the glacial terminations and inceptions with a rather long -roughly 100 kyr -period- stood) that both types of dynamics typically lead to large values of T Donner et al., , 2011bZou et al., 2010), so that this result is consistent with theoretical expectations. During the Early Pleistocene, the signatures at both sites decouple from each other, which could be the result of an enhancement of the atmospheric Walker circulation (Ravelo et al., 2004). For the last about 1.5 Myr, the variations of transitivity become more similar between ODP sites 721/722 and 967, particularly highlighting the Mid Pleistocene transition between 1.2 and 0.7 Myr BP (Fig. 6a), which corresponds to a change in the dominating Milankovich-type periodicity. The results obtained for the average path length L (Fig. 10) are mostly consistent with these findings, also highlighting the Mid Pliocene, Early Pleistocene, and Mid Pleistocene as periods with changes in the long-term dust flux variability. Specifically, L tends to show significant peaks at abrupt change points between regular and more erratic climate variability, as indicated by T (see Marwan et al. (2009) for a theoretical explanation of this behaviour).
The oxygen isotope anomaly obtained from the analysis of benthic foraminifera characterises a distinctively different climatic parameter (i.e. global ice volume) than terrigenous dust flux, so that it can be expected that the variability recorded by this proxy differs from that of the dust flux. An inspection of the RN properties indeed confirms this expectation. Specifically, the transitivity T does not show any systematic maxima at all (Fig. 9a), which is in clear contrast J. F. Donges et al.: Identification of dynamical transitions in marine p  icity) (Fig. 10A). Figures 11 and 12   to the aeolian dust flux. The average path length L shows significant maxima around 2.9 Myr BP (possibly being related to the intensification of Northern hemisphere glaciation at around this time), between 1.8 and 1.3 Myr BP (consistent with the corresponding results for the dust flux records, suggesting a high-latitude mechanism behind the large-scale climatic changes during this time period), and after about 900 kyr BP (possibly resulting from the glacial terminations and inceptions with a rather long -roughly 100 kyr -periodicity) (Fig. 10a).
Figures 11 and 12 additionally show the time variability of the two other RN properties assortativity R and diameter D. Since the latter one is closely related to the average path length L , the variability of both measures is very similar. Moreover, we also find some much weaker similarities between the temporal variability patterns of transitivity T and assortativity R, which are less pronounced, since both properties characterise less obviously related aspects of the network geometry in phase space. Specifically, the time interval of suppressed dust flux in ODP 659 and 721/722 during the Mid Pliocene results not only in an increased transitivity, but also a high assortativity. The latter feature can be explained by the fact that a relatively large cluster of state vectors representing this laminar regime emerges in the network, which is rather densely connected (Fig. 6c).
We conclude that the RN measures are not statistically independent in their time evolution (Table 2). For the ODP site 659 δ 18 O and dust flux records, the correlations between Fig. 11. Evolution of RN assortativity R for (A) the δ 18 O record from ODP site 659, and the dust flux records from ODP sites (B) 659, (C) 721, and (D) 967 during the Plio-Pleistocene. Parameters, significance test, and vertical lines are the same as in Fig. 9. of the two other RN properties assortativity R and diameter D. Since the latter one is closely related to the average path length L , the variability of both measures is very similar. Moreover, we also find some much weaker similarities between the temporal variability patterns of transitivity T and assortativity R, which are less pronounced, since both properties characterise not so obviously related aspects of the network geometry in phase space. Specifically, the time interval of suppressed dust flux in ODP 659 and 721/722 during the Mid Pliocene results not only in an increased transitivity, but also a high assortativity. The latter feature can be explained by the fact that a relatively large cluster of state vectors representing this laminar regime emerges in the network, which is rather densely connected (Fig. 6C).
We conclude that the RN measures are not statistically independent in their time evolution (Tab. 2). For the ODP site 659 δ 18 O and dust flux records, the correlations between transitivity T and assortativity R as well as between average path length L and diameter D as measured by Spearman's ρ are most pronounced, which is consistent with theoretical expectations . Correlations are more clearly developed between all four measures in case of the more highly sampled dust flux records from ODP sites 721/722 and 967 (Tab. 1). However, for all records the four transitivity T and assortativity R as well as between average path length L and diameter D as measured by Spearman's ρ are most pronounced, which is consistent with theoretical expectations . Correlations are more clearly developed between all four measures in case of the more highly sampled dust flux records from ODP sites 721/722 and 967 (Table 1). However, for all records the four measures can be considered sufficiently independent to justify including all of them for a broad and thorough nonlinear time series analysis of oxygen isotope and terrigenous dust flux variability.

Robustness of the results
To assure the reliability and robustness of our results, we systematically study their dependence on the relevant algorithmic parameters of our method, in particular, the widths of the recurrence window (W * ) and the detrending window (W * D ) as well as the embedding delay (τ * ). In Figs. 13-15, the results of the significance test are presented as contours at two prescribed significance levels obtained from the observed measure's quantiles with respect to the corresponding test distribution. Green contours represent the lower prescribed quantile (5 %), while black contours indicate the upper one (95 %). This implies that values of the measure under study enclosed by green contours can be considered as exceptionally low, while those lying within black contours are exceptionally large, recalling the interpretation of the applied null-model given in Sect. 2.6. It is, however, important to recognise that the null-hypothesis of stationarity has been tested pointwise, while physical significance requires the null-hypothesis to be rejected over a certain period of time, i.e. for several subsequent time points (Maraun et al., 2007). Therefore, certain line-like structures, particularly those seen in Fig. 15, are likely to reflect statistical fluctuations rather than physically significant dynamical transitions. In the following, we will only present the results for the ODP site 659 dust flux record.
i. Recurrence window size W * : as for the model systems in Sect. 3, we first discuss the sensitivity of our results to the changing width of the recurrence window W * . The corresponding results for the four chosen RN measures are shown in Fig. 13. We recognise that the most significant features persist under varying W * , although the relevant structures become broader and less significant for larger windows. This is to be expected since more and more data from time intervals not directly affected by the origin of specific network properties (e.g. a laminar phase in the dynamics) contribute to the longer windows. As the window width is increased linearly, conelike structures emerge (which is especially well visible for the Mid Pliocene transitivity maximum as the most  Fig. 13. Dependence of (A) transitivity T , (B) average path length L, (C) assortativity R, and (D) diameter D on the recurrence window size W * for the dust flux record from ODP site 659. The recurrence window step size is fixed to ∆W * = 41 kyr, the detrending window size to W * D = 500 kyr. Green and black contours correspond to the 5% and 95% quantiles with respect to the test distribution obtained from 10,000 realisations of our null-model. Other parameters were: embedding dimension m = 3 and delay τ * = 10 kyr, the threshold was chosen to yield a fixed edge density ρ = 0.05. The white bands indicating "no value" at the left and right margins of each panel appear because we plot the network measures g µ,W * at the mid-points t µ of the windows µ used for RN analysis (Sec. 2.4). As the mid-points of the first (last) window move further into the past (present) for increasing W * , the white bands grow linearly for linearly increasing W * . of time series analysis, the consideration of embedding with properly chosen parameters is necessary in order to obtain feasible results.
In contrast to other techniques, RN analysis does not characterise temporal interrelationships within the analysed records (although time information enters indirectly through embedding parameters, however, mostly on short time scales as typically mτ * N ∆T ), but quantifies geometric properties of the sampled dynamical system in its (reconstructed) phase space. The only implicit assumption is that the available sample of observed state vectors {y µ i } represents the spatial distribution of the true state vectors in the (properly reconstructed) phase space of the underlying dynamical system sufficiently well.
In this respect, our approach is very generally applicable and has comparably moderate requirements in terms of the requested number of data (i.e., windows with O(100) data points are sufficient for a reasonable analysis of nonstationary systems). In case of palaeoclimate records, this complementary way for characterising time series avoids conceptual problems of other approaches due to uncertain age models and non-uniform sampling. E.g., the results of RN analysis {g µ } are invariant to changes in the age model {t i }, only the associated windows' mid-points {t µ } change with variations in {t i } (Eq. (9)). However, the aforementioned problems indirectly persist in terms of the necessary embedding of the data and have to be finally resolved in corresponding future work. While the present work focussed on Fig. 13. Dependence of (A) transitivity T , (B) average path length L, (C) assortativity R, and (D) diameter D on the recurrence window size W * for the dust flux record from ODP site 659. The recurrence window step size is fixed to W * = 41 kyr, the detrending window size to W * D = 500 kyr. Green and black contours correspond to the 5 % and 95 % quantiles with respect to the test distribution obtained from 10 000 realisations of our null-model. Other parameters were: embedding dimension m = 3 and delay τ * = 10 kyr, the threshold was chosen to yield a fixed edge density ρ = 0.05. The white bands indicating "no value" at the left and right margins of each panel appear because we plot the network measures g µ,W * at the mid-points t µ of the windows µ used for RN analysis (Sect. 2.4). As the mid-points of the first (last) window move further into the past (present) for increasing W * , the white bands grow linearly for linearly increasing W * . relevant feature). In general, we observe that for our example the transitivity is most robust with respect to changes of W * , whereas the other network measures may lose significance if this parameter of our analysis method is varied. We note, however, that the periods of interest identified in Sect. 4.1 are robust for a wide range of recurrence window sizes, presenting a trade-off between good localisation of identified features (small windows) and reasonable statistical confidence of the calculated network properties (large windows).
ii. Detrending window size W * D : regarding the dependence of our observations on the choice of the detrending window, Fig. 14 shows that the general temporal variability pattern of the different network measures remains unchanged as W * D is altered, whereas the actual significance levels are more strongly influenced. In general, we can conclude, however, that the most significant time periods persist under variations of W * D , which is particularly well expressed for the transitivity during the Mid Pliocene. Together with the fact that RN analysis of the three dust flux records without detrending produces consistent results  this suggests that trends do not have a significant influence on the outcomes of RN analysis as long as W * N T . However, this should be checked in any particular application by comparing the results for the time series data before and after detrending. Note that the results for undetrended time series are approximated by those displayed in Fig. 14 for W * D ≈ N T , since RN analysis is invariant to nearly uniform translations of the data.
iii. Embedding delay τ * : our results are also seen to be robust with respect to reasonable variations of the embedding delay τ * around the previously chosen delay time τ * = 10 kyr (Fig. 15). However, for the Transitivity Average path length

Assortativity Diameter
Fig. 14. Dependence of (A) transitivity T , (B) average path length L, (C) assortativity R, and (D) diameter D on the detrending window size W * D for the dust flux record from ODP site 659. The recurrence window size is fixed to W * = 410 kyr with a step size of ∆W * = 41 kyr. Green and black contours correspond to the 5% and 95% quantiles with respect to the test distribution obtained from 10,000 realisations of our null-model. Other parameters were: embedding dimension m = 3 and delay τ * = 10 kyr, the threshold was chosen to yield a fixed edge density ρ = 0.05. In regions outside the black dashed lines the results are influenced by boundary effects, since the effective detrending window size WD(t) has to decrease towards the time series' limits (Eq. (3)). The white bands indicating "no value" at the left and right margins of each panel appear because we plot the network measures g µ,W * at the mid-points t µ of the windows µ used for RN analysis (Sec. 2.4). In contrast to Fig. 13 their width does not change as W * is fixed here. the technical aspects of applying RN analysis to palaeoclimate time series, an in-depth discussion of the results obtained for the three dust flux records in the light of additional proxy records and palaeontological evidence is given in Donges et al. (subm.). The recurrence window size is fixed to W * = 410 kyr with a step size of W * = 41 kyr. Green and black contours correspond to the 5 % and 95 % quantiles with respect to the test distribution obtained from 10 000 realisations of our null-model. Other parameters were: embedding dimension m = 3 and delay τ * = 10 kyr, the threshold was chosen to yield a fixed edge density ρ = 0.05. In regions outside the black dashed lines the results are influenced by boundary effects, since the effective detrending window size W D (t) has to decrease towards the time series' limits (Eq. 3). The white bands indicating "no value" at the left and right margins of each panel appear because we plot the network measures g µ,W * at the mid-points t µ of the windows µ used for RN analysis (Sect. 2.4). In contrast to Fig. 13 their width does not change as W * is fixed here. embedding delay exceeding τ * = 20 kyr the results and significance levels change considerably. This is expected as for delays larger than 20 kyr, autocorrelations in the time series do not decrease significantly anymore. In the case of the δ 18 O record from ODP site 659 they even increase again due to pronounced (obliquitydriven) Milankovich cycles with a period around 41 kyr (Fig. 3), so that the autocorrelation criterion for the choice of τ * does not apply here anymore.

Conclusions
We have demonstrated that RN analysis allows detecting dynamical transitions in non-stationary model systems as well as real-world palaeoclimate data. Transitivity and average path length have been previously discussed as appropriate network properties indicating qualitative changes in the dynamics of the underlying system. Here we have provided examples that also other global network measures such as assortativity and network diameter trace qualitative changes in dynamical systems, which, however, do not have a similarly straighforward interpretation in terms of basic system properties as the two other aforementioned quantities.
Our results show that the outcomes of RN analysis are quite robust if the fundamental parameters of the method (detrending and recurrence window sizes, embedding delay) are varied within a reasonable range. Unlike for other methods of time series analysis, the consideration of embedding with properly chosen parameters is necessary in order to obtain feasible results.
In contrast to other techniques, RN analysis does not characterise temporal interrelationships within the analysed records (although time information enters indirectly through embedding parameters, however, mostly on short time scales as typically mτ * N T ), but quantifies geometric Fig. 15. Dependence of (A) transitivity T , (B) average path length L, (C) assortativity R, and (D) diameter D on the embedding delay time τ * for the dust flux record from ODP site 659. The recurrence window size is fixed to W * = 410 kyr with a step size of W * = 41 kyr, the detrending window size to W * D = 500 kyr. Green and black contours correspond to the 5 % and 95 % quantiles with respect to the test distribution obtained from 10 000 realisations of our null-model. Vertical line-shaped contours are likely to correspond to statistical fluctuations rather than physically significant time intervals (see text). Other parameters were: embedding dimension m = 3, the threshold was chosen to yield a fixed edge density ρ = 0.05.
properties of the sampled dynamical system in its (reconstructed) phase space. The only implicit assumption is that the available sample of observed state vectors {y µ i } represents the spatial distribution of the true state vectors in the (properly reconstructed) phase space of the underlying dynamical system sufficiently well.
In this respect, our approach is very generally applicable and has comparably moderate requirements in terms of the requested number of data (i.e. windows with O(100) data points are sufficient for a reasonable analysis of nonstationary systems). In case of palaeoclimate records, this complementary way for characterising time series avoids conceptual problems of other approaches due to uncertain age models and non-uniform sampling. E.g. the results of RN analysis {g µ } are invariant to changes in the age model {t i }, only the associated windows' mid-points {t µ } change with variations in {t i } (Eq. 9). However, the aforementioned problems indirectly persist in terms of the necessary embedding of the data and have to be finally resolved in corresponding future work. While the present work focussed on the technical aspects of applying RN analysis to palaeoclimate time series, an in-depth discussion of the results obtained for the three dust flux records in the light of additional proxy records and palaeontological evidence is given in Donges et al. (2011).