# Bringing statistics to storylines: rare event sampling for sudden, transient
extreme events

###### Abstract

A leading goal for climate science and weather risk management is to
accurately model both the physics and statistics of extreme events. These two
goals are fundamentally at odds: the higher a computational model’s
resolution, the more expensive are the ensembles needed to capture accurate
statistics in the tail of the distribution. Here, we focus on events that are
localized in space and time, such as heavy precipitation events, which can
start suddenly and decay rapidly. We advance a method for sampling such events
more efficiently than straightforward climate model simulation. Our method
combines elements of two recent approaches: adaptive multilevel splitting
(AMS), a rare event algorithm that generates rigorous statistics at reduced
cost, but that does not work well for sudden, transient extreme events; and
“ensemble boosting” which generates physically plausible storylines of these
events but not their statistics. We modify AMS by splitting trajectories well
in advance of the event’s onset following the approach of ensemble boosting,
and this is shown to be critical for amplifying and diversifying simulated
events in tests with the Lorenz-96 model. Early splitting requires a rejection
step that reduces efficiency, but nevertheless we demonstrate improved
sampling of extreme local events by a factor of order 10 relative to direct
sampling in Lorenz-96. Our work makes progress on the challenge posed by fast
dynamical timescales for rare event sampling, and it draws connections with
existing methods in reliability engineering which, we believe, can be further
exploited for weather risk assessment.

Journal of Advances in Modeling Earth Systems (JAMES)

Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts
Institute of Technology

Justin<EMAIL_ADDRESS>
Rare event algorithms may help address the challenge of simulating extreme
weather events and quantifying their probability.

When the event of interest is sudden and transient, perturbed ensembles
diversify too slowly for standard rare event algorithms to work.

Using the Lorenz-96 model as a prototype for midlatitude weather, we use early
perturbation and a rejection step to gain a speedup.

## Plain Language Summary

What is the strongest rainstorm that we can expect in a given thousand-year
period? To augment the available $\sim 100$ years of historical data and to
account for climate change, computer simulations are a useful, but expensive,
tool to answer such questions. A model must run for many millennia to deliver
an answer with statistical confidence. _Rare event algorithms_ provide a
promising alternative simulation protocol, in which an ensemble of short
simulations is biased to produce more extreme events and reweighting is used
to correct for the bias when calculating statistics. However, a classical rare
event algorithm fails when the events of interest are short and “bursty” (like
heavy rainstorms) instead of long and slow-moving (like anomalously hot
summers). We modify the rare event algorithm to make it amenable to
precipitation-like events in an idealized dynamical system with chaotic
traveling waves.

##  1 Introduction

In climate modeling, high spatial resolution is important for realistically
representing localized extreme weather events like cyclones producing extreme
precipitation and winds [O’Brien . (2016), van der Wiel . (2016)]. But given
finite computational resources, high resolution has to be traded off with the
need for ensembles of models and simulations to deal with uncertainty related
to model physics, parameters, initial conditions and boundary conditions
including emissions scenarios. Extreme events are particularly challenging
because they occur infrequently, and hence need large ensemble sizes to have
their small probabilities accurately quantified. The conflict for
computational resources therefore comes to a head in the study of extreme
events.

A variety of shortcuts have developed in the past century to alleviate this
conflict. Leading statistical approaches include extreme value theory
[<]EVT;¿[]Coles2001introduction and large deviation theory [Touchette (2009)],
which respectively describe the behavior of _maxima_ and anomalously large
_running means_ in random processes. In principle, we can use these theories
to fit a parametric family to limited data and then extrapolate to even longer
return periods. EVT has become an important tool in risk assessment and
climate change attribution [Kharin . (2007), Naveau . (2020)], while large
deviation theory succinctly encodes the severity of long-lasting, large-area
events such as persistent heat waves [Gálfi . (2021)]. Statistical theories
help make the most of a fixed dataset, but parameter estimation can be
unstable given the restrictive underlying assumptions and the limited datasets
available [WK. Huang . (2016), Gálfi . (2017)]. For example, EVT only holds in
the limit of large blocks of data or high thresholds for extremity, which
directly conflicts with the requirement of many samples for low-variance
parameter estimation. Moreover, statistical theories don’t provide spatio-
temporal resolved extreme events (e.g., the spatial field of rainfall and
other fields on the day of an extreme event) which are needed to drive impact
models.

Statistical or dynamical downscaling is another way to address the problem of
extremes by reducing the computational cost of obtaining high-resolution
output from long simulations or large ensembles [X. Huang . (2020), Lee .
(2020), Emanuel (2021), Saha Ravela (2022), Krouma . (2022)]. Downscaling
nevertheless has some drawbacks. Dynamical downscaling using regional climate
models faces the challenge of correctly forcing a regional model with output
from a different global model, and the regional model inherits errors in
large-scale fields from the global model [Adachi Tomita (2020)], while
statistical downscaling assumptions can create systematic errors [Schmidli .
(2007)] and may not generalize to different climates.

The focus of this paper is _rare event sampling_ , which is a strategy for
allocating more of the computational effort towards rare events, and less
effort towards the long intervening periods of comparatively mild behavior.
This is usually achieved by _splitting_ methods, which consist of three steps
repeated in a cycle: (1) run an ensemble of simulations forward, (2) identify
the ensemble members making the most progress towards the extreme event, and
(3) clone these most-promising ensemble members (applying small perturbations)
while discarding the less-promising members, resulting in a new ensemble that
is more prone to extremes than was the original ensemble. With repeated rounds
of splitting, one can populate the tail of the probability distribution more
fully, while neglecting the more typical behavior of lesser interest.
Crucially, in statistical analysis of the ensemble, one must compensate for
the bias by weighting each clone with a factor less than one, relying on the
_importance sampling_ formalism. See Bucklew2004introduction for an
introduction to rare event sampling.

This generic procedure has many possible variants, which have been developed
largely in the fields of physics [Kahn Harris (1951), Giardinà . (2006)],
chemistry [Kästner (2011), Zuckerman Chong (2017)], and reliability
engineering [Au Beck (2001)], but have recently started to make an impact on
Earth and planetary sciences. For example, extreme European heat waves were
sampled by Ragone2018computation and Ragone2021rare with genealogical particle
analysis (GPA), and by Yiou2020simulation with empirical importance sampling.
Wouters2023rare sampled extreme European seasonal precipitation accumulations,
also using GPA. Webber2019practical developed a quantile-based variant of GPA
to sample more extreme versions of tropical cyclones. Planetary science
applications include jet nucleation [Bouchet . (2019)] and orbit
destabilization [Abbot . (2021)]. For studies of climate, rare event sampling
can be applied to global models or paired with the dynamical and statistical
downscaling approaches mentioned earlier.

We have elected to use a particular rare event algorithm called _adaptive
multilevel splitting_ (AMS), which was first established by Cerou2007adaptive
and is similar to the earlier RESTART algorithm [Villén-Altamirano . (1991)].
Lestang2018computing successfully applied AMS to the Ornstein-Ulhenbeck
process, while Lucente2022coupling and Baars2021application used AMS to study
regime transitions in idealized climate models. AMS has also been usefully
employed in other diverse fields such as molecular dynamics and air traffic
control (see Cerou2019adaptive for a recent review). The distinguishing
feature of AMS is that it operates on the level of full trajectories over a
fixed time horizon, and applies the small perturbation to trajectories at the
instant that they first cross a threshold of extremity. The “child” trajectory
is identical to its parent up until this time, whereas it diverges from its
parent afterward to give a new realization of the extreme event. All ensemble
members failing to cross the threshold are discarded, and the threshold is
then raised for repeated rounds of splitting and killing.

A related approach, “ensemble boosting”, is a novel technique for generating
“storylines” of unprecedented climate extremes [Gessner . (2021), Gessner
(2022)]. In this approach, one identifies several extreme events from a long
climate simulation, perturbs the antecedent conditions (1-3 weeks ahead of
time), and re-simulates the event to generate alternative realities, which
sometimes turn out even more extreme. While similar to splitting methods,
ensemble boosting does not explicitly quantify statistics. As explained below,
a major goal of this paper is to combine the benefits of ensemble boosting
with that of rare event algorithms, in particular AMS.

Given the successes in using rare event sampling discussed above, it is
desirable to also use it to sample shorter-term extreme weather events, such
as daily precipitation extremes, which have large societal impacts in the
current climate [Wright . (2021), Thompson . (2017)] and are expected to
intensify under climate change [O’Gorman (2015), Pfahl . (2017), Tandon .
(2018), Myhre . (2019)]. However, heavy precipitation events (or high wind
events) have some dynamical characteristics that distinguish them from the
previous applications and pose challenges to existing rare event algorithms.
Unlike continental-scale, seasonally averaged anomalies studied previously
[Ragone . (2018), Wouters . (2023)], heavy precipitation events of interest
are often sudden, transient, and relatively small-scale. Their timescale at a
particular location is often limited by the propagation of the dynamical
feature causing the precipitation such as cyclones and fronts [Dwyer O’Gorman
(2017)]. The strategy used in Ragone2018computation and Wouters2023rare relies
on some slow-moving notion of _progress_ towards the extreme event, naturally
given by the integrated temperature anomaly itself when targeting extreme
seasonal average temperatures, in order to decide which simulations to clone
or kill. In the precipitation study of Wouters2023rare, the extreme event is
again a seasonal total, for which a mid-seasonal total is a reasonable measure
of progress. But for individual precipitation events, if one uses
precipitation itself to measure progress towards the event, and applies
perturbations to a simulation when precipitation picks up, it is too late for
these perturbations to take effect by the time of maximum precipitation. The
event simply comes and goes faster than perturbed simulations diverge.
Lestang2018computing found a similar pathology with AMS when sampling extreme
pressure fluctuations on a body embedded in a turbulent channel flow. There,
the extreme events were caused by vortices sweeping past the body, roughly
analogous to cyclones sweeping past a location on Earth, and the rapidity of
the fluctuation crippled the effectiveness of the standard splitting strategy.

To isolate and solve the problem of applying rare event algorithms to sudden,
transient extremes, we postpone the specific application to precipitation and
first descend the model hierarchy to the Lorenz-96 model [Lorenz (1996)], a
spatiotemporal chaotic system often used as a toy model for the atmosphere.
The model produces extreme events posing the same algorithmic challenges as
precipitation extremes: intermittent, short-lived bursts carried by traveling
waves with unpredictable amplitudes. It has been used in numerous past studies
of extreme event statistics and predictability [Sterk van Kekem (2017), Qi
Majda (2016), Hu . (2019)]. With this cheap but behaviorally rich model, we
have developed a simple modification to AMS, drawing inspiration from ensemble
boosting by simply applying a split in advance of the event’s onset by some
advance split time $\delta$—hence, “trying early” AMS (TEAMS). To make this
statistically rigorous, a rejection step is necessary, which comes at an
efficiency cost, but still enables moderate speedups of $\sim$10 relative to
direct sampling. Fig. 1 displays a schematic diagram for TEAMS, which will be
elaborated in section 3. In fact, TEAMS is a repurposing of a more general
method called _subset simulation_ [Au Beck (2001)] from structural reliability
engineering, a field whose sophisticated rare event algorithms could benefit
the climate risk community.

Figure 1: Schematic of the splitting step in (a) AMS and (b) TEAMS. Black
curves represent an initial ensemble member, or ancestor, which exceeds the
first level $\ell_{1}$ and has been selected for cloning in the first round.
In AMS, the perturbation is applied at the instant $t_{0}(\ell_{1})$ when the
ancestor first exceeds $\ell_{1}$, resulting in a descendant trajectory (blue)
which essentially replicates the extreme event because the separation
timescale is longer than the event itself. On the other hand, in TEAMS (right)
we apply the perturbation in advance, by some margin $\delta>0$. This can
sometimes result in rejection (blue descendant), i.e., failure to cross
$\ell_{1}$. However, when a descendant is accepted (red) it will be more
distinct from the ancestor than the corresponding descendant in AMS and have
the potential to reach a substantially higher peak value.

This paper is organized as follows. In section 2, we present a stochastically
forced Lorenz-96 model and the behavior of its extreme events as a function of
stochastic forcing strength. In section 3, we first introduce the general
framework of subset simulation. In section 3.1, we specialize to AMS, and in
section 3.2 we show that AMS fails in the low-noise forcing regime, which is
often most relevant for weather and climate models. In section 3.3, we modify
AMS to use a “trying early” step with rejection sampling and recover a
substantial speedup. In section 4, we further explore the relationship between
the advance splitting time—a key algorithmic parameter—and classical notions
of predictability timescales. Finally, in section 5 we point out directions
for further development.

##  2 Lorenz-96: a customizable spatiotemporal chaotic system

Table 1: Physical parameters for Lorenz-96 system (upper section), and algorithmic parameters for the TEAMS algorithm (lower section). Symbol | Explanation | Value or range
---|---|---
$K$ | Number of longitude sites | 40
$a$ | Strength of advection term |  $\\{1,0\\}$ (mostly 1)
$F_{0}$ | Constant background forcing | 6
$m$ | Wavenumber for stochastic forcing |  $\\{1,4,7,10\\}$ (mostly 4)
$F_{m}$ | Strength of stochastic forcing at wavenumber $m$ | $\\{3,1,0.5,0.25,0\\}$
$N$ | Number of initial ensemble members | 128
$\kappa$ | Number of members to kill each round | 1
$J$ | Number of rounds of splitting | 896
$T$ | Time horizon | 6
$\delta$ | Advance split time | $[0,2]$

