The identification of cycles in periodic signals is a ubiquitous problem in time series analysis. Many real-world datasets only record a signal as a series of discrete events or symbols. In some cases, only a sequence of (non-equidistant) times can be assessed. Many of these signals are furthermore corrupted by noise and offer a limited number of samples, e.g., cardiac signals, astronomical light curves, stock market data, or extreme weather events. We propose a novel method that provides a power spectral estimate for discrete data. The edit distance is a distance measure that allows us to quantify similarities between non-equidistant event sequences of unequal lengths. However, its potential to quantify the frequency content of discrete signals has so far remained unexplored. We define a measure of serial dependence based on the edit distance, which can be transformed into a power spectral estimate (EDSPEC), analogous to the Wiener–Khinchin theorem for continuous signals. The proposed method is applied to a variety of discrete paradigmatic signals representing random, correlated, chaotic, and periodic occurrences of events. It is effective at detecting periodic cycles even in the presence of noise and for short event series. Finally, we apply the EDSPEC method to a novel catalog of European atmospheric rivers (ARs). ARs are narrow filaments of extensive water vapor transport in the lower troposphere and can cause hazardous extreme precipitation events. Using the EDSPEC method, we conduct the first spectral analysis of European ARs, uncovering seasonal and multi-annual cycles along different spatial domains. The proposed method opens new research avenues in studying of periodic discrete signals in complex real-world systems.

Many dynamical processes exhibit characteristic periodic time scales, which can be assessed by power spectral density estimation. This is the usual and standard procedure of time series analysis, where the data are, in general, equally sampled and follows a more or less continuous nature. However, in specific applications, only events are observable, which could be irregularly occurring activities or extreme events. Identification of cyclical behavior in such kind of data is difficult. We propose a novel but simple tool to estimate the power spectral density in such discrete data. This method can be applied on symbolic sequences, discrete data, or (extreme) event series.

## I. INTRODUCTION

The estimation of the power spectral density (PSD) of cyclical signals is an important field in signal and data analysis.^{1–4} However, the estimation of PSD from discrete data, such as symbolic time series or event series, is difficult. A time series can generally be denoted by a set of ordered pairs ${( t i, x i)}$ of time $ t i$ with $ t i + 1> t i$ and corresponding data value $ x i$; and with sampling index $i$. PSD estimation becomes even more challenging, if instead only the time points of the events are available (i.e., in contrast to the time series definition above, we have only a set of time points ${ t i}$ which indicate the presence of an event; time series analysis tools usually handle only data in the form ${( t i, x i)}$).

Analyzing discrete or event series is of general interest in many scientific fields. For example, heart beats are discrete events, described as the R-waves in an electrocardiogram (ECG). Pathologic states cause significant changes in the beat-to-beat intervals, which can be used to diagnose or even predict life-threatening cardiac events.^{5,6} Specific methods designed for discrete time series analysis are available for classifying cardiological signals, based on symbolic dynamics and ordinal patterns.^{7,8} Brain activity is controlled by the firing of neurons, occurring as sequences of events, also called “spike trains.”^{9} The synchronization between firing neurons is the base of transmitting information through the nervous system but can also cause serious pathological states, e.g., seizures.^{10,11} Symbolic dynamics is a promising tool to study seizure activity using electroencephalographic measurement.^{12} In climatology, a profound understanding of extreme weather events such as droughts, floods, or heavy precipitation is becoming more crucial.^{13,14} Detrimental societal and technological impacts of single events and cascades of compound weather extremes are projected to increase on a warming planet. In this context, atmospheric rivers (ARs) receive increasing interest as they can trigger extreme rainfall.^{15} Frequency and intensity of AR events are controlled by the Madden–Julian oscillation, the Quasi-Biennial Oscillation, the El Niño/Southern Oscillation (ENSO), and the Pacific Decadal Oscillation.^{16,17} Extreme rainfall and river flood event series are studied with event synchronization, edit distance, and complex networks to shed light on their serial dependency and spatial propagation.^{18,19} Discrete data also appear in engineering, e.g., in the investigation of how weather conditions affect car crashes, using integer-valued autoregressive models.^{20} In the analysis of stock markets, identification of regularities in the occurrence of trades in the order book or (ultrafast) extreme events is a common research objective.^{21,22} The temporal structure, heterogeneity, and serial dependencies in financial data are studied using response functions, spectral analysis, and recurrence plots.^{23–25} Communication between human individuals via social media platforms is intricatedly linked to recurrence of periodic social events and is, for instance, studied by means of periodic patterns in tweets on Twitter.^{26} Finally, detection of periodicities in irregular and event-like time series is also one of the main challenges faced in computational astrophysics.^{27,28}

