# Revealing the statistics of extreme events hidden in short weather forecast
data

Justin Finkel<EMAIL_ADDRESS>Department of Earth, Atmospheric, and Planetary
Sciences, Massachusetts Institute of Technology  Edwin P. Gerber  Courant
Institute of Mathematical Sciences, New York University  Dorian S. Abbot
Department of Geophysical Sciences, University of Chicago  Jonathan Weare
Courant Institute of Mathematical Sciences, New York University

###### Abstract

Extreme weather events have significant consequences, dominating the impact of
climate on society. While high-resolution weather models can forecast many
types of extreme events on synoptic timescales, long-term climatological risk
assessment is an altogether different problem. A once-in-a-century event
takes, on average, 100 years of simulation time to appear just once, far
beyond the typical integration length of a weather forecast model. Therefore,
this task is left to cheaper, but less accurate, low-resolution or statistical
models. But there is untapped potential in weather model output: despite being
short in duration, weather forecast ensembles are produced multiple times a
week. Integrations are launched with independent perturbations, causing them
to spread apart over time and broadly sample phase space. Collectively, these
integrations add up to thousands of years of data. We establish methods to
extract climatological information from these short weather simulations. Using
ensemble hindcasts by the European Center for Medium-range Weather Forecasting
(ECMWF) archived in the subseasonal-to-seasonal (S2S) database, we
characterize sudden stratospheric warming (SSW) events with multi-centennial
return times. Consistent results are found between alternative methods,
including basic counting strategies and Markov state modeling. By carefully
combining trajectories together, we obtain estimates of SSW frequencies and
their seasonal distributions that are consistent with reanalysis-derived
estimates for moderately rare events, but with much tighter uncertainty
bounds, and which can be extended to events of unprecedented severity that
have not yet been observed historically. These methods hold potential for
assessing extreme events throughout the climate system, beyond this example of
stratospheric extremes.

## Key points

  1. 1.

Extreme weather risk is inherently difficult to quantify because of data
scarcity.

  2. 2.

Subseasonal weather forecast ensembles are an untapped resource for computing
statistics of extremes, done here by statistical weighting.

  3. 3.

We characterize stratospheric extremes up to once in 500-year events from
46-day hindcast ensembles across 20 winters.

## Plain Language Summary

Weather extremes are a continually recurring threat to human life,
infrastructure, and economies. Yet, we only have sparse datasets of extremes,
both simulated and observed, because by definition they occur rarely. We
introduce an approach to extract reliable extreme event statistics from a non-
traditional data source: short, high-resolution weather simulations. With only
20 years of 46-day weather forecasts, we estimate the magnitudes of once-
in-500-year events.

##  1 Introduction

The atmosphere’s extreme, irregular behavior is, in some ways, more important
to characterize than its typical climatology. A society optimized for average
historical weather patterns is highly exposed to damage from extreme heat and
cold, flooding, and other natural hazards. Moreover, extremes may respond more
sensitively than mean behavior to climate change, an argument supported by
elementary statistics Wigley (2009), empirical observations Coumou & Rahmstorf
(2012); AghaKouchak et al. (2014); O’Gorman (2012); Huntingford et al. (2014);
Naveau et al. (2020) and simulations Pfahl et al. (2017); Myhre et al. (2019).
Recent unprecedented extreme weather events demonstrate the serious human
impacts Mishra & Shah (2018); Van Oldenborgh et al. (2017); Goss et al.
(2020); Fischer et al. (2021). The overall “climate sensitivity” Hansen et al.
(1984), summarized by a change in global-mean temperature, does not do justice
to these consequences, which has led to the development of “event-based
storylines” Shepherd et al. (2018); Sillmann et al. (2021) as a more tangible
expression of climate risk.

The intermittency of extreme events makes precise risk assessment exceedingly
difficult. 100 flips of a biased coin with $\mathbb{P}\\{\text{Heads}\\}=0.01$
is almost as likely to yield zero heads (probability 0.366) as one head
(probability 0.370), and half as likely to yield two heads (probability
0.185). Similarly, in a 100-year climate simulation or historical record, a
once-per-century event will more likely appear either non-existent or twice as
likely as it really is. This difficulty is present for a stationary climate,
but worsens in the presence of time-dependent forcing, anthropogenic or
otherwise. The limited historical record forces us to use numerical models as
approximations, introducing a dilemma: we can run cheap, coarse-resolution
models for long integrations, providing reliable statistics of a biased
system, or expensive, high-resolution models for short integrations, which
have lower bias but provide statistics with higher variance due to under-
sampling. For example, the Integrated Forecast System (IFS) of the European
Center for Medium-Range Weather Forecasts (ECMWF) is one of the most accurate
weather models available today, running at high resolutions of $\sim$16-32 km
ECMWF (2016). The forecasts are skillful, but typically last for a single
season or less—far too short a duration to estimate rare event probabilities
directly.

However, these forecasts are launched multiple times every week in large
parallel ensembles, which can be exploited to bridge the gap from weather to
climate timescales. The key is to include the data from ensemble members in a
statistically principled way. Our main contribution in this paper is to
introduce methods to achieve this, using the ensemble forecasts archived in
the subseasonal-to-seasonal (S2S) project at ECMWF Vitart et al. (2017).

Specifically, in this work we estimate probabilities of sudden stratospheric
warming (SSW) events, in which the winter stratospheric polar vortex rapidly
breaks down from its typical state with a strong cyclonic circulation over the
winter-hemisphere pole. The associated subsidence of air in the polar
stratosphere leads to adiabatic warming, causing lower-stratospheric
temperatures to rise up to 40 K or more over a few days Baldwin et al. (2021).
The breakdown of the stratospheric vortex exerts a “downward influence” on
tropospheric circulation Baldwin & Dunkerton (2001); Baldwin et al. (2003);
Hitchcock & Simpson (2014); Kidston et al. (2015). The midlatitude jet and
storm track shift equatorward, bringing extreme cold spells and other
anomalous weather to nearby regions Kolstad et al. (2010); Kretschmer et al.
(2018). For example, King et al. (2019) documents the impact of an SSW on
extreme winter weather over the British Isles, the so-called “Beast from the
East” in February 2018. SSWs are a demonstrated source of surface weather
predictability on the subseasonal-to-seasonal (S2S) timescale Sigmond et al.
(2013); Butler et al. (2019); Scaife et al. (2022). Pushing this “frontier” of
weather forecasting can improve disaster preparation and resource management
in the face of meteorological extremes White et al. (2017); Bloomfield et al.
(2021). For these reasons, there is keen interest in improving (i) the
prediction of SSW itself beyond the horizon of $\sim$10 days that marks the
current state-of-the-art Tripathi et al. (2016); Domeisen et al. (2020), and
(ii) understanding of the long-term frequency and other climatological
statistics of SSWs Butler et al. (2015); Gerber et al. (2022).

##  2 Data and definitions

Fig. 1(a,b) show the evolution of zonal-mean zonal wind at 10 hPa and 60∘N, a
standard index for the strength of the stratospheric polar vortex Butler et
al. (2019), which we abbreviate $U$. Blue timeseries show $U$ through two
consecutive winters: (a) 1998-1999, in which an extreme SSW occurred as
quantified by the deep drop in $U$ in mid-December, and (b) 2009-2010, when a
more mild SSW occurred in February. Both timeseries are superimposed upon the
1959-2019 ERA-5 climatology.

$U$ is typically positive throughout the winter months, characterizing a
strong circumpolar jet that forms in the stratosphere during the polar night.
The standard definition of an SSW event is that $U$ changes sign Butler et al.
(2015), but it does not capture the range of intensities between events.
Clearly, December 1998 exhibited a much stronger breakdown of the vortex than
February 2010. More intense SSW events have been linked to stronger
tropospheric impacts Karpechko et al. (2017); Baldwin et al. (2021), which
motivates our efforts to distinguish between them. Historical data can provide
reasonably robust estimates of moderately rare events such as February 2010,
in which $U$ barely reversed sign; events of this magnitude occur on average
every two years. On the other hand, extraordinary events like December 1998
have only been observed a few times.

We define an SSW as the first decrease in $U$ below a threshold
$U^{(\mathrm{th})}$ during the “SSW season” of Nov. 1-Feb. 28. We only count
the first event of a winter to exclude the subsequent oscillations of $U$
about $U^{\text{(th)}}$ as separate SSW events, without complicating the
definition with a minimum separation time as in Charlton & Polvani (2007). The
main quantity of interest is the _rate_ : the average number of SSW events per
year, a number between zero and one. Equivalently, the reciprocal of the rate
is called the _return period_ : the expected number of years to wait before an
event of a given severity. For the standard threshold $U^{(\mathrm{th})}=0$,
the rate is approximately $0.6$ Baldwin et al. (2021), but we will consider a
range of severities by varying $U^{(\mathrm{th})}$ down to $-52$ m/s.

One can estimate the rate with reanalysis by counting the fraction of years
with an SSW event. Fig. 2a shows two reanalysis-derived rate estimates as a
function of $U^{\text{(th)}}$: the blue points use 61 years of data
(1959-2019) while the orange points use only 20 years of data (1996-2015). The
corresponding error bars encompass the 50% (thick lines) and 95% (thin lines)
confidence intervals of $(X_{1}+\ldots+X_{n})/n$, where the $X_{i}$’s are
independent Bernoulli random variables with success rate given by the
estimated rate, and a number of trials equal to $20$ (blue) or $61$ (orange).
Fig. 2b-e shows the corresponding seasonal distribution of events at four
selected thresholds, with histograms normalized to have unit area. It may
appear inconsistent that the support of the distribution at
$U^{\mathrm{(th)}}=-8$ is not fully contained in the support of the
distribution at $U^{\mathrm{(th)}}=0$; for example, the third blue and orange
bins of February have positive weight in panel (c), but zero weight in panel
(b). There is no contradiction: although every winter with an SSW at level
$-8$ also must have an SSW at level $0$, the weaker threshold is crossed first
and is sometimes counted in previous weeks.

Figure 1: Climatology of polar vortex and illustration of dataset. Black
curves show the mean seasonal cycle of $\overline{u}$(10 hPa, 60∘N),
abbreviated as $U$, and two gray envelopes show the percentile ranges 25-75
and 0-100, respectively. All statistics are computed with respect to the
61-year ERA-5 dataset between 1959 and 2019. Two individual years are shown in
blue: 1998-1999 (a) and 2009-2010 (b). Two ensembles of S2S hindcasts (red)
are shown each winter, a small sample from the large S2S dataset of two
ensembles _per week_ from the ECMWF IFS. Horizontal dashed lines mark several
different SSW thresholds $U^{\mathrm{(th)}}$ used in this study, including the
standard threshold of 0 m/s and several more extreme ones. The time window
Nov. 1 - Feb. 28 is marked by vertical dashed lines. When $U$ crosses
$U^{\mathrm{(th)}}$ downward for the first time within the time window, an SSW
has occurred. (c): schematic of the flux-counting method. (d): schematic of
the Markov State Model (MSM) method.

The other estimates displayed in Fig. 2 are derived from the S2S dataset,
which consists of hindcast trajectories launched in ten-member perturbed
ensembles (plus a control member that we omit from our analysis). We use only
hindcast data produced by the 2017 version of the ECMWF IFS: that is,
integrations initialized from past initial conditions for the 20 years prior,
in our case from autumn 1996 to spring 2016 (labeled 1996-2015 in the plots).
Each ensemble member has small perturbations applied to its initial
conditions, and is integrated forward with stochastically perturbed tendencies
Buizza et al. (1999); Berner et al. (2009). For details on the the model, see
Vitart et al. (2017) and ECMWF (2016). The dataset is publically accessible at
https://apps.ecmwf.int/datasets/data/s2s/.

The total number of days contained therein is roughly

$\displaystyle 20\text{ years}\times\frac{52\text{ weeks}}{\text{
year}}\times\frac{2\text{ ensembles}}{\text{week}}\times\frac{10\text{
members}}{\text{ensemble}}\times\frac{47\text{ days}}{\text{member}}\approx
2700\text{ years}$ (1)