Lorenz1996predictability introduced a simple dynamical system (L96 hereafter)
meant to capture some crucial aspects of atmospheric dynamics. The model state
space consists of $K$ ($\geq 4$) variables $\\{x_{k}\\}_{k=1}^{K}$ arranged on
a one-dimensional periodic lattice, each $k$ representing a longitude sector
on Earth. $x_{k}$ represents a generic atmospheric variable like wind speed or
vorticity and evolves according to the coupled equations

$\frac{dx_{k}}{dt}=ax_{k-1}(x_{k+1}-x_{k-2})-x_{k}+\mathcal{F}_{k},\hskip
28.45274ptk=0,...,K-1,$ (1)

where $x_{k+K}$ is identified with $x_{k}$. The quadratic terms on the right-
hand side represent advection, like the quadratic nonlinearity in the material
derivative of the Navier-Stokes equations, which on its own conserves “energy”
$\frac{1}{2}\sum_{k}x_{k}^{2}$. The linear term $-x_{k}$ represents damping
due to friction, and the additive term $\mathcal{F}_{k}$ represents external
forcing, like a meridional insolation gradient. The latter two terms destroy
exact energy conservation, but balance out in a time-averaged sense to make
for a statistically steady state. Lorenz1996predictability introduced the
above model with $\mathcal{F}_{k}$ constant in $k$ and also a version in which
$\mathcal{F}_{k}$ is a “subgrid-scale forcing” that is a function of an
additional tier of dynamical variables representing finer scales, and this
version has proven useful for testing stochastic parameterization schemes
[<]e.g.,¿[]Wilks2005effects,Hu2019effects,Gagne2020machine. Here, we also
allow $\mathcal{F}_{k}$ to vary stochastically with longitude ($k$) and time:

$\mathcal{F}_{k}=F_{0}+F_{m}\bigg{[}\eta_{1}\cos\bigg{(}\frac{2\pi
mk}{K}\bigg{)}+\eta_{2}\sin\bigg{(}\frac{2\pi mk}{K}\bigg{)}\bigg{]}$ (2)

where $\eta_{1,2}$ are independent Gaussian white-noise processes, and $m$ is
an integer wavenumber. Formally, Eq. (2) renders Eq. (1) a diffusion process,
using the Itô convention for stochastic integrals [Pavliotis (2014)]. This
simple stochastic forcing is analagous to a stochastic parameterization in a
weather or climate model, and in the AMS framework it allows us to easily
generate new ensemble members by splitting an existing ensemble member at a
certain time. We verify below that for weak amplitudes the stochastic forcing
does not appreciably alter model statistics.

The parameters used here are summarized in the upper section of Table 1. We
set $K=40$, following Lorenz1998optimal. We fix the constant part of the
forcing to be $F_{0}=6.0$, which is sufficient for weak turbulence (a larger
value would be needed with smaller $K$). We choose the stochastic forcing
wavenumber as $m=4$ because that empirically seems to drive ensemble members
apart slightly faster than very small or large wavenumbers (see section 4.2).
Indeed the stochastically perturbed parameterization tendencies (SPPT) method
developed at ECMWF uses noise that is spatially correlated at a $\sim
10^{\circ}$ length scale [Buizza . (1999), Palmer . (2009)]. The amplitude of
$F_{m}(=F_{4})$ will be explored systematically below. One further parameter,
the coefficient $a$, determines the strength of the advection term. $a=1$ is
standard for L96, while $a=0$ gives an array of correlated Ornstein-Uhlenbeck
(OU) processes [Pavliotis (2014)]. Retaining the OU process as a special case
of L96 is useful to provide a reference case on which existing rare event
splitting algorithms excel. Results for $a=0$ are shown in supplementary Figs.
S1 and S2, and all other results presented are for $a=1$.

Fig. 2 displays short numerical integrations of L96 with four different
parameter choices. We used the Euler-Maruyama method with a timestep of 0.001
to integrate Eq. (1), saving out every 0.05 time units. For comparison,
Lorenz1998optimal interpret a single time unit as 5 days. The left column
shows single-site variables $x_{0}(t)$ for each parameter set, while the right
column shows corresponding Hovmöller diagrams. In the standard deterministic
system $F_{4}=0$ in the top row, $x_{0}(t)$ fluctuates with a semi-regular
period of $\sim 2$ time units (10 “days”) but with irregular amplitudes, the
largest of which are precisely the extreme events we choose to study here. The
Hovmöller diagram reveals these fluctuations to arise from a field of
traveling waves, with roughly eight peaks and troughs moving with negative
(“westward”) phase velocity. The waves experience intermittent disturbances,
sometimes getting stuck in place for several turnover times and setting up
favorable conditions for extreme events. Globally, these stagnations manifest
as kinks that propagate in the positive (“eastward”) direction. This is
reminiscent of atmospheric Rossby waves, whose phase and group velocities have
opposite signs (up to a Doppler shift due to the mean flow) [Lorenz Emanuel
(1998)]. Thus, we can loosely think of the peaks and troughs as being like
highs and lows in the midlatitude atmosphere.

Fig. 2 rows 2 and 3 show analogous pictures for moderate ($F_{4}=1$) and
strong $(F_{4}=3)$ stochastic forcing, respectively. As noise increases the
traveling waves transition from unidirectional to zigzagging. The timeseries
become more jagged and more liable to take large excursions from their mean
and hover there for longer durations.

Figure 2: Time evolution of the L96 model expressed as timeseries of
$x_{0}(t)$ (left column) and Hovmöller diagrams (right column) with three
different levels of stochastic forcing. (a,b) have $F_{4}=0$ (the
deterministic system); (c,d) have $F_{4}=1$ (moderate forcing); (e,f) have
$F_{4}=3$ (strong forcing).

Fig. 3a overlays PDFs of the single-site value ($x_{0}$) for all these
parameter regimes, plus two more: $F_{4}=0.5$ and $0.25$. Reducing the noise
roughly preserves the mode but shrinks the tails. The PDF appears basically
converged for $F_{4}\leq 0.5$. Fig. 3b confirms this is true even in the far
tail, with a log-transformed plot of return level vs. return time for
$x_{0}^{2}$. The limiting case $F_{4}=0$ has a bounded tail, which is easy to
see with an energy argument (see also Qi2016predicting): defining
$\overline{x}=\frac{1}{K}\sum_{k=1}^{K}x_{k}$, the energy
$E=\frac{1}{2}\sum_{k}x_{k}^{2}$ evolves as
$\frac{dE}{dt}=-2E+KF\overline{x}$. Since
$|\overline{x}|\leq\sqrt{\overline{x^{2}}}=\sqrt{2E/K}$ by the Cauchy-Schwarz
inequality, the first term dominates for $E$ larger than some critical
$E_{0}$, which must therefore bound the steady-state distribution’s tail.
However, $E_{0}$ would increase with $K$, i.e., higher-dimensional systems can
in principle support heavier tails [<]e.g.¿[ch. 4 discusses general
relationships between the shape parameter and the attractor
dimension]Lucarini2016extremes. This is part of our motivation to set $K$
relatively large.

The return level vs. return period plot (as in Fig. 3b) will be used
throughout the paper, and we calculate it using the “modified block maximum”
method of Lestang2018computing. For a fixed return level $\ell$, the return
period $\tau(\ell)$ is defined as the mean (over initial conditions and noise
realizations) of the waiting time until an exceedance occurs:
$\tau(\ell)=\mathbb{E}[\min\\{t:R(x(t))>\ell\\}$], where $R$ is some
observable of interest for the dynamical system. We take $R(x)=x_{0}^{2}$, the
local energy (times two) at longitude $k=0$. Lestang2018computing approximates
the exceedance times by a Poisson process for high $\ell$ to give

$\tau(\ell)=-\frac{T}{\log\big{[}1-p_{T}(\ell)\big{]}}.$ (3)

where $p_{T}(\ell)$ is the probability of at least one exceedance in a fixed
time $T$. $p_{T}(\ell)$ can be estimated from any collection of length-$T$
blocks of data— _either from a single continuous timeseries or not_. This is
very useful because rare event splitting algorithms generate branching trees
of short trajectories, from which we can estimate block-wise exceedances but
not return times directly.

To produce Fig. 3b, we started with simulations of length $1.28\times 10^{6}$
(after discarding the first 50 for spinup), split them into $B$ blocks of
length $T=6$, and measure the maxima $M_{1},...,M_{B}$ of $x_{0}^{2}$ over
each block. Letting $M_{(b)}$ denote the $b$th largest block maximum, we use
the empirical (complementary) CDF estimator, $\widehat{p}_{T}(M_{(b)})=b/B$.
Hence, the return curve should interpolate the ordered pairs
$(\tau_{b},\ell_{b})=\big{(}-\frac{T}{\log(1-b/B)},M_{b}\big{)}$. Because it
is common to think of $\ell$ as a function of $\tau$, and to consider
logarithmically spaced return periods, we linearly interpolate $M_{(b)}$ over
$\log\tau_{B}$ to get a curve $\hat{\ell}(\tau)$. We bootstrap to estimate
uncertainty, resampling the blocks $1,...,B$ with replacement and repeating
the above procedure 5000 times. Shading indicates the basic bootstrap 95%
confidence interval [Wasserman (2004)], meaning
$\hat{\ell}(\tau)+(\hat{\ell}(\tau)-\ell^{*}_{0.975}(\tau),\hat{\ell}(\tau)-\ell^{*}_{0.025}(\tau))$,
where $\ell^{*}_{\alpha}$ denotes the $\alpha$th quantile of the bootstrap
distribution of $\hat{\ell}$ for each $\tau$. Note that when
$\ell^{*}_{0.025}(\tau)$ is much less than $\hat{\ell}(\tau)$, we get a very
large _upper_ bound on the confidence interval, because it suggests via the
basic bootstrap philosophy that $\hat{\ell}(\tau)$ could be very much less
than the true parameter $\ell(\tau)$. The lowest-noise curves are close to
within uncertainty even in the far tails, demonstrating the convergence of
extreme value statistics for $F_{4}\leq 0.5$. This confirms that stochastic
forcing, when sufficiently weak, does not alter the system’s statistics very
much, which allows us to approximate the deterministic system’s rare events
while remaining within the AMS framework which relies on explicit randomness.

The longest return period estimable by this method of “direct numerical
simulation” (DNS) is $\sim 8\times 10^{5}$, the simulation’s length. Rare
event algorithms can sample physical realizations of extreme events at long
return periods $\tau(\ell)$ with much less computation time than $\tau(\ell)$,
but have not yet been applied to local events in L96 with weak stochastic
forcing. Wouters2016rare did apply rare event algorithms to L96, but their
system parameters differed substantially from ours, with $F_{0}=256$ giving a
much more turbulent regime reminiscent of a stochastic process. Moreover,
their target quantity of interest was a globally averaged energy, whereas we
target local energy at one longitude as a closer analogue to extreme
precipitation or winds hitting a particular location.

Figure 3: Steady-state statistics of the L96 model as a function of noise
strength, calculated from a long simulation of length $1.28\times 10^{6}$. (a)
Histograms of the model variable at one site ($x_{0}$) and (b) return level
vs. return period for (twice) the local energy $x_{0}^{2}$. Shading in (b)
represents 95% bootstrapped confidence intervals from the modified block
maximum method. See text for details.

The parameters $a$ and $F_{4}$ allow us to test the performance of AMS for a
range of problems, from systems on which AMS performs well to more difficult
systems akin to the extreme local precipitation problem. $a=0$ (the OU
process) is an easy setting for AMS; $a=1$ with large noise $F_{4}$ is harder,
but still doable because of the dominance of noise. Shrinking $F_{4}$ further,
towards the system of actual interest, gradually renders standard AMS
ineffective and leads us to a modified version of the algorithm called TEAMS
that allows for early splitting. The next sections present the basic algorithm
and its modification along this parameter path.

##  3 Subset simulation

TEAMS (and the special case AMS) may be viewed as a version of _subset
simulation_ (SS), which we use to frame our overall approach, and which we
believe has considerable potential for application to climate problems. SS was
introduced in Au2001estimation and has been most widely used in structural
reliability engineering [X. Huang . (2016)]. For a short pedagogical
introduction, see Zuev2015subset. The description below will introduce several
tunable algorithmic parameters, which are summarized in the lower section of
Table 1.

The goal is to estimate the probability that a random variable $x$ from a
distribution $\rho$ gives rise to large values of some quantity of interest
$S(x)$,

$p(\ell)=\int\mathbb{I}\\{S(x)>\ell\\}\rho(x)\,dx=\mathbb{E}_{\rho}\big{[}\mathbb{I}\\{S(X)>\ell\\}\big{]},$
(4)

given only the ability to draw samples $X_{1},X_{2},...\sim\rho$.
$\mathbb{I}\\{\cdot\\}$ denotes the indicator function: one if the argument is
true, zero if false. For us, each $X_{i}=\\{X_{i}(t):0\leq t\leq T\\}$ is a
length-$T$ trajectories of L96 (with stochastic forcing); the score function
is a maximum over the interval, $S(X)=\max_{0\leq t<T}R(X(t))$; and $\rho(x)$
is the distribution over _trajectories of length_ $T$ induced by the
stochastically forced L96 system. In structural engineering, $X$ might be the
state of a building or dam, with $\rho(x)$ induced by a probability
distribution over external stresses like wind, earthquakes, or rainfall, while
$S(x)$ would measure the proximity to failure. Because the probabilities of
interest are very small, a set of independent samples $\\{X_{n}\\}_{n=1}^{N}$
from $\rho$ will usually have few if any exceedances, making the “vanilla”
Monte Carlo estimate of $p(\ell)$ (the fraction of exceedances) subject to
high relative uncertainty. The ratio of the estimator’s variance to its mean
is approximately $1/\sqrt{Np(\ell)}$ [Zuev (2015)]. If we want to aim for a
tenfold-longer return period with the same uncertainty, we need to generate
tenfold more samples. Worse, to reduce uncertainty tenfold we would need one
hundredfold more samples, which may be untenable.