Although discrete or (extreme) event time series are ubiquitous in many scientific disciplines, the study of cyclicities within these events is, to our best knowledge, still in its infancy. Potential methods would be Walsh transform^{29–31} and spectral analysis using the Haar wavelet.^{32–34} The Walsh transform provides a specific set of discrete basis functions consisting of binary sequences of different shape and length. Therefore, the interpretation of a plain Walsh spectrum is not straight-forward. Nevertheless, Walsh transforms are recently receiving more interest, e.g., in combination with machine learning tools for classification purposes.^{35} Applications of the Haar wavelet to non-equidistant and event-like data have so far remained sparse.^{36} While not providing a PSD-like framework, another family of methods based on line-search algorithms, phase folding, and maximum-likelihood estimation aims at finding periods in event series.^{37} Finally, some methods characterize event time series based on inter-event times.^{7,38}

The Haar wavelet and Walsh power spectrum require time series in the form of ordered pairs ${( t i, x i)}$. The Walsh power spectrum even needs equidistant time series (i.e., $ t i + 1\u2212 t i=const.$). However, series of (extreme) events, e.g., extracted from data by thresholding and further considering only the times of these events $ t i$, consist only of the set of these time stamps ${ t i}$. Alternatively, if assigning corresponding values to get an (extreme) event time series ${( t i, x i)}$, it would be one which is strongly non-equidistant (i.e., $ t i + 1\u2212 t i\u2260const.$) or, as a further alternative, to ensure equidistant sampling, it would consist mainly of zero values. Most methods of time series analysis are limited for such kind of data. Specific methods considering the extreme nature of discrete data points have their roots in neuroscience.^{39–41} The main purpose of these methods is testing coincidence or synchronization between event series. Among these methods are event synchronization,^{40} edit distance,^{39} and coincidence analysis.^{42} Event synchronization has been applied in climate science to construct complex networks to investigate the spatiotemporal dynamics of extreme rainfall.^{18,43} Edit distance has been applied in different applications and in combination with recurrence analysis, to study serial dependence^{19,24,44} and to allow analyzing non-equidistantly (palaeoclimate) time series.^{45,46}

We propose a novel power spectral analysis for discrete and event series ${ t i}$ based on the edit distance metric (modified Levenshtein metric^{39,47}). This method allows us to estimate a power spectrum directly from the event sequence without interpolation.

In Sec. II, we introduce the edit distance and the power spectral estimation based on this metric. We illustrate this approach in Sec. III on selected paradigmatic examples of random, periodic, and chaotic event sequences. Finally, we apply it on atmospheric rivers in Europe to find typical recurring cycles in Sec. IV.

## II. METHOD

For the estimation of a power spectrum, we propose a combination of a similarity measure for discrete time series (applied as an auto-correlation measure of this discrete time series) and the Wiener–Khinchin theorem. We use the edit distance metric as a similarity measure. Instead of the standard definition of a time series as ordered pairs of time and value ${( t i, x i)}$, we consider to have only the time stamps of the discrete or extreme events ${ t i}$, finally, forming a sequence $ S={ t 1, t 2,\u2026, t N}$ of $N$ events.

### A. Edit distance

The edit distance quantifies the dissimilarity between two sequences $ S a$ and $ S b$ of (discrete) symbols by means of how much it would minimally cost to transform $ S a$ into $ S b$. It was initially introduced in computer science (Levenshtein metric) where the two sequences are usually identified with two strings or words.^{47} The operations that are typically allowed include inserting a symbol, deleting a symbol, and shifting two symbols (the Levenshtein metric would consider addition, deletion, and replacement). The edit distance offers the advantage that it neither requires sequences of equal lengths nor equally displaced symbols. In practice, measuring similarities with the edit distance is, thus, not limited to strings, but it has been applied to marked point processes,^{24,39} extreme events,^{19} spatial trajectories,^{48} and irregularly sampled time series.^{46,49} Recent modifications to the edit distance include accounting for a saturation of shifting costs given a large temporal displacement between events^{19} and a correction method for non-stationary sampling rates.^{50} Another recent study already hints at the potential of the edit distance to reveal scale-dependent variability in the context of recurrence analysis.^{51}