or, only counting Nov-Feb, 900 years. Many of these extra ensemble members
reach farther into the negative-$U$ tails than the reanalysis. Thinking of
these as alternative realities, we can calculate otherwise inaccessible
probabilities.

Figure 2: Rate estimates derived from S2S and reanalysis. Left: SSW rate
(inverse annual probability of SSW) as a function of zonal wind threshold,
$U^{(\mathrm{th})}$, estimated by the four methods described in the text.
Error bars indicate 95% confidence intervals. MSM and flux-counting error bars
are computed by bootstrapping on entire years of data. ERA5 error bars are
computed analytically as the 2.5-97.5 percentile range of the success rate of
a binomial random variable with a success probability given by the estimated
rate and a number of trials given by the number of years in the record (20 or
61). Error bars going off the bottom of the plot include zero (note the log
scale). Right: seasonal distribution of SSW events at four selected
thresholds, according to each of four methods. All histograms have a bin width
of 7 days and are rescaled to have unit area.

##  3 Two estimates of long return times from short trajectories

To take advantage of the S2S data, we have to overcome two complications.
First, not all trajectories are independently sampled: on the contrary, all
members of an ensemble are initialized close to reanalysis, and take several
days to separate. Thus, the effective sample size is smaller than 900 years.
Second, no individual ensemble can directly provide an SSW probability beyond
the 46-day time horizon, which is well short of the 120 days between November
1 and February 28 when SSWs are allowed to happen. We cannot use hindcasts
directly to estimate the rate, because we need to know what would have
unfolded if the 46-day simulation were to continue. The challenge is to make
use of the “hanging” trajectory endpoints, such as the eight members of the
first ensemble shown in Fig. 1 which do not dip below the threshold. Below, we
present two related, but distinct methods: flux-counting and Markov state
modeling.

###  3.1 Flux-counting for direct estimates

The first approach is quite simple, as sketched in Fig. 2c: we compute the
probability of an SSW on each day by calculating the fraction of trajectories
that cross the threshold on that day (avoiding double counting by keeping
track of “active” trajectories as detailed below), and the sum up all the
daily probabilities over the season. Formally, we decompose the winter months
of interest into a sequence of one-day windows, which is the sampling
resolution of S2S:

$\displaystyle T_{0}=\text{Nov 1},T_{0}+1=\text{Nov.
2},\ldots,T_{1}=\text{Feb. 28}$ (2)

and estimate the probability of an SSW separately for each calendar day. By
our definition, an SSW can happen at most once per season, to ensure the
events are disjoint and have additive probabilities:

Rate $\displaystyle=\sum_{t=T_{0}}^{T_{1}}\mathbb{P}\\{\text{SSW on day }t\\}$
(3) $\displaystyle=\sum_{t=T_{0}}^{T_{1}}\mathbb{P}\Big{\\{}\min_{T_{0}\leq
s<t}U(s)>U^{\text{(th)}}\text{ and }U(t)\leq U^{\text{(th)}}\Big{\\}}$ (4)

The summand can be considered a probability per day of crossing the threshold
$U^{(\mathrm{th})}$, i.e., one of the horizontal dashed lines in Fig. 1a,c. It
is estimated by averaging over all hindcast trajectories that are “active” on
calendar day $t$, meaning those launched some day between $t-46$ and $t$. More
precisely, if we enumerate the active trajectories by
$i\in\mathcal{I}(t)=\\{1,...,N(t)\\}$ and denote the $i$’th trajectory’s zonal
wind by $U_{i}$, then the estimate of daily SSW probability is

$\displaystyle\mathbb{P}\\{\text{SSW on day
}t\\}=\frac{1}{N(t)}\sum_{i\in\mathcal{I}(t)}\mathbb{I}\Big{\\{}\min_{T_{0}\leq
s<t}U_{i}(s)>U^{\text{(th)}}\Big{\\}}\mathbb{I}\Big{\\{}U_{i}(t)\leq
U^{\text{(th)}}\Big{\\}}$ (5)

where $\mathbb{I}$ is an indicator function, equal to 1 if the argument is
true and 0 if the argument is false. In words, we count the trajectories that
dip below $U^{\text{(th)}}$ for the first time on day $t$, as a fraction of
all trajectories that are active on that day. The past of ensemble member $i$
before its initialization date is given by the corresponding reanalysis from
which it branched.

Summing up these probabilities from Nov. 1 to Feb. 28, and sweeping over all
thresholds $U^{\text{(th)}}$, we obtain the black curve in Fig. 2a. Error bars
come from a bootstrapping procedure: we apply the estimate (5) to 20 different
random 10-year subsets of $\\{1996,\ldots,2015\\}$, calculate the 2.5th and
97.5th percentiles of rate estimates, and form the pivotal 95% bootstrap
confidence interval (see Wasserman (2004), chapter 8, for a formal account,
although we have modified the procedure by sampling without replacement to
maintain independence of different years.)

The average (5) is a sum of _dependent_ random variables, with all ensemble
members in a given year sharing common history. This increases the variance of
the estimator or, in other words, reduces the effective sample size from 900
years. This situation is common in the Monte Carlo simulation for inverse
problems. But the error bars make clear that flux-counting enjoys a tremendous
advantage over the direct ERA5 estimate. At all thresholds, the flux-counting
error bar overlaps with the ERA5 error bar, but is much smaller. This gives us
confidence to trust the flux-counting estimate farther into the tail where no
ERA5 data are available.

###  3.2 Markov state model

The second method is more intricate, but delivers more insight into the
predictability of SSWs. We construct a _Markov state model_ (MSM) Deuflhard et
al. (1999); Pande et al. (2010); Chodera & Noé (2014) which is sketched in
Fig. 2d. On each day $t$, we partition state space into a disjoint collection
of bins $S_{t,1},S_{t,2},\ldots,S_{t,M_{t}}$ and approximate the transition
probability matrix for each time-step from $t$ to $t+1$,

$\displaystyle P_{t,t+1}(j,k)$ $\displaystyle=\mathbb{P}\\{\mathbf{X}(t+1)\in
S_{t+1,k}|\mathbf{X}(t)\in S_{t,j}\\},$ (6)

by counting the transitions between corresponding boxes. Explicitly,

$\displaystyle P_{t,t+1}(j,k)$
$\displaystyle=\frac{\sum_{i\in\mathcal{I}(t)}\mathbb{I}\\{\mathbf{X}_{i}(t)\in
S_{t,j}\\}\mathbb{I}\\{\mathbf{X}_{i}(t+1)\in
S_{t+1,k}\\}}{\sum_{i\in\mathcal{I}(t)}\mathbb{I}\\{\mathbf{X}_{i}(t)\in
S_{t,j}\\}}$ (7)

The matrices are row-normalized, which corrects for the redundancy and
statistical dependence between ensemble members. This sequence of matrices is
the key ingredient that enables all downstream calculations. Choosing the
partition of state space is a crucial step which involves a tradeoff: too few
clusters will coarsen the dynamics too much, whereas too many clusters will
reduce the number of data points in each cluster and thus increase the
statistical noise involved in estimating $P_{t,t+1}(j,k)$. There is a lack of
general theory on how to construct MSMs, but here we exploit the particular
structure of the dataset to validate our algorithmic choices, as explained in
the supplement. Here, we focus on conveying the general MSM procedure.

We build the sets $S_{t,i}$ using $k$-means clustering of the data using the
scikit-learn package Pedregosa et al. (2011). As input to $k$-means, we use a
vector of feature $\Phi$ consisting of time-delays of $U$:

$\displaystyle\Phi(\mathbf{X}(t))=[U(\mathbf{X}(t)),U(\mathbf{X}(t-1)),\ldots,U(\mathbf{X}(t-\delta+1))]$
(8)

where $\delta$ days is the number of retained time-delays. This time-delay
embedding encodes additional information about the atmospheric state, enabling
a model based just on the zonal mean wind at 10 hPa. Heuristically, the
embedding captures approximate time-derivatives up to order $\delta$-1. The
technique has precedent in climate science Ghil et al. (2002), and a growing
body of theoretical and empirical evidence supports the use of time-delay
coordinates as reliable features for encoding dynamical attractors Takens
(1981); Kamb et al. (2020); Broomhead & King (1986); Giannakis & Majda (2012);
Brunton et al. (2017); Thiede et al. (2019); Strahan et al. (2021). We have
also experimented with richer feature spaces including EOFs of geopotential
height, but found it unnecessary.

We find that any $\delta$ from 2 to 10 and any number of clusters (denoted
$M_{t}$) from 50 to 150 gives similar results. In Fig. 2 we display results of
a single representative choice of $\delta=5$ days and $M_{t}=150$, along with
a shaded 95% confidence interval derived from the pivotal bootstrap procedure
Wasserman (2004) with 20 independent resamplings of the data (but without
replacement). The supplement further explains how we selected these parameters
to simultaneously optimize the MSM’s fidelity and robustness on a simple
performance benchmark. We emphasize that these clusters are not supposed to
identify metastable weather regimes in the tradition of, e.g., Michelangeli et
al. (1995); rather, they are a discretization of state space meant to
represent continuous functions over that space, encoding gradual progress
towards an SSW event.

Given the clusters $\\{S_{t,j}\\}$ and the transition matrices
$\\{P_{t,t+1}\\}$, we can calculate the rate and seasonal distribution of SSW
events with the following procedure.

  1. 1.

Let $B$ denote the set of “weak-vortex” clusters: all $(t,j)$ such that
$T_{0}\leq t\leq T_{1}$ and the majority of data points in $S_{t,j}$ have
$U<U^{\text{(th)}}$. Let $A$ denote the set of “non-winter” clusters: all
$(t,j)$ such that $t<T_{0}$ or $t>T_{1}$. With this setup, an SSW event is a
_transition from $A$ to $B$_.

  2. 2.

Compute the _committor probability_ ,

$\displaystyle q^{+}_{t}(j)$ $\displaystyle=\mathbb{P}\Big{\\{}U(s)\leq
U^{\text{(th)}}\text{ for some }s\in[t,T_{1}]\ \Big{|}\ \mathbf{X}(t)\in
S_{t,j}\Big{\\}},$ (9)

by solving the following terminal/boundary-value problem. By definition,
$q^{+}_{T_{1}+1}(j)=0$ for all clusters $j$ at the end of winter, while
$q^{+}_{t}(j)=1$ for all $(t,j)\in B$. Stepping backward through time, we have
a recursion relation:

$\displaystyle q^{+}_{t}(j)$
$\displaystyle=\sum_{k=1}^{M_{t+1}}P_{t,t+1}(j,k)q^{+}_{t+1}(k).$ (10)

In words, for a vortex that is initially strong today $(t)$ to break down by
Feb. 28 ($T_{1}$), it must break down sometime between tomorrow ($t+1$) and
$T_{1}$. Hence $q^{+}_{t}(j)$ is a weighted combination of $q^{+}_{t+1}(k)$
for all possible scenarios $k$ for tomorrow. This equation is simply the
Kolmogorov Backward Equation in discrete form E et al. (2019). In this light,
viewing Eq. (10) as a discretized partial differential equation, the clusters
$\\{S_{t,j}\\}$ can be seen as members of a finite element basis and
$P_{t,t+1}(i,j)$ as stiffness matrices. Indeed, here we use an MSM as a
“dynamical Galerkin approximation”, a basis expansion approach to computing
forecast quantities like the committor probability from short trajectory data
that was originally developed for chemistry applications Thiede et al. (2019);
Strahan et al. (2021) and has recently been applied to climate dynamics Finkel
et al. (2021, 2022); Jacques-Dumas et al. (2022).

  3. 3.

Estimate an empirical probability distribution over clusters at the beginning
of winter,

$\displaystyle\pi_{T_{0}}(j)$ $\displaystyle=\mathbb{P}\\{\mathbf{X}(T_{0})\in
S_{T_{0},j}\\}$ (11)

However, in practice, the result is not sensitive to the choice of initial
probability distribution. This is because $T_{0}$ is early enough in the
winter season that the distribution of $U$ is still narrow (see Fig. 1) and
the memory of initial conditions is practically erased by the time of the
first SSW. We can also propagate $\pi$ to each day of the season, using the
Kolmogorov Forward equation (a.k.a. the Fokker Planck equation) in discrete
form:

$\displaystyle\pi_{t+1}(k)$
$\displaystyle=\sum_{j=1}^{M_{t}}\pi_{t}(j)P_{t,t+1}(j,k)$ (12)

  4. 4.

Compute the rate as the average of committor probabilities on the first day of
the SSW season, weighted by the probability distribution $\pi_{T_{0}}$:

$\displaystyle R$
$\displaystyle=\sum_{j=1}^{M_{T_{0}}}q^{+}_{T_{0}}(j)\pi_{T_{0}}(j)$ (13)

In words, the probability of SSW in a random year is the sum of probabilities
from every possible initial condition, weighted by the probability of that
initial condition. Fig. 2 shows in purple the rate according to the MSM ,
which matches remarkably well with the flux-counting method. Error bars
indicate the 95% confidence interval, obtained with the same bootstrapping
procedure that we used for flux-counting. In particular, the _entire
clustering_ procedure is repeated for each 10-year subset of data.

  5. 5.

Compute the seasonal distribution by decomposing the rate over all possible
entrance times to $B$, rather than exit points (i.e., initial conditions) from
$A$:

$\displaystyle R$
$\displaystyle=\sum_{t=T_{0}-1}^{T_{1}}\sum_{j=1}^{M_{t}}\sum_{k=1}^{M_{t+1}}q_{t}^{-}(j)P_{t,t+1}(j,k)\mathbb{I}\\{(t+1,k)\in
B\\}$ (14)

where $q_{t}^{-}(j)=\mathbb{P}\\{\text{no SSW has occurred yet between $T_{0}$
and }t|\mathbf{X}(t)=j\\}$ is known as the _backward committor_. The backward
committor obeys a recursion analogous to that of $q_{t}^{+}$, but moving
backward through time and with a time-reversed transition matrix:

$\displaystyle q_{t+1}^{-}(k)$
$\displaystyle=\begin{cases}\sum_{j=1}^{M_{t}}P_{t,t+1}(j,k)\frac{\pi_{t}(j)}{\pi_{t+1}(k)}q_{t}^{-}(j)&(t,j)\notin
A\cup B\\\ 0&(t,j)\in B\\\ 1&(t,j)\in A\end{cases}$ (15)

The purple histogram in Fig. 2b-e is given by the individual summands (in
groups of 7, according to the bin width of 7 days).

The committor, defined in step 2 above, measures probabilistic progress
towards an SSW event (how likely). To measure _temporal_ progress (how soon),
we further define the _hitting time_ as

$\displaystyle\tau^{+}_{t}=\min\\{s\geq 0:(t+s,\mathbf{X}(t+s))\in B\\}$ (16)

This is a random variable that tells you the timing of the SSW, depending on
the realization of $\mathbf{X}$. We compute two summary statistics of this
random variable. First, its cumulative probability mass function
$\mathbb{P}\\{\tau_{t}^{+}<\sigma\\}$ is a _time-limited_ version of the
committor, which we use to validate our choice of MSM parameters (see the
supplement). Second, its average value, conditional on the vortex actually
breaking down the same winter, is called the _expected lead time_ :

$\displaystyle\eta_{t}^{+}=\mathbb{E}\big{[}\tau_{t}^{+}|t+\tau_{t}^{+}\leq
T_{1}\big{]}$ (17)

This is another useful summary statistic to quantify how far away the system
is from an SSW event. We displayed a similar quantity in Finkel et al. (2021,
2022) in the context of an idealized model. The expected lead time can also be
computed by recursion with the MSM, but the formula is slightly more involved
and left to the supplement.

Let us take a brief aside to reference some mathematical context for the
method above. The general framework that we have used to combine committor
probabilities to compute rates and other steady-state statistics of rare
transitions is _transition path theory_ (TPT) Vanden-Eijnden (2014). TPT has
been applied to molecular dynamics Noé et al. (2009); Meng et al. (2016);
Strahan et al. (2021); Antoszewski et al. (2021), atmospheric and oceanic
sciences Finkel et al. (2020, 2022); Miron et al. (2021, 2022) and social
sciences Helfmann et al. (2021). Though TPT is typically formulated in a time-
homogeneous setting, here we have built in explicit time-dependence to deal
with the seasonal cycle, similarly to Helfmann et al. (2020).

Our MSM-based approximation of the committor probability is similar in spirit
to analogue forecasting van den Dool (1989), which is enjoying a renaissance
with novel data-driven techniques, especially for characterizing extreme
weather Chattopadhyay et al. (2020); Lucente et al. (2022). Dynamical Galerkin
approximation (using a basis different than the one used here) and a short
trajectory variant of analogue forecasting are tested on several benchmark
problems in Jacques-Dumas et al. (2022). Formally, the transition operator
encoded by the matrix in (6) is related to linear inverse models Penland &
Sardeshmukh (1995), which have also been used to predict subseasonal extremes
Tseng et al. (2021). Both MSMs and linear inverse models involve finite-
dimensional approximations of the transition operator (or Koopman operator for
deterministic dynamics) Mezić (2013, 2005); Klus et al. (2018).

##  4 Results

###  4.1 Rate estimates

Fig. 2 compares rate estimates from the MSM and flux-counting methods against
the reanalysis rates. We include the ERA-5 estimator based on just 1996-2015
to overlap with the S2S period. For the mild thresholds of
$U^{\mathrm{(th)}}=0,-4$ m/s, corresponding to return times of 2-3 years, the
MSM and flux-counting estimates agree with both short- and long-term
reanalysis estimates. Moving to moderate thresholds of
$U^{\mathrm{(th)}}\sim-28$ m/s, the MSM and flux-counting rates track somewhat
closer with the 61-year estimate (orange), which has slightly lower rates
across the board. The S2S data were initialized from the 20-year time period
corresponding to the blue curve, but the S2S hindcasts recover the longer-term
climatology, despite the (slightly) greater SSW frequency from the period in
which they were initialized.

One can think of these SSW frequencies as the _climatology according to the
Integrated Forecast System_ , given the boundary conditions of the 1996-2015
period. At least in the “model world” of the IFS, it does not appear that a
differences in atmospheric boundary conditions (e.g., sea surface
temperatures) caused a systematic increase in SSW frequency between 1996 and
2015; rather the observed increase in SSWs was luck of the draw. It is
possible, however, that systematic model error could be obscuring the
systematic differences suggested by Reichler et al. (2012) and Dimdore-Miles
et al. (2021).

At all levels of $U^{\mathrm{(th)}}$, but especially in the negative extremes,
the confidence intervals from the two S2S estimates are smaller than those
from the ERA5 estimates, thanks to the large amount of S2S data that the MSM
and flux-counting methods can exploit. The agreement between MSM and ERA-5 on
common events gives us more confidence to trust the MSM on less common events
in the negative-$U^{(\mathrm{th})}$ tail, where ERA-5 data are too sparse to
give a meaningful rate estimate. Both the direct counting and MSM approaches
suggest the potential for events where the vortex becomes so disrupted it
spins -40 m/s (stronger than average, but in the opposite direction), albeit
only once or twice in a millennium.

Several recent studies have performed the same task of filling out a sparse
climate distribution using models Horan & Reichler (2017); Kelder et al.
(2020), but with uninterrupted long runs of a global climate model. The
techniques we have introduced—MSM and flux-counting—offer a novel way to
estimate such quantities from short trajectories only, without access to a
centennial-scale run of the IFS which the standard estimation method would
require. We believe the higher resolution IFS is also more appropriate for
capturing the most extreme SSWs.

###  4.2 Seasonal distribution

Fig. 2b-e illustrate that S2S data also offers an advantage for describing the
seasonal distribution of SSW events, an inherently noisier statistic than the
full-winter rate estimate because one has to split the data into finer
categories. The S2S-derived histograms, in black and purple, are able to bring
out seasonal structure that is ambiguous in the reanalysis data directly. For
$U^{\mathrm{(th)}}=0$ m/s, there is a gradual increase in SSW frequency from
November to January, followed by a plateau in February, consistent with prior
studies of seasonality at monthly resolution Charlton & Polvani (2007) and
supporting the late winter maximum found by Horan & Reichler (2017).

As $U^{\mathrm{(th)}}$ becomes more negative, the reanalysis histograms
dwindle and degenerate into a few isolated spikes, whereas the S2S histograms
become intriguingly bimodal, but retain their smoothness. The S2S histograms
show a persisting SSW occurrence throughout February after the second peak, a
feature that is also faintly present in the longer reanalysis period
(1959-2019), but not at all in the shorter reanalysis period (1996-2015) from
which S2S was initialized. Again, the IFS recovers features of the longer-term
climatology.

The January/February peak is documented in the literature, e.g., by Horan &
Reichler (2017), who diagnosed the peak as a balance between two time-varying
signals: the background strength of the polar vortex, and the vertical flux of
wave activity capable of disturbing the vortex. The bimodal structure seen in
S2S has also been found tentatively in prior studies with both reanalysis and
models Horan & Reichler (2017); Ayarzagüena et al. (2019), and more robustly
in other features of the boreal winter, e.g., the midwinter suppression of
Pacific storm activity Nakamura (1992). We speculate that the early peak
represents Canadian warmings Meriwether & Gerrard (2004), which our result
suggests may deserve a more decisive classification. Seasonal differences are
associated with dynamical differences in SSW events. For example, “Canadian
warmings” shift the Aleutian high and occur earlier in the winter Butler et
al. (2015). Categorizing SSWs by their seasonality may reveal preferred
timings that indicate when and why the polar vortex is most vulnerable Horan &
Reichler (2017).

###  4.3 Statistical predictors of SSWs

Figure 3: Sparse regression results. Heat maps show the importance of each
feature (listed on the vertical axis) for predicting the expected lead time
$\eta^{+}$ at a range of zonal wind thresholds $U^{(\mathrm{th})}$ (listed
along the horizontal axis). The left-hand heat map shows the Gini importances
from random forest regression, and the right panel shows LASSO coefficients.
We also used $\overline{u}$ at lower levels than 10 hPa as input features, but
found none of them to have any importance, and so omitted them from the
figure.

Estimates of long return times alone do not provide physical insight into the
mechanisms driving the event. The committor probability and expected lead time
estimates provided by the MSM approach encode information on the dynamics and
predictability of SSW events, and on extreme events in general. These
quantities cannot be computed by the flux-counting approach. A number of
recent articles have pursued committor probabilities as windows into
transitional dynamics, e.g., Miloshevich et al. (2022) for European heat waves
and Frishman & Grafke (2022) for the spread of turbulence in a pipe. On SSWs
specifically, our own previous studies with a simple SSW model Finkel et al.
(2021, 2022) found through sparse regression that a small set of physical
variables could explain key variability in the committor.

Here we analyze the S2S dataset in a similar way, using sparse regression to
reveal the main determinants of the committor $q_{t}^{+}$ (how likely is an
SSW to occur?) and expected lead time $\eta_{t}^{+}$ (if it does occur, how
soon?), among a large collection of candidate variables including zonal-mean
zonal winds and meridional eddy heat fluxes at various time delays, altitudes,
and wavenumbers (listed on the left of Fig. 3). We have performed two kinds of
regression: linear regression with a sparsity-promoting $L1$ penalty of $0.1$,
also known as LASSO Tibshirani (1996), and random forest regression Hastie et
al. (2009) with 10 trees of depth 3. Both algorithms, as implemented using
scikit-learn, provide not only predictions of the output variable but also
notions of relevance for each input feature: nonzero coefficients in the case
of LASSO, and Gini importances in the case of the random forest Pedregosa et
al. (2011). These relevances are of greater interest to us than the raw skill
of the regression.

We focus on the early part of the SSW season to connect the results with the
rate formula (13). The target variable for regression is $\log(\eta_{t}^{+})$,
which guarantees the predicted $\eta_{t}^{+}$ is positive and also emphasizes
variability in small values of $\eta_{t}^{+}$ (when an SSW is close) rather
than large values (when an SSW is distant). The training data consist of
trajectory snapshots as inputs and MSM-labeled $\eta_{t}^{+}$ values as
outputs. We include only those snapshots between Nov. 1 and Nov. 30, and
strictly outside of sets $A$ and $B$, where $0<q_{t}^{+}<1$. (We also
regressed on $q_{t}^{+}$ and discovered similar but more subtle patterns of
importance; for brevity we show results only for $\eta_{t}^{+}$.)