SS breaks down this task into a sequence of easier tasks by setting up a
series of intermediate levels $\ell_{1}<\ell_{2}<...<\ell_{J}=\ell$ where $J$
is the number of levels, and estimating a sequence of conditional
probabilities
$\mathbb{P}\\{S(X)>\ell_{j+1}|S(X)>\ell_{j}\\}=:p(\ell_{j+1}|\ell_{j})$, which
all have moderate magnitudes and are expected to be easier to estimate. Their
product provides an estimate for the target probability:

$\hat{p}_{\mathrm{SS}}(\ell)=\hat{p}(\ell_{1})\hat{p}(\ell_{2}|\ell_{1})...\hat{p}(\ell_{J}|\ell_{J-1}).$
(5)

The first term can be estimated by vanilla Monte Carlo: generate $N$ samples
$X_{1},...,X_{N}$, and attach unit weights to each: $W_{n}=1$ for $n=1,...,N$.
Rank the samples by $S$ so that $S(X_{(1)})\leq S(X_{(2)})\leq...\leq
S(X_{(N)})$, and let $\hat{p}(\ell_{1})=(N-\kappa_{1})/N$, where $\kappa_{1}$
is chosen so that $S(X_{(\kappa_{1})})\leq\ell_{1}<S(X_{(\kappa_{1}+1)})$. The
parameter $\kappa_{1}$ is the number of trajectories that are “killed” meaning
they don’t appear in the first subset (see below). For the case of AMS,
$\kappa_{1}$ is chosen as a parameter of the algorithm, and $\ell_{1}$ is then
set adaptively as
$\ell_{1}=\frac{1}{2}[S(X_{(\kappa_{1})})+S(X_{(\kappa_{1}+1)})]$.

The second term $\hat{p}(\ell_{2}|\ell_{1})$ is estimated with a splitting
strategy in which we focus in on the “subset” of samples that exceed the first
threshold: $\\{S(X)>\ell_{1}\\}$ containing samples $X_{(i)}$ with
$\kappa_{1}<i\leq N$. To better sample this subset, we spawn additional
samples from it via a “Modified Metropolis algorithm”:

  1. 1.

Initialize a list $\mathbb{X}_{1}=\\{X_{(\kappa_{1}+1)},...,X_{(N)}\\}$, which
will eventually grow to a (user-chosen) size $N_{1}$ as well as a first-in-
first-out queue $\mathbb{Q}$ of the same elements but in a random order: the
“parent queue”.

  2. 2.

Pop $\mathbb{Q}$ to yield the next parent $X$. Apply some small perturbation
to $X$ to generate a new sample $\widetilde{X}$, which itself is drawn from
$\rho$ but correlated to $X$. A general way to do this is with one step of the
Metropolis-Hastings algorithm which involves an accept/reject step, but an
easier approach is available in the particular case of AMS as described in the
next section.

  3. 3.

Evaluate $S(\widetilde{X})$. If it exceeds $\ell_{1}$, we have successfully
generated a new sample from the subset. Accept the new sample, meaning insert
$\widetilde{X}$ into both $\mathbb{Q}$ and $\mathbb{X}_{1}$ and assign it a
weight equal to that of its parent $X$. Otherwise, if
$S(\widetilde{X})\leq\ell_{1}$, reject $\widetilde{X}$. Re-insert $X$ into
$\mathbb{Q}$ and add a copy of $X$ to $\mathbb{X}_{1}$. (In implementation, we
don’t store two copies of the high-dimensional object $X$, but rather we
assign a multiplicity to each member and increment $X$’s multiplicity by one.)

  4. 4.

Repeat steps 2 and 3 until $\mathbb{X}_{1}$ has $N_{1}$ elements (counting
multiplicity).

  5. 5.

Multiply the weights of all members of $\mathbb{X}_{1}$ by a factor
$(N-\kappa_{1})/N_{1}$, which preserves the total weight $N$ of the original
ensemble $\\{X_{n}\\}_{n=1}^{N}$ while spreading that weight over more
members.

Having expanded to $N_{1}$ samples from the subset $\\{S(X)>\ell_{1}\\}$, we
can now proceed to the next level and generate additional samples from the
next subset $\\{S(X)>\ell_{2}\\}$ so that it contains $N_{2}$ samples, where
$\ell_{2}$ can be determined adaptively as an order statistic of
$\mathbb{X}_{1}$, i.e., the average of the $\kappa_{2}$th and the
$(\kappa_{2}+1)$th ranked values. The same procedure is repeated to generate
the next subset $\mathbb{X}_{2}$ (and $\mathbb{Q}$ is initialized with only
unique elements, not counting multiplicity, in order to maintain as much
diversity as possible). $\mathbb{X}_{3},\mathbb{X}_{4},...,\mathbb{X}_{J}$ are
generated in the same fashion, until either a computational budget is reached,
an ultimate target threshold is overcome, or some other halting criterion is
met.

Ultimately we are left with a weighted ensemble
$\\{(X_{1},W_{1}),...,(X_{M},W_{M})\\}$, where
$M=\kappa_{1}+\kappa_{2}+...+\kappa_{J}+N_{J}$. The sampling
$\\{S(X_{m})\\}_{m=1}^{M}$ is over-represented in the tails, but with
correspondingly smaller weights there, and all weights sum to $N$. Any
expectation of an observable $\Phi(x)$ can be estimated as

$\mathbb{E}[\Phi(X)]=\int\Phi(x)\rho(x)\,dx\approx\hat{\Phi}=\frac{1}{N}\sum_{m=1}^{M}\Phi(X_{m})W_{m}.$
(6)

The SS algorithm will generally help to improve this estimate for functions
$\Phi$ most sensitive to the tail region of $S(x)$, rather than its central
bulk. In particular, setting $\Phi(x)=\mathbb{I}\\{S(x)>\ell\\}$, we recover
the estimator $\hat{p}_{\mathrm{SS}}(\ell)$:

$\displaystyle\mathbb{E}[\mathbb{I}\\{S(X)>\ell\\}]=p(\ell)\approx\frac{1}{N}\sum_{m:S(X_{m})>\ell}W_{m}=\hat{p}_{\mathrm{SS}}(\ell).$
(7)

An important set of algorithmic choices are the population parameters
$N,N_{1},...,N_{J}$, the killing numbers
$\kappa_{1},\kappa_{2},...,\kappa_{J}$, as well as the halting criterion which
determines $J$. Cerou2019adaptive reviews theoretical bases for several
different choices, but here for simplicity we opt for the same rule as used in
Lestang2018computing: $\kappa_{j}=\kappa=1$ (the “drop 1” rule) and $N_{j}=N$
for all $j=1,...,J$ (the population is replenished after each new level is
set). Note that with $\kappa_{j}=1$, only a single parent is selected from
$\mathbb{Q}$ at each round before the level is raised and the queue re-
initialized.

###  3.1 Adaptive multilevel splitting (AMS)

AMS (in particular “trajectory AMS (TAMS)” in the nomenclature of
Lestang2018computing) can be seen as a special case of SS where each
$X=\\{X(t):0\leq t\leq T\\}$ is a length-$T$ trajectory of a stochastic
dynamical system, $S(X)=\max_{0\leq t<T}R(X(t))$ for some time-dependent score
function $R$, and with a particular choice for splitting trajectories.
Trajectories are split by constructing a new forcing sequence
$\widetilde{\eta}(t)$ ($\widetilde{\eta}_{1,2}(t)$ for our L96 model) to drive
the child trajectory $\widetilde{X}(t)$ starting from the old forcing sequence
$\eta(t)$ that drove the parent. First, copy the initial condition
$\widetilde{X}(0)=X(0)$. Then, copy $\widetilde{\eta}(t)=\eta(t)$ up until
some _split time_ $t_{\mathrm{sp}}$, which is chosen as first time
$t_{0}(\ell)$ that the parent clears the threshold:

$t_{\mathrm{sp}}=t_{0}(\ell_{1})=\min\\{t\in[0,T]:R(X(t))>\ell_{1}\\}.$ (8)

For following times $t\geq t_{\mathrm{sp}}$, swap in a new and independent
noise forcing sequence for $\widetilde{\eta}(t)$. No Metropolis-style
accept/reject step is needed for step (2) above; each newly sampled Brownian
increment of $\widetilde{\eta}(t)$ is drawn independently from
$\mathcal{N}(0,\Delta t)$, and so $\widetilde{\eta}(t)$ is a proper sample
from the same noise-generating distribution as $\eta(t)$. Furthermore, the
choice of $t_{\mathrm{sp}}=t_{0}(\ell_{1})$ guarantees $\widetilde{X}(t)=X(t)$
for all $t\leq t_{0}(\ell_{1})$, so that $S(\widetilde{X})>\ell_{1}$, and
acceptance is guaranteed in step (3) as well.

The change in forcing for $t\geq t_{\mathrm{sp}}$ will cause the child to
diverge from the parent, producing a new—but correlated—sample (Fig. 1a). How
correlated $\widetilde{X}$ is to its parent $X$ depends on $t_{\mathrm{sp}}$,
with later $t_{\mathrm{sp}}$ implying a longer shared history and less
independence. Applying the split at $t_{\mathrm{sp}}=t_{0}(\ell)$ maximizes
the independence of the child—and ultimately the diversity of the AMS
ensemble—while guaranteeing $S(\widetilde{X})$ exceeds $\ell_{1}$, and
therefore is accepted in the modified Metropolis Algorithm. The same procedure
is carried out for every subsequent level.

We performed a sequence of AMS experiments with the following parameters:

  1. 1.

Physical constants and timescales: $F_{4}\in\\{3,1,0.5,0.25\\}$ for the
default case $a=1$ which gives the stochastically forced L96 model, and
$F_{4}=3$ for the case $a=0$ which gives the OU process (shown in
supplementary Figs. S1 and S2). We fix $F_{0}=6$, and $K=40$ throughout, and
set the time horizon to $T=6$.

  2. 2.

Ensemble sizes and population control: $N=N_{j}=128$ and $\kappa_{j}=1$ for
$j=1,2,...,J=896$ adhering to a fixed computational budget of 1024 time
horizons simulated. One additional halting criterion is imposed: if the
population loses so much diversity that all active ensemble members descend
from the same ancestor, we terminate the algorithm early.

  3. 3.

We repeat the whole procedure $M=56$ times for each parameter set, with
different seeds for pseudo-random number generation. Each repetition will be
called a “run” of AMS. Having multiple runs allows us to assess variance, and
by using pooled estimates from all runs to hedge against stagnation within
local optima of phase space in a particular run.

The initial $N$-member ensemble is generated as a sequence of consecutive
blocks from a moderate initialization simulation of length $N\times T$ ($T=6$
is the time horizon), after discarding the first 50 units as spinup. The
spinup is initialized as $x_{k}(0)=F_{0}+\frac{1}{1000}\sin\big{(}\frac{2\pi
k}{K})$. The random number generator used to create the noise forcing
sequences $\eta_{1,2}(t)$ is seeded with $s\in\\{0,...,M-1\\}$, a different
value for each AMS run with a fixed parameter set. The $N$ initial blocks,
although weakly correlated, comprise a sample from the steady-state
distribution of the stochastic L96 system. Larger $N$ reduces the variability
of the AMS results, but it also means more up-front cost and more rounds of
splitting needed to reach return times long enough to make the algorithm
worthwhile.

We compare our results from AMS to a long DNS simulation of length $1.28\times
10^{6}$ (separate from the initialization), which is then further elongated by
a factor of 40 (concatenating all $K$ timeseries end-to-end) into $5.12\times
10^{7}$, exploiting the statistical equivalence of all $K=40$ sites of L96.
This curve is our best estimate of ground truth. Note that the symmetry is
only exploited to extend the DNS estimate, not the AMS estimate. In a climate
model with zonal inhomogeneities, such as continents, it would be
inappropriate to aggregate different longitudes together.