Here, we consider sequences of binary symbols with no continuous amplitude. These could, e.g., represent spike trains or series of time stamps of the occurrence of extreme events. Hence, we refer to the symbols of a sequence as events in the following, with sequence $ S$ of $N$ events.

^{49}

### B. Power spectral estimate

*continuous*-amplitude) time series ${( t i, x i)}$, the power spectrum is given by the squared absolute value of the signal’s Fourier transform

A technical note needs to be added to the computation of this power spectral estimate: in the analysis of binary event series, events are usually recorded just by their time instance, time stamp ${ t i}$. It is, therefore, sufficient to add the time delay $\tau $ to all time instances in order to shift the event series against itself. This, however, results in non-overlapping segments at the beginning and end of the event series when computing $ R S edit(\tau )$ for some specific time delay $\tau $. In order to not include the transformation costs that solely result from the transformation of the non-overlapping events in the numerical computation, we prune the shifted event series to their largest intersection. The selection of the considered maximum time delay will change the lowest assessible frequency resolved by the power spectral estimate. While a large maximum value of $\tau $ allows for studying longer cycles, the length of the overlap shrinks and introduces undesired effects such as spectral leakage. We thus recommend to carefully select this parameter based on the number of observed events and research question at hand.

## III. PARADIGMATIC EXAMPLES

We illustrate the proposed EDSPEC by power spectral analysis of selected paradigmatic example systems (Fig. 1). In particular, we consider sequences with events of random [Fig. 1(a)], periodic [Fig. 1(b)], stochastic with serial dependence [Fig. 1(c)], and chaotic [Fig. 1(d)] occurrences. For each example, we generate $N=1000$ events over a time span of $T=10000$ arbitrary time units (a.u.). For the examples studied here, the considered maximum time delay $\tau $ is chosen as $\tau =T/10$. For evaluating the significance of the potentially detected periodicities, we perform a surrogate based statistical test. We generate $n=100$ surrogate event series by drawing $N$ time instances between 0 and $T$ with equal probability. The corresponding null-hypothesis is the absence of any non-stochastic/periodic pattern that controls the sequence of events. The $95%$-quantile of spectral power for each frequency from the surrogate ensemble is then used as the respective $95%$-confidence level.

### A. Random sequence

First, we generate a random sequence of events by drawing exponentially distributed waiting times between consecutive events [Fig. 1(a)]. The resulting process is a Poisson process and is, e.g., used as a Markovian model for arrival of costumers in a store.^{52} The EDACF rapidly declines to a value close to zero, reflecting vanishing serial dependence as it is expected for a Poisson process [Fig. 1(a)]. The EDSPEC follows a power law over two orders of magnitude. No significant frequencies exceed the 95%-confidence interval for the studied Poisson process, confirming the absence of any periodicities in this event series.

### B. Periodic sequence

The opposite is expected for a periodic event series. We generate an event series of 1000 equally spaced events with a frequency of $\omega =1/10$. The resulting EDACF is strictly periodic [Fig. 1(b)]. Consequently, the EDSPEC gives the expected result of a sharp peak at a frequency of $f=\omega $, exceeding the 95%-confidence interval. The observed spectral leakage is similar to what is found with standard power spectral density when it is computed for a monochromatic sinusoidal with a rectangular window function. Aliases of the original frequency can be identified at (positive) integer multiples of $\omega $. Spectral power below $\omega $ is constant.

### C. Stochastic sequence with serial dependence