Fig. 3 summarizes the results of regression across all $U^{\mathrm{(th)}}$
thresholds. Random forest importances are always nonnegative and represented
on a yellow-black color scale, while regression coefficients are signed and
represented on a red-blue color scale. Red is associated with a weaker vortex,
meaning a shorter lead time. The correlation coefficient $R^{2}$ remains
between 0.4 and 0.6 for both methods across all thresholds, indicating that
these regressions are imperfect expressions of the expected lead time, but do
explain a significant part of the variance.

The models illuminate some interesting patterns, some obvious and some
surprising. _A priori_ , one expects $\overline{u}$(10 hPa, 60∘N, $t$) itself
[the bottom listed feature, abbreviated $U(t)$] to dominate the regression,
since it defines the event. This is true at a mild threshold of
$U^{(\mathrm{th})}$—stronger zonal wind means longer expected lead time,
according to the positive coefficient in the panel’s lower right corner—but
for more extreme thresholds, it is actually the _time-delayed_ zonal wind
$U(t-1),\ldots,U(t-4)$ that is more relevant. Furthermore, the corresponding
LASSO coefficients are negative, suggesting that the decrease over time of $U$
is more important than its value today. At the most extreme thresholds, it
even appears that _strong_ $U(t)$ portends a sooner vortex collapse,
suggesting that the most extreme SSW events (those reaching the most
_negative_ zonal wind) follow from precursor states with anomalously _strong_
zonal wind.

Another important set of features is the 10-day averaged meridional heat flux
$\overline{vT}$ averaged over 45-75∘N, although LASSO and Random Forest
regressions emphasize different altitudes and wavenumber components. Both
methods agree that the 10 hPa heat flux at wavenumbers 0 and 1 exert strong
and competing (statistical) influences on expected lead time: a stronger
wavenumber-0 component ($\overline{v}\overline{T}$) means vortex collapse is
farther away, while a stronger wavenumber-1 component means vortex collapse is
sooner. At lower levels of the atmosphere, the eddy heat fluxes exert
significant but diminishing influences, although they remain important for the
most extreme SSW events (at least according to LASSO).

Figure 4: Committor and expected lead time. Aggregating all S2S data from Nov.
1 to Nov. 30, this figure displays the committor (left column) and expected
lead time (right column) in shading, as well as the climatological probability
density $\pi$ in black contours, on a logarithmic scale ($\pi$ has arbitrary
units, normalized to have unit integral over state space). Two zonal wind
thresholds are considered: 0 m/s (top), and $-8$ m/s (middle) (bottom). All
results are derived from the MSM. We show $\eta^{+}$ and $q^{+}$ as functions
of two variables only: $\overline{u}$(10 hPa, 60∘N, $t$) and meridional heat
flux averaged over 10 days between 45∘N and 75∘N at 100 hPa and wavenumber 1.
The remaining variables are averaged out and weighted by $\pi$ in this
display.

How can we make sense of all these correlations? One simple method of
visualization is to plot the committor and lead time as approximate functions
of two variables (averaging over remaining variables). The regression results
present us with many possible pairs of important variables. Here we select
just one pair: zonal-mean zonal wind at 10 hPa, and wavenumber 1 meridional
heat flux at 100 hPa, averaged over the preceding 10 days. The latter feature
is assigned high importance by the random forest, though not by LASSO, and is
especially interesting as a signal coming from a lower altitude than 10 hPa,
possibly related to the two-way influence characteristic of coupled
troposphere/stratosphere dynamics. Fig. 4 displays the committor (left column)
and expected lead time (right column) as a function of these two variables at
two thresholds: $U^{(\mathrm{th})}=0$ m/s (top), and $-8$ m/s (bottom).
Contours of the climatological probability density $\pi$ signal which regions
are more frequently visited and which ones are rare. We only average over the
first month of the SSW season, Nov. 1-Nov. 30, to represent a map of possible
“initial conditions” for the winter vortex evolution.

The orientation of contours in phase space reveals a pattern of influence that
would be hard to intuit from the regression coefficients alone. At
$U^{(\mathrm{th})}=0$, the $q^{+}$ contours run almost perpendicular to the
$U$ axis, confirming that the zonal wind itself primarily determines how
likely an SSW is for the coming season. But the $\eta^{+}$ contours tell a
different side of the story: at stronger $U$, the contours progressively tilt
away from vertical towards horizontal, indicating that the time until SSW
depends strongly on heat flux—at least in the regime of strong $U$, from where
an SSW is unlikely to begin with. The influence of heat flux grows more
significant as the threshold $U^{(\mathrm{th})}$ is lowered in row 2 of the
figure, even for the committor. We infer a general pattern: the 10-hPa zonal
wind strength in November determines how likely an SSW is for the coming
winter, but when it is rather unlikely, the lower-stratospheric wavenumber-1
heat flux determines when the SSW will happen.

What the phase space images reveal most of all is that $q^{+}$ and $\eta^{+}$
are nonlinear functions: the influence of a variable depends on the state of
all the other variables. Nonlinear regression methods, such as random forests,
are therefore crucial to uncover a complete description.

##  5 Discussion

By comparing S2S results with reanalysis, we are measuring the composition of
potentially three separate error sources: (i) forecast model error, (ii) non-
stationarity of the climate _with respect to SSW events_ over the reanalysis
period, and (iii) numerical errors, both statistical (from the finite sample
size) and systematic (from the projection of forecast functions onto a finite
basis in the case of the MSM). We briefly address each error source in turn.

The S2S trajectories were realized only in simulation, not in the physical
world. Accordingly, our S2S estimates apply strictly to the climatology of the
2017 IFS, a statistical ensemble that could (at least in principle) be
concretely realized by running the model uninterrupted for millennia, with
external climatic parameters sampled from their variability in the 20-year
time window of 1996-2015. Long, equilibrated simulations have been performed
with coarser models by, e.g., Kelder et al. (2020) to assess UK flood risk
(the so-called “UNSEEN” method), and by Horan & Reichler (2017) to assess SSW
frequencies, but this is not practical given the constraints and mission of
the ECMWF IFS. The S2S dataset is an ensemble of opportunity. It was created
to compare the skill of different forecast systems on S2S timescales, not at
all for the purpose of establishing a climatology of SSWs.

And while the IFS model has proven outstanding in its medium-range forecast
skill Vitart (2014); Kim et al. (2014); Vitart & Robertson (2018), it was
designed for short forecasts. It is not clear how it would behave if allowed
to run for hundreds of years as a climate model, which requires careful
attention to the boundary condition and conservation issues. Even if the
climate were to remain stationary with its 1996-2015 parameters, numerical and
model errors would inject some bias into the equilibrated simulation.
Repeatedly initializing S2S forecasts with reanalysis ensures a realistic
background climatology, and allows us to rely on the IFS strictly for the
short-term integrations that it was designed for. Our method may be used as a
diagnostic tool to compare different models against each other, with specific
attention paid to their rare event rates. A useful extension of this work
would be to repeat the analysis on multiple data streams from all 11
forecasting centers worldwide that contribute to the S2S project, as a
different way to compare different models’ ability to represent extremes.

Boreal SSWs provide an ideal demonstration of our method, providing both
moderately and extremely rare events. A natural and intriguing future
application is the rate of Southern-hemisphere SSW events, in the spirit of
Jucker et al. (2021). which is postponed to future work for the sake of
brevity. The method may be extended to other kinds of extremes as well, though
care must be exercised when defining the event (e.g., sets $A$ and $B$) and
choosing features in which to do clustering (for the MSM approach), especially
in the case of more spatially localized events.

The rate we estimate from the S2S data set is based on 1996-2015 boundary
conditions (sea surface temperatures, CO2), and our MSM method assumes the
climate was stationary over this period. Our results indicate that according
to the 2017 IFS, 1996-2015 conditions were more similar to 1959-2019 than
direct counting of SSW events might suggest. This could mean that the IFS was
missing some key climatological variable during that period Dimdore-Miles et
al. (2021). There is, however, substantial uncertaintuy on the impact of
global warming on SSWs, even under 4xCO2 forcing Ayarzagüena et al. (2020). By
repeating our analysis on different historical periods, or simulations
initialized from climate model integrations under different forcing, one could
discern a more decisive signal of forced changes than would be available from
raw data. Moreover, the expression for the SSW rate (13) as a “dot product”
between a committor and a climatological probability density would allow us to
decompose small changes in SSW frequency as changes in these two components
separately. A changing probability density $\pi$ would reflect changes in the
slow background conditions, whereas a changing committor $q^{+}$ would reflect
a change in the system dynamics.

Error source (iii) is the most open to scrutiny and improvement. We have used
the short S2S hindcasts directly to validate our parameter choices for the MSM
(see the supplement). In a sequence of preceding papers Finkel et al. (2021,
2022), we have benchmarked the performance of the method on a highly idealized
SSW model due to Holton & Mass (1976). Nevertheless, large-scale atmospheric
models are a mostly-unexplored frontier for this class of methods.

Our method exceeds what is possible directly from reanalysis, but we are not
yet fully liberated from observations: every S2S trajectory is initialized
near reanalysis, and it only has 46 days to explore state space before
terminating. This fundamentally limits how far we can explore the tail of the
SSW distribution. In other words, the real climate system sets the _sampling
distribution_ which is a flexible but important component in rare event
estimation problems Thiede et al. (2019); Strahan et al. (2021); Finkel et al.
(2021). With an executable model, we could initialize secondary and tertiary
generations of short trajectories to push into more negative-$U$ territory and
maintain statistical power for increasingly extreme SSW events. This is the
essence of many rare-event sampling algorithms, such as those reviewed in
Bouchet et al. (2019) and Sapsis (2021). For example, a splitting large-
deviation algorithm was used in Ragone et al. (2018) to sample extreme
European heat waves and estimate their return times. Quantile diffusion Monte
Carlo was used in Webber et al. (2019) to simulate intense hurricanes, and in
Abbot et al. (2021) to estimate the probability of extreme orbital variations
of Mercury. A natural extension of these various techniques would combine
elements of active rare event sampling with committor estimation via MSMs.
Early developments of such a coupling procedure are presented in Lucente et
al. (2022).

##  6 Conclusion

Extreme weather events present a fundamental challenge to Earth system
modeling. Very long simulations are needed to generate sufficiently many
extreme events to reduce statistical error, but high-fidelity models are
needed to simulate those events accurately. Conventionally, no single model
can provide both, due to computational costs. Here, we have demonstrated an
alternative approach that leverages _ensembles_ of _short_ , high-fidelity
weather model forecasts to calculate extreme weather statistics, with specific
application to sudden stratospheric warming (SSW). By exploiting the huge
database of forecasts stored in the subseasonal-to-seasonal (S2S) database
Vitart et al. (2017), we have obtained estimates of the rate and seasonal
distribution of extreme SSW events. From just 20 years of data, we obtain
probability estimates of events with a 500 year return time, which are so
extreme that the vortex is as strong in the easterly direction as its typical
westerly climatology. These events have never been observed historically, but
can be pieced together using our analysis method.

Our method uses data to estimate the dynamics on a subspace relevant for SSW,
namely the polar vortex strength as measured by zonal-mean zonal wind. This
single observable, augmented by time-delay embedding, gives a simple set of
coordinates sufficient to estimate rate and seasonal distributions. Our
demonstration opens the door to address many other data-limited questions of
basic physical interest. For instance, a high-resolution model could be used
in ensemble forecast mode, but initialized around a decade at the end of this
century provided by a climate model, to understand the impact of global
warming on extremes.

##  7 Open Research

Our analysis is based on publicly available datasets from the European Center
for Medium-Range Weather Forecasts ECMWF (2022b) and from associated
Copernicus Climate Data Store ECMWF (2022a). Python scripts to download the
necessary data and reproduce the paper’s analysis will be made available in a
public Zenodo repository at publication time.