Figure 4:  Scores for single ancestors and their descendents within the AMS
algorithm (special case of TEAMS with $\delta=0$). For each stochastic forcing
amplitude, 56 independent runs of AMS were carried out (indexed 0-55) with
$N=128$ ensemble members (0-127). (a) Time-dependent score function $R(X(t))$
for the 7th initial ensemble member (ancestor) of run 14 for $F_{4}=3$. A
black circle indicates the scalar score $S(X)=\max_{t}R(X(t))$. $R(X(t)$ and
$S(X)$ are also shown for a single lineage (path down the family tree) in a
sequence of brightening colors, ending with the highest scoring descendant’s
score in red. (b) Scores in gray dots, with the horizontal axis numbering all
descendants from ancestor 7 of run 14 for $F_{4}=3$. Colored circles indicate
those descendants in the lineage from (a). The dashed gray curve indicates the
levels $\ell$ from which each descendant was split. (c,e,g) are the same as
(a), and (d,f,h) are the same as (b), but with stochastic forcing strength
decreasing to $F_{4}=1,0.5$, and 0.25 respectively. In each case, the run and
ancestor were hand-selected among the ancestors with the maximum boosting.
Figure 5: Performance of the AMS algorithm (special case of TEAMS with
$\delta=0$). (a) Return level vs. return period plots for $F_{4}=3$. Blue
lines show estimates from the initial 128 members of each AMS run; red lines
show estimates from the completed AMS runs; black line shows DNS. (b) Return
level vs. return period for a pooled AMS ensemble containing all $56\times
1024$ members. Blue and red envelopes indicate 95% confidence intervals (see
text for details). Gray envelope is a 95% confidence interval based on subsets
of DNS equal in total cost to the 56 AMS runs. Thus, the dashed red line and
shading from AMS is of equal cost to the gray shading from DNS. (c) Unweighted
histogram of scores for AMS initialization (blue), completed AMS (red), and
DNS (black). Following rows are same as first row, but with noise decreasing
to $F_{4}=1,0.5$, and 0.25, respectively. The slight variability in TEAMS
costs listed to the left are due to the early halting criterion of one single
ancestor remaining (see section 3).

Fig. 4a,b illustrates the effect of successive mutations over the course of
the AMS algorithm, on the relatively easy test case with strong stochastic
forcing, $F_{4}=3$ and $a=1$ (the even easier case of $a=0$—the OU process
with no interference from advection—is documented in Lestang2018computing and
included in supplementary Figs. S1 and S2 for completeness). By design, the
levels increase monotonically over the course of generations and the
descendant scores march upward, ultimately mutating the moderate ancestor into
an extreme descendant. Going beyond this successful “anecdote”, Fig. 5(a,b,c)
confirm the benefit of AMS for a _statistically accurate_ sampling of the
distribution’s tails. Fig. 5a shows return period curves calculated with the
modified block maximum method according to three datasets: the full weighted
ensemble from AMS; the initialization (“Init”), consisting of $N$ ensemble
members per AMS run; and the long DNS simulation. The return levels are
interpolated onto a common logarithmically spaced grid of return periods for
easy comparison between the three data sources. Whereas return level estimates
based on the initializations alone (blue) scatter considerably around the
ground truth, AMS provides a tighter range of estimates (red) around the
ground truth, and for $\sim 3$ orders of magnitude-longer return periods, at
only 8 times the cost of initialization (1024 members from an initial 128).
Moreover, each AMS run is $\sim 5000$ times less costly than the DNS run that
gave the ground truth curve; altogether, the 56 AMS runs are $\sim 100$ times
less costly.

Another way of comparing AMS to DNS is by pooling together all members from
the 56 ensembles and considering them as one larger ensemble of size $56\times
1024=57344$. Fig. 5b shows the resulting statistics which have the advantage
of extending to considerably longer return periods than the individual AMS
runs. Here, as in Fig. 3, the error bars are given by the basic bootstrap 95%
confidence interval using 5000 bootstrap samples, but in the case of DNS (gray
error bar), each bootstrap resampling contains only enough blocks to match the
total simulation time used by AMS (including all independent runs). This lets
us compare the uncertainties fairly between the two methods. In the case of
AMS error bars, the members within a single run are not independent of each
other, and so we resample the AMS runs. That is, we sample the numbers
$\\{0,...,55\\}$ 5000 times with replacement, and for each resampling we pool
together all members from the corresponding list of AMS runs, including
repetitions. Fig. 5c shows the unweighted histogram of scores coming from the
three data sources. The difference in shape of the AMS histogram compared to
the DNS histogram demonstrates the main effect of AMS: to undersample the low
end of the distribution and oversample the tail, shifting the computational
burden to where it is more useful for sampling extremes.

We consider AMS to “win” over DNS if either of two criteria are met: (i) the
AMS estimate remains close to the DNS (relative to error bar width) for return
periods well beyond the AMS total simulation time $T_{\mathrm{AMS}}$; (ii) the
AMS error bar is much smaller than the DNS error bar at $T_{\mathrm{AMS}}$.
Under strong stochastic forcing, AMS performs very well by both criteria,
accurately (and confidently) estimating return periods as long as $10^{7}$ in
the pooled estimate using only $3.4\times 10^{5}$ time units of computation.
This aligns with the demonstration in Lestang2018computing for the OU process,
and serves as a departure point for our modification of the algorithm.

###  3.2 Failure of AMS in the regime of weak stochastic forcing

The story gets more complicated when the stochastic forcing is weak and
nonlinear dynamics dominate. In deterministic chaos, perturbations grow
exponentially with a rate inversely proportional to the _Lyapunov timescale_
—at least, so long as the perturbations remain infinitesimal. Only after
several elapsed Lyapunov times—what we call the _divergence timescale_ ,
quantified further in section 4—do perturbations become large enough to be
useful for splitting algorithms, but also at which size nonlinear effects take
over. In contrast to deterministic chaos, white noise realizations diverge
immediately. The stochastic L96 system inherits both behaviors to some extent,
determined by the relative strength of stochastic forcing. Our main thesis is
that when nonlinear dynamics dominate, and divergence time exceeds the
duration of the event of interest, standard AMS is inadequate, but this can be
remedied by adjusting the choice of splitting time $t_{\mathrm{sp}}$ as shown
in the next section.

Fig. 4c-h show ancestors and descendents for AMS, analogous to Fig. 4a,b and
with identical algorithmic parameters, but with decreasing levels of
stochastic forcing: $F_{4}=1,0.5,0.25$. For all four stochastic forcing
strengths, ancestors can spawn more extreme descendants. However, there is a
key difference between the strong- and weak-stochastic forcing regimes. With
strong stochastic forcing $F_{4}=3$ (Fig. 4a,b), each descendant along the
lineage improves upon the _same event_. In other words, the sequence of
maximum scores comes from a peak in the timeseries for $R(X(t))$ that grows
taller and taller, drifting only slightly forward in time. With weaker
stochastic forcing (Fig. 4 c-d, e-f and especially g-h), events tend to see
only modest boosts from generation to generation. The only way for a child
$\widetilde{X}$ to improve _substantially_ over its parent $X$ is by creating
a whole new event—a new peak later in the time horizon—rather than building on
an existing event. This happens because the stochastic forcing is too weak to
open a large gap between $R(\widetilde{X}(t))$ and $R(X(t))$ during the short
interval between the splitting time $t_{0}(\ell)$, when $R(X(t))$ first
exceeds $\ell$, and the peak $\text{argmax}_{t}R(X(t))$. The child ends up
essentially replicating the parent’s peak, which is the same behavior
illustrated schematically in Fig. 1a. The characteristic time scale of the
peak (what we will call the event duration) is set by the zonal propagation of
waves, and this timescale is not long enough compared to the divergence time
for AMS to work well. The same phenomenon was observed in
Lestang2020numerical): extreme spikes in the force on a body in a turbulent
channel flow (see their Fig. 14) could not be boosted via AMS, which was
attributed to the “sweeping” of vortices past the body. Similar reasoning
holds for the zonal propagation of waves in L96 and the passage of midlatitude
cyclones or fronts past a location in the midlatitudes.

Fig. 5 summarizes the performance of AMS for different strengths of stochastic
forcing. The suspicion of failure raised by Fig. 4 is confirmed by the clear
degradation of performance as $F_{4}$ shrinks. In particular, the individual
AMS return level curves tend to fall farther and farther underneath the true
return level curves (left column of Fig. 5). There is a large scatter in the
individual runs, and in the case $F_{4}=0.5$, a lucky few of the 56 runs
salvage the pooled estimate for a decent approximation of the DNS return
levels, but the width and asymmetry of the confidence intervals indicate the
unreliability of this result (Fig. 5h). The problem becomes particularly acute
as $F_{4}$ drops to 0.25, with the individual AMS runs barely improving upon
the initial scores (Fig. 5j) and a large underestimate at longer return
periods for the pooled estimate (Fig. 5k).

It thus appears that standard AMS is dead on arrival for cases where the
divergence timescale is longer than the event duration. In principle, there is
a canonical fix for this problem, namely to use a more intelligent score
function than the quantity of interest $R(X(t))$ itself. The ideal such proxy
is the _committor_ : the probability, given an initial condition $X(t)=x$,
that $R(X(s))$ will exceed $\ell$ at some time $s\in(t,T)$ before the time
horizon ends. By definition, the committor incorporates information about the
model state $X(t)$ that is not available from $R(X(t))=x_{0}^{2}$, for example
the speeds and magnitudes of different wave packets scattered across the
domain that may all soon converge at $k=0$ and result in an extreme burst of
energy. The committor is an _optimal_ score function for AMS in terms of
minimizing the variance for $\hat{p}(\ell)$ [Lestang . (2018), Cérou . (2019),
Lucente, Rolland . (2022)]. Considerable research has recently pursued
approximation strategies for the committor in various climate applications
[<]e.g.,¿[]Tantet2015early,Finkel2021learning,Lucente2022committor,Miloshevich2022probabilistic,JacquesDumas2023data.

Unfortunately, these strategies all require either a high volume of training
data—potentially canceling out the savings of a rare event algorithm, which is
useful precisely in the low-data regime—or very specific knowledge of phase
space geometry, such as a bistable structure, which is not typically available
for realistic climate models. A second, related problem is that the optimality
property only holds true for a single committor with a fixed threshold $\ell$.
What if we seek return periods for a whole range of thresholds? We would have
to sacrifice the accuracy of some return periods in favor of others.
Alternatively, we could use the committor for a single very high threshold
$\ell_{\mathrm{max}}$, but then even less training data would be available.
Although it is interesting and worthwhile to search for committor functions
based on traveling-wave dynamics, we leave that to future work, and in the
next section we describe a simpler strategy to get around the stagnation issue
seen in Fig. 4.

###  3.3 Trying-early adaptive multilevel splitting (TEAMS)

To address the failure of AMS in the nonlinear regime, we adjust
$t_{\mathrm{sp}}=t_{\delta}(\ell)=:t_{0}(\ell)-\delta$ by an _advance split
time_ $\delta>0$, allowing some time for the child $\widetilde{X}$ to drift
farther away from the parent and possibly achieve a higher maximum score.
Indeed, ensemble boosting [Gessner . (2021)] does exactly that, systematically
applying perturbations every day from 19 to 7 days in advance of heat wave
onset, although ensemble boosting does not by itself allow the calculation of
return periods for the boosted events. When splitting early we lose the
guarantee that $R(\widetilde{X}(t))$ clears the current level $\ell$
(decpicted schematically in Fig.1b), which is why we frame our modified
algorithm using subset simulation (see section 3) which includes an
accept/reject step: when a child fails to score higher than $\ell$, it is
discarded from the ensemble and its parent is duplicated instead (in other
words, doubling its statistical weight). The resulting algorithm, which we
call TEAMS (“trying-early adaptive multilevel splitting”), incurs additional
cost due to rejected samples, but also gains back the ability to build
significantly upon ancestral scores. One can interpret $\delta$ as setting the
width of the proposal distribution, a key parameter in Markov chain Monte
Carlo methods. A wider proposal allows the child to explore farther afield
from its parent, but increases the risk of rejection. Proposal width often has
to be tuned carefully, and the sampling community has devoted substantial
efforts to adaptively designing the proposal [Walter R. Gilks Sahu (1998),
Andrieu Thoms (2008)]. Such methods will surely prove useful for complex
climate models, but in our present proof-of-concept study of the algorithm, we
found approximately optimal $\delta$ values by exhaustive grid search for each
noise level. Section 4 explains this procedure and shows that the optimal
$\delta$ can be related to the error saturation timescale, a classical measure
of predictability.

We performed a sequence of TEAMS experiments with
$(F_{4},\delta)\in\\{3,1,0.5,0.25\\}\times\\{0,0.2,0.4,...,2.0\\}$. We adjust
the time horizon $T=6+\delta$ to give each parameter choice the same length of
score to boost. All other parameters are as before for the AMS experiments.

Fig. 6 shows TEAMS in action for the same parameter sets from Fig. 4, but with
(roughly optimal) advance splitting times $\delta=0.0,0.6,1.0$, and $1.4$ for
the decreasing noise levels (at $F_{4}=3$, $\delta=0$ still works best, and
panel (a) is the same as in Fig. 4a)). Note that the score functions $R(X(t))$
are only defined for times $t>\delta$, because if $t_{0}(\ell)<\delta$ then
$t_{\delta}(\ell)<0$, so we cannot apply the split early enough. This is
implemented by setting the early scores to NaN, and lengthening the time
horizon from $T$ to $T+\delta$ as mentioned above. We account for this extra
cost in all the performance calculations to follow, but we omit the first
$\delta$ time units from the plots. For all four stochastic forcing strengths,
we see examples of children building significantly, and directly, upon a
parent’s maximum, without having to discover a new peak farther into the
future. The values of the scores form continuous point clouds in panels
(b,d,f,h), unlike the discrete horizontal bands appearing in Fig. 4(f,h) where
$\delta=0$ and stochastic forcing is weak. The negative side-effect is that
many gray dots fall short of the gray dashed line, indicating a rejected
sample. Clearly, increasing $\delta$ brings both higher risk and higher
reward.

Figure 6: Scores for single ancestors and their descendants generated by the
TEAMS algorithm: the same as Fig. 4 but with advance split times $\delta$
chosen to be approximately optimal for each noise level: $\delta=$ 0, 0.6, 1,
and 1.4 for $F_{4}=3,1,0.5$, and 0.25, respectively. Because $\delta=0$ is
optimal for $F_{4}=3$, (a,b) is the same as Fig. 4a,b. Section 4 explains how
the $\delta$ values were chosen.