For the first 100 delays, the EDACF declines to a value slightly below 0.5. On top of this decline, we observe small oscillations. The corresponding EDSPEC exhibits higher spectral power than expected from a random null-model for the lowest frequencies and for a band of higher frequencies around $f=1/10$. In the scaling region up to $f=1/20$, scaling is more steep than for the random event series. The serial dependence in the serially correlated event series does consequently not result in sharp significant frequencies (as expected) but gives rise to significantly high spectral power in a broader but limited band around the average time interval between consecutive events. For event series with a unimodal distribution of inter-event times, this average sampling interval could effectively be understood as a cycle. This is due to the fact that inter-events times center around this preferential interval, in contrast to the exponentially distributed inter-event times in Fig. 1(a).

### D. Chaotic sequence

^{53}

### E. Influence of noise and number of events

In the analysis of real-world event series, it is common to find considerable uncertainty in the timing of events. This uncertainty can, for example, be represented by jitter around the true time instance of an event. The impact of this effect on the estimation of the proposed edit distance-based power spectrum can be investigated by superimposing observational noise onto an event series. As the proposed method is predominantly designed to identify periodicities in event series, we generate a periodic event series with two distinct frequencies $ \omega 1= 1 20$ and $ \omega 2= 1 12.7$ (i.e., frequencies are not integer multiples of each other). Another challenge in the study of observational event series is that only a low number of events can be recorded, e.g., due to constrains in the measurement procedures. In order to test the sensitivity of our method against the presence of observational noise and shortness of event series, we generate two samples of this periodic event series, one with 130 and the other with only 20 events between 0 and $T=1000$ a.u. Due to the short event series lengths, we set a maximum time delay of $\tau =T/2$ in order to still assess lower frequencies. To establish a comparison with a typical alternative method, we also compute the regular power spectral density (PSD) of the event series using the Fourier periodogram of the binary event series (Fig. 2). Overall, the proposed EDSPEC generates more distinct peaks in presence of noise and deteriorates less rapidly when the number of events is reduced.

In the absence of noise, both the regular PSD and the EDSPEC effectively identify both frequencies as distinct peaks [Figs. 2(a)/2(b), black curves]. Additional peaks can be found at integer multiples of each frequency as well as at integer multiples of the sum and difference of them. If we crop the event series such that there are only $N=20$ events, both peaks are still visible but become harder to delineate [Figs. 2(c)/2(d), black curves].

Next, we superimpose normally distributed noise $ N(0,\sigma \Delta )$ onto the event series with noise strength $\sigma $ the average time interval $\Delta $ between consecutive events. It appears that the EDSPEC gives visually more convincing peaks at the expected frequencies for relatively high noise strength $\sigma =0.5$ (Fig. 2, orange curves). To systematically study the noise sensitivity of our method, we increase $\sigma \u2208[0,1]$ and obtain the regular PSD and EDSPEC for 500 distinct noise realizations that are superimposed onto the periodic event series. As a quantitative measure of peak quality, we compute the prominence of each of the two peaks by computing the vertical distance between the peak height and the height of its baseline. The height of each baseline is estimated by identifying the closest (local) minimum value of spectral power in proximity of the peak. Consequently, we obtain two values of peak prominence for each noise strength ( $ \omega 1$ and $ \omega 2$). We normalize both peak prominence curves separately for the two methods to their respective value in the noise-free case ( $\sigma =0$). We do the same for the shorter event series.

The identification of both frequencies $ \omega 1, \omega 2$ deteriorates considerably more slowly for the proposed EDSPEC than for the standard periodogram (Fig. 3). The prominence of the peak that is associated with the lower frequency $ \omega 1$ is on average more easy to detect than the one for $ \omega 2$. Interestingly, small amounts of noise (up to $30%$) seem to enhance prominence of the peak associated with $ \omega 1$ when we use the EDSPEC method. The validity of this finding will, however, most likely depend on the selected measure of peak prominence. In the study of very short event series ( $N=20$, black dotted lines), peak prominence for detecting both frequencies is reduced by about $1/5$ for the EDSPEC method, i.e., roughly proportional to the reduction in length as compared to $N=130$ (black solid line). While this reduction is certainly significant and also visually hampers detection of peaks especially in the presence of noise [Fig. 2(d)], the traditional PSD reveals considerably stronger deterioration in peak prominence for very short event series (Fig. 3, blue dotted line). Finally, distinct noise realizations can allow considerably better or worse peak detection, as indicated by the shown peak prominence curves for single realizations (gray lines).