## Acknowledgments

We extend special thanks to Andrew Charlton-Perez, who suggested the S2S
dataset as a case study for the methodology, and Simon Lee, who helped
familiarize us with the data. We thank Amy Butler for help with the
reanalysis. Our collaborators at the University of Chicago, including Aaron
Dinner, John Strahan, and Chatipat Lorpaiboon, offered helpful methodological
advice. Computations for this project were performed on the Greene cluster at
New York University.

J.F. was supported by the U.S. DOE, Office of Science, Office of Advanced
Scientific Computing Research, Department of Energy Computational Science
Graduate Fellowship under Award Number DE-SC0019323 at the time of writing. J.
F. also acknowledges continuing support from the MIT Climate Grand Challenge
on Weather and Climate Extremes. E.P.G. acknowledges support from the NSF
through award OAC-2004572. J.W. acknowledges support from the NSF through
award DMS-2054306 and from the Advanced Scientific Computing Research Program
within the DOE Office of Science through award DE-SC0020427.

##  Appendix A Supporting information

Our work relies completely on publicly available datasets of reanalysis and
hindcasts, which we describe in the subsequent section. We then lay out the
numerical procedure to compute rates and seasonal distributions using
transition path theory (TPT). We then present the formulas used to display
results in the main text. Finally, we document the method used to select
parameters. A software implementation will be made available in a public
repository at the time of publication.

Dataset description

We use two different datasets for this study.

  * •

ERA-5: reanalysis product from ECMWF Hersbach et al. (2020), spanning
1959-2019. We used daily averaged temperature, geopotential height, zonal and
meridional wind at 10, 100, 500, and 850 hPa levels at a resolution of 2.5∘.
ERA-5 was downloaded from the Copernicus Data Store
https://cds.climate.copernicus.eu/

  * •

S2S: perturbed reforecast (hindcast) ensembles from the 2017 model version of
the ECMWF IFS, launched every Monday and Thursday, from fall 1996 through
spring 2016, at the same resolution as ERA5 above. Each hindcast ensemble has
10 perturbed members which run for 47 day, including the initialization date.
We used the same time and space resolution as with ERA5. S2S was downloaded
from the ECMWF data portal https://ecmwf.int.

The S2S dataset can be summarized as follows:

$\displaystyle\text{S2S dataset}=\\{\mathbf{X}_{i}(t):i=1,\ldots,N\\}$ (18)

Here, $\mathbf{X}$ denotes the state vector of all relevant meteorological
variables. Associated with each calendar day $t$ is a subset of “active”
indices, $\mathcal{I}(t)$, containing the trajectories launched sometime
between $t-46$ and $t$. The same notation can be used for the zonal-mean zonal
wind itself, i.e., $U_{i}(t)$ as a function of $\mathbf{X}_{i}(t)$.

##  Appendix B MSM calculations

The main text explains the computation of the committor, $q_{t}^{+}$, and
probability distribution, $\pi_{t}$, for a Markov chain. Here we further
explain how to compute statistics of the hitting time,

$\displaystyle\tau^{+}_{t}=\min\\{s\geq 0:(t+s,\mathbf{X}(t+s))\in B\\}$ (19)

####  B.0.1 PMF of hitting time

The probability mass function of $\tau_{t}^{+}$ conditioned on an initial set
$S_{t,j}$ can be computed with a recursion relation similar to that for
$q_{t}^{+}$. As a base case, we observe that only way for $\tau_{t}^{+}=0$ is
for the system to already be in state $B$:

$\displaystyle\mathbb{P}\\{\tau^{+}_{t}=0|\mathbf{X}(t)\in
j\\}=\begin{cases}1&\text{ if }(t,j)\in B\\\ 0&\ \text{ otherwise}\end{cases}$
(20)

For the case $\tau_{t}^{+}\geq 1$, the hitting time is one plus the hitting
time at the next step:

$\displaystyle\mathbb{P}\\{\tau_{t}^{+}=s|\mathbf{X}(t)\in
j\\}=\sum_{k}P_{t,t+1}(j,k)\mathbb{P}\\{\tau_{t+1}^{+}=s-1|\mathbf{X}(t+1)\in
S_{t+1,k}\\}$ (21)

To carry out this recursion, we compute the left-hand side for all $t$ and $j$
with a fixed $s$ before incrementing $s\rightarrow s+1$. After doing this for
$0,1,\ldots,s$, we then have a _time-limited_ committor,

$\displaystyle
q^{+}_{t,\sigma}(j)=\mathbb{P}\\{\tau_{t}^{+}\leq\sigma|\mathbf{X}(t)\in
S_{t,j}\\}=\sum_{s=0}^{\sigma}\mathbb{P}\\{\tau_{t}^{+}=s|\mathbf{X}(t)\in
S_{t,j}\\},$ (22)

which is the preferred version of the committor for some other studies, e.g.,
Lucente et al. (2022). Below we use the time-limited committor for
hyperparameter tuning.

####  B.0.2 Expectation of hitting time

The most intricate computation is that of the expected lead time
$\eta_{t}^{+}=\mathbb{E}[\tau_{t}^{+}|\tau_{t}^{+}+t<T_{1}]$. We do this via
the _moment-generating function_ ,

$\displaystyle
M_{t}(j;\lambda)=\mathbb{E}\big{[}\exp\big{(}\lambda\tau_{t}^{+}\big{)}\mathbb{I}\\{\tau_{t}^{+}+t<T_{1}\\}\big{|}\mathbf{X}(t)\in
S_{j}\big{]}$ (23)

Along a single trajectory, we have $\tau_{t}^{+}=\tau_{t+1}^{+}+1$ and thus
$\exp\big{(}\lambda\tau_{t}^{+}\big{)}=\exp\big{(}\lambda\tau_{t+1}^{+}\big{)}\exp\big{(}\lambda\big{)}$.
Using this fact, there is a recursive relationship between $M_{t}$ and
$M_{t+1}$:

$\displaystyle
M_{t}(j;\lambda)=\begin{cases}\sum_{k}P_{t,t+1}(j,k)e^{\lambda}M_{t+1}(k;\lambda)&(t,j)\notin
A\cup B\\\ 1&(t,j)\in B\\\ 0&(t,j)\in A\end{cases}$ (24)

One can observe that $M(t;\lambda=0)$ is the committor itself, $q_{t}^{+}$.
But to get at $\tau_{t}^{+}$, we now must differentiate with respect to
$\lambda$:

$\displaystyle\frac{\partial}{\partial\lambda}M_{t}(j;\lambda)$
$\displaystyle=\begin{cases}\sum_{k}P_{t,t+1}(j,k)e^{\lambda}\Big{[}M_{t+1}(k,\lambda)+\frac{\partial}{\partial\lambda}M_{t+1}(k,\lambda)\Big{]}&(t,j)\notin
A\cup B\\\ 0&(t,j)\in B\\\ 0&(t,j)\in A\end{cases}$ (25)

The lead time can be expressed

$\displaystyle\eta_{t}^{+}(j)=\frac{\mathbb{E}[e^{\lambda\tau_{t}^{+}}\mathbb{I}\\{\tau_{t}^{+}+t<T_{t}\\}|\mathbf{X}(t)\in
S_{j}]}{\mathbb{P}\\{\tau_{t}^{+}+t<T_{1}|\mathbf{X}(t)\in
S_{j}\\}}=\frac{[\partial
M_{t}(j;\lambda)/\partial\lambda]_{\lambda=0}}{q_{t}^{+}(j)}$ (26)

Therefore, the recursion relation for $\eta_{t}^{+}$ is found by setting
$\lambda=0$ in the recursion relation for $M_{t}(j;\lambda)$:

$\displaystyle\eta_{t}^{+}(t,j)$
$\displaystyle=\begin{cases}\frac{1}{q_{t}^{+}(j)}\sum_{k}P_{t,t+1}(j,k)e^{\lambda}q_{t+1}(k)\big{[}1+\eta_{t+1}(k)\big{]}&(t,j)\notin
A\cup B\\\ 0&(t,j)\in B\\\ \text{undefined}&(t,j)\in A\end{cases}$ (27)

The formulas for $q^{+}$ and $\eta^{+}$ are exact, straightforward to
implement, and fast to compute after a Markov chain is constructed. How to
choose parameters to construct the chain in the first place is the focus of
the next section.

##  Appendix C MSM hyperparameter selection

Figure 5: Time-limited committors for hyperparameter selection. Left two
columns: shading indicates time-limited committor probabilities according to
counting S2S hindcasts (left) and the MSM (right), at two different
thresholds: $U^{\mathrm{(}th)}=0$ m/s (top) and $-16$ m/s (bottom). Black
contours delineate level sets of the climatological probability density,
derived from the same two methods. The MSM was constructed using 5 time delays
and 150 clusters, a choice informed by systematic evaluation of the MSMs
performance. The upper right panel shows the normalized skill of the MSM
averaged over zonal wind thresholds, while the lower right panel plots the
skill against thresholds for a few selected parameter sets.

The MSM procedure uses two key parameters: the number of time delays $\delta$,
and the number of clusters to use in the resulting $(\delta+1)$-dimensional
feature space. We optimize this choice by comparing the time-limited committor
with $\sigma=20$ according to the MSM ($q^{\mathrm{(MSM)}}$, which takes on
continuous values between 0 and 1) with that according to the S2S data
directly ($q^{\mathrm{(S2S)}}$, which takes on discrete values of 0 or 1).
Forgive the temporary abuse of sub- and superscripts on $q$ in this
subsection. For any ensemble member $\mathbf{X}_{i}(t)$ at some time $t$
between 0 and 46, the latter can be assessed directly by asking whether
$\mathbf{X}_{i}$ achieves SSW before a time $\sigma$ has elapsed. Fig. 5 shows
$q^{+}_{t,\sigma}$ as a function of $t$ and $U^{\mathrm{(th)}}$ for
$\sigma=20$ days and two different zonal wind thresholds, 0 m/s (top) and
$-16$ m/s (bottom). Black contours are level sets of the climatological
probability density, $\pi$. The two columns show $q^{+}_{t,\sigma}$ estimated
from S2S directly (left) and the MSM (right), and they match approximately by
eye. To quantify their agreement, we use the log-likelihood commonly used for
logistic regression, which is equivalent to the (negative) cross entropy of
Bernoulli$\big{(}q^{(\mathrm{MSM})}\big{)}$ relative to
Bernoulli$\big{(}q^{(\mathrm{S2S})}\big{)}$:

$\displaystyle\text{LL}\big{(}q^{(\mathrm{S2S})},q^{(\mathrm{MSM})}\big{)}=\text{mean}\Big{\\{}q^{\mathrm{(S2S)}}_{i}\log\Big{(}q^{\mathrm{(MSM)}}_{i}\Big{)}+\Big{(}1-q^{\mathrm{(S2S)}}_{i}\Big{)}\log\Big{(}1-q^{\mathrm{(MSM)}}_{i}\Big{)}\Big{\\}}.$
(28)

The committors are evaluated at day 23 of the member $i$, and the mean is
taken over the ensemble members $i$ such that $q_{i}^{\mathrm{(MSM)}}$ is
strictly between zero and one on the 23rd day. We choose 23 days because (i)
it is beyond the 15-day maximum lag time, before which not all time-lagged
features are defined for $\delta=15$, and (ii) it leaves room for 20 days of
extra lead time for validation purposes. Other choices in this neighborhood do
not affect our results appreciably.

For extremely rare events—say, with 1% probability, so that the sample mean
$\overline{q}^{(\mathrm{S2S})}=0.01$—there is a class imbalance problem:
$q^{(\mathrm{MSM})}$ can achieve a high LL score with a “climatological”
forecast of $q_{i}^{(\mathrm{MSM})}\equiv\overline{q}^{(\mathrm{S2S})}$. LL
can still be used as a relative score to discriminating between parameter
choices at a fixed $\overline{q}^{(\mathrm{S2S})}$, but we wish to choose MSM
parameters that are optimal on average across a wide range of
$U^{(\mathrm{th})}$ and therefore of $\overline{q}^{(\mathrm{S2S})}$. To
compare these scales meaningfully, we convert LL into an absolute scale
following Benedetti (2010); Miloshevich et al. (2022):