Fig. 7 quantitatively confirms the hopeful suggestion of Fig. 6: that
increasing $\delta$ can give TEAMS a speedup over DNS in the weak stochastic
forcing regime. For all cases shown, TEAMS extends the estimated return
period, _accurately_ , well beyond the gray envelope which marks the limit
achievable by an equal-cost run of DNS. The black ground truth curve remains
within the 95% confidence band of TEAMS to return periods of $\sim 10^{7}$
across all forcing levels. Simultaneously, the TEAMS confidence band is
narrower than the DNS band.

Fig. 7 shows TEAMS gives a good estimate of the return values when all runs
are pooled together, but that most individual TEAMS runs underestimate the
true return values while a few overestimate them to allow for a good pooled
estimate. As in Lucente2022coupling, we can attribute this behavior to
_apparent bias_ , which is best explained by analogy: an experiment consisting
of 100 flips of a coin with $p=\mathbb{P}(\text{heads})=0.001$ has a nine in
ten chance of landing no heads, yielding a probability estimate $\hat{p}=0$.
But one experiment out of ten will yield $\hat{p}=0.01$, a gross over-
estimate, and only by pooling these two scenarios together can we see the
estimator’s lack of bias. Unlike the coin-flipping experiment, TEAMS is
designed to preferentially sample extreme values, but a given AMS run for L96
may still get stuck in a local optimum yielding underestimated return values,
especially if the stochastic forcing is too weak to jolt a trajectory out of
it. Thus, pooling over multiple runs is especially crucial in the
deterministic limit.

Figure 7:  Performance of the TEAMS algorithm: the same as Fig. 5 but with
advance split times $\delta$ chosen to be approximately optimal for each noise
level: $\delta=$ 0, 0.6, 1, and 1.4 for $F_{4}=3,1,0.5$, and 0.25,
respectively. Because $\delta=0$ is optimal for $F_{4}=3$, (a-c) are the same
as Fig. 5a-c.

##  4 Optimizing advance split time

In this section, we explain how we determined optimal values of the advance
split time $\delta$ using a simple exhaustive search. We then investigate the
behavior of $\delta$ as a function of stochastic forcing strength as a guide
for choosing $\delta$ prior to running TEAMS on a more expensive model for
which exhaustive search would not be feasible.

###  4.1 Exhaustive search

We selected the “optimal” $\delta$ values based on two simple performance
metrics, which are plotted in Fig. 8.

  1. 1.

Return level RMSE: the root-mean-square difference of return level between a
TEAMS estimate (from a single run) and the DNS-determined ground truth, where
the mean is taken over uniform bins in $\log\tau$ space. This metric is
proportional to the $L^{2}$-norm between a red line and the black line in the
left columns of Figs. 5 and 7. In cases where the red line stops before the
black line, it is extrapolated to longer return periods with a constant given
by its maximum to penalize the algorithm getting stuck at a false upper bound.
We calculate statistics of the return level RMSE across runs, including the
mean and quantiles, which are displayed in Fig. 8(a,c,e,g). Note that these
correspond to _percentile bootstrap_ confidence intervals [Wasserman (2004)],
as opposed to the _basic bootstrap_ confidence intervals shown in Figs. 5 and
7. Here we use the percentile bootstrap as a means of sensitivity analysis, to
show the range of results that might occur due to sampling variability. The
basic bootstrap, by contrast, is intended to bracket the ground truth of some
parameter value. The return level RMSE can also be calculated for the pooled
estimate, and it shows similar but noisier trends.

  2. 2.

Mean family gain: the maximum improvement (difference in scores) from ancestor
to descendant over all $N$ ancestors, averaged over the 56 runs. This does not
measure statistical accuracy, but only the consistent ability to generate
extreme events out of moderate events. Fig. 8 (b,d,f,h) shows mean family
gain. Other metrics of gain, such as the maximum descendant score minus the
maximum ancestral score (not necessarily from the same family tree) yield very
similar trends with $\delta$, albeit different absolute values.

A good choice of $\delta$ should have a small return level RMSE and a large
mean family gain. Based on both performance metrics, we selected optimal
$\delta=0,0.6,1,1.4$ for $F_{4}=3,1,0.5,0.25$, respectively. These optimal
values are marked with vertical gray lines in Fig. 8, and they are used in
Figs. 6 and 7. For $F_{4}=0.5$, the two metrics gave slightly difference
optimal values ( $\delta=1.2$ for return level RMSE or $\delta=1$ for mean
family gain); we chose $\delta=1$ because it gave the better pooled estimate.
We emphasize that the optimal values are only discernible by averaging over
many independent runs. For completeness, we display all 44 return level vs.
return period plots (4 values of $F_{4}$ $\times$ 11 values of $\delta$) in
the supplement. In general, shifting the optimal $\delta$ by $\pm 0.2$ doesn’t
change the results qualitatively, but larger shifts can degrade performance.
The absolute values of errors should not be compared between stochastic
forcing levels, since each has its own statistical steady state. Rather, the
important takeaway is the increase in optimal $\delta$ as the stochastic
forcing weakens. Indeed, in the singular limit of zero stochastic forcing the
advance split time must necessarily go to infinity to have any effect at all,
and initial condition perturbations would be needed to split trajectories.

Figure 8: Performance of TEAMS as a function of advance split time $\delta$
and as measured by (a,c,e,g) return level RMSE and (b,d,f,h) mean family gain
for $F_{4}=$ (a,b) 3, (c,d) 1, (e,f) 0.5, and (g,h) 0.25. Return level RMSE is
computed separately for each run. Thick red lines show the average over runs,
and red envelopes show the quantile ranges 25%-75% (or interquartile range,
IQR) and 2.5%-97.5% across the 56 runs. Mean family gain is maximum gain in
score within a single family averaged over the 56 runs. Vertical gray lines
show the optimal values of $\delta$ used in Figs. 6 and 7.

To summarize, we have found that some choices of $\delta$ can make TEAMS
effective where AMS is not effective, and that the optimal $\delta$ increases
as stochastic forcing magnitude decreases. In the next section we relate this
behavior to the predictability time, which points toward a cheap method to
estimate an optimal—or at least, reasonably performant—$\delta$, without
having to repeatedly run TEAMS.

###  4.2 Relation between optimal advance time and error saturation timescales

Heuristically, we expect the optimal advance time $\delta$ to reflect the
divergence timescale of perturbed trajectories that are introduced in
splitting. Can this be related to classical predictability timescales?
Lyapunov analysis describes perturbation growth by way of Lyapunov exponents
and singular vectors [Cencini Ginelli (2013), Norwood . (2013), Pazo . (2010),
Maiocchi . (2024)], but it applies to the regime of _infinitesimal_
perturbations. The kind of perturbations we strive for in rare event sampling
are finite and nonlinear, turning peaks into substantially larger peaks as in
Figs. 4, 6. “Finite size Lyapunov exponents” (FSLEs) [Boffetta . (1998),
Cencini Vulpiani (2013)] are closer to what we need, generalizing the Lyapunov
exponent to be dependent on an initial error amplitude. Typically, error grows
in two stages: first exponentially, during which the FSLE equals the leading
Lyapunov exponent, and then diffusively (scaling as a power law with time),
during which the FSLE declines. The divergence timescale is bounded below by
this change point, which approaches zero as stochastic forcing becomes
dominant: indeed, the variance of pure Brownian motion grows linearly in $t$
immediately.

On the other hand, the optimal $\delta$ is bounded above by the _error
saturation timescale_ , when perturbed ensemble members become independent and
inhabit totally different regions of the attractor. By then, the root-mean-
square error (RMSE) of the ensemble will equal the root-mean-square distance
(RMSD) between two randomly chosen points on the attractor. In climate models,
the saturation timescale is a convenient and effective measure of
predictability [Sheshadri . (2021)]. Clearly, $\delta$ must be chosen shorter
than the time to saturation, since a child trajectory ought to take advantage
of pre-existing maxima produced by its parent. To investigate this
relationship, the following experiments measure time in terms of fraction of
saturation.

For each $F_{4}$ considered, we ran a moderate-length control simulation
$x(t)$ for $0\leq t\leq 1050$ (discarding the first 50 as spinup), and
measured the RMSD for this simulation. At initialization times
$50,70,90,...,990$ (48 in total) we branched a 16-member ensemble with
identical initial conditions $x(t)$ but independent stochastic forcing
realizations (a convenient feature of stochastic forcing is that errors grow
even from perfect initial conditions, removing dependence on initial
perturbation amplitude). We integrated each member for 15 time units,
calculated RMSE as a function of time (separately for each ensemble), and
inverted to find the times $t_{\epsilon}$ at which the fraction of saturation
$\epsilon=$RMSE/RMSD reached a given value. In other words,
$\mathrm{RMSE}(t_{\epsilon})=\epsilon\times\mathrm{RMSD}$. Finally, we take
the average across initializations to get a single value
$\overline{t_{\epsilon}}$ for each of several $\epsilon$ values. The total
cost of this experiment is $1.2\times 10^{4}$ time units, roughly equal to 1.5
runs of AMS and much cheaper than the 56-run pooled estimate. Moreover,
halving the number of initializations used yields qualitatively similar
results.

Fig. 9 shows timeseries of $x_{0}(t)$ (both control and perturbed) and error
growth for two such ensembles from the high and low stochastic forcing cases.
The time axis is truncated to 10 days past initialization. The early linear
growth of $\epsilon$ vs. $\overline{t_{\epsilon}}$ indicates a steady decline
in relative growth rate, meaning that the perturbations begin to enter the
diffusive (sub-exponential) growth regime quite early. This is thanks to
stochastic forcing, which is visible in the top row as the emergence of red
members from the shadow of the control trajectory. As expected, the error
growth is much faster for the higher value of stochastic forcing.

Figure 9: Growth of perturbations in the experiments described in subsection
4.2 for one representative initialization time $t_{0}=70$ and two values of
the stochastic forcing: (a,c) $F_{4}=3$ and (b,d) $F_{4}=0.25$ . (a,b) show
$x_{0}(t)$ for the control simulation (black) and 16 simulations with the same
initial condition but different white-noise forcing realizations (red). (c,d)
show Euclidean distance between each ensemble member to the control as a
fraction of RMSD versus time (red), and the fraction of saturation RMSE/RMSD
versus the time until each $\epsilon$ value is reached averaged across all
initializations and ensemble members (black), i.e., $\epsilon$ vs.
$\overline{t_{\epsilon}}$. Dots indicate $\epsilon=1/32,1/16,1/8,1/4,3/8,1/2$,
and these same values reflected about $1/2$.

If the optimal $\delta$ could be predicted from the error growth rates alone,
the TEAMS algorithm could be calibrated simply and cheaply before being
deployed. Fig. 10 shows the time $\overline{t_{3/8}}$ when RMSE reaches a
fixed fraction of RMSD ($3/8$) as compared to the optimal $\delta$ values
determined from Fig. 8, as a function of the strength of stochastic forcing.
We include results from forcing at wavenumbers $m=1,4,7,10$. There is an
encouraging similarity between the dependence of optimal $\delta$ and
$\overline{t_{3/8}}$ on stochastic forcing strength, suggesting that the
fractional saturation time might be useful to provide an estimate for
$\delta$.

Another interesting and less obvious feature is the dependence on wavenumber
of error growth (albeit a weak dependence): medium-length wave forcing ($m=4$
and $m=7$) drives error to saturation faster than very short ($m=10$) or long
($m=1$) wave forcing, which informed our choice of $m=4$ throughout the TEAMS
experiments. However, the variability due to initial conditions (indicated by
$\pm 1\sigma$ error bars) tend to exceed systematic differences between
wavenumbers. This variability reflects a distribution of divergence timescales
across the attractor, which was also found be be quite heterogeneous in
Maiocchi2024heterogeneity (there measured by Lyapunov exponents). It also
suggests that the best strategy may be to not fix a single $\delta$, but to
allow the algorithm to adaptively set a $\delta$, or sample from a range, to
account for differing divergence timescales between different initial
conditions, and this could be investigated in future work.

Figure 10: Time $\overline{t_{3/8}}$ until the perturbations described in
subsection 4.2 reach a fixed fraction (3/8) of RMSD as a function of
stochastic forcing strength $F_{m}$ for different wavenumbers $m$. Error bars
are $\pm 1$ standard deviation of the distribution over different initial
conditions. Optimized values of $\delta$ (determined from the performance
metrics in Fig. 8) are shown in the black dashed line for $m=4$.

##  5 Conclusions and Outlook

A vexing challenge in climate science is reliably quantifying the probability
of extreme weather events, which are fundamentally difficult to characterize
because of data scarcity. Among various competing strategies, rare event
algorithms hold several key advantages, chiefly (i) access to dynamical
samples of the events, rather than just return period curves which extreme
value theory might provide, and (ii) more statistical rigor than storyline-
based approaches like “ensemble boosting” [Gessner . (2021)], thanks to
careful re-weighting of cloned trajectories. Inspired by recent successes of
rare event algorithms on long-lasting heat waves [Ragone . (2018)] and
idealized models of regime transitions [Lucente, Rolland . (2022), Jacques-
Dumas . (2023)], we have investigated the ability of a particular algorithm,
adaptive multilevel splitting (AMS) to sample extreme events of a different
character: intermittent, short-lived bursts of energy in the Lorenz-96 model
which have some similar characteristics as extreme daily rain or wind
associated with midlatitude cyclones.