Overall, we find that the proposed EDSPEC is less sensitive to timing jitter/observational noise than the regular PSD when applied to a simple periodic event series with two frequencies. The prominence of peaks deteriorates less strongly when the number of events is reduced. This should be of high relevance for sparse and noisy event series obtained from observational data. The robustness of other important properties beyond peak detection, such as spectral scaling, should be investigated in future studies.

## IV. APPLICATION TO EUROPEAN ATMOSPHERIC RIVERS

Atmospheric rivers (ARs) are narrow filaments of extensive water vapor transport in the lower troposphere. They can transport moisture for thousands of kilometers, predominately from the tropics toward the mid-latitudes.^{15} Their life times can range between a few hours and several days. Contrasting common beliefs, ARs represent the largest “rivers” on the Earth as they on average transport more than twice the amount of freshwater of the Amazon River.^{54} While some recent studies hint at alternative possible origins of ARs, most ARs are associated with the front of an extratropical cyclone.^{55} Low-intensity ARs provide vital amounts of freshwater and rarely result in heavy precipitation events. On the other hand, high-level ARs can cause detrimental impacts, including floods and high economic damage, when they landfall or ascent along an orographic barrier.

In Europe, links between ARs and heavy precipitation events have recently started to attract attention^{56,57} (similar to research of ARs at the American west coast^{58}). The frequency of ARs and related heavy precipitation events in Europe is projected to increase due to enhanced atmospheric moisture caused by recent climate change, whereas considerable uncertainties exist.^{59} Impact by high-category ARs is regionally heterogeneous over Europe and likely linked to regional presence of orographic barriers. Fall and winter have been identified as the seasons where most strong ARs occur, indicating the presence of a seasonal cycle in AR occurrence.^{56} Furthermore, AR frequency in northern and southern Europe has been found to correlate with the phase of the North Atlantic Oscillation (NAO).^{56}

Here, we study periodic cycles in the occurrence of European ARs between 1979 and 2019 by means of the proposed power spectral measure. To this extent, we track European ARs and regard the final (continental) position along their track as an event, as the termination of an AR potentially triggers rainfall (“landfalling” AR). Detection of ARs is based on localizing anomalous atmospheric transport of moisture. Many approaches track ARs by defining a threshold that identifies local anomalies in integrated vapor transport (IVT). Such approaches effectively presume stationary atmospheric moisture levels and often exclude low-level ARs. Here, we instead employ an AR-detection framework (*ARtracks*^{60}) based on global ERA5 reanalysis data that utilizes image-processing techniques (using the IPART algorithm^{61}). We implement the most widely used scale that characterizes European ARs based on their strength and persistence^{57} to distinguish five categories of ARs. While ARs of categories 1 and 2 are predominately beneficial, ARs of categories 4 and 5 often entail hazardous heavy precipitation events; category-3 ARs rank in between.

First, we generate the European AR catalog using the *ARtracks* algorithm. The used parameters can be found in the documentation of *ARtracks* and align with standard calibration of the IPART algorithm.^{61} As an intermediate result, we obtain a catalog of (uncategorized) ARs along with detailed AR properties, including unique AR identifiers and AR strength. To generate the AR catalog, ERA5 reanalysis dataset is re-gridded to a temporal resolution of six hours and a spatial resolution of $ 0.75 \xb0\xd7 0.75 \xb0$. Next, we apply a set of additional selection criteria to this AR catalog and categorize every landfalling AR. A more detailed description can be found in the Appendix. Based on the resulting filtered and categorized AR catalog, $65.5%$ of days between 1979 and 2019 have at least one active AR in the covered spatial region [Figs. 4 and 5(b)]. We identify various hot spots of landfalling ARs along the European west coast, such as the Iberian Peninsula, the UK, and Scandinavia (Fig. 4). In the Atlantic, Greenland, as well as some islands attract landfalling ARs. On the other hand, we also identify landfalling ARs along the Mediterranean coast of Italy.

For $7.1%$ of all days, there is a landfalling AR which we identify with an AR event. For by far the most ARs that are studied here, the final landfall location is mostly covering the European continent (see Appendix). Only few AR tracks terminate predominately over Africa, Asia, or North America (but still exhibit a European contribution). Most European ARs are potentially beneficial as they are assigned to category 1 or 2 (see Appendix).

