Weather persistence on sub-seasonal to seasonal timescales: a methodological review

Persistence is an important concept in meteorology. It refers to surface weather or the atmospheric circulation either remaining in approximately the same state (quasi-stationarity) or repeatedly occupying the same state (recurrence) over some prolonged period of time. Persistence can be found at many different timescales; however, sub-seasonal to seasonal (S2S) timescales are especially relevant in terms of impacts and atmospheric predictability. For these reasons, S2S persistence has been attracting increasing attention from the scientific community. The dynamics responsible for persistence and their potential evolution under climate change are a notable focus of active research. However, one important challenge facing the community is how to define persistence from both a qualitative and quantitative perspective. Despite a general agreement on the concept, many different definitions and perspectives have been proposed over the years, among which it is not always easy to find one's way. The purpose of this review is to present and discuss existing concepts of weather persistence, associated methodologies and physical interpretations. In particular, we call attention to the fact that persistence can be defined as a global or as a local property of a system, with important implications in terms of methods and impacts. We also highlight the importance of timescale and similarity metric selection and illustrate some of the concepts using the example of summertime atmospheric circulation over western Europe.

1 Introduction

Surface weather persistence at sub-seasonal to seasonal (S2S) timescales can have severe impacts on human and natural systems. Long-lasting dry conditions, for instance, can lead to droughts and wildfires and can affect agriculture and energy production. Long-lasting wet spells may cause severe flooding and crop loss. Persistent surface weather can result either from quasi-stationary, long-lived atmospheric circulation conditions (quasi-stationarity) or from repeated, shorter-lived circulation features (recurrence). Recurrence refers to the repeated occurrence of similar large-scale circulation types or weather systems within some (S2S) time interval, usually with brief interruptions. Many recent high-impact weather and climate events were linked to persistent quasi-stationary or recurrent weather conditions. An example for recurrence are the western European floods in July 2021 that occurred at the end of an extreme wet spell in western Europe. The wet spell resulted from repeated atmospheric blocks and Rossby wave breaking episodes (Tuel et al., 2022b). Other examples of recurrence include the floods in the UK during winter 2013–2014 and in Queensland (Australia) in February–April 2022 both of which were caused by sequences of cyclones (Huntingford et al., 2014; Wikipedia, 2022; Floodlist, 2022). An example for quasi-stationarity is the catastrophic flooding in Pakistan in summer 2022 that resulted (in part) from long-lasting and particularly heavy monsoon rains (Mallapaty, 2022). Intense heat waves and associated atmospheric circulations also tend to be persistent (often arising from a combination of recurrence and quasi-stationarity) (Lorenz et al., 2010), as in western Europe in 2003 (Black et al., 2004; García-Herrera et al., 2010), western Russia in 2010 (Drouard and Woollings, 2018; Di Capua et al., 2021), or China (World Meteorological Organization, 2022) and India (Bloomberg, 2022) in 2022.

S2S weather persistence offers the potential for improved predictability at the S2S timescale (Franzke et al., 2011), which is highly relevant for risk preparedness and is attracting increased attention from the research community (e.g., Vitart et al., 2017; Meehl et al., 2021; Domeisen et al., 2022). However, a persistent state is not necessarily highly predictable, and persistent states with low predictability can cause large errors in sub-seasonal weather forecasts (Quandt et al., 2017).

Persistence is also an important aspect of climate model evaluation and climate projections. Whether global climate models are able to correctly simulate persistence is key to the robustness of long-term projections, especially of high-impact weather events – all the more so as climate projections suggest enhanced persistence (Li and Thompson, 2021; Hoffmann et al., 2021; Tuel and Martius, 2021a) in the future.

Characterizing weather persistence is therefore key to our understanding of the atmospheric circulation and its predictability and the associated hazards. Previous studies have focused on weather persistence from varied perspectives. Some assessed the persistence of specific weather systems or features, like atmospheric blocking (e.g., Liu, 1994), Rossby waves (Röthlisberger et al., 2019) or teleconnection patterns (Barnes and Hartmann, 2010). Others analyzed specific episodes of particularly persistent weather conditions (e.g., Black et al., 2004; Di Capua et al., 2021; Tuel et al., 2022b) or characterized the overall tendency of the atmospheric circulation and surface weather to exhibit persistence (e.g., MacDonald, 1992; Li and Thompson, 2021; Hoffmann et al., 2021). However, while previous studies generally agree on what persistence means conceptually, past work on this topic has involved many different definitions, often causing confusion and leading to different interpretations of persistence. Many case studies have also described observed situations as persistent based on subjective analyses rather than quantitative metrics. A further source of confusion is that “recurrence” is used in the literature to refer not only to successive occurrences of the same weather pattern at close intervals – what we will focus on in this review – but also to the states of the atmosphere with the highest probability of occurrence (Michelangeli et al., 1995).

It is difficult to give a unique definition of weather persistence. Besides, it may not even be desirable as different interpretations are possible and useful, depending on the system and timescale of analysis and on the motivations and goals of the study. Our goal here is therefore to review existing concepts of weather persistence, associated methodologies and physical interpretations. We present and structure a wide variety of approaches, definitions, techniques and metrics that have been used to analyze these concepts and that allow for answering one or several of the following questions.

• Is there persistence in the data?

• What are the persistent timescales in the data?

• In which specific periods does persistence occur?

• What are the persistent locations in the state space?

While persistence occurs at many different timescales, we specifically focus on the S2S timescale – often highly impact relevant and certainly important for predictability – but most methods and arguments in principle also apply to longer timescales. An overview of the methods is given in Fig. 1. We begin in Sect. 2 by introducing the two aspects of persistence: quasi-stationarity and recurrence. Section 3 then discusses several different perspectives on persistence to consider when choosing an analysis methodology. Finally, we present a detailed list of methods to detect or quantify persistence in Sects. 4 and 5. We illustrate most of the methods with examples taken from the literature or with our own analyses of summertime atmospheric circulation over Europe (details about the data we use are given in the data availability section). We keep the interpretation of the results to a minimum, the point being to illustrate the methods and not to analyze European summer circulation persistence in detail.

Figure 1 Overview of the persistence methods discussed in this paper. Section numbers relative to each method are indicated in bold between brackets.

2 Persistence: quasi-stationarity or recurrence?

Persistence in a dynamical system (climate variable, atmospheric circulation field, etc.) arises from the repeated occurrence of the same value(s) or pattern(s) over a period of time. Successive occurrences can follow each other continuously – a situation we refer to as “quasi-stationary” – or in an interrupted sequence – in which case we speak of recurrence. Persistence therefore comes in two flavors, quasi-stationarity and recurrence, which we illustrate with the example of extreme warm conditions in a region north of the Black Sea in Fig. 2. Here, extreme warmth is defined as a daily mean temperature exceeding 1.5 standard deviations above its annual cycle. In summer 2010, the region experienced persistent (quasi-stationary) extreme warmth, with temperatures remaining continuously above the threshold for 41 d straight (Fig. 2a). By contrast, in spring 1979, extreme warm temperatures occurred frequently throughout May, but the three extreme warm episodes within 1 month were separated by at least 5 d with near-average temperatures (Fig. 2c). The corresponding evolution of the atmospheric circulation can be conceptually described as a slow-moving trajectory for quasi-stationary conditions (Fig. 2b) and a more rapidly evolving trajectory that revisits the same point repetitively (Fig. 2d).

The literature often uses persistence as a synonym for quasi-stationarity (e.g., Dole and Gordon, 1983; Barnston and Livezey, 1987; Franzke et al., 2011; Di Lorenzo and Mantua, 2016; Fereday, 2017; Liu et al., 2018; Du et al., 2019; Francis et al., 2020; Hoffmann et al., 2021; Li and Thompson, 2021), i.e., persistence is associated with long-lived flow anomalies and little change in atmospheric circulation and surface weather. Recurrence has by contrast attracted less scientific attention. Most studies on recurrence have focused on extreme or impactful events (e.g., Mailier et al., 2006; Barton et al., 2016; Dacre and Pinto, 2020; Tuel and Martius, 2022a) and Rossby waves (e.g., Röthlisberger et al., 2019; Ali et al., 2021). Yet, from a physical perspective, quasi-stationarity and recurrence are intimately related. Recurring weather systems typically result from quasi-stationary favorable large-scale conditions, like sea surface temperature (SST) anomalies or the location of extratropical jets (e.g., Tuel and Martius, 2022b). Temporal dependence, or memory, in a system can thus translate in practice into both quasi-stationarity and recurrence (Franzke, 2013). Additionally, from the impacts perspective, it makes sense to look at quasi-stationarity and recurrence together, since both can cause prolonged, impactful surface weather. Recurrent Rossby waves, for instance, modulate the persistence of surface temperature and precipitation anomalies (Röthlisberger et al., 2019; Ali et al., 2021) while long-lived SST anomaly patterns like El Niño–Southern Oscillation (ENSO) can remotely trigger recurrent extreme weather (Gershunov and Barnett, 1998).

Note that recurrence, as we define it here, is sometimes referred to as “temporal” or “serial clustering” (Franzke, 2013), for instance in the case of recurrent cyclones (Mailier et al., 2006) or heavy precipitation events (Barton et al., 2016). Note also that recurrence is also frequently used in meteorology to refer to preferred patterns that repeatedly occur in a time series but not necessarily at close intervals (e.g., Vigaud et al., 2018; Kornhuber et al., 2019; Son et al., 2021) – in the case of atmospheric circulation patterns, one speaks of circulation or weather “regimes” (Michelangeli et al., 1995). This is different from our definition, in which recurrence specifically relates to the repeated occurrence of the same patterns over S2S timescales. Such patterns may be rare in the full dataset and would therefore not be considered “recurrent” in the regime approach but can be highly relevant from an impacts perspective. Hannachi et al. (2017) gave a comprehensive overview of the weather regime approach and associated methodologies. We also discuss the relevance of the weather regime perspective for quasi-stationarity analysis in Sect. 4.2.1.

Figure 2 Illustrating quasi-stationarity and recurrence on S2S timescales. (a, c) Example of daily mean 2 m temperature series averaged over the 50–57^∘ N, 33–43^∘ E region illustrating (a) quasi-stationarity (summer 2010) and (c) recurrence (spring 1979) in extreme warm temperatures on S2S timescales (black: observations; gray: mean annual cycle; light gray: +1.5 standard deviation from the mean). Data are from the ERA5 reanalysis (Hersbach et al., 2020). (b, d) Idealized system trajectories (thick black lines) in the phase space corresponding to quasi-stationarity (b) and recurrence (d). The background PDF is shown in light contours. Panels (b) and (d) reproduced from Hannachi et al. (2017). © 2017 American Geophysical Union, all rights reserved.

Because quasi-stationarity and recurrence are two faces of the same coin, distinguishing one from the other may not be evident nor necessarily relevant.

First, the distinction often depends on the variable of interest. Recurrent weather systems can indeed result in quasi-stationary surface weather anomalies and vice versa. For example, the long-lived heat waves of the 2010 and 2021 summers in western Russia (Fig. 2a) and the Baltic were linked to recurrent atmospheric blocks (Drouard and Woollings, 2018; Tuel et al., 2022b). Likewise, prolonged droughts or wet spells may result from recurrent Rossby wave activity (Röthlisberger et al., 2019; Ali et al., 2021, 2022). Quasi-stationary surface warm and humid conditions can also trigger recurrent thunderstorm activity (Mohr et al., 2020). Conversely, recurrent extreme precipitation events (e.g., Barton et al., 2022) or extratropical cyclones (e.g., Dacre and Pinto, 2020) are often linked to quasi-stationary jet states that last for much longer than the lifetime of individual weather systems.

Second, the longer the timescale of analysis, the less obvious the difference between quasi-stationarity and recurrence becomes. Impacts, for example, often depend on anomalies of surface temperature or precipitation averaged or accumulated over several weeks to months (e.g., droughts). Weekly or monthly values may thus sometimes be preferred to daily ones, in which case synoptic-scale quasi-stationarity and recurrence would both result in large weekly or monthly anomalies that would result from simple “persistence”.

Nevertheless, the distinction between quasi-stationarity and recurrence remains highly relevant for several reasons: from a methodological perspective (Sects. 4 and 5) but also for forecasting and process understanding at the synoptic timescale and for some stakeholders like insurers (for whom it matters whether impacts resulted from a single or multiple events). For process understanding, taking a weather systems perspective is often very relevant. In such a case, distinct weather systems (like cyclones) can in principle be separated from one another, and long-lived single systems be distinguished from multiple short-lived systems occurring in close succession. It therefore matters whether persistence is driven by recurrence or quasi-stationarity. The distinction is also important to assess whether numerical models simulate persistence for the right physical reasons.

3 Considerations and methodological approaches for persistence analysis in weather and climate data

Before reviewing how persistence can be apprehended and quantified, let us begin with some basic notation and definitions. In the following, we denote the dynamical system under analysis with x(t)∈ℝ^m.

