Detecting dynamical anomalies in time series from different palaeoclimate proxy archives using windowed recurrence network analysis

Analysing palaeoclimate proxy time series using windowed recurrence network analysis (wRNA) has been shown to provide valuable information on past climate variability. In turn, it has also been found that the robustness of the obtained results differs among proxies from different palaeoclimate archives. To systematically test the suitability of wRNA for studying different types of palaeoclimate proxy time series, we use the framework of forward proxy modelling. For this, we create artificial input time series with different properties and compare the areawise significant anomalies detected using wRNA 5 of the input and the model output time series. Also, taking into account results for general filtering of different time series, we find that the variability of the network transitivity is altered for stochastic input time series while being rather robust for deterministic input. In terms of significant anomalies of the network transitivity, we observe that these anomalies may be missed by proxies from tree and lake archives after the non-linear filtering by the corresponding proxy system models. For proxies from speleothems, we additionally observe falsely identified significant anomalies that are not present in the input 10 time series. Finally, for proxies from ice cores, the wRNA results show the best correspondence with those for the input data. Our results contribute to improve the interpretation of windowed recurrence network analysis results obtained from real-world palaeoclimate time series.

Windowed recurrence network analysis (wRNA) has already been successfully used to detect dynamical anomalies in time series from different palaeoclimate archives (Donges et al., 2011a, b;Eroglu et al., 2016;Lekscha and Donner, 2018). But it has also been observed that this method often yields high numbers of false positive significant points, and, that not all palaeoclimate 25 archives give equally robust results. In order to improve the interpretation of results obtained with wRNA for real-world time series, we here systematically test the suitability of the approach for analysing palaeoclimate data from different types of archives by employing proxy system models.
Proxy system models are forward models that simulate the formation of a palaeoclimate proxy based on the systemic understanding of that proxy . That is, given the climate input variables, the proxy system model outputs a proxy 30 time series. We here use intermediate complexity models for tree ring width, branched glycerol dialkyl glycerol tetraethers in lake sediments, and the isotopic composition of speleothems and ice cores (Tolwinski-Ward et al., 2011;Dee et al., 2015Dee et al., , 2018. We first create artificial climate input time series with different properties. In particular, we consider two stochastic processes, Gaussian white noise and an autoregressive process of order 1, and two non-stationary time series from the paradigmatic Rössler and Lorenz systems. Additionally, we use climate input from the last millennium reanalysis project as an estimate of 35 realistic past climate variability (Hakim et al., 2016;Tardif et al., 2019). We then compare the input and the model output with respect to the properties of the time series and the results of wRNA. To quantify significant dynamical anomalies, we use an areawise significance test that was recently introduced for wRNA (Lekscha and Donner, 2019).
With this work, we want to contribute to a better understanding and an improved interpretation of results obtained with wRNA for palaeoclimate applications. We first introduce our analysis framework in Sec. 2, the proxy system models in Sec. 3,40 and the input time series in Sec. 4. Then, we show and discuss the results of the wRNA for the different input and model output time series in Sec. 5 before we present our main conclusions in Sec. 6. Additional information and figures are included in the Supplementary Material accompanying this paper.
2 Analysis framework 2.1 Phase space reconstruction 45 The first step when analysing data using recurrence based approaches is to reconstruct the higher-dimensional phase space of the system from the measured univariate time series x(t) with observations made at times {t i } N i=1 . For this, we use uniform time delay embedding (Takens, 1980;Packard et al., 1980) which relates the higher-dimensional coordinates of the system's phase space to delayed versions of the measured, univariate time series (1) space can give information on the original system dynamics. Still, this theorem only provides a sufficient condition for the 55 embedding dimension, while also smaller-dimensional embeddings may give useful information about the system's dynamics.
In practical applications, however, the box-counting dimension of the original attractor is usually unknown and the data are finite and subject to noise such that the choice of the embedding parameters plays an important role for the quality of the phase space reconstruction. In particular, the embedding dimension can be estimated using the method of false nearest neighbours (Kennel et al., 1992), while the first zero of the autocorrelation function or the first minimum of the mutual infor-60 mation can serve as an estimate for a suitable embedding delay (Abarbanel et al., 1993;Fraser and Swinney, 1986;Kantz and Schreiber, 2004). We here use an embedding dimension of m = 3 as a compromise between the number of available data points and the minimum number of independent coordinates as indicated by the false nearest neighbour criterion. For estimating the embedding delay, we use the first zero of the global lagged autocorrelation function.

65
We analyse the embedded time series x(t) by means of windowed recurrence network analysis (wRNA), using a sliding window approach and dividing the embedded time series into windows of width W with mutual offset dW . For each of those N = (N − (m − 1)τ − W )/dW windows, we construct a network from the time series by identifying the different as nodes of the network and drawing a link between two nodes x i and x j if they are mutually closer in phase space than some threshold . For this purpose, we measure the distances in phase space using the maximum norm, which has, due to its particularly simple form when calculated in Euclidean space, been widely used in recurrence-based analyses. Then, the resulting recurrence network can be described by its adjacency matrix A with entries where θ(·) is the Heaviside function and the delta function δ i,j excludes self-loops in the network. Here, we select the threshold 75 adaptively by choosing a fixed edge density ρ = 0.05. This means, for every window, we choose the threshold such that a fraction ρ of all possible links in the network are realised (corresponding to the 100ρ% mutually closest pairs of state vectors).
From the adjacency matrix, we can estimate various network properties. In the course of this work, we will restrict ourselves to the network transitivity 80 which gives the probability that two randomly chosen neighbours of a randomly chosen node are mutually connected. Network transitivity is particularly well suited for detecting dynamical anomalies in non-stationary time series, since this network measure has been shown to be closely related to the dimensionality of the system dynamics, in particular, when using the maximum norm to measure distances in phase space (Donner et al., 2010a(Donner et al., , 2011b. Specifically, the quantity log T / log(3/4) then corresponds to a generalised fractal dimension (Donner et al., 2011a the offset of the windowed analysis to be always dW = 1.

Significance testing and confidence levels
In order to identify dynamical anomalies from the resulting network transitivity, we first perform a pointwise significance test using random shuffling surrogates. That is, for every window width W , we randomly draw N s = 1, 000 times W embedded state vectors from x(t) and calculate the network transitivity of this set of W state vectors. We then take the empirical 2.5th 95 and 97.5th percentiles from the N s realisations to obtain a two-sided confidence interval of s pw = 95%. All transitivity values of the original analysis results that fall outside this interval are considered to show pointwise significant anomalies, that is, the null hypothesis of the corresponding values of transitivity resulting from random data with the same amplitude distribution as the original data is rejected. The result of this pointwise significance test can be summarised in the binary significance matrix S pw of the same size as the matrix of the transitivity results Q. The entries of S pw equal one if the corresponding value of the 100 transitivity has been found to be pointwise significant, and zero otherwise.
As intrinsic correlations of the analysis results in both, the time domain (due to the short offset dW ) and the window width domain can lead to patches of false positives in the pointwise significance matrix, we additionally perform an areawise significance test (Lekscha and Donner, 2019;Maraun et al., 2007). A pointwise significant point is defined to be areawise significant, if all neighbouring points within a given rectangle are also pointwise significant, i. e., if the pointwise significant 105 point lies within a patch of pointwise significant points that is larger than this rectangle. The side lengths of the rectangle depend on the intrinsic correlations of a chosen null model that are estimated as described in detail in Lekscha and Donner (2019).
Here, we employ a data-adaptive null model using iterative amplitude-adjusted Fourier transform surrogates of the original time series (Schreiber and Schmitz, 2000). That is, we test whether the same analysis results could have been obtained from 110 data with the same amplitude distribution and linear correlation structure as the original data but that are otherwise random, i. e., originate from linear stochastic processes with prescribed correlations.

Proxy system models
Forward modelling of palaeoclimate proxies offers the possibility to gain insights into the underlying processes that influence the sensitivity of a given proxy to climate variations and can thus be used to investigate characteristic properties of time series 115 of different palaeoclimate archives and their implications for further analyses. We here use four models for typical proxies from tree rings, lake sediments, ice cores and speleothems, respectively, in order to test how well dynamical anomalies can be identified when applying wRNA to time series originating from those archives.
Generally, a proxy system model can be divided into an environment, a sensor, an archive and an observation sub-model . The environment model can be used to model the environmental factors that the archive is sensitive to using the 120 climatic input variables. The sensor model then describes how the archive reacts to the environment, and the archive model accounts for how this reaction is written into the archive. Finally, an observation model can be used to simulate uncertainties in dating or measurements. Here, depending on the archive, we use different combinations of the environment, sensor and archive sub-models, while neglecting possible dating and measurement uncertainties. In the following, we will give a brief overview over the model structures and parameters. A full description of the models can be found in the corresponding references. . This is a reduced version of the full Vaganov-Shashkin model (Vaganov et al., 2006) and requires monthly input data of temperature T and 130 either precipitation P or soil moisture M . Additional model parameters are thresholds for temperature and soil moisture below which growth is not possible and above which growth is optimal (T 1 , M 1 , T 2 , M 2 ), the latitude of the site Φ, and integration start and end months I 0 and I f that set the period over which climate is responsible for growth in a given year.
If precipitation is given, the Leaky Bucket model (Huang et al., 1996) is used as environment model to calculate the soil moisture M (t) based on the water balance in soil. This model requires additional parameters such as the initial moisture 135 content of the soil M 0 , minimum and maximum soil moisture content M min,max , runoff parameters m, α and µ and the root depth d r . The sensor model for the tree ring width then basically relies on the principle of limiting factors (Fritts, 1976), that is, it assumes that tree ring growth is limited by the resource that is the scarcest for optimal growing conditions, i. e., in this case, the temperature or the soil moisture. A temperature-based growth response is calculated as 140 and similarly, a growth response g M (t) is calculated for soil moisture. In addition, a third insolation-based growth response g E (t) is calculated based on the mean of the normalised lengths of the day for each month. The total growth response g(t) of the tree to the climatic input is then given as the minimum of the temperature and moisture based growth responses modulated by the insolation based growth response, To obtain the annual growth response from those monthly data, the model integrates the growth response g(t) over those months that are specified as the start and end months I 0 and I f . Finally, the annually resolved time series of tree ring width anomalies is obtained by normalising the annual growth response to zero mean and unit variance.
It should be noted that this model does not take into account juvenile tree growth. Real tree ring width data are of course subject to juvenile growth and the effect is usually subtracted from the measured data. Problems arising from this are thus 150 disregarded in the model which we will further address when discussing the results for the model.
To set the model parameters to realistic values, we use an exemplary real-world data set of a local tree ring width index chronology from eastern Canada (54.2 • N, 70.3 • W) that was previously used for regional summer temperature reconstruction (Gennaretti et al., 2014). The quality of this data set with respect to data homogeneity, sample replication, growth coherence, chronology development, and climate signal has been ranked high (Esper et al., 2016). Regional average annual 155 temperature is about −3.8 • C with average maximum temperatures of 16 • C in July and minimum temperatures of −23 • C in January. The average annual precipitation sum is 693 mm with a monthly minimum of about 30 mm in February/March and a maximum of about 100 mm in September. This information has been used in order to choose the mean and standard deviations of the climatic input variables and their annual cycles. An overview of the climate input and also the model parameters is given in Table 1. To determine the threshold parameters, we used the Bayesian parameter estimation as suggested in Tolwinski-Ward

Lake sediments
Records from lake sediments are available from many regions worldwide and can provide information about past temperatures and precipitation, depending on the regional boundary conditions and measured proxy (Cohen, 2003). We here model branched glycerol dialkyl glycerol tetraethers (brGDGTs) from lacustrine archives that have been related to air temperatures by using 165 one of the sensor models provided in the PRYSM v2.0 framework (Dee et al., 2018). For this, as model input, a time series of mean annual air temperature T is required.
BrGDGTs are produced by bacteria and their degree of cyclisation and methylation has been related to soil temperatures, lake pH, and also to mean annual air temperatures (Weijers et al., 2007;De Jonge et al., 2014;Russell et al., 2018). The degree of methylation is quantified by a methylation of branched tetraether (MBT) index. The sensor model employs the MBT 5ME 170 index that only uses 5-methyl isomers. In particular, the calibration of Russell et al. (2018) is used, in which the mean annual air temperature is related to the MBT 5ME index via the equation The archive model then accounts for bioturbation and mixing of the sediments using the TURBO2 model (Trauth, 2013).
TURBO2 models the benthic mixing effects on individual sediment particles by assuming that in a mixed layer of a specified To tune the climatic input variables, we use the climatic setup corresponding to the one used for the tree ring archive in eastern Canada, while for the model parameters, we use typical values for lake sediments that are taken as default in the PRYSM implementation of the lake archive model (Dee et al., 2018). In particular, it should be noted that for the abundances 185 of the input species and the mixed layer thickness, we use constant values over time. The climatic and model input parameters are specified in Table 2.

Speleothems
Oxygen isotope fractions of speleothems have been shown to provide valuable insights into past climate variability (Wong and Breecker, 2015). We here model the isotopic composition of speleothems by using the speleothem model presented within 190 the PRYSM framework (Dee et al., 2015). This intermediate-complexity model is based on the model in Partin et al. (2013) and requires the mean annual temperature T and the mean of the precipitation-weighted annual isotopic composition of the precipitation δ 18 O P as input. Additionally, the ground water residence time τ gw has to be specified. The sensor model covers processes in the karst and the cave, while processes in the soil such as evapotranspiration are neglected. The model filters the δ 18 O P signal by applying an aquifer recharge model to simulate the storage and thus the mixing 195 of water of different ages in the karst. This process is parameterised by the mean transit time τ gw . The isotopic composition of the cave water is then given as the convolution ( * ) of δ 18 O w with the impulse response of the aquifer recharge model with the temperature T a being the decadal average of T that is calculated using a Butterworth filter (Zumbahlen, 2008).
The parameters for the speleothem δ 18 O model are tuned by using the data of stalagmite DA from Dongge cave (Wang et al., 2005) which is one of the most studied speleothem data sets. Dongge cave is located in Southern China ( There is no information available about the mean transit time, i. e., the time the water has spent inside the karst before entering the cave. We here choose to use an average transit time of five years, which is slightly larger than the average 210 sampling rate of the data (4.2 years). An overview of the input variables and model parameters is given in Table 3.

Ice cores
Proxy time series from ice cores have been used in a variety of contexts to study past climate variability on different time scales (Jouzel, 2013;Thompson et al., 2005). As for the speleothem δ 18 O, we use the model for ice core δ 18 O that is implemented and presented within the PRYSM framework (Dee et al., 2015). The model requires the precipitation-weighted mean T m , the altitude of the glacier z, the mean surface pressure p, the mean accumulation rate at the site A, and the total depth of the core h max which is given by the time span of the observations times the average accumulation rate.
The sensor model corrects the isotopic composition of the precipitation for the altitude of the glacier by using the relation 220 with a describing the altitude effect. The archive model then accounts for compaction and diffusion within the ice core. First, the density of the core has to be calculated as a function of the depth of the core. For this, an adapted version of the firn densification model by Herron and Langway is used (Herron and Langway, 1980). From the density and the mean temperature T m , we can then compute the diffusion length σ within the core as a function of the depth h (Johnsen et al., 2000). This, in turn, is used to calculate the final proxy time series δ 18 O d by convolving the isotopic signal of the ice δ 18 O ice with a Gaussian 225 kernel function with standard deviation equalling the diffusion length at a given depth, where the convolution is again denoted by the asterisk ( * ).
To tune the climatic input variables and model parameters of the ice core δ 18 O model, we use an exemplary real world data set of the Quelccaya ice cap  which is one of the most studied ice core data sets outside the polar  Table 4.
4 Input data 235 We now introduce the data sets that we use as input for the proxy system models. We first consider two stationary stochastic processes, namely Gaussian white noise (GWN) and an autoregressive process of order 1 (AR(1) process) to evaluate whether such input can lead to the detection of dynamical anomalies in the proxy time series. Then, we consider non-stationary versions of the two well-known Rössler (ROS) and Lorenz (LOR) systems. For all those processes, time series of length N = 1, 000 are series are normalised to zero mean and unit standard deviation, and for each model, the mean is adjusted to the corresponding climatic boundary conditions.

Gaussian white noise
For the case of GWN, we draw N data points independently at random from the probability distribution

250
where σ = 1 is the standard deviation and µ = 0 is the mean of the distribution. We do not expect to detect significant dynamical anomalies from time series created in this way as the process is stationary.

Autoregressive process of order 1
For the AR(1) process, we create a time series of length N by using the relation with t a Gaussian random variable with zero mean and constant standard deviation σ = 0.5. The scaling factor is given as α = 0.7 corresponding to the approximate value that we obtained when fitting an AR(1) process to the tree ring width data from eastern Canada (Gennaretti et al., 2014). As initial condition, we use x(0) = 0.3. The resulting time series is normalised to zero mean and unit standard deviation. As for GWN, we do not expect to detect significant dynamical anomalies in the corresponding time series.

Non-stationary Rössler system
In a next step, we use data from a non-stationary version of the Rössler system which exhibits non-trivial and rich cascades of bifurcations despite its rather simple attractor topology. The Rössler system is defined by the set of ordinary differential equations (ODEs) (Rössler, 1976) x It should be noted that this model is not meant to simulate real world climate dynamics.
We From the Feigenbaum diagram, it becomes clear that we expect to detect alternating periods of lower and higher dimensional 275 dynamics in the time series. In particular, we stress that we do not expect to detect the bifurcation points but periods of outstandingly high-or low-dimensional dynamics in between them as we use random shuffling surrogates for the pointwise significance test.

Non-stationary Lorenz system
The Lorenz system shows a more complicated attractor topology than the Rössler system and has been originally introduced 280 as a simple toy model for atmospheric convection processes (Lorenz, 1963). It is given by the following set of ODEs (Lorenz, 1963):  The stationary Lorenz system has been found to exhibit a shift from periodic to chaotic dynamics at b = 166.0 (Barrio and Serrano, 2007), while for the transient system as described above, transitions could be detected at b = 161.0 and b = 166.5 (Donges et al., 2011a). Thus, we expect to detect regimes of more periodic dynamics for b < 166.5 and of more chaotic dynamics for b ≥ 166.5 in terms of significantly high and low values of the network transitivity, respectively.

Last millennium reanalysis data 295
Finally, we consider reconstructed temperature and precipitation data of the years 501 − 2000 AD from the last millennium reanalysis project version 2 (Hakim et al., 2016;Tardif et al., 2019), which combines information from general circulation models and from proxy measurements using palaeoclimate data assimilation. As no information about the isotopic composition of the precipitation is available, we only consider the models for tree ring width and lake sediments with this input. From the available global gridded data, we use the standardised, thus, unit-less, ensemble average time series of temperature and 300 precipitation from the grid points with the coordinates closest to those for which we have calibrated the tree and the lake model, that is, at 54 • N, 70 • W.

Results
Before studying the different proxy system models, we note that applying some general filtering techniques such as moving average filtering or exponential smoothing to the described stochastic input time series on the one hand and the deterministic  the processing through the model causes the initially significant (false positive) points to be missed and creates an additional patch of areawise significant network transitivity in the case of the speleothem model. In the case of the AR(1) process (Fig. 2), we observe an areawise significant patch of anomalously high values of network transitivity for the input temperature that is not apparent in the tree ring width and lake sediment model output. For the isotope input, we do not find any areawise significant points, while the speleothem model shows two large and two small falsely 325 identified significant patches of network transitivity. That is, when recalling the results for the filtering of the AR(1) process, we find that the speleothem model seems to particularly involve a way of filtering that produces patches of areawise significant high values of the network transitivity while the other models rather change the variability of the stochastic input.
For the non-stationary Rössler system (Fig. 3) The input time series of the non-stationary Lorenz system (Fig. 4)   has often be attributed to more stable climate conditions while the LIA has often been associated with more variable climate conditions even though this imprint has varied locally and has been mainly discussed for Europe. Thus, given the theoretical relation between network transitivity and the dimensionality of the system's dynamics (Donner et al., 2011a), we tentatively conclude that the higher values of the network transitivity during the MCA and the lower values of the network transitivity during the LIA reflect a lower/higher dimensional dynamics of the system at this particular location in terms of the complexity 360 of temporal variations rather than just a change in variance. In terms of the time series properties for the different periods, we indeed additionally observe an increase in variance of the time series from the MCA to the LIA, which is very likely also reflected in the different recurrence networks and, thus, the resulting network transitivity. We note that this non-stationarity in variance along with the MCA-LIA transition in the European/North Atlantic sector has also been reported as being reflected in other non-linear characteristics, which have been previously interpreted as a hallmark of some dynamical anomaly (Franke 365 and Donner, 2017;Schleussner et al., 2015).
In order to quantify the effect of the proxy system models as non-linear filters of the input signal on the detection of areawise significant points, Table 5 displays the fraction of falsely identified and missed significant points in the different proxy models for the different input time series. We observe that for the tree ring width model, the lacustrine sediment model, and the ice core model missed significant points are more common than falsely identified significant points, with the exception of the Lorenz 370 input for the tree ring width model. For the speleothem model, both falsely identified and missed significant points occur. In proxies and/or archives from the same region to obtain further climatological knowledge from such kind of analysis (Donges et al., 2015a;Franke and Donner, 2017).

Discussion and conclusions
In this paper, we have studied the suitability of windowed recurrence network analysis (wRNA) for detecting dynamical 380 anomalies in time series from different proxy archives. For this, we used proxy system models that simulate the formation of proxy archives, such as for example tree rings, lake sediments, speleothems, and ice cores, given some climatic input variables like temperature and precipitation. We created artificial input time series with different properties and additionally used temperature and precipitation data from the last millennium reanalysis project (Hakim et al., 2016;Tardif et al., 2019). We then processed the input time series through the different proxy models and contrasted the results of wRNA for the different 385 input and model output time series.
We first compared the results of wRNA for stochastic and deterministic input to corresponding results when filtering the time series using moving average filtering and exponential smoothing (results not shown). We found that filtering alters the variability of the network transitivity with a bias towards additional and extended patches of areawise significant high transitivity values for stochastic input. The network transitivity of deterministic input seems to be rather robust under such filtering. When 390 processing the input time series through the different proxy system models, these differences for stochastic and deterministic input were also apparent. In terms of areawise significant anomalies of the network transitivity, we found that time series of tree ring width and brGDGTs in lake sediments have problems with missing areawise significant points, while the isotopic composition of speleothems also exhibits falsely identified significant points of high values of the network transitivity probably related to the bias towards higher values of the network transitivity due to the particular filtering in the model. Time series of 395 the isotopic composition of ice yield comparable results to the corresponding input but also sometimes miss significant points.
Taken together, our results show the need of further studying the effects of different filtering mechanisms on the results of the wRNA in order to interpret the results and draw reliable conclusions when analysing real-world data. This particularly concerns data from speleothems as they have been studied quite often using windowed recurrence analysis (e. g., Donges et al. (2015a); Eroglu et al. (2016)) and filtering effects seem to have a non-negligible influence on the results for this type of archive.

400
Thus, the role of the model parameters that control the filtering (here, the mean aquifer transit time) should be studied in more detail and in the context of the real-world applications, should be taken into account for estimates of these parameters for the particular data to better interpret corresponding results of the network transitivity.
Still, we want to stress that even then, providing a general recipe for interpreting the resulting network transitivity is hardly possible in the (palaeo)climate context. Climate-related interpretations always vary depending on the location and thus, local 405 boundary conditions have to be taken into account. As motivated in Section 2, the network transitivity has been related to the dynamical regularity of the variations in the analysed time series (e.g., Donges et al. (2011b)) with higher values of the network transitivity corresponding to less irregular variability and vice versa. This is in accordance with the interpretation of the network transitivity as an indicator of the dimensionality of the system's dynamics. In this regard, detected anomalies in the network transitivity could be related to some tipping point, but do not have to be.

410
Future work should also include the study of alternative proxy system models within this framework. Results of proxy system models for both, the same proxies (but with more detailed systemic understanding of the formation of the proxies, as for example a tree ring width model accounting for juvenile growth of the trees) and different proxy variables (in particular, for other lacustrine proxy variables), will complement the improved understanding of the suitability of wRNA for these types of time series and will advance the interpretation of the corresponding results. Also, sensitivity studies for the different model pa-415 rameters are of interest to better interpret results obtained with wRNA for a given real-world data set. This concerns particularly the mean aquifer transit time of the speleothem model.
For the stochastic input time series as for example the isotope input of GWN, we found some areawise significant artefacts in single realisations. To improve the reliability of the results for the these processes, more realisations should be considered to confirm the results and to exclude the influence of random artefacts. As we applied an areawise significance test to identify of non-linear observability might give an interesting new perspective on this as the filtering can be seen as creating a new observable, and the choice of observable has already been shown to influence results of recurrence quantification analysis and recurrence network analysis (Portes et al., 2014(Portes et al., , 2019. Moreover, further systematically studying the relation between 430 the autocorrelation of a time series and the resulting network properties might yield additional information on the role of the different archives for the results of wRNA. Code availability. Exemplary Python code for the windowed recurrence network analysis including the areawise significance test is available at https://gitlab.pik-potsdam.de/lekscha/awsig. A comprehensive implementation of recurrence network analysis can be found in the opensource Python package pyunicorn (Donges et al., 2015b), which can be found at https://github.com/pik-copan/pyunicorn.