We generate event series from the filtered and categorized AR catalog by both temporally and spatially integrating ARs. ARs are integrated onto a monthly time axis by selecting the highest category AR for each month. For the spatial integration, we largely follow the spatial domains oriented along three reference meridians as defined in Ref. 62, we distinguish between the (1) Iberian Peninsula, (2) west coast of France, (3) the UK, (4) southern Scandinavia and the Netherlands, and (5) northern Scandinavia [Fig. 5(b)]. We only extend the latter domain (northern Scandinavia, longitudes: $ 4.25 \xb0$ to $ 14.25 \xb0$) to include more AR events for sufficient spectral estimates. Each monthly and domain-specific time series of AR categories can be regarded as a symbolic sequence with integer values (symbols) between 1 and 5 [Fig. 5(a)]. Most (integrated) landfalling AR events are identified along the Iberian Peninsula and northern Scandinavia (see Appendix). In particular, the studied domains of France, UK, and southern Scandinavia exhibit an overproportionally high number of category-5 ARs, whereas northern Scandinavia exhibits the most regular distribution of AR categories.

For the spectral analysis, we only distinguish between low-category (category $<3$) and (moderate to) high-category (category $\u22653$) ARs. For each of these two groups, we obtain a single event series, indicating the occurrence of a low- or high-category AR at times ${ t i}$. We chose the considered maximum time delay as $\tau =30$ years to include decadal frequencies of AR occurrence.

We find evidence for a clear seasonal cycle in the occurrence of low-category ARs in all domains except the UK, corroborating findings from other studies on European ARs^{56,57} [Figs. 8(a)–8(e)]. This showcases the ability of the proposed method to identify known cycles in complex observational discrete data. While we do not analyze the phases of the identified cycles, it is likely that they align with the known enhanced frequency of AR occurrence in boreal autumn/winter. In the UK domain, the (significant) seasonal cycle is still noticeable, but it is suppressed compared to the other domains. This suggests a reduced seasonal predictability of ARs along parts of the British coast. Tested against random event series ( $99.99%$ confidence level from 1000 surrogates), broad bands of frequencies stand out as significant. This hints at a certain degree of predictability of ARs based on non-random cyclicities that could be exploited in forecasting approaches.

Even for high-category ARs, a broad range of significant cycles exist [Figs. 8(\,f\,)–8(j)]. High-category ARs only exhibit a pronounced seasonal cycle along the Iberian peninsula and northern Scandinavia. Considering both low- and high-category ARs, seasonality is expressed least between $ 0 \xb0$ and $ 5 \xb0$W longitude (France and UK). For high-category ARs, spectral power at lower frequencies is enhanced as compared to low-category ARs. These multi-annual/decadal cycles could be linked to the NAO: A negative (positive) NAO is likely linked to enhanced AR frequency in southern (northern) Europe.^{56}

## V. CONCLUSION

In many disciplines, research focuses on the identification of cycles in periodic discrete data (e.g., symbolic time series). In particularly challenging cases, only a sequence of (non-equidistant) times is recorded. Standard techniques, such as the Fourier PSD are not designed to quantify the spectral content of such signals. Noise contamination and limited signal length are additional obstacles in the spectral analysis of such systems.

In this work, we proposed a novel power spectral estimate for discrete data. It is based on the edit distance, a powerful distance measure that quantifies the similarity of two non-equidistant event sequences of unequal lengths by means of their transformation costs. The fundamental idea proposed here is to, first, normalize the edit distance as an auto-correlation measure of a discrete time series and, second, transform this measure into a power spectral estimate using the Fourier transform (analogously to the Wiener–Khinchin theorem). The proposed spectrum is referred to as EDSPEC. It enables us to identify cycles in non-equidistant discrete ordered pairs of time and values ${( t i, x i)}$ as well as for the more challenging case when only time stamps of the discrete events ${ t i}$ can be assessed. Here, we focus on the latter case.