Even in this simple model, we have elucidated some key obstacles that hinder
rare event splitting algorithms on sudden, short-lived events, and furthermore
taken some steps to overcome them. AMS sets up a sequence of thresholds for an
observable of interest and estimates conditional exceedance probabilities in
stages by cloning and perturbing “successful” ensemble members when they cross
a threshold, to generate new “successful” samples. This simple prescription
suffers a fatal problem when the events are short-lived compared to the
divergence timescale (how long it takes a perturbation to grow appreciably): a
perturbed ensemble member essentially replicates its parent’s success, and
doesn’t develop its own history until after the event is over. Neither the
magnitude nor the diversity of rare event samples is enhanced. To fix this
problem, we have drawn inspiration from ensemble boosting to apply a
perturbation _in advance_ of the rare event by some lead time $\delta$. But we
have also retained rigorous statistics for these “storylines” by exploiting a
more general rare event algorithm, subset simulation (SS), of which AMS is
only a special case. We name the resulting algorithm “trying-early AMS”
(TEAMS) and demonstrate its success in sampling the tails of the rare event
distribution more efficiently than direct numerical simulation can do, despite
an extra computational cost due to rejected samples.

Our study is a proof of concept that suggests a path forward, but with some
open questions and directions for improvement, which we summarize here:

  * •

The most crucial algorithmic parameter is the advance split time, $\delta$,
which is equivalent to a proposal distribution width. Our grid search for
optimal $\delta$, though not a scalable solution, demonstrates a relationship
with the time over which perturbations grow to a fraction of saturation. An
important goal for future work is to assess this result for other underlying
models such as general circulation models or for other error growth metrics.
Given the localized nature of our observable ($x_{0}^{2}$ is the energy at a
single longitude site), it is interesting that a _global_ Euclidean metric
correlates with the optimal $\delta$. Weighting the metric more heavily for
grid points near $k=0$ might further improve this relationship.

  * •

The weak stochastic forcing limit $F_{m}\to 0$ is important to confront for
climate models, which may be more practical to perturb just at the splitting
time rather than continuously at every time step, especially if the climate
model is not already equipped with a stochastic subgrid parameterization. In
the TEAMS framework, this would translate to perturbing a simulation at a lead
time $\delta$ ahead of the event, but not at all following times. Perturbing
at just one time makes a given perturbation magnitude less powerful—but also
opens up interesting possibilities such as the use of deterministic
optimization strategies to more efficiently discover the most extreme event
possible from a given initial condition. For example, some directions of
perturbation (singular vectors) grow much faster than others, a fact which has
informed ensemble design in operational weather forecasting [Palmer Zanna
(2013)], and could be used to further improve the algorithm. Methods such as
conditional nonlinear optimal perturbation [<]¿[and references
therein]Wang2020useful, generalized stability theory [Farrell Ioannou (1996)],
and large deviation theory [Dematteis . (2018), Dematteis . (2019), Schorlepp
. (2023)] may prove useful for this task.

  * •

Related to the previous point, it is desirable to have greater efficiency with
samples in order to deploy rare event algorithms at scale. For example, we
should not simply discard rejected samples, but rather learn from their
“mistakes” to design better perturbations. Frameworks like Bayesian
optimization and adaptive importance sampling based on model reduction have
been developed for this task, and have been used in safety assessment for
reliability engineering
[<]e.g.,¿Cousins2014quantification,Huang2016assessing,Mohamad2018sequential,Sapsis2020output,Uribe2021cross,Zhang2022koopman.

Rare event algorithms represent a new way to allocate computational resources
to where they matter most. To realize their considerable potential for
efficiency gains, we have taken one of the necessary steps to make them
flexible enough to target intermittent, localized, transient events that
characterize phenomena such as heavy precipitation in complex global climate
models. The Lorenz-96 model is an invaluable prototype as a cheap system that
poses similar algorithmic challenges. Forthcoming papers will use the insight
gained here as a stepping stone to more complex and realistic models.

## Data availability statement

The software to simulate and sample extreme events in Lorenz-96 using TEAMS is
available in a public Zenodo repository at
https://zenodo.org/doi/10.5281/zenodo.10608187. Interested readers are
encouraged to try out the algorithm on other systems of interest, and should
not hesitate to contact J. F. for assistance.

###### Acknowledgements.

This research is part of the MIT Climate Grand Challenge on Weather and
Climate Extremes. It received support by the generosity of Eric and Wendy
Schmidt by recommendation of Schmidt Sciences as part of its Virtual Earth
System Research Institute (VESRI). Computations for this research were carried
out on the MIT Engaging cluster.