x(t) could, for example, be the weather evolution in a region, the time evolution of one variable at one location or the general circulation. (x(t))_t evolves within a state space 𝒮 that consists of the set of all possible system states ${\left\{s_i\right\}}_i$ that the system can occupy. A single system value and a single system state corresponds to each time t, x(t). The same system state can be attained at multiple different time steps. The series of successive system values constitutes the system trajectory. In practice, only a finite number of observations are available $\left(\mathit{x}\right(t_k){)}_{\mathrm{1}\le k\le N}$ to characterize the persistence of a system.

Characterizing a system as “persistent” can mean different things. Consequently, it is important to always specify the perspective that is taken to avoid confusion. First, there are different flavors to persistence (global, state or episodic persistence; see Sect. 3.1). Second, persistence can be studied from a Lagrangian or an Eulerian perspective (Sect. 3.2). Finally, persistence is linked to a similarity metric (Sect. 3.3) and a notion of timescale (Sect. 3.4), for which different approaches are possible.

3.1 Global, state and episodic persistence

Persistence (whether quasi-stationarity or recurrence) is a broad concept that covers different kinds of behavior in dynamical systems. We might say, for instance, that temperature is more persistent than precipitation because temperature evolves, on average, over longer timescales than precipitation. In this sense, persistence quantifies the system's inertia. However, if we qualify last summer's weather as particularly persistent, we mean something different: namely, last summer's weather varied much less than what one could have reasonably expected in a normal summer. Here, persistence refers to some unusual behavior of the system over a particular period. Saying that zonal jets or atmospheric blocks are persistent again means something else, i.e., that these particular states of the circulation tend to be more long-lived or recurrent than other states.

This leads us to make the distinction between three types of persistence: global, state and episodic persistence.

1. Global persistence. This characterizes the tendency to quasi-stationarity or recurrence across the whole trajectory of the dynamical system. We choose this term because persistence in this sense is a “global” property of the system, meaning that it is not restricted to any particular system state or time period. Global quasi-stationarity translates into the tendency for the system to change little at small timescales (successive values being close to each other). In mathematical terms, this translates into $<\left|\frac{\mathrm{d}\mathit{x}\left(t\right)}{\mathrm{d}t}\right|>\ll \frac{{\mathit{\sigma }}_x}{T}$, where $<\cdot >$ is a time average, σ_x is the standard deviation of the series x(t) and T its typical timescale of evolution. Global persistence can sometimes be impact relevant (e.g., trends in global persistence under climate change can be important for impacts). However, it does not make distinctions between the different system states. Consequently, global persistence is not suited to characterize persistent system states or specific time intervals with persistent system behavior and is generally not the best-suited approach for risk assessment. Global persistence is, however, strongly related to intrinsic system predictability, since present values of the system largely constrain its future values. It can thus yield important information for numerical forecasting, including at the S2S timescale. Most global persistence methods focus on quasi-stationarity – like the autocorrelation coefficient (Sect. 4.1.1) or the Hurst exponent (Sect. 4.1.2) – but some exist for recurrence as well (Sect. 5.2.2). Many studies look at global persistence, and we may cite, for instance, MacDonald (1992); Pfleiderer and Coumou (2018); Pfleiderer et al. (2019) and Li and Thompson (2021), who analyze the persistence of temperature series, and Hoffmann et al. (2021), who consider 10 d atmospheric flow persistence.

2. State persistence. This relates to the persistent behavior of specific system states. State space persistence analysis consists of either identifying persistent system states – e.g., with the quasi-stationary (Sect. 4.2.2) or optimally persistent pattern (Sect. 4.2.3) methods – or characterizing the persistence of given states with methods like residence times (Sect. 4.2.5) for quasi-stationarity or Ripley's K for recurrence (Sect. 5.1.3). State persistence is highly relevant from the impact, forecasting and process understanding perspectives. Knowing which states are persistent and to what extent is useful (i) to make the link to surface impacts and (ii) to know which weather patterns or sequences may be more predictable than others. It also helps to determine the physical processes that support or are responsible for the persistence. The state approach to persistence has been used to shed light on quasi-stationary states of the North Atlantic circulation and their predictability (Faranda et al., 2017b), to characterize the quasi-stationarity of continental-scale weather patterns (Francis et al., 2018, 2020) or drought (Ford and Labosier, 2014), or to quantify recurrence in extra-tropical cyclones (Mailier et al., 2006) and extreme precipitation (Tuel and Martius, 2021a).

3. Episodic persistence. This is tied to specific time intervals (or “episodes”) during which the system exhibits quasi-stationary or recurrent behavior. It is in that sense a purely local property that characterizes the anomalous behavior of the system over a limited time period. Importantly, episodic persistence can occur simply by chance in any dynamical system, even in systems that exhibit no global persistence. Episodic persistence is therefore not necessarily relevant for predictability and process understanding. Persistent periods can still always be analyzed to look for potential drivers that discriminate pure “statistical flukes” from possibly predictable events. Episodic persistence is, however, well suited for impacts analysis because it can make a direct link between periods of persistent weather and impacts. Relevant methods include running window techniques (Sect. 4.3) for quasi-stationarity and window counts (Sect. 5.1.1) for recurrence. Hoffmann et al. (2021) investigated for example quasi-stationarity in 10 d sequences of atmospheric circulation, while Bevacqua et al. (2020) and Kopp et al. (2021) looked at sub-seasonal periods with recurrent cyclones and extreme precipitation events.

While it helps to capture the various interpretations of persistence, this classification is not perfect, and there is some overlap between categories. In practice, the state or episodic perspectives can also be used to quantify global persistence (by averaging persistence metrics across systems states or time intervals) and global persistence metrics can be computed on subsets of the data to identify persistent periods. Some methods, like recurrence plots (Sect. 5.2.2), can even deal with all three types of persistence.

The three types of persistence are also not independent from one another. Global persistence, for instance, can emerge from repeated occurrences in one or a handful of system states only, while the rest of the trajectory, if analyzed separately, may not be qualified as persistent. There are also strong relationships between state and episodic persistence. The most common system states to occur during persistent periods are indeed likely to be persistent states. Correspondingly, persistent states, when they occur, are likely to be associated with persistent periods. Persistent states can thus be uncovered from the knowledge of persistent periods, for instance with pattern recognition or clustering algorithms applied to system values during persistent periods, or simply by averaging system values during persistent periods (e.g., Faranda et al., 2017b; Hoffmann et al., 2021).

However, state persistence only characterizes the average behavior of system states – it does not imply that all occurrences of a persistent state will necessarily be persistent. Consequently, there is no one-to-one relationship between persistent states and persistent periods. Some system states can behave in a persistent way under certain conditions but not under others. The lifetime and travel speed of atmospheric blocks, for instance, is affected by land–atmosphere feedbacks or upstream latent heating (Steinfeld et al., 2020). Similarly, extratropical cyclones may occur in sequences but also as single events (Dacre and Pinto, 2020). Additionally, persistent periods may be associated with a variety of system states. While we expect persistent states to be the most frequent ones during persistent periods, some states which on average are not persistent can still at times behave persistently by pure chance.

However, this classification is useful to illustrate the different methodological ways that weather persistence can be tackled, and we rely on it to structure the description of methods in Sects. 4 and 5.

3.2 Lagrangian and Eulerian perspectives

Weather persistence is most often analyzed from an Eulerian perspective, i.e., persistence of the same weather pattern or conditions at a fixed location in space. By contrast, in the Lagrangian perspective the focus is on the persistence of a given weather pattern in time. In the Eulerian framework, the system x(t) represents a time series over a domain fixed with time, whereas in the Lagrangian one x(t) follows individual weather patterns or the background atmospheric flow. The Eulerian stationarity of a quantity ψ(x,t) translates into $\partial \mathit{\psi }/\partial t=\mathrm{0}$, while the Lagrangian stationarity implies $\partial \mathit{\psi }/\partial t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }=\mathrm{0}$, where u is the background flow. Similarly, quasi-stationarity translates into $\partial \mathit{\psi }/\partial t\approx \mathrm{0}$ in the Eulerian framework and $\partial \mathit{\psi }/\partial t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }\approx \mathrm{0}$ in the Lagrangian one. For example, temperature anomalies can be tracked over time (e.g., Kornhuber and Tamarin-Brodsky, 2021) or analyzed at a fixed location (e.g., Pfleiderer and Coumou, 2018; Li and Thompson, 2021). Similarly, weather systems such as blocking, cyclones or vortices (Bray and Cavallo, 2022) can be analyzed at a fixed location or following the weather systems. For example, Økland and Lejenäs (1987) contrast the persistence of blocking at fixed longitudes against the persistence of individual blocking episodes. Another example is Kossin (2018), who characterize the persistence in tropical cyclones by their translation speed.

Both the Eulerian and the Lagrangian perspectives are relevant for impact and risk assessment. The former links persistence to impacts at a given location, and the latter highlights impacts along the trajectory of a weather pattern (system) during its lifetime. Indeed, the same weather system can produce hazardous weather over large areas, putting strain on the resources of insurance companies or governments. The two perspectives can be brought together by considering the translation speed and lifetime of the tracked weather systems. Systems with long lifetimes and low translation speed lead to both Eulerian and Lagrangian persistence. By contrast, long-lived systems that travel quickly are persistent from a Lagrangian perspective only. Likewise, slow-moving but short-lived systems are not Lagrangian persistent, but Eulerian persistence can still be detected if several such systems occur over the same area in close succession.

3.3 Quantifying similarity

Assessing persistence typically requires quantifying the self-similarity of system values x(t) with a metric (e.g., Wharton et al., 2008; Zerzucha and Walczak, 2012; Ali et al., 2020; Ontañón, 2020). By self-similarity, we refer to the tendency for successive values of x(t) to be similar to each other according to some metric. This is not to be confused with the concept of geometric self-similarity in fractal geometry. Metric selection is an important step that should be done with care because it conditions how persistence is quantified and some metrics are not suited to certain kinds of data (e.g., heavily skewed). There are two main classes of similarity metrics.

1. Categorical metrics. With these metrics, system values are either similar (if they belong to the same category) or not (if they do not).

2. Continuous metrics. These metrics measure the degree of similarity between two system values in a continuous way.

Categorical metrics focus on specific features of the system, such as the presence of a given weather pattern or the occurrence of a specific event. They require the set of all x(t) values to be classified into distinct categories (usually from two to a few dozen). Categories can represent predefined system states of interest (e.g., a warm or cold anomaly, a given phase of a teleconnection mode, or the occurrence of a specific weather pattern, like a block) (e.g., Pinto et al., 2014; Drouard and Woollings, 2018; Pfleiderer and Coumou, 2018; Ali et al., 2021; Kopp et al., 2021) but can also be obtained objectively with dimension reduction or pattern recognition methods. Examples include principal component analysis (or empirical orthogonal functions, EOFs) (Fereday, 2017, e.g.,), optimally interpolated patterns (OPPs) (Hannachi, 2008), self-organizing maps (SOMs) (e.g., Francis et al., 2018; Weiland et al., 2021; Rousi et al., 2022b) or clustering algorithms (probabilistic, hierarchical or non-hierarchical) (e.g., Demuzere et al., 2011; Hannachi et al., 2012; Grams et al., 2017). When x(t) represents 2D circulation data (like sea level pressure or geopotential height), the ensemble of system states is often referred to as “weather regimes” (Michelangeli et al., 1995; Grams et al., 2017; Francis et al., 2018) (see Sect. 4.2.1). In EOF analysis, the distance metric is the L² norm (Euclidean distance), but most clustering methods can work with any distance metric. The challenge with dimension reduction methods is choosing the number of categories to retain (EOFs, clusters, SOM nodes, etc.). A high number can capture rare states of the system, but this comes at the cost of making persistence more difficult to assess (since sequences when the system remains in the exact same category will become less frequent).

Continuous metrics measure the degree of similarity between any pair of system values in a continuous way. Common examples of continuous metrics include the Euclidean distance (e.g., Faranda et al., 2017b), pattern correlation(e.g., Mo and Ghil, 1987), or more complex similarity indices like the Teweles–Wobus score (e.g., Horton et al., 2017; Blanchet et al., 2018) or the image structural similarity index (SSIM) (e.g., Hoffmann et al., 2021).

In comparison to categorical metrics, continuous metrics offer the advantage that they do not require specifying features or events of interest beforehand. They are also more flexible insofar as similarity can be quantified with respect to any system value as reference and not just representative values for each category. Since they do not require simplifying the state space, continuous metrics may also be able to pick up rare persistent patterns that are missed by dimension reduction methods that focus on the most common patterns in a series. These advantages come at a cost: working with continuous metrics, especially complex ones, can be more computationally intensive and sometimes more difficult to interpret physically.

3.4 Persistence timescales

Persistence is linked to a notion of timescale during which the system continuously remains in the same state (for quasi-stationarity) or occupies that same state repeatedly (for recurrence). There are three common ways to approach persistence timescales.