We applied the proposed EDSPEC to a range of discrete paradigmatic signals, both of stochastic and deterministic origin. A simple hypothesis test based on random event surrogates is provided to test for the presence of significant cycles. The EDSPEC proves effective in identifying known cycles and broad bands of significant frequencies, also in the presence of relatively high noise intensity (or timing jitter). Even for a very low number of events, the EDSPEC performs comparably well.

In a challenging real-world application, we used the EDSPEC to investigate potential cycles in a novel catalog of European ARs. With the increasing moisture capacity of a warming atmosphere, ARs are expected to become a more and more vital source of atmospheric water supply, including heavy rainfall events. While seasonality in European ARs is undisputed and has been inferred from other AR catalogs (e.g., based on composites), our study confirms the presence of a seasonal cycle in AR occurrence based on our more suitable spectral analysis framework. This warrants further credibility to (spatially heterogeneous) seasonality in the landfall of European ARs. Furthermore, the proposed framework allows to separately study the presence of seasonality as well as multi-annual cycles for different AR categories. We found that hazardous high-category ARs show less seasonal regularity but offer predictability based on a broad range of significant cycles, potentially linked to the phase of the NAO.

In future studies, additional features of the proposed EDSPEC should be studied, e.g., spectral scaling and leakage effects. A multitude of complex discrete signals from real-world systems calls for application of the proposed measure, potentially revealing unknown periodicities.

## ACKNOWLEDGMENTS

This study was supported by BMBF Project climXtreme (No. 01LP1902J) “Spatial synchronization patterns of heavy precipitation events.” We thank Sara M. Vallejo-Bernal for fruitful discussions.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Norbert Marwan:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Tobias Braun:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Code used for the creation of event series and EDSPEC calculation of the paradigmatic examples, as well as code for the EDSPEC calculation of the ARs is available in Zenodo at 10.5281/zenodo.7555049. AR events are extracted from the ARtracks Atmospheric River Catalog available in Zenodo at https://doi.org/10.5281/zenodo.7018725.

### APPENDIX: SPECTRAL ANALYSIS OF AR CATALOG

We filter the AR tracks as returned by *ARtracks* with respect to several criteria in order to ensure the obtained European ARs are physically meaningful with respect to our research objective. First, the information on which proportion of the AR is located over Europe is used to ensure that only ARs are considered that are located over Europe for most track instances [Fig. 4(a)]. This still preserves ARs that overall have mixed continental contributions, e.g., over Africa. Next, an AR is fully discarded if it (a) exhibits IVT less than 250 kg/ms for any instance along its track, (b) does not persist longer than 6 h, (c) terminates over the ocean, and (d) is solely continental. We further only consider AR tracks that can be identified with a continuous track, i.e., no time 6-hourly instance is missing. Track instances where an AR is fully located over a different continent than Europe are discarded (while not discarding the full track). Finally, we categorize ARs adhering to the scale proposed in Ref. 57. If an AR exceeds the scale, either in terms of persistence of strength, we assign the most suitable category at the respective edge of the scale (e.g., an AR with maximum IVT magnitude between 750 and 1000 kg/ms and persistence of more than 72 h is still ranked as category 4). Categories are computed exclusively from all instances of a track where the AR landfalls, i.e., some fraction of its area is located over Europe. For the strength-based classification of ARs, we choose an AR’s category based on which category would be correct for its strength along most instances of its track. The resulting spatial distribution of ARs (based on the final instances of their centroids) is strongly oriented along the European west coast (Fig. 4). There are various hot spots of landfalling ARs along the European west coast, such as the Iberian Peninsula, the UK, and Scandinavia. In the Atlantic, Greenland as well as some islands attract landfalling ARs. On the other hand, we also identify landfalling ARs along the Mediterranean coast of Italy. For by far the most ARs that are studied here, the final landfall location is mostly covering the European continent (Fig. 6, gray bars). Only few AR tracks terminate predominately over Africa, Asia, or North America (but still exhibit a European contribution). Most European ARs are potentially beneficial as they are assigned to category 1 or 2 (Fig. 6, colored bars).

A more detailed assessment of the robustness of the obtained spatial distribution and statistics of AR properties beyond the findings mentioned here is investigated in a study that is currently prepared.

## REFERENCES

*Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing*(Association for Computational Linguistics, 2013), pp. 977–988.