## References

  * Abbot . (2021) Abbot2021rareAbbot, DS., Webber, RJ., Hadden, S., Seligman, D. Weare, J.  2021Dec.  Rare Event Sampling Improves Mercury Instability Statistics Rare Event Sampling Improves Mercury Instability Statistics.  The Astrophysical Journal9232236.  https://dx.doi.org/10.3847/1538-4357/ac2fa8 10.3847/1538-4357/ac2fa8
  * Adachi Tomita (2020) Adachi2020methodologyAdachi, SA. Tomita, H.  2020\.  Methodology of the Constraint Condition in Dynamical Downscaling for Regional Climate Evaluation: A Review Methodology of the Constraint Condition in Dynamical Downscaling for Regional Climate Evaluation: A Review.  Journal of Geophysical Research: Atmospheres12511e2019JD032166.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD032166 e2019JD032166 10.1029/2019JD032166  https://doi.org/10.1029/2019JD032166
  * Andrieu Thoms (2008) Andrieu2008tutorialAndrieu, C. Thoms, J.  2008Dec.  A tutorial on adaptive MCMC A tutorial on adaptive MCMC.  Statistics and Computing184343-373.  https://doi.org/10.1007/s11222-008-9110-y 10.1007/s11222-008-9110-y
  * Au Beck (2001) Au2001estimationAu, SK. Beck, JL.  2001\.  Estimation of small failure probabilities in high dimensions by subset simulation Estimation of small failure probabilities in high dimensions by subset simulation.  Probabilistic Engineering Mechanics164263-277.  https://www.sciencedirect.com/science/article/pii/S0266892001000194 https://doi.org/10.1016/S0266-8920(01)00019-4
  * Baars . (2021) Baars2021applicationBaars, S., Castellana, D., Wubs, F. Dijkstra, H.  2021\.  Application of adaptive multilevel splitting to high-dimensional dynamical systems Application of adaptive multilevel splitting to high-dimensional dynamical systems.  Journal of Computational Physics424109876.  https://www.sciencedirect.com/science/article/pii/S0021999120306501 https://doi.org/10.1016/j.jcp.2020.109876
  * Boffetta . (1998) Boffetta1998extensionBoffetta, G., Giuliani, P., Paladin, G. Vulpiani, A.  1998\.  An Extension of the Lyapunov Analysis for the Predictability Problem An Extension of the Lyapunov Analysis for the Predictability Problem.  Journal of the Atmospheric Sciences55233409 - 3416\.  https://journals.ametsoc.org/view/journals/atsc/55/23/1520-0469_1998_055_3409_aeotla_2.0.co_2.xml https://doi.org/10.1175/1520-0469(1998)055¡3409:AEOTLA¿2.0.CO;2
  * Bouchet . (2019) Bouchet2019rare2Bouchet, F., Rolland, J. Simonnet, E.  2019Feb.  Rare Event Algorithm Links Transitions in Turbulent Flows with Activated Nucleations Rare Event Algorithm Links Transitions in Turbulent Flows with Activated Nucleations.  Phys. Rev. Lett.122074502.  https://link.aps.org/doi/10.1103/PhysRevLett.122.074502 10.1103/PhysRevLett.122.074502
  * Bucklew (2004) Bucklew2004introductionBucklew, JA.  2004\.  Introduction to Rare Event Simulation Introduction to Rare Event Simulation (1 ).  Springer New York, NY.  https://doi.org/10.1007/978-1-4757-4078-3
  * Buizza . (1999) Buizza1999stochasticBuizza, R., Milleer, M. Palmer, TN.  1999\.  Stochastic representation of model uncertainties in the ECMWF ensemble prediction system Stochastic representation of model uncertainties in the ECMWF ensemble prediction system.  Quarterly Journal of the Royal Meteorological Society1255602887-2908.  https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.49712556006 https://doi.org/10.1002/qj.49712556006
  * Cencini Ginelli (2013) Cencini2013lyapunovCencini, M. Ginelli, F.  2013Jun.  Lyapunov analysis: from dynamical systems theory to applications Lyapunov analysis: from dynamical systems theory to applications.  Journal of Physics A: Mathematical and Theoretical4625250301.  https://dx.doi.org/10.1088/1751-8113/46/25/250301 10.1088/1751-8113/46/25/250301
  * Cencini Vulpiani (2013) Cencini2013finiteCencini, M. Vulpiani, A.  2013Jun.  Finite size Lyapunov exponent: review on applications Finite size Lyapunov exponent: review on applications.  Journal of Physics A: Mathematical and Theoretical4625254019.  https://dx.doi.org/10.1088/1751-8113/46/25/254019 10.1088/1751-8113/46/25/254019
  * Coles (2001) Coles2001introductionColes, S.  2001\.  An introduction to statistical modeling of extreme values An introduction to statistical modeling of extreme values (1 ).  Springer.  10.1007/978-1-4471-3675-0
  * Cousins Sapsis (2014) Cousins2014quantificationCousins, W. Sapsis, TP.  2014\.  Quantification and prediction of extreme events in a one-dimensional nonlinear dispersive wave model Quantification and prediction of extreme events in a one-dimensional nonlinear dispersive wave model.  Physica D: Nonlinear Phenomena280-28148-58.  https://www.sciencedirect.com/science/article/pii/S016727891400092X https://doi.org/10.1016/j.physd.2014.04.012
  * Cérou Guyader (2007) Cerou2007adaptiveCérou, F. Guyader, A.  2007\.  Adaptive Multilevel Splitting for Rare Event Analysis Adaptive Multilevel Splitting for Rare Event Analysis.  Stochastic Analysis and Applications252417-443.  https://doi.org/10.1080/07362990601139628 10.1080/07362990601139628
  * Cérou . (2019) Cerou2019adaptiveCérou, F., Guyader, A. Rousset, M.  2019\.  Adaptive multilevel splitting: Historical perspective and recent results Adaptive multilevel splitting: Historical perspective and recent results.  Chaos: An Interdisciplinary Journal of Nonlinear Science294043108.  https://doi.org/10.1063/1.5082247 10.1063/1.5082247
  * Dematteis . (2018) Dematteis2017rogueDematteis, G., Grafke, T. Vanden-Eijnden, E.  2018\.  Rogue waves and large deviations in deep sea Rogue waves and large deviations in deep sea.  Proceedings of the National Academy of Sciences1155855-860.  https://www.pnas.org/doi/abs/10.1073/pnas.1710670115 10.1073/pnas.1710670115
  * Dematteis . (2019) Dematteis2019extremeDematteis, G., Grafke, T. Vanden-Eijnden, E.  2019\.  Extreme Event Quantification in Dynamical Systems with Random Components Extreme Event Quantification in Dynamical Systems with Random Components.  SIAM/ASA Journal on Uncertainty Quantification731029-1059.  https://doi.org/10.1137/18M1211003 10.1137/18M1211003
  * Dwyer O’Gorman (2017) Dwyer2017changingDwyer, JG. O’Gorman, PA.  2017\.  Changing duration and spatial extent of midlatitude precipitation extremes across different climates Changing duration and spatial extent of midlatitude precipitation extremes across different climates.  Geophysical Research Letters44115863-5871.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL072855 https://doi.org/10.1002/2017GL072855
  * Emanuel (2021) Emanuel2021responseEmanuel, K.  2021\.  Response of Global Tropical Cyclone Activity to Increasing CO2: Results from Downscaling CMIP6 Models Response of Global Tropical Cyclone Activity to Increasing CO2: Results from Downscaling CMIP6 Models.  Journal of Climate34157 - 70.  https://journals.ametsoc.org/view/journals/clim/34/1/jcliD200367.xml https://doi.org/10.1175/JCLI-D-20-0367.1
  * Farrell Ioannou (1996) Farrell1996generalizedFarrell, BF. Ioannou, PJ.  1996\.  Generalized Stability Theory. Part I: Autonomous Operators Generalized Stability Theory. Part I: Autonomous Operators.  Journal of Atmospheric Sciences53142025 - 2040.  https://journals.ametsoc.org/view/journals/atsc/53/14/1520-0469_1996_053_2025_gstpia_2_0_co_2.xml 10.1175/1520-0469(1996)053¡2025:GSTPIA¿2.0.CO;2
  * Finkel . (2021) Finkel2021learningFinkel, J., Webber, RJ., Gerber, EP., Abbot, DS. Weare, J.  2021\.  Learning Forecasts of Rare Stratospheric Transitions from Short Simulations Learning Forecasts of Rare Stratospheric Transitions from Short Simulations.  Monthly Weather Review149113647 - 3669.  https://journals.ametsoc.org/view/journals/mwre/149/11/MWR-D-21-0024.1.xml 10.1175/MWR-D-21-0024.1
  * Gagne II . (2020) Gagne2020machineGagne II, DJ., Christensen, HM., Subramanian, AC. Monahan, AH.  2020\.  Machine Learning for Stochastic Parameterization: Generative Adversarial Networks in the Lorenz ’96 Model Machine Learning for Stochastic Parameterization: Generative Adversarial Networks in the Lorenz ’96 Model.  Journal of Advances in Modeling Earth Systems123e2019MS001896.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001896 e2019MS001896 10.1029/2019MS001896  https://doi.org/10.1029/2019MS001896
  * Gálfi . (2017) Galfi2017convergenceGálfi, VM., Bódai, T. Lucarini, V.  2017Sep06.  Convergence of Extreme Value Statistics in a Two-Layer Quasi-Geostrophic Atmospheric Model Convergence of Extreme Value Statistics in a Two-Layer Quasi-Geostrophic Atmospheric Model.  Complexity20175340858.  https://doi.org/10.1155/2017/5340858 10.1155/2017/5340858
  * Gálfi . (2021) Galfi2021applicationsGálfi, VM., Lucarini, V., Ragone, F. Wouters, J.  2021Jun01.  Applications of large deviation theory in geophysical fluid dynamics and climate science Applications of large deviation theory in geophysical fluid dynamics and climate science.  La Rivista del Nuovo Cimento446291-363.  https://doi.org/10.1007/s40766-021-00020-z 10.1007/s40766-021-00020-z
  * Gessner (2022) Gessner2022physicalGessner, C.  2022\.  Physical storylines for very rare climate extremes Physical storylines for very rare climate extremes .  ETH Zurich.
  * Gessner . (2021) Gessner2021veryGessner, C., Fischer, EM., Beyerle, U. Knutti, R.  2021\.  Very Rare Heat Extremes: Quantifying and Understanding Using Ensemble Reinitialization Very Rare Heat Extremes: Quantifying and Understanding Using Ensemble Reinitialization.  Journal of Climate34166619 - 6634.  https://journals.ametsoc.org/view/journals/clim/34/16/JCLI-D-20-0916.1.xml 10.1175/JCLI-D-20-0916.1
  * Giardinà . (2006) Giardina2006directGiardinà, C., Kurchan, J. Peliti, L.  2006Mar.  Direct Evaluation of Large-Deviation Functions Direct Evaluation of Large-Deviation Functions.  Phys. Rev. Lett.96120603.  https://link.aps.org/doi/10.1103/PhysRevLett.96.120603 10.1103/PhysRevLett.96.120603
  * Hu . (2019) Hu2019effectsHu, G., Bódai, T. Lucarini, V.  2019\.  Effects of stochastic parametrization on extreme value statistics Effects of stochastic parametrization on extreme value statistics.  Chaos: An Interdisciplinary Journal of Nonlinear Science298083102.  https://doi.org/10.1063/1.5095756 10.1063/1.5095756
  * WK. Huang . (2016) Huang2016estimatingHuang, WK., Stein, ML., McInerney, DJ., Sun, S. Moyer, EJ.  2016\.  Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (GEV) distributions Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (GEV) distributions.  Advances in Statistical Climatology, Meteorology and Oceanography2179–103.  https://ascmo.copernicus.org/articles/2/79/2016/ 10.5194/ascmo-2-79-2016
  * X. Huang . (2016) Huang2016assessingHuang, X., Chen, J. Zhu, H.  2016\.  Assessing small failure probabilities by AK–SS: An active learning method combining Kriging and Subset Simulation Assessing small failure probabilities by AK–SS: An active learning method combining Kriging and Subset Simulation.  Structural Safety5986-95.  https://www.sciencedirect.com/science/article/pii/S0167473016000035 https://doi.org/10.1016/j.strusafe.2015.12.003
  * X. Huang . (2020) Huang2020futureHuang, X., Swain, DL. Hall, AD.  2020\.  Future precipitation increase from very high resolution ensemble downscaling of extreme atmospheric river storms in California Future precipitation increase from very high resolution ensemble downscaling of extreme atmospheric river storms in California.  Science Advances629eaba1323.  https://www.science.org/doi/abs/10.1126/sciadv.aba1323 10.1126/sciadv.aba1323
  * Jacques-Dumas . (2023) JacquesDumas2023dataJacques-Dumas, V., van Westen, RM., Bouchet, F. Dijkstra, HA.  2023\.  Data-driven methods to estimate the committor function in conceptual ocean models Data-driven methods to estimate the committor function in conceptual ocean models.  Nonlinear Processes in Geophysics302195–216.  https://npg.copernicus.org/articles/30/195/2023/ 10.5194/npg-30-195-2023
  * Kahn Harris (1951) Kahn1951estimationKahn, H. Harris, TE.  1951\.  Estimation of particle transmission by random sampling Estimation of particle transmission by random sampling.  National Bureau of Standards applied mathematics series1227–30.
  * Kharin . (2007) Kharin2007changesKharin, VV., Zwiers, FW., Zhang, X. Hegerl, GC.  2007\.  Changes in Temperature and Precipitation Extremes in the IPCC Ensemble of Global Coupled Model Simulations Changes in Temperature and Precipitation Extremes in the IPCC Ensemble of Global Coupled Model Simulations.  Journal of Climate2081419 - 1444.  https://journals.ametsoc.org/view/journals/clim/20/8/jcli4066.1.xml https://doi.org/10.1175/JCLI4066.1
  * Krouma . (2022) Krouma2022assessmentKrouma, M., Yiou, P., Déandreis, C. Thao, S.  2022\.  Assessment of stochastic weather forecast of precipitation near European cities, based on analogs of circulation Assessment of stochastic weather forecast of precipitation near European cities, based on analogs of circulation.  Geoscientific Model Development15124941–4958.  https://gmd.copernicus.org/articles/15/4941/2022/ 10.5194/gmd-15-4941-2022
  * Kästner (2011) Kastner2011umbrellaKästner, J.  2011\.  Umbrella sampling Umbrella sampling.  WIREs Computational Molecular Science16932-942.  https://wires.onlinelibrary.wiley.com/doi/abs/10.1002/wcms.66 https://doi.org/10.1002/wcms.66
  * Lee . (2020) Lee2020statisticalLee, CY., CaMargo, SJ., Sobel, AH. Tippett, MK.  2020\.  Statistical–Dynamical Downscaling Projections of Tropical Cyclone Activity in a Warming Climate: Two Diverging Genesis Scenarios Statistical–Dynamical Downscaling Projections of Tropical Cyclone Activity in a Warming Climate: Two Diverging Genesis Scenarios.  Journal of Climate33114815 - 4834.  https://journals.ametsoc.org/view/journals/clim/33/11/jcli-d-19-0452.1.xml 10.1175/JCLI-D-19-0452.1
  * Lestang . (2020) Lestang2020numericalLestang, T., Bouchet, F. Lévêque, E.  2020\.  Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques.  Journal of Fluid Mechanics895A19.  10.1017/jfm.2020.293
  * Lestang . (2018) Lestang2018computingLestang, T., Ragone, F., Bréhier, CE., Herbert, C. Bouchet, F.  2018Apr.  Computing return times or return periods with rare event algorithms Computing return times or return periods with rare event algorithms.  Journal of Statistical Mechanics: Theory and Experiment20184043213.  https://doi.org/10.1088/1742-5468/aab856 10.1088/1742-5468/aab856
  * Lorenz (1996) Lorenz1996predictabilityLorenz, EN.  1996\.  Predictability: A problem partly solved Predictability: A problem partly solved.  Proc. Seminar on predictability Proc. Seminar on predictability ( 1).  https://www.cambridge.org/core/books/abs/predictability-of-weather-and-climate/predictability-a-problem-partly-solved/3221BDE379DEB669BA52C66263AF3206
  * Lorenz Emanuel (1998) Lorenz1998optimalLorenz, EN. Emanuel, KA.  1998\.  Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model.  Journal of the Atmospheric Sciences553399 - 414.  https://journals.ametsoc.org/view/journals/atsc/55/3/1520-0469_1998_055_0399_osfswo_2.0.co_2.xml 10.1175/1520-0469(1998)055¡0399:OSFSWO¿2.0.CO;2
  * Lucarini . (2016) Lucarini2016extremesLucarini, V., Faranda, D., de Freitas, JMM., Holland, M., Kuna, T., Nicol, M.others  2016\.  Extremes and recurrence in dynamical systems Extremes and recurrence in dynamical systems.  John Wiley & Sons.
  * Lucente, Herbert Bouchet (2022) Lucente2022committorLucente, D., Herbert, C. Bouchet, F.  2022\.  Committor Functions for Climate Phenomena at the Predictability Margin: The Example of El Niño Southern Oscillation in the Jin and Timmermann Model Committor Functions for Climate Phenomena at the Predictability Margin: The Example of El Niño Southern Oscillation in the Jin and Timmermann Model.  Journal of the Atmospheric Sciences.  https://journals.ametsoc.org/view/journals/atsc/aop/JAS-D-22-0038.1/JAS-D-22-0038.1.xml 10.1175/JAS-D-22-0038.1
  * Lucente, Rolland . (2022) Lucente2022couplingLucente, D., Rolland, J., Herbert, C. Bouchet, F.  2022Aug.  Coupling rare event algorithms with data-based learned committor functions using the analogue Markov chain Coupling rare event algorithms with data-based learned committor functions using the analogue Markov chain.  Journal of Statistical Mechanics: Theory and Experiment20228083201.  https://dx.doi.org/10.1088/1742-5468/ac7aa7 10.1088/1742-5468/ac7aa7
  * Maiocchi . (2024) Maiocchi2024heterogeneityMaiocchi, CC., Lucarini, V., Gritsun, A. Sato, Y.  2024\.  Heterogeneity of the attractor of the Lorenz ’96 model: Lyapunov analysis, unstable periodic orbits, and shadowing properties Heterogeneity of the attractor of the Lorenz ’96 model: Lyapunov analysis, unstable periodic orbits, and shadowing properties.  Physica D: Nonlinear Phenomena457133970.  https://www.sciencedirect.com/science/article/pii/S016727892300324X https://doi.org/10.1016/j.physd.2023.133970
  * Miloshevich . (2023) Miloshevich2022probabilisticMiloshevich, G., Cozian, B., Abry, P., Borgnat, P. Bouchet, F.  2023Apr.  Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data ( 8).  American Physical Society.  https://link.aps.org/doi/10.1103/PhysRevFluids.8.040501 10.1103/PhysRevFluids.8.040501
  * Mohamad Sapsis (2018) Mohamad2018sequentialMohamad, MA. Sapsis, TP.  2018\.  Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems.  Proceedings of the National Academy of Sciences1154411138-11143.  https://www.pnas.org/doi/abs/10.1073/pnas.1813263115 10.1073/pnas.1813263115
  * Myhre . (2019) Myhre2019frequencyMyhre, G., Alterskjær, K., Stjern, CW., Hodnebrog, Ø., Marelle, L., Samset, BH.Stohl, A.  2019Nov05.  Frequency of extreme precipitation increases extensively with event rareness under global warming Frequency of extreme precipitation increases extensively with event rareness under global warming.  Scientific Reports9116063.
  * Naveau . (2020) Naveau2020statisticalNaveau, P., Hannart, A. Ribes, A.  2020\.  Statistical Methods for Extreme Event Attribution in Climate Science Statistical Methods for Extreme Event Attribution in Climate Science.  Annual Review of Statistics and Its Application7189-110.  https://doi.org/10.1146/annurev-statistics-031219-041314 10.1146/annurev-statistics-031219-041314
  * Norwood . (2013) Norwood2013lyapunovNorwood, A., Kalnay, E., Ide, K., Yang, SC. Wolfe, C.  2013Jun.  Lyapunov, singular and bred vectors in a multi-scale system: an empirical exploration of vectors related to instabilities Lyapunov, singular and bred vectors in a multi-scale system: an empirical exploration of vectors related to instabilities.  Journal of Physics A: Mathematical and Theoretical4625254021.  https://dx.doi.org/10.1088/1751-8113/46/25/254021 10.1088/1751-8113/46/25/254021
  * O’Brien . (2016) OBrien2016resolutionO’Brien, TA., Collins, WD., Kashinath, K., Rübel, O., Byna, S., Gu, J.Ullrich, PA.  2016\.  Resolution dependence of precipitation statistical fidelity in hindcast simulations Resolution dependence of precipitation statistical fidelity in hindcast simulations.  Journal of Advances in Modeling Earth Systems82976-990.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016MS000671 https://doi.org/10.1002/2016MS000671
  * O’Gorman (2015) OGorman2015precipitationO’Gorman, PA.  2015Jun01.  Precipitation Extremes Under Climate Change Precipitation Extremes Under Climate Change.  Current Climate Change Reports1249-59.  https://doi.org/10.1007/s40641-015-0009-3 10.1007/s40641-015-0009-3
  * Palmer . (2009) Palmer2009stochasticPalmer, TN., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, GJ.Weisheimer, A.  2009\.  Stochastic parametrization and model uncertainty Stochastic parametrization and model uncertainty.  ECMWF Technical Memoranda.
  * Palmer Zanna (2013) Palmer2013singularPalmer, TN. Zanna, L.  2013Jun.  Singular vectors, predictability and ensemble forecasting for weather and climate Singular vectors, predictability and ensemble forecasting for weather and climate.  Journal of Physics A: Mathematical and Theoretical4625254018.  https://dx.doi.org/10.1088/1751-8113/46/25/254018 10.1088/1751-8113/46/25/254018
  * Pavliotis (2014) Pavliotis2014stochasticPavliotis, GA.  2014\.  Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations ( 60).  Springer.
  * Pazo . (2010) Pazo2010spatiotemporalPazo, D., Rodriguez, MA. Lopez, JM.  2010\.  Spatio-temporal evolution of perturbations in ensembles initialized by bred, Lyapunov and singular vectors Spatio-temporal evolution of perturbations in ensembles initialized by bred, Lyapunov and singular vectors.  Tellus A62110-23.  https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1600-0870.2009.00419.x https://doi.org/10.1111/j.1600-0870.2009.00419.x
  * Pfahl . (2017) Pfahl2017understandingPfahl, S., O’Gorman, PA. Fischer, EM.  2017Jun01.  Understanding the regional pattern of projected future changes in extreme precipitation Understanding the regional pattern of projected future changes in extreme precipitation.  Nature Climate Change76423-427.  https://doi.org/10.1038/nclimate3287 10.1038/nclimate3287
  * Qi Majda (2016) Qi2016predictingQi, D. Majda, AJ.  2016\.  Predicting fat-tailed intermittent probability distributions in passive scalar turbulence with imperfect models through empirical information theory Predicting fat-tailed intermittent probability distributions in passive scalar turbulence with imperfect models through empirical information theory.  Communications in Mathematical Sciences1461687–1722.
  * Ragone Bouchet (2021) Ragone2021rareRagone, F. Bouchet, F.  2021\.  Rare Event Algorithm Study of Extreme Warm Summers and Heatwaves Over Europe Rare Event Algorithm Study of Extreme Warm Summers and Heatwaves Over Europe.  Geophysical Research Letters4812e2020GL091197.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020GL091197 e2020GL091197 2020GL091197  https://doi.org/10.1029/2020GL091197
  * Ragone . (2018) Ragone2018computationRagone, F., Wouters, J. Bouchet, F.  2018\.  Computation of extreme heat waves in climate models using a large deviation algorithm Computation of extreme heat waves in climate models using a large deviation algorithm.  Proceedings of the National Academy of Sciences115124–29.  https://www.pnas.org/content/115/1/24 10.1073/pnas.1712645115
  * Saha Ravela (2022) Saha2022downscalingSaha, A. Ravela, S.  2022\.  Downscaling Extreme Rainfall Using Physical-Statistical Generative Adversarial Learning. Downscaling Extreme Rainfall Using Physical-Statistical Generative Adversarial Learning.  https://arxiv.org/abs/2212.01446
  * Sapsis (2020) Sapsis2020outputSapsis, TP.  2020\.  Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples Output-weighted optimal sampling for Bayesian regression and rare event statistics using few samples.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences476223420190834.  https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2019.0834 10.1098/rspa.2019.0834
  * Schmidli . (2007) Schmidli2007statisticalSchmidli, J., Goodess, CM., Frei, C., Haylock, MR., HunDecha, Y., Ribalaygua, J. Schmith, T.  2007\.  Statistical and dynamical downscaling of precipitation: An evaluation and comparison of scenarios for the European Alps Statistical and dynamical downscaling of precipitation: An evaluation and comparison of scenarios for the European Alps.  Journal of Geophysical Research: Atmospheres112D4.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2005JD007026 https://doi.org/10.1029/2005JD007026
  * Schorlepp . (2023) Schorlepp2023scalableSchorlepp, T., Tong, S., Grafke, T. Stadler, G.  202310\.  Scalable methods for computing sharp extreme event probabilities in infinite-dimensional stochastic systems Scalable methods for computing sharp extreme event probabilities in infinite-dimensional stochastic systems.  Statistics and Computing336.  http://dx.doi.org/10.1007/s11222-023-10307-2 10.1007/s11222-023-10307-2
  * Sheshadri . (2021) Sheshadri2021midlatitudeSheshadri, A., Borrus, M., Yoder, M. Robinson, T.  2021\.  Midlatitude Error Growth in Atmospheric GCMs: The Role of Eddy Growth Rate Midlatitude Error Growth in Atmospheric GCMs: The Role of Eddy Growth Rate.  Geophysical Research Letters4823e2021GL096126.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL096126 e2021GL096126 2021GL096126  https://doi.org/10.1029/2021GL096126
  * Sterk van Kekem (2017) Sterk2017predictabilitySterk, AE. van Kekem, DL.  2017Sep24.  Predictability of Extreme Waves in the Lorenz-96 Model Near Intermittency and Quasi-Periodicity Predictability of Extreme Waves in the Lorenz-96 Model Near Intermittency and Quasi-Periodicity.  Complexity20179419024.  https://doi.org/10.1155/2017/9419024 10.1155/2017/9419024
  * Tandon . (2018) Tandon2018understandingTandon, NF., Zhang, X. Sobel, AH.  2018\.  Understanding the Dynamics of Future Changes in Extreme Precipitation Intensity Understanding the Dynamics of Future Changes in Extreme Precipitation Intensity.  Geophysical Research Letters4562870-2878.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL076361 https://doi.org/10.1002/2017GL076361
  * Tantet . (2015) Tantet2015earlyTantet, A., van der Burgt, FR. Dijkstra, HA.  2015\.  An early warning indicator for atmospheric blocking events using transfer operators An early warning indicator for atmospheric blocking events using transfer operators.  Chaos: An Interdisciplinary Journal of Nonlinear Science253036406.  https://doi.org/10.1063/1.4908174 10.1063/1.4908174
  * Thompson . (2017) Thompson2017unseenThompson, V., Dunstone, NJ., Scaife, AA., Smith, DM., Slingo, JM., Brown, S. Belcher, SE.  2017Jul24.  High risk of unprecedented UK rainfall in the current climate High risk of unprecedented UK rainfall in the current climate.  Nature Communications81107.  https://doi.org/10.1038/s41467-017-00275-3 10.1038/s41467-017-00275-3
  * Touchette (2009) Touchette2009largeTouchette, H.  2009\.  The large deviation approach to statistical mechanics The large deviation approach to statistical mechanics.  Physics Reports4781-31–69.
  * Uribe . (2021) Uribe2021crossUribe, F., Papaioannou, I., Marzouk, YM. Straub, D.  2021\.  Cross-Entropy-Based Importance Sampling with Failure-Informed Dimension Reduction for Rare Event Simulation Cross-Entropy-Based Importance Sampling with Failure-Informed Dimension Reduction for Rare Event Simulation.  SIAM/ASA Journal on Uncertainty Quantification92818-847.  https://doi.org/10.1137/20M1344585 10.1137/20M1344585
  * van der Wiel . (2016) vanderWiel2016resolutionvan der Wiel, K., Kapnick, SB., Vecchi, GA., Cooke, WF., Delworth, TL., Jia, L.Zeng, F.  2016\.  The Resolution Dependence of Contiguous U.S. Precipitation Extremes in Response to CO2 Forcing The Resolution Dependence of Contiguous U.S. Precipitation Extremes in Response to CO2 Forcing.  Journal of Climate29227991 - 8012.  https://journals.ametsoc.org/view/journals/clim/29/22/jcli-d-16-0307.1.xml 10.1175/JCLI-D-16-0307.1
  * Villén-Altamirano . (1991) Villen1991restartVillén-Altamirano, M., Villén-Altamirano, J. .  1991\.  RESTART: a method for accelerating rare event simulations RESTART: a method for accelerating rare event simulations.  Queueing, Performance and Control in ATM (ITC-13)71–76.
  * Walter R. Gilks Sahu (1998) Gilks1998adaptiveWalter R. Gilks, GOR. Sahu, SK.  1998\.  Adaptive Markov Chain Monte Carlo through Regeneration Adaptive Markov Chain Monte Carlo through Regeneration.  Journal of the American Statistical Association934431045-1054.  https://doi.org/10.1080/01621459.1998.10473766 10.1080/01621459.1998.10473766
  * Wang . (2020) Wang2020usefulWang, Q., Mu, M. Sun, G.  2020Jan01.  A useful approach to sensitivity and predictability studies in geophysical fluid dynamics: conditional non-linear optimal perturbation A useful approach to sensitivity and predictability studies in geophysical fluid dynamics: conditional non-linear optimal perturbation.  National Science Review71214-223.  https://doi.org/10.1093/nsr/nwz039 10.1093/nsr/nwz039
  * Wasserman (2004) Wasserman2004allWasserman, L.  2004\.  All of statistics All of statistics.  New YorkSpringer New York, NY.  10.1007/978-0-387-21736-9
  * Webber . (2019) Webber2019practicalWebber, RJ., Plotkin, DA., O’Neill, ME., Abbot, DS. Weare, J.  2019\.  Practical rare event sampling for extreme mesoscale weather Practical rare event sampling for extreme mesoscale weather.  Chaos: An Interdisciplinary Journal of Nonlinear Science295053109.  https://doi.org/10.1063/1.5081461 10.1063/1.5081461
  * Wilks (2005) Wilks2005effectsWilks, DS.  2005\.  Effects of stochastic parametrizations in the Lorenz ’96 system Effects of stochastic parametrizations in the Lorenz ’96 system.  Quarterly Journal of the Royal Meteorological Society131606389-407.  https://rmets.onlinelibrary.wiley.com/doi/abs/10.1256/qj.04.03 https://doi.org/10.1256/qj.04.03
  * Wouters Bouchet (2016) Wouters2016rareWouters, J. Bouchet, F.  2016Aug.  Rare event computation in deterministic chaotic systems using genealogical particle analysis Rare event computation in deterministic chaotic systems using genealogical particle analysis.  Journal of Physics A: Mathematical and Theoretical4937374002.  https://dx.doi.org/10.1088/1751-8113/49/37/374002 10.1088/1751-8113/49/37/374002
  * Wouters . (2023) Wouters2023rareWouters, J., Schiemann, RKH. Shaffrey, LC.  2023\.  Rare Event Simulation of Extreme European Winter Rainfall in an Intermediate Complexity Climate Model Rare Event Simulation of Extreme European Winter Rainfall in an Intermediate Complexity Climate Model.  Journal of Advances in Modeling Earth Systems154e2022MS003537.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022MS003537 e2022MS003537 2022MS003537  https://doi.org/10.1029/2022MS003537
  * Wright . (2021) Wright2021resilienceWright, DB., SaMaras, C. Lopez-Cantu, T.  2021\.  Resilience to Extreme Rainfall Starts with Science Resilience to Extreme Rainfall Starts with Science.  Bulletin of the American Meteorological Society1024E808 - E813.  https://journals.ametsoc.org/view/journals/bams/102/4/BAMS-D-20-0267.1.xml 10.1175/BAMS-D-20-0267.1
  * Yiou Jezequel (2020) Yiou2020simulationYiou, P. Jezequel, A.  2020\.  Simulation of extreme heat waves with empirical importance sampling Simulation of extreme heat waves with empirical importance sampling.  Geoscientific Model Development132763–781.  https://gmd.copernicus.org/articles/13/763/2020/ 10.5194/gmd-13-763-2020
  * Zhang . (2022) Zhang2022koopmanZhang, BJ., Sahai, T. Marzouk, YM.  2022\.  A Koopman framework for rare event simulation in stochastic differential equations A Koopman framework for rare event simulation in stochastic differential equations.  Journal of Computational Physics456111025.  https://www.sciencedirect.com/science/article/pii/S0021999122000870 https://doi.org/10.1016/j.jcp.2022.111025
  * Zuckerman Chong (2017) Zuckerman2017weightedZuckerman, DM. Chong, LT.  2017\.  Weighted Ensemble Simulation: Review of Methodology, Applications, and Software Weighted Ensemble Simulation: Review of Methodology, Applications, and Software.  Annual Review of Biophysics46143-57.  https://doi.org/10.1146/annurev-biophys-070816-033834 PMID: 28301772  10.1146/annurev-biophys-070816-033834
  * Zuev (2015) Zuev2015subsetZuev, K.  2015\.  Subset Simulation Method for Rare Event Estimation: An Introduction. Subset Simulation Method for Rare Event Estimation: An Introduction.  https://arxiv.org/abs/1505.03506