The first option is to choose a single, fixed timescale for analysis. This choice can be guided by impact and forecasting considerations, by physical knowledge of the underlying system, or by observations of persistent events (e.g., Huntingford et al., 2014; Lawrence et al., 2020; Overland and Wang, 2021; Rakovec et al., 2022; Rousi et al., 2022a). This is the most common approach to analyze episodic persistence, but it is also applicable to global persistence. Quasi-stationarity can be quantified by the average similarity between n successive states (for continuous similarity metrics see, e.g., Kolstad et al., 2017; Hoffmann et al., 2021) or by the variety of system states during an n-step window (for categorical metrics see, e.g., Fereday, 2017; Richardson et al., 2019). Quasi-stationarity can also be inferred from extreme anomalies of circulation, temperature or precipitation during n-step windows. For instance, Gálfi et al. (2019) and Tuel and Martius (2023) identify persistent warm and cold spells by averaging temperature anomalies over 1–3 weeks.

Recurrence can similarly be assessed by calculating the number of times that a particular system state or event occurs during n-step windows (e.g., Mailier et al., 2006; Pinto et al., 2014; Kopp et al., 2021; Tuel and Martius, 2022a). A fixed timescale can also be used as a threshold to separate quasi-stationary from non-stationary events, by requiring quasi-stationary events to last at least n steps. For instance, Francis et al. (2018) and Francis et al. (2020) define persistent periods by requiring the circulation pattern in a region to remain in the same state for at least 4 consecutive days. Mann et al. (2018) similarly define persistent resonant wave events as those lasting at least 10 d. Dole and Gordon (1983) also take this approach for the persistence of point-wise geopotential anomalies. Note that this fixed timescale can also consist of a single time step. In this case, persistence characterizes how much past system values determine future ones. Li and Thompson (2021), for instance, characterize persistence with the lag-1 autocorrelation coefficient, Röthlisberger and Martius (2019) look at 1 d transition probabilities between different system states, and Vautard (1990); Michelangeli et al. (1995) and Hannachi et al. (2017) calculate time derivatives of geopotential fields to identify quasi-stationary states.

The second option is to select an analysis method that explores a range of timescales and pinpoints the relevant persistence timescales, sometimes accompanied by some notion of statistical significance. This approach only works for global persistence. Autocorrelation analysis, for instance, detects the timescales at which the system exhibits significant lagged memory (Sect. 4.1.1). Spectral analysis can likewise highlight important timescales of variability that can be linked to quasi-stationarity (Sect. 4.1.2). For recurrence, methods like Ripley's K function indicate the timescales at which recurrence is statistically significant (Sect. 5.1.3).

Finally, the third option is to characterize persistence not by a single timescale but by a distribution of timescales. To assess quasi-stationarity, one can typically work with the distribution of persistent event durations. Persistent events are periods during which the system satisfies a persistence criterion: for quasi-stationarity, successive system values must be similar, and for recurrence the same event must occur multiple times, with each occurrence being separated by at most n time steps from the previous one. For recurrence, it is also possible to consider the distribution of inter-event times (e.g., Altmann and Kantz, 2005). The distribution of event lengths can be further summarized by considering the average or maximum event length (Faranda et al., 2017b; Rousi et al., 2022b) or by modeling it with an exponential or power law distribution (Sect. 4.2.5). This approach has been applied to numerous cases: heat waves (Lorenz et al., 2010), droughts (Meng et al., 2017; Moon et al., 2018), wet spells (Ali et al., 2021), atmospheric blocks (Liu, 1994), circulation patterns (Huguenin et al., 2020), and midlatitude cyclone clustering (Bevacqua et al., 2020).

4 Quasi-stationarity

The diversity of perspectives on persistence translates into a wide range of methods, of which we give an overview in the following two sections dedicated, respectively, to quasi-stationarity and recurrence. Following the distinctions introduced in Sect. 3, we separate methods that quantify global, state, and episodic persistence (though some methods can be used for more than one). We also specify whether methods can only be used with a single timescale or whether they quantify persistence across timescales. An overview of methods is shown in Table 1.

4.1 Global quasi-stationarity

We begin with several methods that quantify global quasi-stationarity in one-dimensional time series. They characterize quasi-stationarity in the series as a whole but are generally unable to identify quasi-stationary states or periods. However, they often allow the user to characterize the timescales of variability and persistence in the data and are hence relevant for system predictability and process understanding.

4.1.1 Autocorrelation

Autocorrelation is a frequently used measure of quasi-stationarity in weather and climate science. If X_t is a continuous, one-dimensional process of mean μ and variance σ², its Pearson autocorrelation coefficient at lag k is defined as follows:

$$\begin{array}{}\text{(1)}& \mathit{\rho }\left(k\right)={\displaystyle \frac{\mathbb{E}\left[\left(X_t-\mathit{\mu }\right)\left(X_{t+k}-\mathit{\mu }\right)\right]}{{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

where 𝔼 denotes the expectancy with respect to the distribution of X_t. Potential cycles and long-term trends should be removed from the data before analysis (Weiss and Weiss, 1999). In Eq. (1), X_t can also be replaced by its rank, which yields the alternative Spearman autocorrelation that is more robust to nonlinear behavior in the data. The definition of Eq. (1) can also be extended to third-order statistics to capture interactions coming from nonlinear correlation, yielding the bi-correlation $\mathit{\rho }(k_{\mathrm{1}},k_{\mathrm{2}})=\frac{\mathbb{E}\left[\left(X_t-\mathit{\mu }\right)\left(X_{t+k_{\mathrm{1}}}-\mathit{\mu }\right)\left(X_{t+k_{\mathrm{2}}}-\mathit{\mu }\right)\right]}{{\mathit{\sigma }}^{\mathrm{3}}}$. (Pires and Hannachi, 2021).

Several summary metrics for autocorrelation exist, like the autocorrelation timescale:

$$\begin{array}{}\text{(2)}& T_{\mathit{\rho }}^{\mathrm{e}}=\underset{k>\mathrm{0}}{min}\left\{\mathit{\rho }\left(k\right)\le \exp \left(\mathrm{1}\right)\right\};\end{array}$$

the decorrelation time (Hannachi, 2021)

$$\begin{array}{}\text{(3)}& T_{\mathit{\rho }}^{\mathrm{d}}=\mathrm{1}+\mathrm{2}\sum _{k=\mathrm{1}}^N\left|\mathit{\rho }\right(k\left)\right|,\end{array}$$

where N is the number of time steps (if ρ is integrable); or the characteristic time $T_{\mathit{\rho }}^{\mathrm{c}}$ (Trenberth, 1984) (Fig. 3a):

$$\begin{array}{}\text{(4)}& T_{\mathit{\rho }}^{\mathrm{c}}=\mathrm{1}+\sum _{k=\mathrm{1}}^N\left(\mathrm{1}-{\displaystyle \frac{k}{N}}\right)\left|\mathit{\rho }\right(k\left)\right|.\end{array}$$

The absolute value on ρ(k) in Eqs. (3) and (4) is not required but can help avoid underestimated timescales in case ρ exhibits oscillations (Franzke et al., 2005).

Both $T_{\mathit{\rho }}^{\mathrm{d}}$ and $T_{\mathit{\rho }}^{\mathrm{c}}$ roughly approximate the average time between successive independent values. The lag-1 autocorrelation ρ(1) can also be used as a quasi-stationarity metric (Li and Thompson, 2021), which is equivalent to fitting a linear regression model between successive system values: $\mathit{x}(t+\mathrm{1})=\mathit{\beta }\mathit{x}\left(t\right)+\mathit{ϵ}\left(t\right)$, where ϵ(t) is a white-noise process. If x(t) is normalized, the regression slope β is equal to ρ(1). Note that more complex models are possible: linear models for continuous series can be extended by the large class of autoregressive-moving average (ARMA) models.

Autocorrelation has many advantages: it has a simple definition, the ease of interpretation (ρ(k) of which is related to the linear regression coefficient of X_t+k against X_t); it has high flexibility; it is already implemented in common programming languages; it requires no subjective threshold choice; and its results can be easily reproduced. By varying the lag k, it can measure quasi-stationarity at all timescales. It also comes with a notion of statistical significance: given a confidence level, it is possible to say whether the obtained autocorrelation is significant (indicating a link at lag k) or not, pointing to the relevant quasi-stationarity timescales in the series. On the downside, autocorrelation only works for one-dimensional data and requires a large number of values as input. Therefore, it is best suited to measuring quasi-stationarity globally. One can compute autocorrelation on a subset of the data only (summer values or a specific time interval, for instance) if the number of data points is large enough, in which case autocorrelation may be used to characterize quasi-stationarity locally in time. Additionally, autocorrelation only highlights linear relationships (though bi-correlation can help capture nonlinearity). Finally, autocorrelation measures the strength of the connection between lagged system values but not how far apart they might be in the state space.

Example applications include Horel (1985a) and Barnston and Livezey (1987), who calculated lag autocorrelation on principal component time series of monthly Northern Hemisphere geopotential fields. MacDonald (1992) used autocorrelation to detect quasi-stationarity in monthly temperature series, as did Weiss and Weiss (1999) to assess quasi-stationarity in ENSO. Weatherhead et al. (2010) characterized daily temperature quasi-stationarity in the Arctic based on lag-1 autocorrelation. Degenhardt and Ólafsson (2019) and Kolstad et al. (2015) calculated lag-1 autocorrelation to highlight intra-seasonal quasi-stationarity of monthly mean temperatures in Iceland and Europe, respectively, and Li and Thompson (2021) applied autocorrelation to daily temperature series and found it was strongly related to the average length of warm and cold episodes across the world. Kolstad et al. (2017) even used autocorrelation analysis in a causal discovery framework by regressing temperature values against previous ones and including potential covariates.

We show the characteristic time for daily rescaled 500 hPa geopotential height over Europe during summer in Fig. 3a. While the results say nothing about the role of individual weather systems, the larger values over the British Isles and western Russia seem consistent with more frequent persistent atmospheric blocks over these regions.

Figure 3 (a) Autocorrelation characteristic time ($T_{\mathit{\rho }}^{\mathrm{c}}$; Eq. (4)) of June–July–August (JJA) daily rescaled Z500. (b) Long-range dependence parameter d of unfiltered surface wind speeds. Only those values that are significant at the 5 % level are displayed. Reproduced from Franzke (2013). © 2013 The Author(s) Published by the Royal Society. All rights reserved.

4.1.2 Asymptotic methods to characterize variability across timescales

Persistence in a system is associated with memory effects that lead to variability being concentrated at long rather than short timescales. How variability in the series is distributed across timescales is therefore an important indicator of global persistence and can point to relevant persistence timescales. Specifically, persistence often translates into scaling laws: in the time series variability as a function of frequency or timescale (for quasi-stationarity) (Bunde et al., 2002) but also in the distribution of inter-event times (for recurrence, see Sect. 5.1.4). The slope of these laws indicates the degree of persistence. Time series also often exhibit different scaling laws over different frequency intervals, highlighting how persistence may differ across timescales in the data.

A first scaling law can be obtained by considering how autocorrelation decreases as a function of time lag. The more quasi-stationary a series, the less rapidly its autocorrelation should decrease. Conceptually, time series can be broadly divided between short-range and long-range dependence series. For a short-range dependent series (like autoregressive models), the autocorrelation decreases rapidly with the time lag, eventually reaching 0 after a certain lag or decaying exponentially.

$$\begin{array}{}\text{(5)}& \mathit{\rho }\left(k\right)\sim {\mathit{\alpha }}^k,\end{array}$$

as k→∞ and with $\mathrm{0}<\left|\mathit{\alpha }\right|<\mathrm{1}$ (with α<0 corresponding to anti-persistence). By contrast, in a long-range dependent series, the autocorrelation decays following a power law:

$$\begin{array}{}\text{(6)}& \mathit{\rho }\left(k\right)\sim k^{\mathrm{2}d-\mathrm{1}},\end{array}$$

where d is called the dependence parameter ($\mathrm{0}<d<\frac{\mathrm{1}}{\mathrm{2}}$) (Beran, 2017; Franzke, 2013). A lower d is associated with a slower decay of the autocorrelation function and hence with more quasi-stationarity systems. White noise has d=0. Witt and Malamud (2013) and Franzke (2013) discuss several ways how d can be estimated.

Like autocorrelation, α and d (the dependence parameter) measure global quasi-stationarity in continuous time series because they are based on asymptotic relationships. However, they cannot highlight specific quasi-stationarity timescales in a series as autocorrelation does. Applications to atmospheric time series show that temperature exhibits long-term dependence (e.g., Yuan et al., 2010; Koscielny-Bunde et al., 1998; Eichner et al., 2003), while precipitation exhibits either short- or long-term dependence (Potter, 1979; Hannachi, 2014; Yang and Fu, 2019). Franzke (2013) computed d and found long-range dependence in North Atlantic winds (Fig. 3b).

The Hurst coefficient (or exponent) H (Hurst, 1951) is a common measure of memory in a time series which is obtained from a second scaling law. Specifically, H characterizes how a time series (X_t)_t fluctuates relative to its mean. Noting

$$\begin{array}{}\text{(7)}& Z_i^n=\sum _{t=(i-\mathrm{1})n+\mathrm{1}}^{in}\left(X_t-\mathbb{E}\left[X_t\right]\right),\end{array}$$

the cumulative fluctuations in (X_t)_t relative to its mean, calculated over intervals of size n, Hurst argued, empirically from geophysical time series, that these fluctuations could exhibit scale-invariant properties. In other words, cumulative fluctuations over two different timescales n and m are related through

$$\begin{array}{}\text{(8)}& Z_i^n\textcolor{red}{}{\left({\displaystyle \frac{n}{m}}\right)}^H{Z}_j^m\forall (i,j,k,l)>\mathrm{0},\end{array}$$

where stands for equality in distribution. H ranges from 0 to 1. $H<\mathrm{1}/\mathrm{2}$ indicates anti-persistent behavior, such that successive increments of X_t relative to its mean (X_t−𝔼[X_t]) tend to be negatively correlated, and the time series fluctuates substantially at short timescales. By contrast, $H>\mathrm{1}/\mathrm{2}$ points to persistent behavior, in which successive increments are positively correlated, and variability is concentrated at long timescales. $H=\mathrm{1}/\mathrm{2}$ corresponds to white noise (no temporal correlation). Thus, the higher H is, the smoother the time series. H is also theoretically related to the dependency parameter (as $H=\mathrm{2}d+\mathrm{1}$) and the power spectrum exponent (see below) (Franzke et al., 2020).

H can be estimated in many ways (see, e.g., Koutsoyiannis, 2003; De las Nieves López García and Requena, 2019; Franzke et al., 2020). In the literature, H has mainly been used to characterize long-term dependence or (quasi-)stationarity (Mandelbrot and Wallis, 1969), for instance in series of monthly or annual mean temperature (MacDonald, 1992; Kumar et al., 2013), precipitation (Bunde et al., 2013), or drought indices (Tatli, 2015). However, some studies also calculated it for daily time series (e.g., Rehman and Siddiqi, 2009; Velásquez Valle et al., 2013).

Like the Hurst exponent, spectral analysis also characterizes how a time series' variability is distributed across timescales. The power spectrum is commonly defined as the Fourier transform of the autocorrelation function ρ:

$$\begin{array}{}\text{(9)}& S\left(f\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathit{\rho }\left(k\right)e^{-\mathrm{2}i\mathit{\pi }fk}\mathrm{d}k,\end{array}$$

where f is the frequency and k the time lag. In a pure white-noise series, the variability is distributed equally across frequencies. Hence, the power spectral density is a constant. If temporal dependence is present, however, the power spectrum typically exhibits a power law decrease with frequency f:

$$\begin{array}{}\text{(10)}& S\left(f\right)\sim f^{-\mathit{\beta }},\end{array}$$

where β, called the power spectrum exponent, indicates the degree of quasi-stationarity in the time series ($\mathrm{0}<\mathit{\beta }<\mathrm{1}$). For statistically stationary series, one can also show that $\mathit{\beta }=\mathrm{2}H-\mathrm{1}$, where H is the Hurst coefficient (Parzen, 1986).

The larger β is, the more variance is concentrated at low frequencies. This implies a memory effect at low frequencies that relates to quasi-stationarity in the time series. Such scaling behavior is very common in climatological series, and the spectrum is often divided into distinct scaling regimes corresponding to specific frequency intervals. Fraedrich and Larnder (1993) discuss the example of precipitation in continental Europe and relate the different regimes to specific physical processes and timescales (from individual storms at high frequencies to climate fluctuations at low frequencies) (Fig. 4). Yang and Fu (2019) obtain similar results using hourly and daily precipitation data for the United States. Pelletier and Turcotte (1997) also used power spectra to quantify quasi-stationarity in various monthly climatic series, as did Ault et al. (2014) for drought quasi-stationarity and MacDonald (1992) for temperature quasi-stationarity. Telesca et al. (2016) analyzed quasi-stationary regimes in 10 min wind series across Switzerland.

Figure 4 Schematic diagram of the scaling regimes of continental European rainfall obtained by spectral analysis, along with the hypothesized meteorological interpretations of the various regimes. Reproduced from Fraedrich and Larnder (1993) under the terms of the Creative Commons CC BY license.

Note that, as with the autocorrelation, the definition of the power spectrum can be extended to third-order statistics to yield the bispectrum (the Fourier transform of the bi-correlation function) (Pires and Hannachi, 2021). While the power spectrum characterizes how the signal's variance is distributed across timescales, the bispectrum provides the contribution of each pair of frequency to the signal's skewness. Pires and Hannachi (2021), for instance, calculated the bispectrum of the 3-monthly ENSO series to better understand its predictability.

4.2 State quasi-stationarity

We now turn to methods that focus on the quasi-stationarity of specific system states. Unlike global methods, which take one-dimensional time series as input, several of the following methods are directly applicable to multidimensional data. We begin with methods that identify quasi-stationary states from the system trajectory (Sects. 4.2.1–4.2.4), before discussing methods that quantify the average quasi-stationarity of a system state (Sects. 4.2.5–4.2.6).

4.2.1 Weather regimes

We begin with an identification method for quasi-stationary states based on the concept of weather regimes (Michelangeli et al., 1995). This concept emerges from the realization that the extra-tropical atmospheric circulation evolves mainly as a succession of a handful of large-scale circulation patterns (Hannachi et al., 2017). These preferred sub-seasonal flow patterns, or weather regimes, tend to be quasi-stationary over timescales of a few days to a few weeks and are therefore strongly related to weather persistence. They account for much of the low-frequency atmospheric variability at intra-seasonal timescales (Pandolfo, 1993; Hannachi et al., 2017). The existence of weather regimes in the midlatitudes has long been recognized, such as the concept of Grosswetterlagen (Baur, 1951) or atmospheric blocking (e.g., Namias, 1964). They have attracted considerable attention, in particular because of the potential long-range predictability they offer (e.g., Ghil and Robertson, 2002; Büeler et al., 2021) and their link to surface impacts (e.g., Grams et al., 2017).

There are many ways to calculate weather regimes (see Huth et al., 2008, and Hannachi et al., 2017, for detailed overviews). The most common methods include pattern recognition and dimensionality reduction techniques, applied to a proxy variable for the atmospheric circulation like 500 hPa geopotential height or sea level pressure. Examples include orthogonal pattern decomposition techniques like principal component analysis (in the temporal or spectral domain) (Fereday, 2017; Grams et al., 2017) or optimally interpolated patterns (OPPs) (Hannachi, 2008), self-organizing maps (Huth et al., 2008; Francis et al., 2020; Weiland et al., 2021), archetypal analysis (Hannachi and Trendafilov, 2017; Chapman et al., 2022), and clustering algorithms (Fig. 5). The latter can be probabilistic (e.g., Gaussian mixture models; Woollings et al., 2010), hierarchical (e.g., Ward clustering; Hannachi et al., 2012) and non-hierarchical (e.g., k means or partitioning around medoids; Grams et al., 2017). Various statistics (gap statistic, silhouette coefficient, etc.) can help objectively select an optimal number of clusters, many of which are available from the R package clusterCrit (Desgraupes, 2018). Note that, in practice, the input data should be normalized to remove long-term and seasonal trends to focus exclusively on intra-seasonal variability (Grams et al., 2017).

Regimes can also be identified as local maxima in the (multidimensional) probability distribution function (PDF) of the target field, obtained empirically through e.g., kernel smoothing (Kimoto and Ghil, 1993; Woollings et al., 2010), or with more complex tools of network theory (Mukhin et al., 2022) and topology (Strommen et al., 2022). Finally, Franzke et al. (2011) identify regimes in the North Atlantic jet position with a hidden Markov model (HMM). HMMs are a powerful tool that brings together Markov models and Gaussian mixture models. Given N unknown (hidden) states, the HMM models the distribution of the observed series conditionally on each hidden state, with the sequence of hidden states following a first-order Markov process (Franzke et al., 2008). The transition matrix, hidden states and conditional distributions can be estimated simultaneously.

The main drawback of the weather regime approach is that it is biased toward preferred or frequent states. Quasi-stationary but rare flow patterns that fall outside the range of the major regimes may therefore be missed. Additionally, certain patterns may fall in between regimes and are not classified or are misclassified. Despite this, weather regimes are very useful because they transform complex multidimensional systems into categorized, one-dimensional series (according to which regime the system is closest to at each time step).

Figure 5 Summer weather regimes computed by k-means clustering from daily Z500 fields over the North Atlantic–European sector from 1950 to 2003 (data from NCEP/NCAR reanalysis). To eliminate transient and ambiguous episodes, only sequences of 5 d or more occupied by the same regime are retained (hence, total regime frequency sums to ≈ 70 %). Reproduced from Cassou et al. (2005). © American Meteorological Society. Used with permission.

4.2.2 Quasi-stationary states

One major disadvantage of weather regimes is that they strongly simplify the state space. They also focus on the patterns that account for most of the variability in the data, regardless of their persistence. They may thus miss rare yet impact-relevant persistent patterns. Several other techniques exist to directly extract quasi-stationary patterns from the system trajectory.

Quasi-stationary circulation patterns can be seen as mathematically quasi-stationary solutions of the atmosphere's equations of evolution (Mo and Ghil, 1987). Such solutions are characterized by average system time derivatives close to zero, meaning that the system tends to remain in their vicinity longer than elsewhere in the state space – hence their link to quasi-stationary states. Quasi-stationary states can therefore be directly identified from the system's dynamics by looking for states whose time derivative is close to zero. If we know the system's exact evolution equations (in the case of simplified models, for example), strictly stationary solutions can be directly computed (e.g., Charney and DeVore, 1979; Legras and Ghil, 1985; Mo and Ghil, 1987). In practice, however, this is rarely the case, and quasi-stationary (or “metastable”) states are defined in a statistical sense only as those for which system tendencies (i.e., time derivatives) approach zero (Vautard, 1990):

$$\begin{array}{}\text{(11)}& \left\{{\mathit{x}}^{*}\mathrm{such}\mathrm{that}\mathcal{T}\left({\mathit{x}}^{*}\right)\approx \mathrm{0}\right\}.\end{array}$$

𝒯(x^*) is the composite tendency at x^*, defined as the average (or area-average for multidimensional fields) of instantaneous tendencies at all occurrences of x^* – in practice at all times when the system is in the neighborhood of x^*:

$$\begin{array}{}\text{(12)}& \mathcal{T}\left({\mathit{x}}^{*}\right)=\sum _t{\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}}{|}_{\mathit{x}\left(t\right)}\mathbb{1}\mathit{\left\{}\mathrm{d}\right({\mathit{x}}^{*},\mathit{x}\left(t\right))\le d_{\mathrm{0}}\mathit{\}},\end{array}$$

where $\mathrm{d}(\cdot ,\cdot )$ is a similarity metric and d₀ some small threshold. It is also possible to weigh the terms in Eq. (12) according to $\mathrm{d}({\mathit{x}}^{*},\mathit{x}(t\left)\right)$ (Vautard, 1990; Michelangeli et al., 1995). When dealing with complex multidimensional systems, it is simpler to calculate tendencies on the time series of the system's leading principal components. Equation (12) can then more easily be solved by minimizing $\left|\mathcal{T}\right(\mathit{x}{)}^{\mathrm{2}}|$ in x, with distinct solutions corresponding to different quasi-stationary states. In practice, instantaneous tendencies at any x^* can exhibit a large variance, and time series are commonly low-pass filtered to remove the short-term noise that would complicate solving for Eq. (12). Additionally, $\left|\mathcal{T}\right(\mathit{x}{)}^{\mathrm{2}}|$ can be difficult to minimize as it is piecewise constant (due to the finite number of observed x values). Only approximate solutions are achievable. It may thus be difficult to know precisely how many zeros of $\left|\mathcal{T}\right(\mathit{x}{)}^{\mathrm{2}}|$ exist and whether two approximate solutions correspond to the same minimum. Vautard (1990) presents a method to select relevant solutions.

It is important to note that even if $T_i\left({\mathit{x}}^{*}\right)\approx \mathrm{0}$, instantaneous tendencies $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}{|}_{\mathit{x}\left(t\right)}$ in Eq. (12) may be far from zero (Vautard, 1990). Indeed, the time derivative of x(t) usually depends not only on x but also on other time-dependent variables that may evolve independently from x.

Haines and Hannachi (1995) and Hannachi (1997) estimated quasi-stationary states over the North Pacific in the output from a global climate forced by perpetual January conditions, by projecting simplified dynamics (e.g., quasi-geostrophy) onto the leading modes of variability in the global climate model simulation. Vautard (1990) found four main quasi-stationary patterns in wintertime daily 700 hPa geopotential height fields over the Atlantic and analyzed their quasi-stationarity and onset and break characteristics. Michelangeli et al. (1995) also looked at 700 hPa geopotential fields and compared quasi-stationary states with the leading EOFs over the Atlantic and Pacific oceans during winter. Mo and Ghil (1987) did a similar comparison for the Southern Hemisphere circulation during the austral winter.

Figure 6 shows five quasi-stationary Z500 anomaly patterns over Europe in summer. We obtained them from the 10 leading EOFs of the 5 d averaged Z500 fields, following Vautard (1990). The patterns were obtained by clustering the resulting approximate solutions with the highest number of clusters for which the resulting patterns were subjectively different. The solutions include blocking-type patterns over Scandinavia and western Russia (Fig. 6a, b) and zonal patterns (Fig. 6c, d) that are similar to the results of Figs. 5 and 11 and to the canonical patterns of variability over the European–Atlantic sector (Grams et al., 2017).

Figure 6 Five quasi-stationary states of JJA 5 d averaged rescaled Z500.

4.2.3 Optimally persistent patterns

The technique of optimally persistent patterns (OPPs) was introduced by DelSole (2001). OPP analysis is conceptually similar to principal component analysis; however, instead of patterns that maximize the variance in the observed series, OPPs are defined as the patterns whose time components (i.e., the projection of the observed series onto the OPP) are the most persistent. Persistence is here measured by the decorrelation time (Eq. 3) or, alternatively, by the squared decorrelation time $T_{\mathrm{2}}=\mathrm{1}+\mathrm{2}{\sum }_{k=\mathrm{1}}^N\mathit{\rho }(k{)}^{\mathrm{2}}$ DelSole (2001) (with ρ(k) the autocorrelation function and N the number of observations). If the dynamical system x(t) is m-dimensional, then for an m-dimensional pattern u, the time component of u is defined as y(t)=u^Tx(t). If ρ_y(k) is its associated autocorrelation function, finding OPPs then consists of maximizing $\mathrm{1}+\mathrm{2}{\sum }_{k=\mathrm{1}}^N{\mathit{\rho }}_y\left(k\right)$ or $\mathrm{1}+\mathrm{2}{\sum }_{k=\mathrm{1}}^N{\mathit{\rho }}_y^{\mathrm{2}}\left(k\right)$. In the first case (maximizing the decorrelation time), the optimization reduces to an eigenvalue problem, details of which can be found in Hannachi (2021). The leading eigenvector maximizes the decorrelation time. OPPs can then be obtained by projecting the observed data x(t) onto the time series associated with each eigenvector. The second case (maximizing the squared decorrelation time) leads to a more complex nonlinear optimization problem that can be solved iteratively (see details in DelSole, 2001). Note that, as in Sect. 4.2.2, the system should first be embedded into a lower-dimensional space, e.g., by projecting it onto the set of leading EOFs. Figure 7 shows the leading OPP obtained by DelSole (2001) by maximizing the squared decorrelation time for daily Northern Hemisphere Z500 fields.

OPPs are especially relevant for forecasting since they correspond to the patterns with the most low-frequency variability. However, like EOFs and quasi-stationary states, OPPs are mathematical objects that do not necessarily correspond to “real” patterns ever attained by the system.

Figure 7 The pattern associated with the leading eigenvector maximizing the squared decorrelation time for daily anomaly fields of 500 hPa geopotential height for the 1950–1999 period (data from NCEP/NCAR reanalysis; units in m). Adapted from DelSole (2001). © American Meteorological Society. Used with permission.

Other methods have been proposed to identify persistent patterns in a spatiotemporal field. One consists of minimizing the one-step-ahead forecast error, where the forecast is obtained by projecting the observed field x(t) onto a set of temporal patterns. These patterns, named “predictive oscillation patterns” by Kooperberg and O'sullivan (1996), can thus be considered the “most predictable” (i.e., persistent) patterns of a spatiotemporal field. Another possibility is to minimize the interpolation error variance between x(t) and its estimate $\widehat{\mathit{x}}\left(t\right)$ obtained from all $\left(\mathit{x}\right(t^{\prime }){)}_{t^{\prime }\ne t}$ values. In this case one speaks of optimally interpolated patterns (Hannachi, 2008). Further details for both methods can be found in Hannachi (2021).

4.2.4 Extreme values and dynamical systems theory

The quasi-stationary state method only focuses on the most quasi-stationary states of a system, and OPPs characterize quasi-stationarity for a handful of possible patterns only. By contrast, dynamical systems theory provides a convenient framework to describe the quasi-stationarity of any system state (Lucarini et al., 2016). In this framework, quasi-stationarity is defined for any point x₀ of the state space as the inverse of the average residence time of trajectories around x₀. The residence time is calculated based on the distance between successive system values, with two successive values deemed similar if their distance is below some small threshold. For any state x₀, the probability that x(t) will remain in a close neighborhood around x₀ (a ball of radius ϵ) can be approximated by Faranda et al. (2017a) and Faranda et al. (2017b)

$$\begin{array}{}\text{(13)}& \mathbb{P}\left(\mathrm{d}\left(\mathit{x}\right(t),{\mathit{x}}_{\mathrm{0}})\le \mathit{ϵ}\right)\simeq \exp \left\{-\mathit{\theta }\left({\mathit{x}}_{\mathrm{0}}\right){\displaystyle \frac{\mathrm{d}\left(\mathit{x}\right(t),{\mathit{x}}_{\mathrm{0}})-\mathit{\mu }\left({\mathit{x}}_{\mathrm{0}}\right)}{\mathit{\sigma }\left({\mathit{x}}_{\mathrm{0}}\right)}}\right\},\end{array}$$

where $\mathrm{d}(\cdot ,\cdot )$ is a distance function (Euclidean distance in Faranda et al., 2017a, b) and θ is called the extremal index. ϵ is usually chosen to be the second percentile of the d(x(t),x₀) values. θ can be estimated in various ways (Hamidieh et al., 2009; Holešovský and Fusek, 2022), with two common ones being the intervals estimator of Ferro and Segers (2003) and the gap estimator of Süveges (Faranda et al., 2017a).

In extreme value statistics, θ is the inverse average duration of consecutive sequences of extreme events and is used to cluster events (Ferro and Segers, 2003). A large θ therefore indicates that event occurrences tend to be isolated, while a low θ indicates that events occur as part of a larger group. In dynamical systems theory, quasi-stationarity is then defined as $\mathrm{1}/\mathit{\theta }$ (some studies directly use θ as a measure of quasi-stationarity e.g., Franzke, 2013). A perfectly stationary point of the system (where $\frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}=\mathrm{0}$) has infinite stationarity. By contrast, if trajectories immediately leave the neighborhood of x₀, the stationarity is equal to 1. In practice, long trajectories of x(t) are required to explore all possible states of the state space (also called the set of attractor points or simply the attractor) and to best approximates sequences of states on the attractor. Note that the choice of ϵ constrains the values that θ can reach. In particular, very small values (necessary for Eq. 13 to hold) usually lead to $\mathrm{1}/\mathit{\theta }$ being around 1–2 d (e.g., Faranda et al., 2017b, a; Holmberg et al., 2023), meaning that the resulting persistence metric is highly local in time and not necessarily relevant for S2S timescales.

The advantage of the method is that it is grounded in mathematical theory and, unlike approaches based on, e.g., transition probabilities, it does not require categorizing the data. It also provides an easily interpretable index of quasi-stationarity for each system state (the average residence time of the system's trajectory around that state). Holmberg et al. (2023) have, for instance, taken this approach to analyze the link between atmospheric circulation persistence and warm spells in Europe. Faranda et al. (2019) averaged θ values across all time steps to quantify global persistence in the North Atlantic flow. Going back to our example of the European summer circulation, the three major persistent Z500 anomaly patterns are obtained by clustering the top 3 % (115) most persistent daily patterns (Fig. 8). We chose five cluster centers empirically, after testing for 2–15 centers and selecting the highest number of clusters for which the resulting patterns were subjectively different. As in Fig. 6, persistent states include blocking regimes (over Scandinavia, western Russia and the North Sea; Fig. 8a, b, c), a zonal regime (Fig. 8d) and a European ridge regime (Fig. 8e).

Figure 8 Most persistent Z500 anomaly patterns during JJA, obtained by clustering the top 3 % (115) daily patterns with the highest quasi-stationarity index $\mathrm{1}/\mathit{\theta }$. The clustering algorithm is partitioning around medoids (PAM). Also shown are the associated 200 hPa zonal wind anomalies (contours; same levels as Z500). The frequency of each pattern among the 115 most persistent ones is indicated in the top left-hand corner of each panel.

4.2.5 Residence times

A common way to quantify quasi-stationarity in time series relies on the concept of “residence time”. Residence times extend the dynamical systems approach of the previous section. For a continuous series, the residence time R at time t is defined as the time during which the system remains similar to its value in t:

$$\begin{array}{}\text{(14)}& R\left(t\right)=\underset{k\ge \mathrm{0}}{max}\left\{\mathrm{d}\left(\mathit{x}\left(t\right),\mathit{x}(t+k^{\prime })\right)\le \mathit{ϵ}\forall \mathrm{0}\le k^{\prime }\le k\right\},\end{array}$$

where $\mathrm{d}(\cdot ,\cdot )$ is a similarity metric and ϵ is a small threshold. For a categorical series, R is similarly defined as the time during which the system remains in its state x(t) before transitioning to another state:

$$\begin{array}{}\text{(15)}& R\left(t\right)=\underset{k\ge \mathrm{0}}{max}\left\{\mathit{x}(t+k^{\prime })=\mathit{x}\left(t\right)\forall \mathrm{0}\le k^{\prime }\le k\right\}.\end{array}$$

The definition can be relaxed to allow for brief interruptions in a sequence of similar states (one “average” day between two sequences of warm or wet days, for instance) (e.g., Ali et al., 2021). This definition can be extended to sets $\stackrel{\mathrm{̃}}{\mathcal{S}}\subset \mathcal{S}$ of several system states: $R\left(t\right)={max}_{k\ge \mathrm{0}}\left\{\mathit{x}(t+k^{\prime })\in \stackrel{\mathrm{̃}}{\mathcal{S}}\forall \mathrm{0}\le k^{\prime }\le k\right\}$ (Richardson et al., 2019). Note that Eqs. (14)–(15) take an Eulerian perspective; their parallel in the Lagrangian perspective is the concept of “survival time”, i.e., the duration of a specific weather system or pattern along its trajectory (Liu et al., 2018; von Lindheim et al., 2021).

The residence time approach can characterize episodic, state and global quasi-stationarity and is particularly useful for predictability and risk assessment (De Luca et al., 2019; Francis et al., 2020; Berkovic and Raveh-Rubin, 2022). From the time perspective, large values of R(t) indicate the most quasi-stationary periods. A minimum threshold is often defined to separate quasi-stationary from non-stationary periods: for instance, 2 d for weather regimes (De Luca et al., 2019; Francis et al., 2020) and extreme precipitation (Du et al., 2022), 5–6 d for warm spells (Berkovic and Raveh-Rubin, 2022; Rousi et al., 2022b), 5–25 d for geopotential anomalies (Dole and Gordon, 1983), or two seasons for droughts (Ford and Labosier, 2014).

In the state space, the quasi-stationarity of any state x₀ can be described from the distribution of its residence times $R_{{\mathit{x}}_{\mathrm{0}}}\sim \mathbb{P}\left\{R\left(t\right)\mathrm{s}.\mathrm{t}.,\mathit{x}\left(t\right)={\mathit{x}}_{\mathrm{0}}\mathrm{and}\mathit{x}(t-\mathrm{1})\ne {\mathit{x}}_{\mathrm{0}}\right\}$ (Fig. 9). The easiest is to calculate the mean (Kyselý and Domonkos, 2006; Kučerová et al., 2017; Richardson et al., 2019), maximum (Rousi et al., 2022b) or some extreme percentile (Zolina et al., 2013) of $R_{{\mathit{x}}_{\mathrm{0}}}$. Like autocorrelation, it is also possible to characterize quasi-stationarity by the dependence of $\mathbb{P}\left(R_{{\mathit{x}}_{\mathrm{0}}}=n\right)$ on n (Sharma and Panu, 2014). Light distribution tails indicate an exponential decay and short-term memory (see Sect. 4.2.6), while heavy tails (e.g., power law (Bunde et al., 2013) or q-exponential (Weber et al., 2019) distributions) point to long-range dependence in the system. For example, Pfleiderer and Coumou (2018) fitted exponential models to the distribution of consecutive warm and cold spells to quantify temperature persistence in the Northern Hemisphere, and Huguenin et al. (2020) did the same for weather types over Central Europe. Many other types of distribution can also be fitted (e.g., Deni et al., 2010; Zolina et al., 2013). Residence times can also be modeled as functions of covariates. For instance, Röthlisberger et al. (2019) and Ali et al. (2021) apply a Weibull regression model to relate the duration of dry and wet spells to Rossby wave activity in the midlatitudes.

Finally, averaging R(t) over time t, or $R_{{\mathit{x}}_{\mathrm{0}}}$ across all possible system states, provides a convenient and easily interpretable global quasi-stationarity metric (e.g., Pfleiderer et al., 2019).

Figure 9 PDF of the North Atlantic jet latitude index (solid) together with the weighted Gaussian PDFs from the HMM: the southern (dashed), northern (dotted) and central (dashed–dotted) regimes. Regime duration curves (southern (black), northern (red) and central (blue)), expressed as the frequency of occurrences lasting at least n days. Dashed curves show the corresponding 2.5 % and 97.5 % confidence levels. Reproduced from Franzke et al. (2011). © American Meteorological Society. Used with permission.

4.2.6 Transition probabilities

Residence times characterize how long the system remains in the same state x₀. By contrast, transition probabilities P describe the likelihood that the state of the system does not change at the next time step:

$$\begin{array}{}\text{(16)}& P\left({\mathit{x}}_{\mathrm{0}}\right)=\mathbb{P}\left(\mathit{x}(t+\mathrm{1})={\mathit{x}}_{\mathrm{0}}\left|\mathit{x}\right(t)={\mathit{x}}_{\mathrm{0}}\right).\end{array}$$

Transition probabilities can be further divided by marginal probabilities to yield probability ratios (e.g., Kolstad et al., 2015):

$$\begin{array}{}\text{(17)}& PR\left({\mathit{x}}_{\mathrm{0}}\right)={\displaystyle \frac{P\left({\mathit{x}}_{\mathrm{0}}\right)}{\mathbb{P}\left(\mathit{x}\left(t\right)={\mathit{x}}_{\mathrm{0}}\right)}}.\end{array}$$

A high transition probability P(x₀) or PR(x₀) means that the system likely remains in x₀, which indicates persistence (Fig. 10a). Like residence times, transition probabilities can be computed on specific system states only or averaged across all system states (according to their frequency) to yield a measure of global quasi-stationarity. Their statistical significance can be estimated by comparing observed transition probabilities to those obtained from a large set of randomly shuffled series. However, unlike residence times, transition probabilities cannot identify quasi-stationary periods: they only work in the state space.

Økland and Lejenäs (1987) applied this method to blocking quasi-stationarity; Ford and Labosier (2014), Wilby et al. (2016), and Moon et al. (2018) to droughts; Fereday (2017) to large-scale sea-level pressure patterns; Guilbert et al. (2015) to daily precipitation; and Röthlisberger and Martius (2019) to dry and warm periods (Fig. 10a).

Figure 10 (a) Climatological May–October (MJJASO) hot-spell survival probabilities (shaded) in ERA-Interim reanalysis data. Red contours show MJJASO cyclone frequencies of 30 % and 40 %. Reproduced from Röthlisberger and Martius (2019). © 2019 American Geophysical Union. All Rights Reserved. (b) Change in the odds ratio of drought occurrence in spring given a single-unit increase in the winter standardized precipitation index estimated from logistic regression. Thick black contours indicate statistical significance with 95 % confidence. Adapted from Ford and Labosier (2014). © 2013 Royal Meteorological Society.

It is often convenient to assume that the distribution of system states at t+1 only depends on the system state at t, in other words to approximate the system with a first-order Markov chain model (Sericola, 2013). Under this assumption, transition probabilities are directly related to the distribution of residence times. Indeed, if the probability of remaining in state x₀ does not depend on how long the system has previously been in state x₀: $\mathbb{P}\left(\mathit{x}(t+\mathrm{1})={\mathit{x}}_{\mathrm{0}}\left|\mathit{x}\right(t)={\mathit{x}}_{\mathrm{0}}\right)=\mathit{\alpha }$, then the probability that the residence time of x₀, $R_{{\mathit{x}}_{\mathrm{0}}}$ exceeds n is equal to the probability of finding a sequence of at least n consecutive x₀ values, which specifically leads to the following equation:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathbb{P}\left(R_{{\mathit{x}}_{\mathrm{0}}}\ge n\right)& ={\mathbb{P}}_t\left(\mathit{x}\right(t-\mathrm{1})\ne {\mathit{x}}_{\mathrm{0}};\mathit{x}(t)={\mathit{x}}_{\mathrm{0}};\mathit{x}(t+\mathrm{1})\\ & ={\mathit{x}}_{\mathrm{0}};\mathrm{\dots };\mathit{x}(t+n-\mathrm{1})={\mathit{x}}_{\mathrm{0}}),\end{array}\end{array}$$

where the subscript in ℙ_t specifies that the probability is taken with respect to the time variable, whereas in $\mathbb{P}\left(T_{{\mathit{x}}_{\mathrm{0}}}\ge n\right)$ it is taken with respect to the ensemble of residence times for x₀. Thanks to the first-order Markov assumption, Eq. (18) factors as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\mathbb{P}\left(R_{{\mathit{x}}_{\mathrm{0}}}\ge n\right)& ={\mathbb{P}}_t\left(\mathit{x}(t-\mathrm{1})\ne {\mathit{x}}_{\mathrm{0}};\mathit{x}\left(t\right)={\mathit{x}}_{\mathrm{0}}\right)\\ & \times \prod _{k=\mathrm{0}}^{n-\mathrm{2}}{\mathbb{P}}_t\left(\mathit{x}(t+k+\mathrm{1})={\mathit{x}}_{\mathrm{0}}\left|\mathit{x}\right(t+k)={\mathit{x}}_{\mathrm{0}}\right).\end{array}\end{array}$$

Thus,

$$\begin{array}{}\text{(20)}& \mathbb{P}\left(R_{{\mathit{x}}_{\mathrm{0}}}\ge n\right)=\mathit{\beta }\times {\mathit{\alpha }}^{n-\mathrm{1}},\end{array}$$

where $\mathit{\beta }={\mathbb{P}}_t\left(\mathit{x}(t-\mathrm{1})\ne {\mathit{x}}_{\mathrm{0}};\mathit{x}\left(t\right)={\mathit{x}}_{\mathrm{0}}\right)$ does not depend on t for a (statistically) stationary system, and

$$\begin{array}{}\text{(21)}& \mathbb{P}\left(R_{{\mathit{x}}_{\mathrm{0}}}=n\right)=\mathit{\beta }{\mathit{\alpha }}^{n-\mathrm{1}}(\mathrm{1}-\mathit{\alpha }).\end{array}$$

The logarithm of the distribution of residence times is therefore a linear function of n, with the slope (log (α)) equal to the transition probability. The Markov assumption is thus consistent with a short-term dependent system in which residence time probabilities decay exponentially with n (see Sect. 4.2.5 and Bunde et al., 2013). For instance, Huguenin et al. (2020) fit such exponential laws and use their slope as measures of weather regime quasi-stationarity over Central Europe in current and future climates. This approach is well suited to comparing two different series (e.g., two different periods or locations) because it separates changes in the marginal frequency of state x₀ (only affecting β) from changes in quasi-stationarity (only affecting α).

4.2.7 Time series modeling

Another possibility to investigate state quasi-stationarity is to fit a statistical model to the data that links successive system values with each other:

$$\begin{array}{}\text{(22)}& \mathit{x}(t+\mathrm{1})\textcolor{red}{}\mathcal{F}({\left\{\mathit{x}\left(k\right)\right\}}_{k=\mathrm{1},\mathrm{\dots },t},{\left\{\mathit{\theta }\left(t\right)\right\}}_{k=\mathrm{1},\mathrm{\dots },t}),\end{array}$$

where stands for equality in distribution and {θ(t)}_k=1…t are covariates. Transition probabilities (Sect. 4.2.6), for example, are equivalent to modeling a system's evolution by a transition matrix between all possible system states. Hidden Markov models (Sect. 4.2.1) likewise estimate transition matrices between hidden states. Logistic regression can be used to assess quasi-stationarity in a given state. Several studies used it to quantify quasi-stationarity in droughts, using previous precipitation or soil moisture anomalies as predictor variables (θ(t) in Eq. 22) (Ford and Labosier, 2014; Meng et al., 2017) (Fig. 10b). The link to drought occurrence at the previous time step is implicitly taken into account with previous precipitation or soil moisture anomalies. If previous system values are explicitly included as covariates to the model, one speaks of auto-logistic regression (Wolters, 2017).

4.3 Running window methods for episodic persistence

The running window approach identifies persistent periods of a given fixed length in continuous or categorical data. It assesses the degree of quasi-stationarity over fixed time intervals by means of a “similarity index”. Given a time interval $\mathcal{T}=\left\{t_{\mathrm{1}},\mathrm{\dots },t_n\right\}$ of length n and a similarity metric $d(\cdot ,\cdot )$, the similarity index S_n(𝒯) is defined as the average similarity between system values in 𝒯. This similarity can be computed with reference to one time step x(t_k) in particular (typically the first one):

$$\begin{array}{}\text{(23)}& S_n\left(\mathcal{T}\right)={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{\mathrm{1}\le k^{\prime }\le n{k}^{\prime }\ne k}\mathrm{d}\left(\mathit{x}\left(t_{k^{\prime }}\right),\mathit{x}\left(t_k\right)\right),\end{array}$$

or it can measure the average similarity across all pairs of values:

$$\begin{array}{}\text{(24)}& S_n\left(\mathcal{T}\right)={\displaystyle \frac{\mathrm{2}}{n(n-\mathrm{1})}}\sum _{\mathrm{1}\le k,k^{\prime }\le nk\ne k^{\prime }}\mathrm{d}\left(\mathit{x}\left(t_k\right),\mathit{x}\left(t_k^{\prime }\right)\right).\end{array}$$

Different time intervals can then be ranked according to their degree of quasi-stationarity. Mo and Ghil (1987), for instance, use pattern correlation as similarity metric and set a threshold of 0.5 to separate quasi-stationary from non-stationary periods. Another possibility is to keep the N periods with the largest S_n(𝒯) (i.e., to use a percentile-based threshold). Quasi-stationary states can then be detected with, e.g., pattern recognition techniques applied to the most quasi-stationary periods (Mo and Ghil, 1987). We show an example for summertime European circulation in Fig. 11. The most common are blocking-type patterns over northern Europe and western Russia (Fig. 11a, b, c), but we also find two pronounced zonal patterns with a southward- or northward-shifted jet (Fig. 11d, e).

Averaging S_n(𝒯) values across multiple periods of length n can additionally provide for a measure of global quasi-stationarity, as Hoffmann et al. (2021) did to analyze 10 d persistence in summer atmospheric circulation. When working with categorical data, similarity indices can be averaged for each system state to highlight the most quasi-stationary ones (Horel, 1985b).

The running window approach is well suited to compute temporal trends in quasi-stationarity, since the similarity index can be defined at all time steps (Hoffmann et al., 2021). One limitation, however, is that results for the same system but calculated for different period lengths are not directly comparable since the marginal distribution of S_n(𝒯) depends on n. Furthermore, S_n(𝒯) only measures the average similarity between successive values. A high S_n(𝒯) does not guarantee a high degree of quasi-stationarity because there could be breaks in between sequences of similar values. This is especially relevant for long periods during which strict persistence is unlikely to occur.

While we try to give a comprehensive view of common methods used in quasi-stationarity analysis, we could not include all the methods that exist in the literature. Kornhuber and Tamarin-Brodsky (2021), for instance, define quasi-stationarity in a Lagrangian context as the zonal velocity of individual weather systems (in their case, localized temperature anomalies), while Hoskins and Woollings (2015) use Rossby wave phase speed as a proxy for weather quasi-stationarity in the midlatitudes. Finally, in the framework of dynamical systems theory, quasi-stationarity can also be equated to the system remaining for some period of time in a small subset of the state space. During that period, the system may explore different configurations, but fewer configurations are explored than if it had been able to evolve freely across the whole state space (Fereday, 2017). In this sense, persistent states are states from which the system takes a (statistically) long time to reach the rest of its state space. This perspective relates to the intransitivity theory of Lorenz (1990), which postulates that the evolution of the atmospheric circulation on S2S timescales can, under certain conditions, be governed by a few well-separated attractors, each with their own preferred states (Weiland et al., 2021).

Table 1 Overview of quasi-stationarity methods: method name, corresponding section, type of persistence it applies to, a brief definition, main features and main limitations.

Figure 11 (a–e) Most persistent 10 d averaged Z500 anomaly patterns during JJA, obtained by clustering the top 10 % (40) 10 d patterns with the highest similarity index (Eq. 24). Similarity is measured with the SSIM index. The clustering algorithm is partitioning around medoids (PAM). The associated 200 hPa zonal wind anomalies are shown by contour lines (same levels as Z500). The frequency of each pattern is also indicated. (f–j) 10 d averaged surface temperature anomalies associated with each persistent Z500 pattern.

5 Recurrence

We now discuss methods to capture recurrence in atmospheric data (see Table 2 for an overview). In contrast to quasi-stationarity, most methods here only quantify recurrence in a specific system state. They typically require binary time series as input (where “1” represents an occurrence of the system state or event of interest) and count event occurrences over time. Recurrence is therefore always defined for a given state or event, usually with the potential for large impacts like tropical (Mumby et al., 2011) and extra-tropical cyclones (Mailier et al., 2006; Vitolo et al., 2009), windstorms (Khare et al., 2015), or extreme precipitation episodes (Kopp et al., 2021). For this reason, we did not explicitly divide this section into global, state and episodic methods. Common methods are available to characterize episodic persistence (Sect. 5.1.1) and the persistence of given states (Sect. 5.1.2–5.2.1). However, few definitions of “global recurrence” exist. Additionally, no study has, to our knowledge, tried to analyze S2S recurrence in a “non-supervised” way, i.e., that would allow us to objectively identify recurrent states. In Sect. 5.2.2, we discuss the potential of recurrence plots to fill that gap.

Because most recurrence methods work directly with binary time series, point processes are a convenient theoretical tool in recurrence analysis. Point processes are a class of probability models for the random occurrence of points in a space (one- or multi-dimensional space). The most simple hypothesis that can be made is that of complete serial randomness, i.e., points occur completely independently of one another. In this case, recurrence occurs by chance only. The homogeneous Poisson point process is a simple model without memory for binary series in which events occur without correlation and with a constant intensity. In a homogeneous Poisson process, the number of points in individual time intervals are independent, and the number of points in an interval of length τ is Poisson distributed with the rate λτ, where λ is the process intensity rate. This model is useful to simulate binary series with complete serial randomness, and therefore provides an empirical basis to assess the statistical significance of various recurrence metrics in observations.

Separating distinct occurrences of a state, pattern or event is critical to distinguish recurrence from quasi-stationarity. Over short periods, like a season, it is in principle possible to visually separate multiple systems occurring in close succession, for instance with Hovmöller diagrams (e.g., Tuel et al., 2022a; Rousi et al., 2022a). Most of the time, however, separation cannot be done by hand, and objective criteria are required. One option is to use physically based detection and tracking algorithms – e.g., for cyclones (Hodges, 1995; Wernli and Schwierz, 2006) or atmospheric blocks (Schwierz et al., 2004; Steinfeld, 2021) – that can separate distinct weather systems. A second option is to use a minimum duration between events to classify them as distinct, with the choice depending on the underlying physical system or on impact considerations. Barton et al. (2016), Kopp et al. (2021) and Tuel and Martius (2021a), for instance, require a minimum of 2 d with non-extreme precipitation between extreme precipitation episodes (a typical timescale for extratropical cyclones). In Pfleiderer et al. (2019), a single day is used to separate distinct dry or warm periods.

5.1 Diagnostic methods

5.1.1 Window counts

A simple way to identify recurrent periods is to look for periods with high event counts. One possibility is to set a time window and require the number of event counts during this window to exceed some threshold (two or higher; Fig. 12a). Kopp et al. (2021) and Tuel and Martius (2022a) take this approach to analyze recurrent extreme precipitation episodes at the 2- and 3-week timescales. Pinto et al. (2014) likewise identify cyclone clusters over the British Isles as periods with four or more consecutive cyclones within 7 d. Another possibility is to look for clusters of events in which each event is separated by at most n time steps from the preceding one. For instance, Bevacqua et al. (2020) define cyclone clusters as sequences of N≥2 consecutive cyclones transiting within a given area and separated by 24 h or less. This method is largely impact driven and does not aim at modeling the binary series, assessing the significance of the clustering or calculating return periods for large event counts (Dacre and Pinto, 2020). For this, further steps must be taken.

Figure 12 (a) Example of a time series of extreme precipitation events (black dots) and 21 d periods containing more than two events (gray shading). Adapted from Kopp et al. (2021) under the terms of the Creative Commons Attribution 4.0 License. (b) Map of dispersion statistic ψ (Eq. 25) of December–January–February (DJF) cyclones over 1975–2004, averaged across an ensemble of CMIP5 models. Reproduced from Bevacqua et al. (2020) under the terms of the Creative Commons Attribution 4.0 license.

5.1.2 Dispersion metrics

Dispersion metrics characterize the distribution of event occurrences in a time series. In a series with no tendency to clustering, high event counts over a given duration should be much less frequent than in a series where clustering is prevalent. If events occur with the same average probability in the two series, then in the clustered series, sequences with no events should mechanically be more frequent, and the variance in event counts will be higher. In a homogeneous Poisson process (see above) N_k(τ), event counts in disjoint intervals of length τ (indexed by k) follow a Poisson distribution with parameter λτ. Hence, the expected value μ(τ) and variance σ(τ)² of N_k(τ) are both equal to λτ. Two common statistical dispersion metrics quantify deviations from this homogeneous behavior. The dispersion statistic is a measure of temporal correlation over different timescales. It is defined as Mailier et al. (2006)

$$\begin{array}{}\text{(25)}& \mathit{\psi }\left(\mathit{\tau }\right)={\displaystyle \frac{\mathit{\sigma }(\mathit{\tau }{)}^{\mathrm{2}}}{\mathit{\mu }\left(\mathit{\tau }\right)}}-\mathrm{1}.\end{array}$$

ψ is related to the Fano factor ($\mathrm{FF}\left(\mathit{\tau }\right)=\mathit{\psi }\left(\mathit{\tau }\right)+\mathrm{1}$) (Telesca, 2007), also called index of dispersion (Mailier et al., 2006). A positive ψ(τ) (σ(τ)²>μ(τ)) indicates over-dispersion, meaning that events tend to occur in clusters. A negative ψ(τ) points to under-dispersion, meaning that events occur more regularly than in a random process. Near-zero values are consistent with a homogeneous Poisson process (Fig. 12b).

Similarly, the Allan factor (AF) is defined as the variance of successive event counts over an interval of length τ divided by twice the average event count in τ steps (Telesca, 2007):

$$\begin{array}{}\text{(26)}& \mathrm{AF}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathbb{E}\left[\right(N_{k+\mathrm{1}}\left(\mathit{\tau }\right)-N_k\left(\mathit{\tau }\right){)}^{\mathrm{2}}]}{\mathrm{2}\mathit{\mu }\left(\mathit{\tau }\right)}}.\end{array}$$

For a homogeneous Poisson process, AF(τ)=1, while AF(τ)>1 indicates clustering behavior. ψ and AF can be evaluated for a single window τ^*, but assessing how they evolve with τ helps better detect recurrent behavior in the series (Serinaldi and Kilsby, 2013). In the absence of clustering, ψ(τ) and AF(τ) remain constant with τ. In clustered processes, however, both metrics scale with τ, often with a power law behavior over some range [τ₁,τ₂], which can be described by an exponent (Telesca, 2007) similar to the spectral analysis of persistence (Sect. 4.1.2). One advantage of ψ and AF is that they are insensitive to the underlying event frequency, which makes them convenient to compare different locations or event identification algorithms with each other (Dacre and Pinto, 2020). However, the corresponding null hypothesis – that the absence of recurrence is equivalent to exponentially distributed event return times under the Poisson assumption – is not always relevant, and other approaches may be better suited to detect deviations from memoryless processes (Blender et al., 2015).

The statistical significance of ψ(τ) and AF(τ) can be assessed empirically through Monte Carlo sampling of random series without memory or analytically. Mailier et al. (2006) discuss a chi-squared test for ψ(τ) (for a fixed window τ), while Serinaldi and Kilsby (2013) introduce the sampling distribution of the empirical AF estimator. Note that seasonality in event counts can artificially inflate the dispersion statistic, meaning that care must be taken either to remove the seasonality (e.g., Tuel and Martius, 2021a) or to simulate random processes with the same seasonality as the observed signal (λ function of the month, for instance). Note also that ψ(τ) and AF(τ) both involve μ(τ) in the denominator; in other words, these metrics are normalized by the marginal event frequency. The same ψ(τ) can be obtained for series with different event frequencies, which makes dispersion metrics less relevant for impacts than other methods (e.g., window counts; Sect. 5.1.1).

Since its introduction by Mailier et al. (2006), the dispersion statistic has been frequently applied in recurrence analysis, namely to detect temporal or serial clustering in cyclones, both tropical (Mumby et al., 2011; Wolff et al., 2016) and extratropical (Vitolo et al., 2009; Pinto et al., 2013, 2014; Economou et al., 2015; Pinto et al., 2016). The Allan factor has been applied to daily precipitation series (Serinaldi and Kilsby, 2013), wave storms (Besio et al., 2017) and to extreme wind speeds in Switzerland (Telesca et al., 2020).

Alternative metrics that characterize the distribution of event counts have been proposed in the literature. For example, Kopp et al. (2021), working with extreme precipitation counts over 2- to 4-week timescales, proposed a new metric, S_cl(τ) – which they named the “clustering metric” – based on a weighted sum of window event counts. Like the dispersion statistic and Allan factor, S_cl(τ) requires choosing a window size τ. Clustering periods of length τ are then identified and ranked by applying a moving window sum to the binary series. The period with the most events is selected first; the corresponding events are removed from the series (to avoid any intersection between clustering periods); and the procedure continues until a minimum number of events is obtained or until no more clustering periods (with at least two events) are found. This yields a set of K periods with decreasing event counts n_k(τ): $n_{\mathrm{1}}\left(\mathit{\tau }\right)\ge n_{\mathrm{2}}\left(\mathit{\tau }\right)\ge \mathrm{\dots }\ge n_K\left(\mathit{\tau }\right)$. Kopp et al. (2021) then define their clustering metric by

$$\begin{array}{}\text{(27)}& S_{\mathrm{cl}}\left(\mathit{\tau }\right)=\sum _{\mathrm{1}\ge k\ge K}{n}_k\left(\mathit{\tau }\right)q_k,\end{array}$$

where q_k are weights that can be defined in various ways, as long as they verify (i) that q_k decreases with k (periods with fewer events are given less weight) and (ii) that $q\left(i\right)-q(i+\mathrm{1})>q(i+\mathrm{1})-q(i+\mathrm{2})$. In the end, this approach is equivalent to a scoring system for clustering periods, and the resulting metric S_cl(τ) correlates well with the dispersion statistic.

5.1.3 Ripley's K function

Given a system state x₀, Ripley's K function measures the average number of occurrences of x₀ within a given neighborhood around a random occurrence of x₀. It was originally developed for spatial data, but it can be applied in a temporal setting to characterize recurrence in x₀. Let the binary variable Y_t=𝟙{X_t is in state x₀}. Ripley's K function for a neighborhood of size n is then equal to Ripley (1981) and Dixon (2014):

$$\begin{array}{}\text{(28)}& K\left(n\right)={\mathit{\lambda }}^{-\mathrm{1}}\left(\mathbb{E}\left[\sum _{k=-n}^n{Y}_{t+k}|Y_t=\mathrm{1}\right]-\mathrm{1}\right),\end{array}$$

where λ=𝔼[Y_t] (frequency of x₀ in X_t). The statistical significance of K(n) can then be assessed by comparing it to values obtained for a homogeneous Poisson process of same length and with the same average rate of occurrence as Y_t (Fig. 13). In Eq. (28), λ⁻¹ can also be discarded, meaning that K(n) represents an average event count and is more directly related to impacts. Ripley's K was used in meteorology to characterize S2S recurrence in extreme precipitation events by Barton et al. (2016), Tuel and Martius (2021a), and Tuel and Martius (2021b). Stephenson et al. (2004) and Hannachi (2010) also used Ripley's K to characterize the clustering of climate modes in the state space.

Figure 13 (a) Assessing the significance of clustering with Ripley's K function: conceptual example of Ripley's K values on an observed series (black crosses) and the corresponding 95 % uncertainty range obtained from random simulations of series with no clustering (blue shading). In this case, observations show significant clustering at the 10–20 d timescale. Reproduced from Tuel and Martius (2021a) under the terms of the Creative Commons CC BY-NC-ND 4.0 license. (b) Map of 20 d Ripley's K values for catchment-averaged extreme precipitation during fall in Switzerland. Hatching indicates statistical significance at the 95 % level. Adapted from Tuel and Martius (2021b) under the terms of the Creative Commons Attribution 4.0 license.

Figure 14 Effect of European–Atlantic zonal regime occurrence (corresponding 500 hPa geopotential mean and anomalies are shown in the left panel) on the temporal clustering of extreme precipitation at the 3-week timescale during winter (measured as a probability multiplicative factor). Reproduced from Barton et al. (2022) under the terms of the Creative Commons Attribution 4.0 license.

5.1.4 Distribution of inter-event times

Quasi-stationarity in a given system state x₀ can be described by the corresponding distribution of residence times (Sect. 4.2.5). Similarly, the distribution of inter-event times (also called “recurrence times” or “waiting times”) characterizes the recurrence in x₀ from a global perspective (Altmann and Kantz, 2005). The recurrence time T_i is the duration between the ith and (i+1)th occurrences of state x₀ in a time series. We can easily see that both short and long recurrence times will be more frequent in a series that exhibits clustering than in a series where events occur independently of one another. The presence of a heavy tail in the distribution of T_i therefore points to a tendency for temporal clustering in the series. If x₀ is rare enough, the distribution of T_i may converge to well-known distributions like the exponential (a sign of short-term dependence, as with homogeneous Poisson processes), stretched exponential or power law distribution (Eichner et al., 2007; Corral, 2015). When events correspond to extremes above a given threshold, the exponents that characterize these distribution are even related to the autocorrelation coefficient of the original series (Santhanam and Kantz, 2008; Kalra and Santhanam, 2021).

5.2 Stochastic modeling of recurrence

5.2.1 Event series

We have so far discussed diagnostic methods that characterize recurrence in binary time series. A different approach consists of fitting point processes to the series. A point process is a random process representing the occurrence times of specific events. It can be characterized by event times T(n) (time of the nth event in the series) or by the distribution of event counts N(t₁,t₂) for any interval [t₁,t₂]. We can define its intensity or hazard function λ(t) by

$$\begin{array}{}\text{(29)}& \mathit{\lambda }\left(t\right)=\underset{{\mathrm{\Delta }}_t\to {\mathrm{0}}^{+}}{lim}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\mathbb{P}\left(N\right(t,t+\mathrm{\Delta }t)=\mathrm{1}),\end{array}$$

and, for any t₁<t₂, its intensity measure Λ(t₁,t₂) as

$$\begin{array}{}\text{(30)}& \mathrm{\Lambda }(t_{\mathrm{1}},t_{\mathrm{2}})=\underset{t_{\mathrm{1}}}{\overset{t_{\mathrm{2}}}{\int }}\mathit{\lambda }\left(t\right)\mathrm{d}t.\end{array}$$

Many different point processes have been described in the statistics literature (e.g., Cox and Isham, 1980), but clustering analyses in weather and climate studies have mainly relied on Poisson processes. In a Poisson process, λ(t) is a deterministic function of time. If λ(t) is constant, the process is said to be homogeneous; otherwise, the process is said to be inhomogeneous. For any t₁<t₂, N(t₁,t₂) follows a Poisson distribution with rate Λ(t₁,t₂):

$$\begin{array}{}\text{(31)}& N(t_{\mathrm{1}},t_{\mathrm{2}})\textcolor{red}{}\mathrm{Poisson}\left(\mathrm{\Lambda }\right(t_{\mathrm{1}},t_{\mathrm{2}}\left)\right),\end{array}$$

and event numbers in disjoint intervals are independent from each other.

The point process approach has mainly been used to relate recurrence in extreme events to atmospheric or climate variability. To that end, λ(t) is assumed to depend on given time covariates X(t). A linear dependence corresponds to Poisson generalized linear models (GLMs):

$$\begin{array}{}\text{(32)}& \mathit{\lambda }\left(t\right)=\mathit{\beta }\cdot X\left(t\right).\end{array}$$

For instance, Villarini et al. (2011), Tuel and Martius (2022a), and Tuel and Martius (2022b) used Poisson GLMs to relate temporal clustering in heavy precipitation to large-scale modes of climate and atmospheric variability. Mailier et al. (2006) did the same for extratropical cyclones in the North Atlantic. The linear constraint can be lifted by using general additive models (GAMs) in which the relationship between λ(t) and the covariates is specified by flexible smooth functions (Hastie, 1992). Villarini et al. (2013) used GAMs to relate flood recurrence in the US Midwest to modes of climate variability like ENSO, and Barton et al. (2022) used GAMs to link temporal clustering in extreme precipitation across Europe to weather regimes in the North Atlantic. Another extension of the Poisson GLM is the Cox regression model (Smith and Karr, 1986), in which

$$\begin{array}{}\text{(33)}& \mathit{\lambda }\left(t\right)={\mathit{\lambda }}_{\mathrm{0}}\left(t\right)\exp (\mathit{\beta }\cdot X(t\left)\right),\end{array}$$

where λ₀(t) is a baseline intensity or hazard function and exp (β⋅X(t)) accounts for the time-varying effects of the selected covariates. Examples include Villarini et al. (2013), Mallakpour et al. (2017), and Yang and Villarini (2019), who applied Cox regression to recurrence in floods and extreme precipitation. More complex models are also possible; for instance, Khare et al. (2015) introduced a Poisson mixture model for windstorm clustering, in which λ(t) is expressed as a stochastic function of time modulated by a gamma distribution (Khare et al., 2015). Note that models can be fitted directly by maximizing the likelihood function (Yang and Villarini, 2019) or indirectly by working with window event counts (Mailier et al., 2006; Tuel and Martius, 2022a; Barton et al., 2022).

5.2.2 Recurrence plots

All previously discussed methods share the same limitation: they require binary time series as input, meaning that they are designed to characterize recurrence only in a given system state. It is not possible with such methods to explore recurrence in an “unsupervised” way, i.e., to automatically detect recurrent system states or sequences of states and associated recurrent periods. Recurrence plots (RPs) may offer a solution to this problem.

An RP is a graphical representation of a system's self-similarity with time (Marwan et al., 2007). Mathematically, it is defined as a two-dimensional binary matrix $\mathbf{R}\in \mathit{\{}\mathrm{0},\mathrm{1}{\mathit{\}}}^{N\times N}$, with time on both axes, that consists of pairwise comparisons of system values:

$$\begin{array}{}\text{(34)}& \mathbf{R}(i,j)=\mathrm{\Theta }\left(\mathit{ϵ}-d\left(\mathit{x}\right(t_i),\mathit{x}(t_j)\right),\end{array}$$

with $d(\cdot ,\cdot )$ a distance or similarity metric (by convention, $d\left(\mathit{x}\right(t_i),\mathit{x}(t_i)=\mathrm{0}$), Θ the Heaviside step function, and ϵ a similarity threshold. If $d\left(\mathit{x}\right(t_i),\mathit{x}(t_j\left)\right)<\mathit{ϵ}$, then x(t_j) and x(t_j) are similar, and $\mathbf{R}(i,j)=\mathrm{1}$ represents a recurrence of the system. By definition, R is symmetric and $\forall \mathrm{1}\le i\le N,\mathbf{R}(i,i)=\mathrm{1}$. A conceptual example is shown in Fig. 15a.

Figure 15 Detecting recurrence with recurrence plots. (a) Conceptual example of a recurrence plot. The red rectangle highlights a period during which the system exhibits recurrence with three successive similar 4 d sequences occurring over a 12 d period (highlighted in red). (b) Zoom over the recurrence period with the three separate 4 d sequences highlighted by dashed green lines.

Figure 16 Two major 3 d recurring sequences of atmospheric circulation (Z500) over Europe during summer obtained by clustering the ≈ 30 21 d windows with ${\mathrm{RI}}_{\mathrm{3},\mathrm{21}}\left(t\right)\ge \mathrm{2}$ with the PAM algorithm. Stippling (hatching) indicates atmospheric blocking (cyclone) frequency anomalies larger than 30 %. Similarity is measured with the SSIM index, with a similarity threshold ϵ in Eq. (34) of 0.25 (corresponding to a recurrence rate of ≈ 5 %). Note that successive panels (a–f) are 1 d apart.

Table 2 The same as Table 1 but for recurrence methods.

Before calculating R, the system x can first be reduced to a lower-dimensional space (as in Sect. 4.2.2), for example by projecting x onto a set of principal components (Marwan et al., 2007; Mukhin et al., 2022). Note that R depends on the selected similarity metric and on ϵ. For a given system x, R is thus not uniquely defined. Choosing a value for ϵ should be done carefully. ϵ should be small enough to account for real recurrences in the data but not too small to have enough recurrences to analyze. Several rules of thumb have been proposed in the literature (see Marwan et al., 2007, and references therein). For instance, ϵ can be selected by requiring an average recurrence rate of 1 %–10 %. For weather and climate data, ϵ could also be selected based on physical considerations and expert judgment.

RPs are a convenient way to visualize a system's trajectory, especially for multi-dimensional systems (like circulation fields). RPs capture persistent behavior: vertical lines indicate quasi-stationarity, while diagonal lines (outside of the main diagonal) indicate recurrence. RPs can thus highlight periods when the trajectory of a system roughly visits the same sequence of states or parts of the state space.

Several measures have been proposed to quantify the presence of specific patterns in RPs. They are known collectively as “recurrence quantification analysis” (Marwan, 2008). Common ones include the recurrence rate:

$$\begin{array}{}\text{(35)}& \mathrm{RR}={\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}\sum _{\mathrm{1}\le i,j\le N}\mathbf{R}(i,j);\end{array}$$

the determinism, which quantifies the fraction of recurrence points forming diagonal lines of minimum length l_min:

$$\begin{array}{}\text{(36)}& \mathrm{DET}={\displaystyle \frac{{\sum }_{l\ge l_{\min }}lD\left(l\right)}{{\sum }_{l\ge \mathrm{1}}lD\left(l\right)}},\end{array}$$

where D(l) is the number of diagonal lines of length l; and the laminarity, which quantifies the fraction of recurrence points forming vertical lines of minimum length l_min:

$$\begin{array}{}\text{(37)}& \mathrm{LAM}={\displaystyle \frac{{\sum }_{l\ge l_{\min }}lV\left(l\right)}{{\sum }_{l\ge \mathrm{1}}lV\left(l\right)}},\end{array}$$

where V(l) is the number of vertical lines of length l. All three can be used as global measures of persistence (e.g., Ramirez-Amaro and Figueroa-Nazuno, 2006). However, these indices do not a priori discriminate between quasi-stationarity and persistence (high values of RR, for instance, could be due to either). To do so would require making sure that similar time steps are separated by “0” in the RP (Fig. 15b).

RPs can also detect local recurrence in time. With RPs, one can extend the counting approach introduced in Sect. 5.1.1 to multiple-day sequences of complex multi-dimensional data. Specifically, to search for sequences of p steps that recur over periods of q steps ($\mathrm{0}<p<q$), we can simply define a recurrence index RI_p,q(t) as the number of diagonal lines of length p in the sub-matrix of R indexed by $[t,t+p]\times [t,t+q]$ (Fig. 15b). To make sure that distinct occurrences of the same sequence are separated, we can require diagonals to be separated by a minimum number of days (2 in Fig. 15b). It is also possible to allow for short breaks (e.g., one step) in the diagonals.

RPs have recently gained attention in a range of disciplines that deal with complex systems (Goswami, 2019) but only rarely in environmental sciences. Ramirez-Amaro and Figueroa-Nazuno (2006) investigated the recurrence properties of major teleconnection patterns with RPs. Yiou et al. (2018) indirectly relied on RPs to identify the most recurrent states of North Atlantic circulation at the intra-seasonal timescale. Adeniji et al. (2018) used RPs to analyze the recurrence characteristics of hourly wind speed in Nigeria, and Ray et al. (2019) did so to investigate daily temperature and humidity data across India. Recently, Mukhin et al. (2022) used RPs to detect weather regimes in the Northern Hemisphere. In Fig. 16, we show results obtained for 3 d sequences of similar upper-level circulation patterns over Europe in summer. We calculate RI_p,q(t) with p=3 d and q=21 d and identify the two most common sequences among the ≈30 21 d windows with ${\mathrm{RI}}_{p,q}\left(t\right)\ge \mathrm{2}$. The method identifies recurring blocks and cyclones over Scandinavia as the most frequent recurring sequences of circulation patterns. The results appear physically meaningful since counts of individual blocks and cyclones during the corresponding windows range from 1–2 (for blocks) and 4–11 (for cyclones). More accurate patterns could possibly be obtained by choosing more than two clusters, but the results of Fig. 16 are only meant to illustrate the potential of the method.

6 Summary and outlook

Weather persistence on S2S timescales has been a topic of research since the early days of meteorology. Quasi-stationary or recurrent behavior are common features of weather dynamics and are strongly related to fundamental physical processes, weather predictability and surface weather impacts. Studying weather persistence is therefore important for theoretical and practical reasons. One challenge is that persistence remains a very broad concept that relates to different behaviors in dynamical systems. We propose a typology/structure for the broad concepts related to persistence. Namely that persistence is used to describe the average behavior of a system across its whole trajectory (global persistence), to refer to specific segments of this trajectory (episodic persistence) or to qualify the behavior of particular system states (state persistence).

A wide range of methods has been introduced in the literature to describe persistence in weather and climate series, and several exist for each type of persistence. They offer many distinct and often complementary perspectives on persistence. Some methods quantify persistence in a statistical framework, which can be very relevant for weather forecasting. Others focus instead on persistent periods, including statistical flukes; this is a useful approach for risk assessment, but it is one that says nothing about the overall behavior of the series. Other methods still aim to identify which weather patterns tend to be more persistent than others.

The diversity of existing methods presented in this review reflects the fact that persistence is a multi-faceted concept. While we can agree on a general definition of the concept, many options exist when it comes to actually quantify persistence in real-world data. Though different methods may be related, each sheds light on a particular aspect of persistence. What is meant by persistence cannot be dissociated from the metric used to quantify it. The choice of method should be guided by the end goal, whether it is process understanding, risk and impact assessment, or predictability. Future research should nevertheless perhaps consider testing the robustness of results to the choice of persistence metric more systematically. This could be particularly important for better characterizing potential trends in persistence under climate change and their associated impacts. As a final note, while we centered our review on S2S persistence, most of the methods we discussed apply in principle to other timescales, be they sub-daily or multi-decadal. What really differs is how to define the system to analyze (daily time steps may not be relevant for inter-annual variability, for instance) and interpret persistence.