Normalized skill
$\displaystyle=\frac{\mathrm{LL}\big{(}q^{(\mathrm{S2S})},q^{(\mathrm{MSM})}\big{)}-\mathrm{LL}\big{(}q^{(\mathrm{S2S})},\overline{q}^{(\mathrm{S2S})}\big{)}}{\mathrm{LL}\big{(}q^{(\mathrm{S2S})},q^{(\mathrm{S2S})}\big{)}-\mathrm{LL}\big{(}q^{(\mathrm{S2S})},\overline{q}^{(\mathrm{S2S})}\big{)}}$
(29)

The numerator is the improvement relative to the climatological forecast, and
the denominator is the maximum possible improvement, when each
$q^{(\text{S2S})}_{i}$ is predicted exactly. Hence the normalized skill is
always less than one, and typically above zero. The right-hand panels of Fig.
5 display the average normalized skill across thresholds (top) and the
normalized skill as a function of threshold (bottom) for a few parameter
choices. Increasing the number of clusters seems to help regardless of the
number of delays. On the other hand, more delays don’t always help; beyond
$\delta=5$ days, the dimensionality may be introducing degeneracies. Among the
parameter choices considered, the top and bottom panel both suggest that
$\delta=5$ delays and $M_{t}=150$ clusters achieves the best performance among
all the MSM choices shown.

Fig. 6 displays further confirmation that this choice is reasonable and
robust. We have repeated the entire MSM pipeline with 20 different random
subsets of the years 1996-2015 (sampled without replacement), and plotted the
corresponding rates as light purple curves in Fig. 6a. The resampled estimates
cluster about the dashed curve, which uses all 20 years simultaneously, for
all but the most extreme events. On the other hand, in Fig. 6b, we used 15
delays and 200 clusters, which gave a slightly better normalized skill score
(not shown). However, the resampled rates are systematically biased away from
the full 20-year estimate, which raises the suspicion of overfitting. For this
reason, we stick to the more modest choice of 5 delays and 150 clusters.
However, all the qualitative results shown—the committors, expected lead
times, and rate curves—remain remarkably consistent with different choices of
clusters and delays. This gives us confidence in the MSM results shown in the
main text.

Figure 6: Bootstrap-resampled rates. Left: reproduction of Fig. 2, plus 20
additional rate curves obtained by resampling the 20 years without into random
subsets of size 10, without replacement. Right: same for 15 delays and 200
clusters. The displacement of the bootstrapped curves from the main curve is a
symptom of overfitting.

## References

  * Abbot et al. (2021) Abbot, D. S., Webber, R. J., Hadden, S., Seligman, D., & Weare, J. (2021).  Rare event sampling improves mercury instability statistics.  The Astrophysical Journal, 923(2), 236.
URL https://dx.doi.org/10.3847/1538-4357/ac2fa8

  * AghaKouchak et al. (2014) AghaKouchak, A., Cheng, L., Mazdiyasni, O., & Farahmand, A. (2014).  Global warming and changes in risk of concurrent climate extremes: Insights from the 2014 california drought.  Geophysical Research Letters, 41(24), 8847–8852.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2014GL062308

  * Antoszewski et al. (2021) Antoszewski, A., Lorpaiboon, C., Strahan, J., & Dinner, A. R. (2021).  Kinetics of phenol escape from the insulin r6 hexamer.  The Journal of Physical Chemistry B, 125(42), 11637–11649.  PMID: 34648712.
URL https://doi.org/10.1021/acs.jpcb.1c06544

  * Ayarzagüena et al. (2019) Ayarzagüena, B., Palmeiro, F. M., Barriopedro, D., Calvo, N., Langematz, U., & Shibata, K. (2019).  On the representation of major stratospheric warmings in reanalyses.  Atmospheric Chemistry and Physics, 19(14), 9469–9484.
URL https://acp.copernicus.org/articles/19/9469/2019/

  * Ayarzagüena et al. (2020) Ayarzagüena, B., Charlton-Perez, A., Butler, A., Hitchcock, P., Simpson, I., Polvani, L., Butchart, N., Gerber, E., Gray, L., Hassler, B., Lin, P., Lott, F., Manzini, E., Mizuta, R., Orbe, C., Osprey, S., Saint-Martin, D., Sigmond, M., Taguchi, M., Volodin, E., & Watanabe, S. (2020).  Uncertainty in the response of sudden stratospheric warmings and stratosphere-troposphere coupling to quadrupled co2 concentrations in cmip6 models.  Journal of Geophysical Research: Atmospheres, 125(6), e2019JD032345.  E2019JD032345 2019JD032345.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD032345

  * Baldwin et al. (2021) Baldwin, M. P., Ayarzagüena, B., Birner, T., Butchart, N., Butler, A. H., Charlton-Perez, A. J., Domeisen, D. I. V., Garfinkel, C. I., Garny, H., Gerber, E. P., Hegglin, M. I., Langematz, U., & Pedatella, N. M. (2021).  Sudden stratospheric warmings.  Reviews of Geophysics, 59(1), e2020RG000708.  E2020RG000708 10.1029/2020RG000708.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020RG000708

  * Baldwin & Dunkerton (2001) Baldwin, M. P., & Dunkerton, T. J. (2001).  Stratospheric harbingers of anomalous weather regimes.  Science, 294(5542), 581–584.
URL https://www.science.org/doi/abs/10.1126/science.1063315

  * Baldwin et al. (2003) Baldwin, M. P., Stephenson, D. B., Thompson, D. W. J., Dunkerton, T. J., Charlton, A. J., & O’Neill, A. (2003).  Stratospheric memory and skill of extended-range weather forecasts.  Science, 301(5633), 636–640.
URL https://www.science.org/doi/abs/10.1126/science.1087143

  * Benedetti (2010) Benedetti, R. (2010).  Scoring rules for forecast verification.  Monthly Weather Review, 138(1), 203 – 211.
URL https://journals.ametsoc.org/view/journals/mwre/138/1/2009mwr2945.1.xml

  * Berner et al. (2009) Berner, J., Shutts, G. J., Leutbecher, M., & Palmer, T. N. (2009).  A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ecmwf ensemble prediction system.  Journal of the Atmospheric Sciences, 66(3), 603 – 626\.
URL https://journals.ametsoc.org/view/journals/atsc/66/3/2008jas2677.1.xml

  * Bloomfield et al. (2021) Bloomfield, H. C., Brayshaw, D. J., Gonzalez, P. L. M., & Charlton-Perez, A. (2021).  Sub-seasonal forecasts of demand and wind power and solar power generation for 28 european countries.  Earth System Science Data, 13(5), 2259–2274.
URL https://essd.copernicus.org/articles/13/2259/2021/

  * Bouchet et al. (2019) Bouchet, F., Rolland, J., & Wouters, J. (2019).  Rare event sampling methods.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(8), 080402.
URL https://doi.org/10.1063/1.5120509

  * Broomhead & King (1986) Broomhead, D., & King, G. P. (1986).  Extracting qualitative dynamics from experimental data.  Physica D: Nonlinear Phenomena, 20(2), 217–236.
URL https://www.sciencedirect.com/science/article/pii/016727898690031X

  * Brunton et al. (2017) Brunton, S. L., Brunton, B. W., Proctor, J. L., Kaiser, E., & Kutz, J. N. (2017).  Chaos as an intermittently forced linear system.  Nature Communications, 8(1), 19.
URL https://doi.org/10.1038/s41467-017-00030-8

  * Buizza et al. (1999) Buizza, R., Milleer, M., & Palmer, T. N. (1999).  Stochastic representation of model uncertainties in the ecmwf ensemble prediction system.  Quarterly Journal of the Royal Meteorological Society, 125(560), 2887–2908.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.49712556006

  * Butler et al. (2019) Butler, A. H., Charlton-Perez, A., Domeisen, D. I., Garfinkel, C., Gerber, E. P., Hitchcock, P., Karpechko, A. Y., Maycock, A. C., Sigmond, M., Simpson, I., & Son, S.-W. (2019).  Chapter 11 - sub-seasonal predictability and the stratosphere.  In A. W. Robertson, & F. Vitart (Eds.) Sub-Seasonal to Seasonal Prediction, (pp. 223–241). Elsevier.
URL https://www.sciencedirect.com/science/article/pii/B9780128117149000115

  * Butler et al. (2015) Butler, A. H., Seidel, D. J., Hardiman, S. C., Butchart, N., Birner, T., & Match, A. (2015).  Defining sudden stratospheric warmings.  Bulletin of the American Meteorological Society, 96(11), 1913 – 1928.
URL
https://journals.ametsoc.org/view/journals/bams/96/11/bams-d-13-00173.1.xml

  * Charlton & Polvani (2007) Charlton, A. J., & Polvani, L. M. (2007).  A new look at stratospheric sudden warmings. part i: Climatology and modeling benchmarks.  Journal of Climate, 20(3), 449 – 469.
URL https://journals.ametsoc.org/view/journals/clim/20/3/jcli3996.1.xml

  * Chattopadhyay et al. (2020) Chattopadhyay, A., Nabizadeh, E., & Hassanzadeh, P. (2020).  Analog forecasting of extreme-causing weather patterns using deep learning.  Journal of Advances in Modeling Earth Systems, 12(2), e2019MS001958.  E2019MS001958 10.1029/2019MS001958.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001958

  * Chodera & Noé (2014) Chodera, J. D., & Noé, F. (2014).  Markov state models of biomolecular conformational dynamics.  Current Opinion in Structural Biology, 25, 135–144.  Theory and simulation / Macromolecular machines.
URL https://www.sciencedirect.com/science/article/pii/S0959440X14000426

  * Coumou & Rahmstorf (2012) Coumou, D., & Rahmstorf, S. (2012).  A decade of weather extremes.  Nature Climate Change, 2(7), 491–496.
URL https://doi.org/10.1038/nclimate1452

  * Deuflhard et al. (1999) Deuflhard, P., Dellnitz, M., Junge, O., & Schütte, C. (1999).  Computation of essential molecular dynamics by subdivision techniques.  In P. Deuflhard, J. Hermans, B. Leimkuhler, A. E. Mark, S. Reich, & R. D. Skeel (Eds.) Computational Molecular Dynamics: Challenges, Methods, Ideas, (pp. 98–115). Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Dimdore-Miles et al. (2021) Dimdore-Miles, O., Gray, L., & Osprey, S. (2021).  Origins of multi-decadal variability in sudden stratospheric warmings.  Weather and Climate Dynamics, 2(1), 205–231.
URL https://wcd.copernicus.org/articles/2/205/2021/

  * Domeisen et al. (2020) Domeisen, D. I. V., Butler, A. H., Charlton-Perez, A. J., Ayarzagüena, B., Baldwin, M. P., Dunn-Sigouin, E., Furtado, J. C., Garfinkel, C. I., Hitchcock, P., Karpechko, A. Y., Kim, H., Knight, J., Lang, A. L., Lim, E.-P., Marshall, A., Roff, G., Schwartz, C., Simpson, I. R., Son, S.-W., & Taguchi, M. (2020).  The role of the stratosphere in subseasonal to seasonal prediction: 2\. predictability arising from stratosphere-troposphere coupling.  Journal of Geophysical Research: Atmospheres, 125(2), e2019JD030923.  E2019JD030923 10.1029/2019JD030923.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD030923

  * E et al. (2019) E, W., Li, T., & Vanden-Eijnden, E. (2019).  Applied stochastic analysis, vol. 199.  American Mathematical Soc.
  * ECMWF (2016) ECMWF (2016).  IFS Documentation CY43R1 - Part V: Ensemble Prediction System.  No. 5 in IFS Documentation. ECMWF.
URL https://www.ecmwf.int/node/17118

  * ECMWF (2022a) ECMWF (2022a).  Copernicus climate data store.