## Supporting Information for “Bringing statistics to storylines: rare event
sampling for sudden, transient extreme events”

Contents

  1. 1.

Figures S1 to S13

Introduction Figs. S1 and S2 display results for AMS applied to the stochastic
L96 model with $a=0$ and $F_{4}=3$, which is really just an array of
correlated OU processes with no advection. Figs. S3-S13 display return level
vs. return period plots for all combinations of stochastic forcing level
$F_{4}\in\\{3,1,0.5,0.25\\}$ and the advance splitting time
$\delta\in\\{0,0.2,0.4,...,2\\}$, only a subset of which are shown in the main
text.

Figure S1: Example of a single lineage generated by AMS applied to the the OU
process (L95 with $a=0,F_{4}=3$), formatted the same as Fig. 4a,b of the main
text. Figure S2:  Statistical results of AMS applied to the OU process (L96
with $a=0,F_{4}=3$) with $N=128$ initial ensemble members and $M=56$ runs.
Format is the same as Fig. 5a,b,c of the main text.  Figure S3: TEAMS
algorithm performance at all four noise levels with $\delta=0.0$ Figure S4:
TEAMS algorithm performance at all four noise levels with $\delta=0.2$ Figure
S5: TEAMS algorithm performance at all four noise levels with $\delta=0.4$
Figure S6: TEAMS algorithm performance at all four noise levels with
$\delta=0.6$ Figure S7: TEAMS algorithm performance at all four noise levels
with $\delta=0.8$ Figure S8: TEAMS algorithm performance at all four noise
levels with $\delta=1.0$ Figure S9: TEAMS algorithm performance at all four
noise levels with $\delta=1.2$ Figure S10: TEAMS algorithm performance at all
four noise levels with $\delta=1.4$ Figure S11: TEAMS algorithm performance at
all four noise levels with $\delta=1.6$ Figure S12: TEAMS algorithm
performance at all four noise levels with $\delta=1.8$ Figure S13: TEAMS
algorithm performance at all four noise levels with $\delta=2.0$