URL https://cds.climate.copernicus.eu

  * ECMWF (2022b) ECMWF (2022b).  Ecmwf public datasets.  This service is based on data and products of the European Centre for Medium-Range Weather Forecasts (ECMWF). This ECMWF data is published under a Creative Commons Attribution 4.0 International (CC BY 4.0). https://creativecommons.org/licenses/by/4.0/. ECMWF does not accept any liability whatsoever for any error or omission in the data, their availability, or for any loss or damage arising from their use.
URL www.ecmwf.int

  * Finkel et al. (2020) Finkel, J., Abbot, D. S., & Weare, J. (2020).  Path properties of atmospheric transitions: Illustration with a low-order sudden stratospheric warming model.  Journal of the Atmospheric Sciences, 77(7), 2327 – 2347\.
URL https://journals.ametsoc.org/view/journals/atsc/77/7/jasD190278.xml

  * Finkel et al. (2021) Finkel, J., Webber, R. J., Gerber, E. P., Abbot, D. S., & Weare, J. (2021).  Learning forecasts of rare stratospheric transitions from short simulations.  Monthly Weather Review, 149(11), 3647 – 3669.
URL https://journals.ametsoc.org/view/journals/mwre/149/11/MWR-D-21-0024.1.xml

  * Finkel et al. (2022) Finkel, J., Webber, R. J., Gerber, E. P., Abbot, D. S., & Weare, J. (2022).  Data-driven transition path analysis yields a statistical understanding of sudden stratospheric warming events in an idealized model.  Journal of the Atmospheric Sciences.
URL
https://journals.ametsoc.org/view/journals/atsc/aop/JAS-D-21-0213.1/JAS-D-21-0213.1.xml

  * Fischer et al. (2021) Fischer, E. M., Sippel, S., & Knutti, R. (2021).  Increasing probability of record-shattering climate extremes.  Nature Climate Change, 11(8), 689–695.
URL https://doi.org/10.1038/s41558-021-01092-9

  * Frishman & Grafke (2022) Frishman, A., & Grafke, T. (2022).  Dynamical landscape of transitional pipe flow.  Phys. Rev. E, 105, 045108.
URL https://link.aps.org/doi/10.1103/PhysRevE.105.045108

  * Gerber et al. (2022) Gerber, E. P., Martineau, P., Ayarzaguena, B., Barriopedro, D., Bracegirdle, T. J., Butler, A. H., N. Calvo, S. C. H., Hitchcock, P., M. Iza, U. L., Lua, H., Marshall, G., Orr, A., Palmeiro, F. M., Son, S.-W., & Taguchi, M. (2022).  Extratropical stratosphere–troposphere coupling, in Stratosphere-troposphere Processes and their Role in Climate (SPARC) Reanalysis Intercomparison Project (S-RIP) (chap. 6), M. Fujiwara, G. L. Manney, L. Gray, and J. S. Wright, Eds.  SPARC Report No. 10, WCRP-6/2021, (pp. 221–264).
  * Ghil et al. (2002) Ghil, M., Allen, M. R., Dettinger, M. D., Ide, K., Kondrashov, D., Mann, M. E., Robertson, A. W., Saunders, A., Tian, Y., Varadi, F., & Yiou, P. (2002).  Advanced spectral methods for climatic time series.  Reviews of Geophysics, 40(1), 3–1–3–41.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2000RG000092

  * Giannakis & Majda (2012) Giannakis, D., & Majda, A. J. (2012).  Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability.  Proceedings of the National Academy of Sciences, 109(7), 2222–2227.
URL https://www.pnas.org/doi/abs/10.1073/pnas.1118984109

  * Goss et al. (2020) Goss, M., Swain, D. L., Abatzoglou, J. T., Sarhadi, A., Kolden, C. A., Williams, A. P., & Diffenbaugh, N. S. (2020).  Climate change is increasing the likelihood of extreme autumn wildfire conditions across california.  Environmental Research Letters, 15(9), 094016.
URL https://doi.org/10.1088/1748-9326/ab83a7

  * Hansen et al. (1984) Hansen, J., Lacis, A., Rind, D., Russell, G., Stone, P., Fung, I., Ruedy, R., & Lerner, J. (1984).  Climate sensitivity: Analysis of feedback mechanisms.  feedback, 1, 1–3.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., Friedman, J. H., & Friedman, J. H. (2009).  The elements of statistical learning: data mining, inference, and prediction, vol. 2.  Springer.
  * Helfmann et al. (2021) Helfmann, L., Heitzig, J., Koltai, P., Kurths, J., & Schütte, C. (2021).  Statistical analysis of tipping pathways in agent-based models.  The European Physical Journal Special Topics, 230(16), 3249–3271.
URL https://doi.org/10.1140/epjs/s11734-021-00191-0

  * Helfmann et al. (2020) Helfmann, L., Ribera Borrell, E., Schütte, C., & Koltai, P. (2020).  Extending transition path theory: Periodically driven and finite-time dynamics.  Journal of Nonlinear Science, 30(6), 3321–3366.
URL https://doi.org/10.1007/s00332-020-09652-7

  * Hersbach et al. (2020) Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., & Thépaut, J.-N. (2020).  The era5 global reanalysis.  Quarterly Journal of the Royal Meteorological Society, 146(730), 1999–2049.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3803

  * Hitchcock & Simpson (2014) Hitchcock, P., & Simpson, I. R. (2014).  The downward influence of stratospheric sudden warmings.  Journal of the Atmospheric Sciences, 71(10), 3856 – 3876\.
URL https://journals.ametsoc.org/view/journals/atsc/71/10/jas-d-14-0012.1.xml

  * Holton & Mass (1976) Holton, J. R., & Mass, C. (1976).  Stratospheric vacillation cycles.  Journal of Atmospheric Sciences, 33(11), 2218 – 2225\.
URL
https://journals.ametsoc.org/view/journals/atsc/33/11/1520-0469_1976_033_2218_svc_2_0_co_2.xml

  * Horan & Reichler (2017) Horan, M. F., & Reichler, T. (2017).  Modeling seasonal sudden stratospheric warming climatology based on polar vortex statistics.  Journal of Climate, 30(24), 10101 – 10116.
URL https://journals.ametsoc.org/view/journals/clim/30/24/jcli-d-17-0257.1.xml

  * Huntingford et al. (2014) Huntingford, C., Marsh, T., Scaife, A. A., Kendon, E. J., Hannaford, J., Kay, A. L., Lockwood, M., Prudhomme, C., Reynard, N. S., Parry, S., Lowe, J. A., Screen, J. A., Ward, H. C., Roberts, M., Stott, P. A., Bell, V. A., Bailey, M., Jenkins, A., Legg, T., Otto, F. E. L., Massey, N., Schaller, N., Slingo, J., & Allen, M. R. (2014).  Potential influences on the united kingdom’s floods of winter 2013/14.  Nature Climate Change, 4(9), 769–777.
URL https://doi.org/10.1038/nclimate2314

  * Jacques-Dumas et al. (2022) Jacques-Dumas, V., van Westen, R. M., Bouchet, F., & Dijkstra, H. A. (2022).  Data-driven methods to estimate the committor function in conceptual ocean models.  EGUsphere, 2022, 1–35.
URL https://egusphere.copernicus.org/preprints/egusphere-2022-1362/

  * Jucker et al. (2021) Jucker, M., Reichler, T., & Waugh, D. W. (2021).  How frequent are antarctic sudden stratospheric warmings in present and future climate?  Geophysical Research Letters, 48(11), e2021GL093215.  E2021GL093215 2021GL093215.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL093215

  * Kamb et al. (2020) Kamb, M., Kaiser, E., Brunton, S. L., & Kutz, J. N. (2020).  Time-delay observables for koopman: Theory and applications.  SIAM Journal on Applied Dynamical Systems, 19(2), 886–917.
URL https://doi.org/10.1137/18M1216572

  * Karpechko et al. (2017) Karpechko, A. Y., Hitchcock, P., Peters, D. H. W., & Schneidereit, A. (2017).  Predictability of downward propagation of major sudden stratospheric warmings.  Quarterly Journal of the Royal Meteorological Society, 143(704), 1459–1470.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3017

  * Kelder et al. (2020) Kelder, T., Müller, M., Slater, L. J., Marjoribanks, T. I., Wilby, R. L., Prudhomme, C., Bohlinger, P., Ferranti, L., & Nipen, T. (2020).  Using unseen trends to detect decadal changes in 100-year precipitation extremes.  npj Climate and Atmospheric Science, 3(1), 47.
URL https://doi.org/10.1038/s41612-020-00149-4

  * Kidston et al. (2015) Kidston, J., Scaife, A. A., Hardiman, S. C., Mitchell, D. M., Butchart, N., Baldwin, M. P., & Gray, L. J. (2015).  Stratospheric influence on tropospheric jet streams, storm tracks and surface weather.  Nature Geoscience, 8(6), 433–440.
URL https://doi.org/10.1038/ngeo2424

  * Kim et al. (2014) Kim, H.-M., Webster, P. J., Toma, V. E., & Kim, D. (2014).  Predictability and prediction skill of the mjo in two operational forecasting systems.  Journal of Climate, 27(14), 5364 – 5378.
URL
https://journals.ametsoc.org/view/journals/clim/27/14/jcli-d-13-00480.1.xml

  * King et al. (2019) King, A. D., Butler, A. H., Jucker, M., Earl, N. O., & Rudeva, I. (2019).  Observed relationships between sudden stratospheric warmings and european climate extremes.  Journal of Geophysical Research: Atmospheres, 124(24), 13943–13961.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD030480

  * Klus et al. (2018) Klus, S., Nüske, F., Koltai, P., Wu, H., Kevrekidis, I., Schütte, C., & Noé, F. (2018).  Data-driven model reduction and transfer operator approximation.  Journal of Nonlinear Science, 28(3), 985–1010.
URL https://doi.org/10.1007/s00332-017-9437-7

  * Kolstad et al. (2010) Kolstad, E. W., Breiteig, T., & Scaife, A. A. (2010).  The association between stratospheric weak polar vortex events and cold air outbreaks in the northern hemisphere.  Quarterly Journal of the Royal Meteorological Society, 136(649), 886–893.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.620

  * Kretschmer et al. (2018) Kretschmer, M., Cohen, J., Matthias, V., Runge, J., & Coumou, D. (2018).  The different stratospheric influence on cold-extremes in eurasia and north america.  npj Climate and Atmospheric Science, 1(1), 44.
URL https://doi.org/10.1038/s41612-018-0054-4

  * Lucente et al. (2022) Lucente, D., Rolland, J., Herbert, C., & Bouchet, F. (2022).  Coupling rare event algorithms with data-based learned committor functions using the analogue markov chain.  Journal of Statistical Mechanics: Theory and Experiment, 2022(8), 083201.
  * Meng et al. (2016) Meng, Y., Shukla, D., Pande, V. S., & Roux, B. (2016).  Transition path theory analysis of c-src kinase activation.  Proceedings of the National Academy of Sciences, 113(33), 9193–9198.
URL https://www.pnas.org/doi/abs/10.1073/pnas.1602790113

  * Meriwether & Gerrard (2004) Meriwether, J. W., & Gerrard, A. J. (2004).  Mesosphere inversion layers and stratosphere temperature enhancements.  Reviews of Geophysics, 42(3).
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2003RG000133

  * Mezić (2005) Mezić, I. (2005).  Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dynamics, 41(1), 309–325.
URL https://doi.org/10.1007/s11071-005-2824-x

  * Mezić (2013) Mezić, I. (2013).  Analysis of fluid flows via spectral properties of the koopman operator.  Annual Review of Fluid Mechanics, 45(1), 357–378.
URL https://doi.org/10.1146/annurev-fluid-011212-140652

  * Michelangeli et al. (1995) Michelangeli, P.-A., Vautard, R., & Legras, B. (1995).  Weather regimes: Recurrence and quasi stationarity.  Journal of Atmospheric Sciences, 52(8), 1237 – 1256.
URL
https://journals.ametsoc.org/view/journals/atsc/52/8/1520-0469_1995_052_1237_wrraqs_2_0_co_2.xml

  * Miloshevich et al. (2022) Miloshevich, G., Cozian, B., Abry, P., Borgnat, P., & Bouchet, F. (2022).  Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data.  arXiv preprint arXiv:2208.00971.
  * Miron et al. (2021) Miron, P., Beron-Vera, F. J., Helfmann, L., & Koltai, P. (2021).  Transition paths of marine debris and the stability of the garbage patches.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(3), 033101.
URL https://doi.org/10.1063/5.0030535

  * Miron et al. (2022) Miron, P., Beron-Vera, F. J., & Olascoaga, M. J. (2022).  Transition paths of north atlantic deep water.  Journal of Atmospheric and Oceanic Technology, 39(7), 959 – 971.
URL https://journals.ametsoc.org/view/journals/atot/39/7/JTECH-D-22-0022.1.xml

  * Mishra & Shah (2018) Mishra, V., & Shah, H. L. (2018).  Hydroclimatological perspective of the kerala flood of 2018.  Journal of the Geological Society of India, 92(5), 645–650.
URL https://doi.org/10.1007/s12594-018-1079-3

  * Myhre et al. (2019) Myhre, G., Alterskjær, K., Stjern, C. W., Hodnebrog, Ø., Marelle, L., Samset, B. H., Sillmann, J., Schaller, N., Fischer, E., Schulz, M., & Stohl, A. (2019).  Frequency of extreme precipitation increases extensively with event rareness under global warming.  Scientific Reports, 9(1), 16063.
  * Nakamura (1992) Nakamura, H. (1992).  Midwinter suppression of baroclinic wave activity in the pacific.  Journal of Atmospheric Sciences, 49(17), 1629 – 1642\.
URL
https://journals.ametsoc.org/view/journals/atsc/49/17/1520-0469_1992_049_1629_msobwa_2_0_co_2.xml

  * Naveau et al. (2020) Naveau, P., Hannart, A., & Ribes, A. (2020).  Statistical methods for extreme event attribution in climate science.  Annual Review of Statistics and Its Application, 7(1), 89–110.
URL https://doi.org/10.1146/annurev-statistics-031219-041314

  * Noé et al. (2009) Noé, F., Schütte, C., Vanden-Eijnden, E., Reich, L., & Weikl, T. R. (2009).  Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations.  Proceedings of the National Academy of Sciences, 106(45), 19011–19016.
URL https://www.pnas.org/doi/abs/10.1073/pnas.0905466106

  * O’Gorman (2012) O’Gorman, P. A. (2012).  Sensitivity of tropical precipitation extremes to climate change.  Nature Geoscience, 5(10), 697–700.
URL https://doi.org/10.1038/ngeo1568

  * Pande et al. (2010) Pande, V. S., Beauchamp, K., & Bowman, G. R. (2010).  Everything you wanted to know about markov state models but were afraid to ask.  Methods, 52(1), 99–105.  Protein Folding.
URL https://www.sciencedirect.com/science/article/pii/S1046202310001568

  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. (2011).  Scikit-learn: Machine learning in python.  the Journal of machine Learning research, 12, 2825–2830.
  * Penland & Sardeshmukh (1995) Penland, C., & Sardeshmukh, P. D. (1995).  The optimal growth of tropical sea surface temperature anomalies.  Journal of Climate, 8(8), 1999 – 2024.
URL
https://journals.ametsoc.org/view/journals/clim/8/8/1520-0442_1995_008_1999_togots_2_0_co_2.xml

  * Pfahl et al. (2017) Pfahl, S., O’Gorman, P. A., & Fischer, E. M. (2017).  Understanding the regional pattern of projected future changes in extreme precipitation.  Nature Climate Change, 7(6), 423–427.
URL https://doi.org/10.1038/nclimate3287

  * Ragone et al. (2018) Ragone, F., Wouters, J., & Bouchet, F. (2018).  Computation of extreme heat waves in climate models using a large deviation algorithm.  Proceedings of the National Academy of Sciences, 115(1), 24–29.
URL https://www.pnas.org/content/115/1/24

  * Reichler et al. (2012) Reichler, T., Kim, J., Manzini, E., & Kröger, J. (2012).  A stratospheric connection to atlantic climate variability.  Nature Geoscience, 5(11), 783–787.
URL https://doi.org/10.1038/ngeo1586

  * Sapsis (2021) Sapsis, T. P. (2021).  Statistics of extreme events in fluid flows and waves.  Annual Review of Fluid Mechanics, 53(1), 85–111.
URL https://doi.org/10.1146/annurev-fluid-030420-032810

  * Scaife et al. (2022) Scaife, A. A., Baldwin, M. P., Butler, A. H., Charlton-Perez, A. J., Domeisen, D. I. V., Garfinkel, C. I., Hardiman, S. C., Haynes, P., Karpechko, A. Y., Lim, E.-P., Noguchi, S., Perlwitz, J., Polvani, L., Richter, J. H., Scinocca, J., Sigmond, M., Shepherd, T. G., Son, S.-W., & Thompson, D. W. J. (2022).  Long-range prediction and the stratosphere.  Atmospheric Chemistry and Physics, 22(4), 2601–2623.
URL https://acp.copernicus.org/articles/22/2601/2022/

  * Shepherd et al. (2018) Shepherd, T. G., Boyd, E., Calel, R. A., Chapman, S. C., Dessai, S., Dima-West, I. M., Fowler, H. J., James, R., Maraun, D., Martius, O., Senior, C. A., Sobel, A. H., Stainforth, D. A., Tett, S. F. B., Trenberth, K. E., van den Hurk, B. J. J. M., Watkins, N. W., Wilby, R. L., & Zenghelis, D. A. (2018).  Storylines: an alternative approach to representing uncertainty in physical aspects of climate change.  Climatic Change, 151(3), 555–571.
URL https://doi.org/10.1007/s10584-018-2317-9

  * Sigmond et al. (2013) Sigmond, M., Scinocca, J. F., Kharin, V. V., & Shepherd, T. G. (2013).  Enhanced seasonal forecast skill following stratospheric sudden warmings.  Nature Geoscience, 6(2), 98–102.
URL https://doi.org/10.1038/ngeo1698

  * Sillmann et al. (2021) Sillmann, J., Shepherd, T. G., van den Hurk, B., Hazeleger, W., Martius, O., Slingo, J., & Zscheischler, J. (2021).  Event-based storylines to address climate risk.  Earth’s Future, 9(2), e2020EF001783.  E2020EF001783 2020EF001783.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020EF001783

  * Strahan et al. (2021) Strahan, J., Antoszewski, A., Lorpaiboon, C., Vani, B. P., Weare, J., & Dinner, A. R. (2021).  Long-time-scale predictions from short-trajectory data: A benchmark analysis of the trp-cage miniprotein.  Journal of Chemical Theory and Computation, 17(5), 2948–2963.  PMID: 33908762.
URL https://doi.org/10.1021/acs.jctc.0c00933

  * Takens (1981) Takens, F. (1981).  Detecting strange attractors in turbulence.  In Dynamical systems and turbulence, Warwick 1980, (pp. 366–381). Springer.
  * Thiede et al. (2019) Thiede, E. H., Giannakis, D., Dinner, A. R., & Weare, J. (2019).  Galerkin approximation of dynamical quantities using trajectory data.  The Journal of Chemical Physics, 150(24), 244111.
URL https://doi.org/10.1063/1.5063730

  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288.
URL
https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.2517-6161.1996.tb02080.x

  * Tripathi et al. (2016) Tripathi, O. P., Baldwin, M., Charlton-Perez, A., Charron, M., Cheung, J. C. H., Eckermann, S. D., Gerber, E., Jackson, D. R., Kuroda, Y., Lang, A., McLay, J., Mizuta, R., Reynolds, C., Roff, G., Sigmond, M., Son, S.-W., & Stockdale, T. (2016).  Examining the predictability of the stratospheric sudden warming of january 2013 using multiple nwp systems.  Monthly Weather Review, 144(5), 1935 – 1960.
URL https://journals.ametsoc.org/view/journals/mwre/144/5/mwr-d-15-0010.1.xml

  * Tseng et al. (2021) Tseng, K.-C., Johnson, N. C., Maloney, E. D., Barnes, E. A., & Kapnick, S. B. (2021).  Mapping large-scale climate variability to hydrological extremes: An application of the linear inverse model to subseasonal prediction.  Journal of Climate, 34(11), 4207 – 4225.
URL https://journals.ametsoc.org/view/journals/clim/34/11/JCLI-D-20-0502.1.xml

  * van den Dool (1989) van den Dool, H. M. (1989).  A new look at weather forecasting through analogues.  Monthly Weather Review, 117(10), 2230 – 2247.
URL
https://journals.ametsoc.org/view/journals/mwre/117/10/1520-0493_1989_117_2230_anlawf_2_0_co_2.xml

  * Van Oldenborgh et al. (2017) Van Oldenborgh, G. J., Van Der Wiel, K., Sebastian, A., Singh, R., Arrighi, J., Otto, F., Haustein, K., Li, S., Vecchi, G., & Cullen, H. (2017).  Attribution of extreme rainfall from hurricane harvey, august 2017.  Environmental Research Letters, 12(12), 124009.
  * Vanden-Eijnden (2014) Vanden-Eijnden, E. (2014).  Transition path theory, (pp. 91–100).  Advances in Experimental Medicine and Biology. Springer New York LLC.
  * Vitart (2014) Vitart, F. (2014).  Evolution of ecmwf sub-seasonal forecast skill scores.  Quarterly Journal of the Royal Meteorological Society, 140(683), 1889–1899.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.2256

  * Vitart et al. (2017) Vitart, F., Ardilouze, C., Bonet, A., Brookshaw, A., Chen, M., Codorean, C., Déqué, M., Ferranti, L., Fucile, E., Fuentes, M., Hendon, H., Hodgson, J., Kang, H.-S., Kumar, A., Lin, H., Liu, G., Liu, X., Malguzzi, P., Mallas, I., Manoussakis, M., Mastrangelo, D., MacLachlan, C., McLean, P., Minami, A., Mladek, R., Nakazawa, T., Najm, S., Nie, Y., Rixen, M., Robertson, A. W., Ruti, P., Sun, C., Takaya, Y., Tolstykh, M., Venuti, F., Waliser, D., Woolnough, S., Wu, T., Won, D.-J., Xiao, H., Zaripov, R., & Zhang, L. (2017).  The subseasonal to seasonal (s2s) prediction project database.  Bulletin of the American Meteorological Society, 98(1), 163 – 173.
URL https://journals.ametsoc.org/view/journals/bams/98/1/bams-d-16-0017.1.xml

  * Vitart & Robertson (2018) Vitart, F., & Robertson, A. W. (2018).  The sub-seasonal to seasonal prediction project (s2s) and the prediction of extreme events.  npj Climate and Atmospheric Science, 1(1), 3.
URL https://doi.org/10.1038/s41612-018-0013-0

  * Wasserman (2004) Wasserman, L. (2004).  All of statistics.  New York: Springer.
  * Webber et al. (2019) Webber, R. J., Plotkin, D. A., O’Neill, M. E., Abbot, D. S., & Weare, J. (2019).  Practical rare event sampling for extreme mesoscale weather.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(5), 053109.
URL https://doi.org/10.1063/1.5081461

  * White et al. (2017) White, C. J., Carlsen, H., Robertson, A. W., Klein, R. J., Lazo, J. K., Kumar, A., Vitart, F., Coughlan de Perez, E., Ray, A. J., Murray, V., Bharwani, S., MacLeod, D., James, R., Fleming, L., Morse, A. P., Eggen, B., Graham, R., Kjellström, E., Becker, E., Pegion, K. V., Holbrook, N. J., McEvoy, D., Depledge, M., Perkins-Kirkpatrick, S., Brown, T. J., Street, R., Jones, L., Remenyi, T. A., Hodgson-Johnston, I., Buontempo, C., Lamb, R., Meinke, H., Arheimer, B., & Zebiak, S. E. (2017).  Potential applications of subseasonal-to-seasonal (s2s) predictions.  Meteorological Applications, 24(3), 315–325.
URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/met.1654

  * Wigley (2009) Wigley, T. M. L. (2009).  The effect of changing climate on the frequency of absolute extreme events.  Climatic Change, 97(1), 67.
URL https://doi.org/10.1007/s10584-009-9654-7