# Identifying the onset and decay of
quasi-stationary families of almost-invariant sets
with an application to atmospheric blocking events

Aleksandar Badza and Gary Froyland
School of Mathematics and Statistics
UNSW Sydney NSW 2052, Australia

###### Abstract

Macroscopic features of dynamical systems such as almost-invariant sets and
coherent sets provide crucial high-level information on how the dynamics
organises phase space. We introduce a method to identify time-parameterised
families of almost-invariant sets in time-dependent dynamical systems, as well
as the families’ emergence and disappearance. In contrast to coherent sets,
which may freely move about in phase space over time, our technique focuses on
families of metastable sets that are quasi-stationary in space. Our
straightforward approach extends successful transfer operator methods for
almost-invariant sets to time-dependent dynamics and utilises the Ulam scheme
for the generator of the transfer operator on a time-expanded domain. The new
methodology is illustrated with an idealised fluid flow and with atmospheric
velocity data. We identify atmospheric blocking events in the 2003 European
heatwave and compare our technique to existing geophysical methods of blocking
diagnosis.

##### Significance statement:

Diagnosing the emergence of persistent structures in turbulent nonlinear flows
is important for many scientific and engineering applications. Computational
and theoretical challenges include detecting the location of almost-invariant
structures at specific times, piecing together the diagnosed objects across
time while simultaneously verifying they are approximately fixed in space, and
determining when these structures appear and cease to exist. Our novel
spacetime approach based on an inflated generator overcomes these challenges
in a single eigencomputation, and relies only on the underlying time-dependent
velocity field. We also outline a simpler averaged generator to provide a
first-cut diagnosis of persistence in space.

##  1 Introduction

Lagrangian approaches to analysing complex flow systems have had great success
over the last two decades. Very broadly speaking these approaches may be
categorised as (i) those that linearise the nonlinear dynamics directly [47,
43, 27, 28, 55, 26, 38] in phase space, through the study of spatial
derivative matrices or the Cauchy–Green tensor, and (ii) those that linearise
the nonlinear dynamics by moving to a space of scalar-valued functions defined
on the phase space, and studying the transfer operator (or cocycles thereof)
[9, 52, 48, 24, 21, 38, 23]. The present work takes the second route and
extends one of the earliest constructions in this approach, namely almost-
invariant sets [9, 52, 15, 22, 18]. As the name suggests, almost-invariant
sets are those sets that remain approximately fixed in space under the
dynamics. They have been an important forerunner to the concept of coherent
sets [24, 16], which share the approximate trajectory-trapping property of
almost-invariant sets, but are able to rapidly move about in phase space as
the dynamics evolve. Much progress has been made on coherent set analysis over
the last decade and one of the aims of this work is to update the concept of
almost-invariant sets by incorporating some of these advances. This aim is not
purely academic. There are real-world situations, such as the atmospheric
blocking events we study here, where one is interested in identifying time-
parameterised families of metastable regions that remain approximately
stationary in space over a significant time duration. This latter criterion is
a restriction that is not present in the definition of coherent sets.

Let $M\subset\mathbb{R}^{d}$ be a smooth full-dimensional compact submanifold
of $\mathbb{R}^{d}$ and let $m$ represent $d$-dimensional volume on $M$. We
assume that the dynamics on $M$ is governed by a smooth time-dependent vector
field $v:[0,\tau]\times M\to\mathbb{R}^{d}$, i.e. that trajectories satisfy
$\dot{x}(t)=v(t,x)$ for $t\in[0,\tau],x\in M$. At each $t\in[0,\tau]$ we seek
one or more almost-invariant sets $A_{t}\subset M$, where the notion of
almost-invariance is made precise in the following section. Collectively,
across time, these individual almost-invariant sets form one or more time-
parameterised families $\\{A_{t}\\}_{t\in[0,\tau]}$, which we require to be
quasi-stationary in the sense that the elements of this family do not vary
much across time; this will be made precise in Section 3. We introduce a
theoretical and numerical methodology for identifying quasi-stationary
families of almost-invariant sets within a long time interval $[0,\tau]$. Of
particular interest is that the family may only exist for times in a strict
subset of $[0,\tau]$. That is, our technique also identifies when these quasi-
stationary families appear and disappear.

Our approach is based on generators of the Perron–Frobenius operators for
flows governed by the vector fields $v(t,x)$ for various fixed $t\in[0,\tau]$.
For steady flows, the leading eigenfunctions of single generators effectively
identify almost-invariant sets [18]. In our time-dependent setting, we connect
different generators across time using ideas from recent work on dynamic
Laplacians on time-expanded domains [19, 2]. This leads to a single inflated
generator acting on spacetime, whose leading eigenfunctions encode the onset
and decay of quasi-stationary families of almost-invariant sets. Our main
contributions include creating and formalising the notion of quasi-stationary
families of almost-invariant sets, demonstrating quantitative relationships
between the eigenspectrum of the inflated generator and almost-invariance
properties of quasi-stationary families, and designing an efficient
computational approach.

A further contribution is our provision of a new approach to diagnosing
atmospheric blocking events, which solely relies on the use of wind fields
rather than other geophysical proxies. A blocking event, or block for short,
is a quasi-stationary, anti-cyclonic air mass formed from the splitting of
westerlies at mid-latitude levels in the Northern Hemisphere. In the Southern
Hemisphere, blocks are analogously defined as cyclonic rather than anti-
cyclonic air masses. This trapping of air parcels inhibits regular atmospheric
flow and can result in increased and prolonged extreme heat [11, 4, 12, 59],
extreme cold [11, 6, 41, 7], or reduced air quality through increased
concentration of pollutants [3, 63, 7]. Despite this phenomenon being
conceptualised decades ago [45, 46], much is still not known about the
underlying physics behind blocks. A better understanding of blocking physics
is crucial for developing clearer links between extreme weather and blocks,
and will lead to better forecasting of and preparedness for these extreme
weather events [37, 57, 13].

Prior research into blocks in the geophysics literature has often focused on
statistics; that is, blocking frequencies at particular lines of latitude.
These frequencies are estimated using blocking proxies, such as geopotential
height gradients [36, 60, 61, 50, 44, 59], geopotential height anomalies [11,
49, 4, 44, 59] or potential vorticity anomalies [53, 50, 44, 59]. This
frequency analysis largely prohibits further examination of the key blocking
physics and impacts on weather patterns. More recently, post-analysis methods
investigating the dynamics of blocks that have been determined from
geophysical proxies have appeared [31, 57, 29].

A range of dynamical systems methods for block detection have also been
developed and implemented on observational data [14, 54, 40, 13, 62] or
simulated data obtained from barotropic vorticity models [56, 58, 51]. These
methods either (i) capture the two-dimensional extent of blocking events at a
single instance of time (in particular without capturing the full evolution of
a block over its lifespan) [56, 54, 51, 13], (ii) identify longitude levels
(but not two-dimensional regions) which remain blocked over a period of time
[40], (iii) detect the timing of the transition of an atmospheric flow regime
to a blocked state and don’t provide information on the spatial location of
the block [14, 58], or (iv) ignore the quasi-stationarity property of blocks
[56, 54, 13, 62].

The aim of our new approach is to provide a full description of a blocking
event over its lifetime, by pinpointing the block’s onset and decay in time,
and tracking its spatial extent throughout its lifetime. Our technique
directly uses wind velocities and does not require other geophysical proxies.
We explicitly incorporate the crucial quasi-stationarity property of blocks
[45, 46, 56, 31, 57, 33, 29], a property which is often of peripheral
importance in frequency analysis [45, 36, 60, 50, 4, 12] but increasingly used
in other dynamical systems based approaches [14, 58, 51, 40, 62]. To
demonstrate this new approach, we study two blocking events from the Northern
Hemisphere summer of 2003, during what is regarded as some of the worst
extreme heat events recorded in Europe [5, 13, 59].

The outline of this paper is as follows. In Section 2, we recap relevant
theory concerning generators of steady flows from [18]. This section ends with
a bridging construction of an averaged generator. In Section 3 we construct
the inflated generator, inspired by the constructions in [18, 19], and detail
how this operator is used to identify quasi-stationary families of almost-
invariant sets. In Section 4, we describe the numerical implementation of the
inflated generator method, including the discretisation of the generator,
selection of key parameters for the generator, and the extraction of time-
parameterised families of almost-invariant sets from its spacetime
eigenfunctions. In Section 5, we illustrate our theory on an idealised fluid
flow that possesses well-defined quasi-stationary almost-invariant sets. In
Section 6, we apply the inflated generator method to European Centre for
Medium-Range Weather Forecasts (ECMWF) atmospheric velocity data to identify
two key blocking events that occurred in the European summer of 2003 [5, 59].
We demonstrate that our inflated generator method produces a clear, detailed
visualisation of the evolution of a block from its onset right up until its
decay. We show that identical results are not obtained from three frequently-
used geophysical block detection methods, partly because of their reliance on
blocking proxies defined through scalar fields on individual time slices,
without explicit consideration of quasi-stationarity of the block.

##  2 Almost-invariant sets and generators of steady flows

As a first step toward our main constructions in the following section, we
suppose for the moment that we have a steady flow given by the vector field
$v(t,\cdot)$ at some fixed time $t$. To simplify the mathematical theory of
spectra of our generators, we will consider the stochastic flow given by

$dx_{s}=v(t,x)\,ds+\epsilon\,dw_{s},$ (1)

where $\epsilon$ is a small positive number representing the magnitude of
Brownian motion $\\{w_{s}\\}$. The corresponding forward propagator of a
function $f:M\to\mathbb{R}$ at time $t$ for duration $s$ under the fixed drift
$v(t,\cdot)$ is given by the Perron–Frobenius operator $P_{t}^{s}.$

Following [18], given a probability measure $\nu$ on $M$ we say that a set
$A\subset M$ is almost-invariant with respect to $\nu$ at time $t$ over a
short time duration $s>0$ if

$\rho_{t,\nu}^{s}(A):=\frac{{\mathrm{Prob}}_{\nu}(x_{t}\in A,x_{t+s}\in
A)}{{\mathrm{Prob}}_{\nu}(x_{t}\in A)}\approx 1,$

where ${\mathrm{Prob}}_{\nu}$ refers to the distribution of points in space at
time $t$ according to $\nu$. If we wish to find almost-invariant sets for this
steady flow in $M$, following the approach in [18] we construct the generator
of the Perron–Frobenius operator $P_{t}^{s}$ for this steady flow at this
fixed time $t$:

$G_{t}f(x):=-\nabla\cdot(v(t,x)f(x))+\frac{\epsilon^{2}}{2}\Delta f(x).$ (2)

To connect the generator and $\rho^{s}_{t,\nu}$, we define the functional

$G_{t,A}f:=\lim_{s\to 0}\int_{A}\frac{P_{t}^{s}f-f}{s}\ dm,$ (3)

which is well defined for $f\in D_{A}(G_{t})$ (denoting the linear subspace of
$L^{1}$ where this limit exists). For example, denoting the indicator on $A$
by $\mathbf{1}_{A}$, Proposition 3.2 [18] states

$\rho^{s}_{t,m}(A)=1+\frac{G_{t,A}\mathbf{1}_{A}}{m(A)}\cdot s+o(s),$ (4)

with $G_{t,A}\mathbf{1}_{A}<0$. Thus, almost-invariant sets at time $t$ are
those sets $A$ for which $G_{t,A}\mathbf{1}_{A}/m(A)\lessapprox 0$. Our
strategy is to use the leading eigenfunctions of $G_{t}$ to find almost-
invariant sets at time $t$. We now briefly motivate this approach recapping
another result from [18].

First, applying zero Neumann (reflecting) boundary conditions at the boundary
of $M$, we have a well-posed eigenvalue problem:

$G_{t}f=\lambda f,\qquad\nabla f(x)\cdot n(x)=0\mbox{ for }x\in\partial M,$
(5)

where $n(x)$ is a vector normal to $\partial M$ at $x\in\partial M$. By
standard theory, e.g. Theorem 8.4.4 [39], the spectrum of $G_{t}$ (considered
as an operator from $H^{2}(M)\cap H_{0}^{1}(M)$ to $L^{2}(M)$) consists of
isolated eigenvalues $\\{\lambda_{k}\\}_{k=1}^{\infty}$ of finite multiplicity
and no finite accumulation points. We are primarily concerned with the
eigenvalues of largest (negative) real part, which corresponds to the
eigenfunctions with the slowest rates of decay under the action of the
generator $G_{t}$.

######  Theorem 1 (Theorem 3.3 [18]).

Let $\lambda<0$ be a real eigenvalue for $G_{t}$ with eigenfunction $f$,
normalised so that $\int_{M}|f|\ dm=1$. Then $\int_{M}f\ dm=0$. Define the
signed measure $\nu(A)=\int_{A}f\ dm$, subsets $A^{+}=\\{f\geq 0\\}$,
$A^{-}=\\{f<0\\}$, and suppose $f\mathbf{1}_{A^{\pm}}\in D_{A}(G_{t})$. Then

$G_{t,A^{+}}(f\mathbf{1}_{A^{+}})/\nu(A^{+})+G_{t,A^{-}}(f\mathbf{1}_{A^{-}})/\nu(A^{-})=\lambda.$
(6)

If we choose $\lambda<0$ to be a real eigenvalue close to 0 then the sets
$A^{\pm}\subset M$ will be almost-invariant with respect to the probability
measure with density $|f|$. Thus, one strategy for finding almost-invariant
sets is to look at positive and negative parts of dominant eigenfunctions of
$G_{t}$. This idea goes back to [9] and has been refined over the years, by
way of optimising thresholds to consider $A^{+}=\\{f\geq c\\}$,
$A^{-}=\\{f<c\\}$ [15], considering multiple eigenvectors [17, 15], and the
subsequent extraction of multiple almost-invariant sets by various means, e.g.
$c$-means clustering [15], constrained nonlinear optimisation [10], and sparse
optimisation [23]. We will use the latter approach of sparse eigenbasis
approximation (SEBA) [23], which is robust and well suited to extracting many
features, particularly when the features do not fully partition the phase
space.

##### An averaged generator for unsteady flows.

We have just discussed connections between generators of steady vector fields
$v(t,\cdot)$ for fixed $t$ and almost invariant sets for the steady flows. As
a partial bridge to the next section, we now consider the full time-dependent
vector field $v(t,x)$ for $0\leq t\leq\tau$, and ask for sets in $M$ that are
almost-invariant in a time-averaged sense across $[0,\tau]$. That is, we seek
$A\subset M$ for which the average residence fraction

$\bar{\rho}^{s}_{m}(A):=\int_{0}^{\tau}\rho^{s}_{t,m}(A)\ dt$

is close to 1. Define $\bar{v}(x)=\frac{1}{\tau}\int_{0}^{\tau}v(t,x)\ dt$ to
be the averaged vector field at a fixed location $x$ across the time interval
$[0,\tau]$. We now introduce the time-averaged generator:

$\bar{G}f(x):=\frac{1}{\tau}\int_{0}^{\tau}G_{t}f(x)\
dt=-\nabla\cdot(\bar{v}(x)f(x))+\frac{\epsilon^{2}}{2}\Delta f(x).$ (7)

From (4) we immediately obtain

$\bar{\rho}^{s}_{m}(A)=1+\frac{\bar{G}_{A}\mathbf{1}_{A}}{m(A)}\cdot s+o(s),$
(8)

with

$\bar{G}_{A}f:=\lim_{s\to 0}\int_{A}\frac{\bar{P}^{s}f-f}{s}\ dm,$ (9)

and $\bar{P}^{s}$ denoting the propagator (Perron–Frobenius operator) of $f$
from time $0$ to $s$ under the SDE (1) with $v(t,x)$ replaced by $\bar{v}$. An
analogous version of Theorem 1 holds with $G_{t},G_{t,A^{+}},G_{t,A^{-}}$
replaced by $\bar{G},\bar{G}_{A^{+}},\bar{G}_{A^{-}}$, respectively. In
summary, to find subsets that are persistently almost invariant across the
full finite time interval $[0,\tau]$ one can use eigenfunctions of $\bar{G}$
in the algorithms that follow. Our ultimate aim is more ambitious, namely
finding quasi-stationary families of almost-invariant sets, including their
emergence and disappearance. The averaged generator may be viewed as a
simpler, less computationally intensive, first-cut analysis tool.

##  3 Quasi-stationary families of almost-invariant sets and an inflated
generator

In this section we explain how to extend the above ideas to find quasi-
stationary families of almost-invariant sets in a time-dependent flow over a
finite time duration $[0,\tau]$, which is one of the main contributions of
this work. We define a time-expanded domain $\mathbb{M}=[0,\tau]\times M$ and
construct a stochastic process on this domain as follows.

  1. 1.

In the one-dimensional time coordinate we have diffusion in time with
intensity $a>0$ and zero drift,

  2. 2.

In the original space coordinates, on $M$ at time $t$ (i.e. on $\\{t\\}\times
M\subset\mathbb{M}$) we have diffusion with intensity $\epsilon$ and drift
given by $v(t,\cdot)$.

More precisely, we define points in $\mathbb{M}$ as $X=(t,x)$, a vector field
$V:\mathbb{M}\to\mathbb{R}^{d+1}$ as $V(t,x)=\begin{pmatrix}0\\\
v(t,x)\end{pmatrix}$ and a $(d+1)\times(d+1)$ diffusion matrix
$\Sigma=\begin{pmatrix}a&0\\\ 0&\epsilon I_{d}\end{pmatrix}$, and consider the
stochastic flow given by

$dX_{s}=V(X_{s})\,ds+\Sigma\,dW_{s},$ (10)

where $dW_{s}$ is Brownian motion in $\mathbb{M}$. Note that the SDE (10) is
homogeneous in time due to the time extension. The idea behind this
construction is as follows. Let $\mathbbm{m}$ denote volume on $\mathbb{M}$.
For a set $\mathbb{A}\subset\mathbb{M}$, denote the $t^{\rm th}$ time fibre
$A_{t}=\mathbb{A}\cap(\\{t\\}\times M)$, which we consider as a subset of $M$;
see Figure 1.

Figure 1:  A schematic of a space-time set
$\mathbb{A}\subset\mathbb{M}=[0,\tau]\times M$ for one-dimensional $M$. We
seek $\mathbb{A}$ that are almost-invariant under the spacetime SDE (10). A
time fibre $A_{t}$ of $\mathbb{A}$ is shown in grey. Trajectories may exit
$\mathbb{A}$ through $\partial\mathbb{A}$ via drift $v$ in space and/or
diffusion in space of strength $\epsilon$, indicated by the solid and dotted
blue arrows respectively. Therefore, the boundary of $\mathbb{A}$ should
attempt to trace through spacetime so that the component of $v(t,x)$ normal to
$\partial\mathbb{A}$ is small, to minimise trajectory loss through drift; this
will ensure that each $A_{t}$ is almost-invariant. Trajectories may also exit
$\mathbb{A}$ through $\partial\mathbb{A}$ via diffusion in the temporal
direction of strength $a$, indicated by the dotted red arrows. If the time
fibres vary slowly in time, the loss of trajectories due to temporal diffusion
will be small; this will ensure that the family $\\{A_{t}\\}_{t=0}^{\tau}$ is
quasi-stationary. Dark blue and red arrows indicate trajectory loss from
$\mathbb{A}$, while pale blue and red arrows illustrate motion internal to
$\mathbb{A}$.

We may now write the space-time set $\mathbb{A}$ as the disjoint union of the
sets $A_{t}$ placed on their individual time-$t$ fibres:
$\mathbb{A}=\bigcup_{t\in[0,\tau]}(\\{t\\}\times A_{t})$. Note that the
dynamics restricted to the time-fibre $A_{t}$ are exactly the spatial dynamics
we discussed in the previous section. For $s>0$ consider the ratio

$\rho_{\mathbbm{m}}^{s}(\mathbb{A}):=\frac{{\mathrm{Prob}}_{\mathbbm{m}}(X_{0}\in\mathbb{A},X_{s}\in\mathbb{A})}{{\mathrm{Prob}}_{\mathbbm{m}}(X_{s}\in\mathbb{A})}.$

If $\rho_{\mathbbm{m}}^{s}(\mathbb{A})$ is close to 1, this means that on
average for a flow duration $s$:

  1. 1.

Spatial almost-invariance: the dynamics of (10) restricted to each time fibre
$A_{t}$ (given by (1)) is such that trajectories rarely leave $A_{t}$.

  2. 2.

Quasi-stationary family in time: the variation of $A_{t}$ with $t$ is small in
order that the pure diffusion in the time direction in (10) allows only rare
exits from the family $\\{A_{t}\\}_{t\in[0,\tau]}$.

It is important to note that the quasi-stationary family
$\\{A_{t}\\}_{t=0}^{\tau}$ may consist of nonempty $A_{t}$ only for a subset
of $[0,\tau]$, for example, for $t\in[t_{1},t_{2}]$, where
$0<t_{1}<t_{2}<\tau$, as illustrated in the example in Figure 1. Thus by
considering almost-invariant sets $\mathbb{A}$ of the SDE (10) we achieve a
third property, in addition to the two above, namely we see when a quasi-
stationary family of almost-invariant sets is born (when the $A_{t}$ become
nonempty) and when this family dissipates (when the $A_{t}$ become empty
again).

We now define a generator on $\mathbb{M}$ in order to identify such quasi-
stationary families of almost-invariant sets from its eigenfunctions. We
define the inflated generator of the Perron–Frobenius operator $\mathcal{P}$
for the flow generated by the SDE (10) by

$\mathcal{G}F(t,x)=-\nabla\cdot(V(t,x)F(t,x))+\frac{\Sigma^{2}}{2}\Delta
F(t,x),$ (11)

where $F:\mathbb{M}\to\mathbb{C}$. In analogy with (3), define the functional

$\mathcal{G}_{\mathbb{A}}F:=\lim_{s\to
0}\int_{\mathbb{A}}\frac{\mathcal{P}^{s}F-F}{s}\ d\mathbbm{m}.$ (12)

Because $\mathcal{G}$ is a generator of a homogeneous process we immediately
obtain (again from Proposition 3.2 [18])

$\rho_{\mathbbm{m}}^{s}(\mathbb{A})=1+\frac{\mathcal{G}_{\mathbb{A}}\mathbf{1}_{\mathbb{A}}}{\mathbbm{m}(\mathbb{A})}\cdot
s+o(s).$ (13)

Thus, almost-invariant sets $\mathbb{A}$ satisfy
$\mathcal{G}_{\mathbb{A}}\mathbf{1}_{\mathbb{A}}/\mathbbm{m}(\mathbb{A})\lessapprox
0$. We apply zero Neumann boundary conditions at the boundary of $\mathbb{M}$
to obtain the eigenvalue problem:

$\mathcal{G}F=\Lambda F,\qquad\nabla F(x)\cdot N(x)=0\mbox{ for
}X\in\partial\mathbb{M},$ (14)

where $N(X)$ is a vector normal to $\partial\mathbb{M}$ at
$X\in\partial\mathbb{M}$. As in the previous subsection, the spectrum of
$\mathcal{G}$ consists of isolated eigenvalues
$\\{\Lambda_{k}\\}_{k=1}^{\infty}$ of finite multiplicity and no finite
accumulation points, and we are interested in the eigenvalues with largest
(negative) real part. It is not hard to see that an analogous version of
Theorem 1 holds.

In practice we will work with the $L^{2}$-adjoint of the inflated generator:

$\mathcal{G}^{*}F(t,x)=V(t,x)\cdot\nabla F(t,x)+\frac{\Sigma^{2}}{2}\Delta
F(t,x).$ (15)

This is consistent with the use of right eigenvectors of Markov chain
approximations of Perron–Frobenius operators to identify almost-invariant sets
[15, 22]. In [15, 22], the transition matrix $P$ of a Markov chain is row-
stochastic, requiring vector multiplication on the left to evolve forward in
time. Thus, right eigenvectors111Strictly speaking, in [15, 22] a symmetrised
transition matrix $R$ was constructed from the non-symmetric transition matrix
$P$ and the right eigenvectors of $R$ were used. Nevertheless the same ideas
apply to the right eigenvectors of the transition matrix $P$ itself, without
the optimality properties afforded by the symmetrisation. of $P$ correspond to
left eigenvectors of $P^{*}$, and therefore to eigenfunctions of
$\mathcal{G}^{*}$ in our current context.

A second reason to consider $\mathcal{G}^{*}$ is to produce robust results
when flows have regions of convergence in $M$. If the vector field $v(t,x)$
has stable approximately fixed points (or more general attracting sets) then
the eigenfunctions of $G_{t}$ will tend to localise about these. On the other
hand, the leading eigenfunction of $G_{t}^{*}$ – corresponding to the
eigenvalue 0 – is always $f_{1}=\mathbf{1}_{M}$, and the remaining leading
eigenfunctions of $G_{t}^{*}$ tend to have a strong level set structure, where
each approximate level set indicates an approximate basin of attraction of a
given stable nearly-fixed point. This is explained in more detail in section V
[25], in the context of finite-state Markov chains with multiple absorbing
regions and applied to map plastic garbage patches in the ocean. There the
transition matrices $P$ of the Markov chains are again row-stochastic and
right eigenvectors are used, which corresponds to left eigenvectors of $P^{*}$
(and therefore to eigenvectors of $G_{t}^{*}$ in our current context). See
also [35] for a related approach based on classical formulae for absorption
probabilities in the presence of a single absorbing state.

The spectrum of $\mathcal{G}^{*}$ also consists of isolated eigenvalues
$\\{\Lambda_{k}\\}_{k=1}^{\infty}$ of finite multiplicity and no finite
accumulation points. The eigenvalue $\Lambda_{1}=0$ has the corresponding
constant eigenfunction $F_{1}^{*}=\mathbf{1}_{\mathbb{M}}$. Because
$\mathcal{G}^{*}$ is the adjoint of the generator of a stochastic process, the
same considerations in [15, 22, 25] apply, and leading eigenfunctions
$F_{1}^{*},F_{2}^{*},\ldots$ of $\mathcal{G}^{*}$ will tend to have level set
structures on almost-invariant sets222When attracting sets exist, we expect
level set structures on basins of attraction of almost-invariant sets. in
$\mathbb{M}$.

The space-time block structure of $V$ and $\Sigma$ allows us to immediately
identify some of the eigenfunctions of $\mathcal{G}^{*}$. Following the
arguments in Section 3.4.2 [19], we note that the constant-in-space functions
$F_{k}^{*}(t,\cdot):=\cos(k\pi t/\tau)$ are eigenfunctions of
$\mathcal{G}^{*}$ with eigenvalues $\Lambda_{k}^{\rm
temp}:=-a^{2}\pi^{2}k^{2}/2\tau^{2}$. We call the $F_{k}^{*}$ temporal
eigenfunctions because they only vary in time and carry no spatial
information. We typically want to choose the temporal diffusion strength $a>0$
large enough so that there is only around one temporal eigenvalue appearing
before we begin to see the non-temporal (spatial) part of the spectrum, whose
eigenvalues contain useful spatiotemporal information relevant to the dynamics
of $\mathbb{M}$.

Finally, we note that if $F\in L^{2}(\mathbb{M})$ is an eigenfunction of
$\mathcal{G}$ corresponding to a non-temporal (spatial) eigenfunction of
$\mathcal{G}^{*}$ then $F$ is orthogonal to all $F_{k}^{*}$. Let
$W_{0}=\\{f\in L^{2}([0,\tau]):\int_{0}^{\tau}f(t)\
dt=0\\}=\mathbf{1}^{\perp}_{[0,\tau]}$. Since the temporal eigenfunctions are
dense in $\\{f\mathbf{1}_{M}:f\in W_{0}\\}\subset L^{2}(\mathbb{M})$, we may
follow the argument of Section 4.2 [19] to see that $\int_{M}F(\cdot,x)\
d\ell(x)=\mathrm{const}\ \mbox{a.e.\ on $[0,\tau]$}$. This provides a
numerical test for whether $F_{k}^{*}$ is a temporal eigenfunction which
involves computing the variation of the spatial integrals of the corresponding
$F_{k}$. In practice it is often clear by inspection which eigenfunctions
$F_{k}^{*}$ are temporal.

##  4 Discretisation of generators

This section concerns the numerical discretisation of $\mathcal{G}$. Because
of the partially separable structure of $\mathcal{G}$ in space and time
evident in (18), we discretise time and space separately. The time interval
$[0,\tau]$ is sampled with equispaced times $0=t_{0}<t_{1}<\cdots<t_{T}=\tau$
with $t_{l}=l\tau/T$, $l=0,\ldots,T$. The manifold $M$ is discretised into
boxes $B_{1},\ldots,B_{N}\subset\mathbb{R}^{d}$ with common side lengths
$\ell_{1},\ldots,\ell_{d}$. For each $t_{l}$, $l=0,\ldots,T$, we discretise
the copy $\\{t_{l}\\}\times M$ of the manifold $M$. This $t_{l}$-fibre is
simply a copy of $M$ on which the vector field is $v(t_{l},\cdot)$.

###  4.1 Discretisation of spatial generators on each time fibre

For simplicity, we use the “Ulam’s method for the generator” scheme elaborated
in [18]. The main construction is given in [18, Equation (4.4)], which
involves a simple spatial integration on the faces of the cubes discretising
$M$ with respect to $d-1$-dimensional volume, denoted $m_{d-1}$. For the
computation on a single time-fibre $\\{t_{l}\\}\times M$, the integrand on a
particular cube face $B_{i}\cap B_{j}$ is the dot product of the vector field
$v(t_{l},\cdot)$ with the outward-pointing unit normal vector $n_{ij}$,
representing the rate of flux out of the face. With $N$ cubes
($\ell_{1}=\ell_{2}=\cdots=\ell_{d}$), one obtains a sparse $N\times N$ matrix
$\mathbf{G}^{t_{l}}$ where $\mathbf{G}^{t_{l}}_{ij}$ is the proportional rate
of outflow from cube $B_{i}$ into $B_{j}$.

$\mathbf{G}^{t_{l}}_{ij}=\begin{cases}\frac{1}{m(B_{j})}\left(\int_{B_{i}\cap
B_{j}}(\max\\{v(t_{l},x)\cdot n_{ij},0\\})\
dm_{d-1}(x)\right)+\epsilon^{2}/2\ell^{2},&\quad i\neq j,\\\ -\sum_{j\neq
i}\mathbf{G}^{t_{l}}_{ij},&\quad\mbox{otherwise.}\end{cases}$ (16)

The term $\epsilon^{2}/2\ell^{2}$ in (16), not present in [18], is a central-
difference estimate of $\epsilon^{2}/2$-strength diffusion that we add to
ensure a simple leading eigenvalue. The small spatial diffusion parameter
$\epsilon$ is set to be

$\epsilon=\sqrt{0.1\bar{v}\ell}$ (17)

where333Since $\epsilon^{2}/2\ell^{2}=0.05\bar{v}/\ell$, one may view this
diffusion term as adding a drift of 5% of the median velocity norm $\bar{v}$
to the integrand in (16). $\bar{v}$ is the median of $\|v(t,x)\|_{2}$ over
$\mathbb{M}$. The only nonzero entries of $\mathbf{G}$ are those where $B_{i}$
and $B_{j}$ share a face. If $\mathbf{g}\in\mathbb{R}^{N}$ is a vector whose
entries represent the mass contained in each cube, then the matrix-vector
product $\mathbf{g}\mathbf{G}$ is a vector containing the rates of change of
the masses in each cube. In our experiments we use julia’s `QuadGK.jl` package
[32] to estimate the integrals in (16) with Gauss-Kronrod quadrature of order
1 with absolute and relative tolerance of $10^{-2}$. In section 6.2 we
generalise (16) to the case of non-identical non-cubic partitions.

###  4.2 Discretisation of the inflated generator

We begin by noting that we can rewrite (11) as

$\mathcal{G}F(t,x)=G_{t}F(t,x)+\frac{a^{2}}{2}\partial_{t}^{2}F(t,x).$ (18)

Each $G_{t}$ is estimated as in the previous subsection as
$\mathbf{G}^{t_{l}}$, $l=0,\ldots,T$. For the second term on the RHS of (18),
we estimate the one-dimensional Laplacian on $[0,\tau]$ by second-order
central differencing. The separation of time and space in the matrix
constructions below follows ideas in [2]. Let $h=t_{l+1}-t_{l}=\tau/T$ be the
discrete spacing on the one-dimensional temporal grid in $[0,\tau]$ and $a$ be
the temporal diffusion strength parameter. The indexing we use for the
discretised manifold $\mathbb{M}$ is: the indices of all cubes in
$\\{t_{0}\\}\times M$, the indices of all cubes in $\\{t_{1}\\}\times M$, and
so on, ending with the indices of all cubes in $\\{t_{T}\\}\times M$.
Including the discretised first components of (18) (namely the
$\mathbf{G}^{t}$) across time and suitably including the discretised second
component of (18) (namely a central difference matrix for one-dimensional
diffusion) across time, we obtain discretised representations
$\mathbf{G}_{\mathrm{spacetime}}$ and
$(a^{2}/2)\mathbf{L}_{\mathrm{spacetime}}$ as below. Each of these matrices is
$TN\times TN$ in size, with $T^{2}$ blocks that are each of size $N\times N$:

$\mathbf{G}_{\mathrm{spacetime}}:=\begin{pmatrix}\mathbf{G}^{t_{0}}&0&0&0&\cdots&0\\\
0&\mathbf{G}^{t_{1}}&0&0&\cdots&0\\\ 0&0&\mathbf{G}^{t_{2}}&0&&0\\\
0&0&0&\ddots&&\vdots&\\\ \vdots&&&\ddots&&\vdots&\\\
0&0&\cdots&0&0&\mathbf{G}^{t_{T}}\end{pmatrix}$
$\mathbf{L}_{\mathrm{spacetime}}:=\begin{pmatrix}-I_{N}&I_{N}&0&0&\cdots&0\\\
I_{N}&-2I_{N}&I_{N}&0&\cdots&0\\\ 0&I_{N}&-2I_{N}&I_{N}&&\vdots\\\
0&0&I_{N}&\ddots&\ddots&0\\\ \vdots&\vdots&&\ddots&-2I_{N}&I_{N}\\\
0&0&\cdots&0&I_{N}&-I_{N}\end{pmatrix}/h^{2}.$

As per (18) these two matrices are combined to form the discretised inflated
generator

$\mathbf{G}_{a}=\mathbf{G}_{\mathrm{spacetime}}+(a^{2}/2)\mathbf{L}_{\mathrm{spacetime}}.$
(19)

If the individual generator matrices $\mathbf{G}^{t_{i}}$, $i=1,\ldots,T$ are
stored in a $T$-vector of $N\times N$ matrices called `Gvector`, the matrix
$\mathbf{G}_{a}$ (denoted `Ga` below) is obtained in julia code with:

⬇

using LinearAlgebra, SparseArrays

Gspacetime = blockdiag(Gvector…) #block-diagonal array of spatial generators

Ltemporal = Tridiagonal(ones(T-1), -2*ones(T), ones(T-1)) #1D central
difference

Ltemporal[1,1] = -1; Ltemporal[T,T] = -1 #handle endpoints

Lspacetime = kron(Ltemporal/h^2, I(N)) #kronecker product with N×N identity

Ga = Gspacetime + (a^2/2)*Lspacetime #combine generator and diffusion
components

####  4.2.1 Selecting the temporal diffusion strength $a$

If $a$ is too small, the lack of significant temporal diffusion means that
distinct time fibres of a spacetime eigenfunction will be poorly linked
together, leading to potentially rapid variation between fibres, which is at
odds with our search for quasi-stationary features. Moreover, temporal
eigenvalues will populate the dominant part of the inflated generator
spectrum, which is unhelpful for distinguishing spatial features. If $a$ is
too large, then only very slight temporal variation between time fibres in a
spacetime eigenfunction will be possible, as small variations will incur high
diffusive mixing in the time coordinate due to strong temporal diffusion. This
forced almost-constancy in time may mask quasi-stationary almost-invariant
behaviour present in the dynamics. In fact, as $a\to\infty$, one expects the
eigenspectrum of $\mathcal{G}$ corresponding to non-temporal eignvalues to
approach the eigenspectrum of the averaged generator $\bar{G}$, introduced at
the end of Section 2.

We now discuss a heuristic for choosing $a$, by matching estimates of leading
temporal and spatial eigenvalues. Recall from the discussion toward the end of
Section 3 that the leading temporal eigenvalue of the inflated generator is
$-a^{2}\pi^{2}/2\tau^{2}$. The real parts of the spatial eigenvalues of
$\mathbf{G}_{a}$ are primarily determined by (i) the spatial diffusion term in
(11) containing $\epsilon$ and (ii) the effective numerical diffusion due to
the Ulam scheme. The former is chosen so that
$\epsilon^{2}/2=0.05\bar{v}\ell$. For the latter, we estimate effective
numerical diffusion, denoted $\epsilon_{\rm eff}$ by $\epsilon_{\rm
eff}^{2}/2=\bar{v}\ell/2$ (see the discussion in Remark 4.6 [18]). We define a
total diffusion $\epsilon_{\rm tot}^{2}/2:=\epsilon^{2}/2+\epsilon_{\rm
eff}^{2}/2=0.55\bar{v}\ell$. Assuming $M$ is a $d$-dimensional hyperrectangle,
the leading nontrivial spatial eigenvalue of $(\epsilon_{\rm
tot}^{2}/2)\Delta$ on $M$ is therefore $-\epsilon_{\rm tot}^{2}\pi^{2}/2L_{\rm
max}^{2}=-0.55\bar{v}\ell\pi^{2}/L_{\max}^{2}$, where $L_{\rm max}$ is the
longest side length of $M$. As discussed in Section 3, it is desirable to
equate diffusion in space and time, so equating these two eigenvalues leads to

$a=\tau\sqrt{1.1\bar{v}\ell}/L_{\rm max}.$ (20)

The value of $a$ can be initialised with this heuristic and adjusted based on
subsequent eigenspectrum computations to better match the leading temporal and
spatial eigenvalues. Related discussions on the parameter $a$ for the inflated
dynamic Laplacian may be found in Section 4.1 [19], where the estimate of the
spatial eigenvalue in [19] arises from a Laplace operator, without any
considerations of a spatial diffusion parameter $\epsilon$.

####  4.2.2 Extracting individual quasi-stationary families with SEBA

Spectral analysis of dynamical systems to determine almost-invariant sets with
transfer operators or generators, and finite-time coherent sets with transfer
operators or dynamic Laplacians, sometimes involves several eigenvectors which
encode several of these sets in a nontrivial way. Individual almost-invariant
or coherent sets need to be disentangled from the leading eigenvectors and a
straightforward way of doing this is with the SEBA algorithm [23]. The input
to SEBA is a collection of eigenvectors and the output is a collection of
sparse vectors (whose supports correspond to individual sets) with
approximately the same span as the input eigenvectors. SEBA is also very
effective in situations where the eigenvectors capture spacetime information,
for example separating individual coherent sets [20] and semi-material
coherent sets [19, 2] across both space and time simultaneously. In our
present setting, if the quasi-stationary family of almost-invariant sets we
seek are not easily identifiable from a _single_ spacetime eigenvector of the
inflated generator, we will apply SEBA to several leading spatial (i.e. non-
temporal) spacetime eigenvectors to extract several individual quasi-
stationary families.

Our overall numerical approach to find quasi-stationary families of almost-
invariant sets from a time-dependent velocity field is summarised below.

######  Algorithm 1 (Identify quasi-stationary families of almost-invariant
sets from vector fields).

Inputs: state space $M$, time-varying velocity field $v$ on $M$, time duration
$[0,\tau]$

Outputs: estimates of families of quasi-stationary almost-invariant sets at
discrete time points in $[0,\tau]$

  1. 1.

Discretise the space space $M$ into cubes $B_{1},\ldots,B_{N}$, and the time
interval $[0,\tau]$ into time nodes $t_{l}$, $l=0,\ldots,T$ as described at
the beginning of Section 4.

  2. 2.

Select a value for $\epsilon$ using equation (17) and compute the discrete
generator matrices $\mathbf{G}^{t_{l}}$ for $l=0,\ldots,T$ as described in
Section 4.1 using (16).

  3. 3.

Combine the $\mathbf{G}^{t_{l}}$ to form $\mathbf{G}_{a}$ as described in
Section 4.2 using (19), where the temporal strength parameter $a$ is initially
selected using (20).

  4. 4.

Compute the real eigenvalues $0=\Lambda_{1}>\Lambda_{2}\geq\cdots$ of the
sparse matrix $\mathbf{G}_{a}^{*}$ that are closest to 0, along with their
corresponding spacetime eigenvectors
$\mathbf{F}_{1}^{*},\mathbf{F}_{2}^{*},\ldots$. Note that $\mathbf{F}_{1}^{*}$
should be a constant vector $\mathbf{1}$.

  5. 5.

Remove the constant-in-space temporal eigenvectors $\mathbf{F}_{k}^{*}$ by
inspection or numerical interrogation; see the discussion at the end of
Section 3 for their functional form. The discretised temporal eigenvectors
will be constant on blocks of indices $lN+1,\ldots,(l+1)N$, $l=0,\ldots,T$.

  6. 6.

If the desired quasi-stationary family of almost-invariant sets is not clearly
apparent in any single eigenvector, but there are signatures of the family
(i.e. large positive and negative values indicate relevant spatial structures)
in an eigenvector $\mathbf{F}_{K}^{*}$ for some $K\geq 2$, apply SEBA to the
$K^{\prime}$ real-valued spatial eigenvectors with index $K$ and below to
obtain spacetime SEBA vectors $\mathbf{S}_{1},\ldots,\mathbf{S}_{K^{\prime}}$,
with $K^{\prime}\leq K$. One of the $\mathbf{S}_{k}$, $k=1,\ldots,K^{\prime}$
should be supported on the desired quasi-stationary family.

##  5 An idealised flow with sharp variation in quasi-stationary structures

In this section we illustrate our constructions on an idealised fluid flow. We
use the “switching double gyre” flow considered in Section 6.1 [2], which is a
modification of the standard double gyre [55].

Figure 2: Some level sets of the stream function for the switching double gyre
flow taken at $t=0$ (left) and $t=1$ (right).

The domain is $M=\left[0,3\right]\times\left[0,2\right]$, and the main almost-
invariant features of the flow are two gyres of unequal size. Initially the
left gyre $A_{0}^{\rm left}$ is half the size of the right gyre $A_{0}^{\rm
right}$, and as the flow progresses from $t=0$ to $t=1$, the sizes change
sharply at times between approximately $t=0.4$ and $t=0.6$, so that at the end
of the flow duration the relative sizes have switched over: $A_{1}^{\rm left}$
is double the size of $A_{1}^{\rm right}$; see Figure 2 for streamline plots
at the initial and final times.

The paper [2] studied the coherence properties of the gyres, attempting to
track the evolution of the gyres in space as their volume changed
dramatically. In the present paper we are instead concerned with determining
the existence of a quasi-stationary family of almost-invariant objects through
time. In contrast to coherent sets, which may freely move about in space while
remaining coherent, a quasi-stationary family of almost-invariant sets are
approximately fixed in space. Thus, the dynamical question we pose is
different, and particularly relevant for e.g. atmospheric blocks, which remain
approximately fixed in space.

For our generator calculations in this experiment, we follow Algorithm 1,
discretising $M$ into a $75\times 50$ grid of square boxes
$B_{1},\ldots,B_{3750}$ of side length $\ell=0.04$. We discretise the time
interval $[0,1]$ into 21 nodes spaced $h=0.05$ time units apart. Following
steps 2 and 3 of Algorithm 1, we calculate $\epsilon\approx 0.0536$ from
$\bar{v}\approx 0.7185$ and $\ell$. We set $a=0.1$ to better match the leading
spatial and temporal eigenvalues, see Figure 3. An initial $a$ value of
approximately 0.0593 was calculated using the heuristic in (20). The sparse
$3750\times 3750$ matrices $\mathbf{G}^{t_{0}},\ldots,\mathbf{G}^{t_{20}}$ are
combined to form the sparse $78750\times 78750$ matrix $\mathbf{G}_{a}$.
Figure 3 shows the leading ten eigenvalues
$0=\Lambda_{1}>\Lambda_{2}>\cdots>\Lambda_{10}$; $\Lambda_{3}$ is a temporal
eigenvalue and is shown as a red cross.

Figure 3: The leading ten eigenvalues of the inflated generator approximation
$\mathbf{G}_{a}$ for the switching double gyre flow.

Figure 4(left) plots the estimates of the second eigenfunction
$\mathbf{F}_{2}^{*}$ at time instances $0,0.1,0.2,\ldots,1$.

Figure 4: Inflated generator construction of quasi-invariant families of
almost-invariant sets in the switching double gyre. _Left:_ Selected time
fibres of the second eigenvector $\mathbf{F}^{*}_{2}$ of $\mathbf{G}_{a}^{*}$.
A cutoff of 0.33 has been applied. _Right:_ Superposition of two SEBA vectors
$\mathbf{S}_{1}$ and $\mathbf{S}_{2}$. A cutoff of 0.2 has been applied.

It will be useful for the reader to refer to Figures 1 and 4 during the
following discussion. At any time $t\in[0,1]$, there is no trajectory escape
from either rectangular gyre $A_{t}^{\rm left}$ or $A_{t}^{\rm right}$ in $M$
under advective from the steady vector field $v(t,\cdot)$. The only
possibility of trajectory escape from $A_{t}^{\rm left}$ or $A_{t}^{\rm
right}$ (restricted to $\\{t\\}\times M$) is by small $\epsilon$-diffusion. On
the other hand, the families $\\{A_{t}^{\rm left}\\}_{t\in[0,1]}$ and
$\\{A_{t}^{\rm right}\\}_{t\in[0,1]}$ undergo a sharp change in spatial
position and area for $t\in[0.4,0.6]$. Therefore if we simply defined
$\mathbb{A}^{\rm left}=\cup_{t\in[0,1]}\\{t\\}\times A_{t}^{\rm
left}\subset\mathbb{M}$ (and defined $\mathbb{A}^{\rm right}$ in a similar
fashion), the resulting spacetime set would have a large loss of trajectories
due to the $a$-diffusion in the time coordinate, along which the set varies
rapidly between times 0.4 and 0.6. This naïve choice of $\mathbb{A}$ would be
optimal in terms of almost-invariance in space, but rather suboptimal in terms
of almost-invariance in time.

Turning to our inflated generator results, in Figure 4(left) we plot time
slices of $\mathbf{F}_{2}^{*}$ across time. One can interpret the red sets as
time fibres of one candidate almost-invariant set $\mathbb{A}^{(1)}$ and the
blue as time fibres of another candidate almost-invariant set
$\mathbb{A}^{(2)}$. Note that parts of $\mathbb{A}^{(1)}$ are cut away
(compared to $\mathbb{A}^{\rm left}$) at times from about $t=0.4$ onwards,
when the velocity field $v(t,\cdot)$ first undergoes rapid change, moving the
separatrix in the vector field to increase the gyre area. This cutting away is
designed to improve the spacetime almost-invariance of the spacetime set
$\mathbb{A}^{(1)}$, and therefore the quasi-stationarity of the family
$\\{A^{(1)}_{t}\\}_{t\in[0,1]}$. This comes at the cost of a slight decrease
in the spatial almost-invariance of $A_{t}^{(1)}$ (compared to $A_{t}^{\rm
left}$) for some $t$, but overall our inflated generator construction
identifies an improved balance between almost-invariance in space and time.
The same interpretation can be given to $\mathbb{A}^{(2)}$, where some parts
are cut away, relative to $\mathbb{A}^{\rm right}$, from time $t=0$ up to
about $t=0.6$.

Inserting $\mathbf{F}^{*}_{1}$ and $\mathbf{F}^{*}_{2}$ into SEBA yields two
SEBA vectors $\mathbf{S}_{1},\mathbf{S}_{2}$, highlighting the individual
spacetime sets $\mathbb{A}^{(1)},\mathbb{A}^{(2)}$, respectively; see Figure
4(right). There is a variation in the depth of colour in the SEBA vectors. The
deepest red indicates those times where one is most certain that points at
those spacetime locations belong to the quasi-stationary family of almost-
invariant sets. We see in $\mathbb{A}^{(1)}$ that the time fibres for $t$ near
the initial time $t=0$ consist of points that strongly belong to this family.
This is because these points belong to a small initial gyre, which will very
likely remain in the family of almost-invariant sets as this family
approximately expands about it in space. On the other hand, points in the
final time fibres of $\mathbb{A}^{(1)}$ (those near $t=1$) are paler and less
certain to remain in the family; points nearer to the gyre “core” are deepest
red and most certain to remain.

##  6 Blocking over Europe

One important practical application of this work involves identifying the
shape and tracing the movement of atmospheric blocking events. We will
identify blocking events by applying the inflated generator method to
atmospheric wind velocity data, and then compare these results to those
obtained from three well-documented block detection methods. We also compare
these results to wind streamlines to identify anticyclonic wind circulation
patterns expected within blocks, which is achieved by seeking closed clockwise
rotating [1] streamline loops.

We use ERA5 atmospheric wind velocity, geopotential height and potential
vorticity reanalysis data obtained from the Copernicus Climate Change Service
Climate Data Store [30]; over the spatial region
$M_{\mathrm{full}}=\left[20^{\circ}\mathrm{W},60^{\circ}\mathrm{E}\right]\times\left[15^{\circ}\mathrm{N},90^{\circ}\mathrm{N}\right]$
and the time period ranging from 00:00 01/01/1979 to 18:00 31/12/2019 (UTC)
inclusive. The potential vorticity data is taken on all available vertical
pressure levels ranging from 150 to 500 hPa (inclusive), while the wind
velocity and geopotential height data are taken at the 500 hPa pressure level
only. The spatial resolution of this data is $0.25^{\circ}$ in both longitude
and latitude, and the temporal resolution for the data is 6 hours. We linearly
interpolate velocities between the spatial data points with julia’s
`Interpolations.jl` package [34]. We use the full data set for calculations
involving geopotential height and potential vorticity, but for our inflated
generator calculations we do not use the entire spatial domain nor the full
spatial resolution available.

###  6.1 Existing Block Detection Methods

We begin by briefly describing three existing geophysical block detection
methods. The reader is referred to the cited literature for more information
on each method.

##### Geopotential height gradients.

Gradients of $\mathrm{Z}_{500}$, which denotes geopotential height taken at a
500 hPa pressure level. While this method was originally used to identify
blocked longitude levels through changes in the $\mathrm{Z}_{500}$ gradient
over a single line of latitude [36, 60], we instead consider a
“two–dimensional” version of this method [50, 6, 44, 59] used to identify
blocked regions at a particular moment in time. Let $\lambda$ and $\phi$
represent longitude and latitude respectively. Using “northern” and “southern”
$\mathrm{Z}_{500}$ gradients, as defined in [50] and denoted here by
$\gamma_{N}$ and $\gamma_{S}$ respectively, a closed set $\Gamma_{t}\subset M$
at an instantaneous moment in time $t$, where

$\Gamma_{t}=\left\\{\left(\lambda,\phi\right)\in M\ :\
\gamma_{N}\left(t,\lambda,\phi\right)<-10,\
\gamma_{S}\left(t,\lambda,\phi\right)>0\right\\}$

is classified as a blocked region.

##### Geopotential height anomalies.

This method also utilises $\mathrm{Z}_{500}$ through the computation of the
long–term anomalies of this quantity, which we denote by
$\mathrm{Z}_{500}^{*}$. We compute $\mathrm{Z}_{500}^{*}$ using the procedure
detailed in [49], and identify blocking events as closed contours of
$\mathrm{Z}_{500}^{*}$ that are above a threshold of +100 geopotential metres
(gpm) [11, 8, 44, 59].

##### Vertically averaged potential vorticity anomalies.

Anomalies of potential vorticity averaged over a range of vertical pressure
levels, $\overline{\mathrm{VPV}}$, which we denote by
$\overline{\mathrm{VPV}}^{*}$. These anomalies are computed using the
methodology described in [44], with the omission of the Fourier transform
harmonics step used to smoothen the “Long Term Daily Means” of
$\overline{\mathrm{VPV}}$. We computed $\overline{\mathrm{VPV}}^{*}$ without
this extra step and found the results to be very similar to those in [59]. We
then identify blocks similarly to the method of the previous paragraph, except
we seek contours of $\overline{\mathrm{VPV}}^{*}$ below a threshold of -1.2
potential vorticity units [53, 44, 59].

Table 1 lists two blocking events occurring over Europe between July and
August 2003 identified in [59] using at least one of these three methods and a
“Ground Truth Dataset” used to verify all identified blocks. Rows 2, 3 and 4
of Figures 7 and 10 show instantaneous snapshots of regions $\Gamma_{t}$
satisfying the $\mathrm{Z}_{500}$ gradient criteria and scalar fields of
$\mathrm{Z}_{500}^{*}$ and $\overline{\mathrm{VPV}}^{*}$ respectively, which
visually identify blocks. The blocking criteria described above should be
maintained for a fixed amount of time (usually between four [60, 42, 13] to
ten [46, 11] days) for a region to be considered blocked, but in practice this
condition is often relaxed.

Block | Date Range | Long. Range | Lat. Range | Identified Using
---|---|---|---|---
East | 29/07/03 – 03/08/03 |  $25^{\circ}$ E – $55^{\circ}$ E |  $55^{\circ}$ N – $70^{\circ}$ N |  $\Gamma_{t}$, $\mathrm{Z}_{500}^{*}$, $\overline{\mathrm{VPV}}^{*}$
West | 02/08/03 – 12/08/03 |  $10^{\circ}$ W – $15^{\circ}$ E |  $40^{\circ}$ N – $60^{\circ}$ N | $\mathrm{Z}_{500}^{*}$
Table 1: The “east” and “west” blocking events as identified in [59]. The
longitudinal and latitudinal ranges of these blocks are estimated from Figure
5 of [59].

###  6.2 Results of block identification using the inflated generator

We wish to apply the inflated generator method to atmospheric velocity data to
identify the onset, evolution, and dissipation of the two blocking events
listed in Table 1, documented [13, 59] as having played a key role in the 2003
European heatwaves. We will take $M$ to be a fixed rectangle restricted to a
surface corresponding to a pressure level of 500 hPa [30] and $v$ to be the
time-varying wind field on this surface. The vector field $v$ need not have
zero normal velocity at the boundary of $M$. Following the approach in section
4.1 we discretise the 500 hPa pressure level into coordinates in degrees
longitude and latitude, which results in rectangular boxes
$B_{1},\ldots,B_{N}$ with areas that depend on the latitude. We denote the
longitudinal and latitudinal extent (in metres) of box $B_{i}$ by
$\ell_{i,{\rm lon}}$ and $\ell_{i,{\rm lat}}$, respectively. Equation (16)
becomes:

$\mathbf{G}^{t_{l}}_{ij}=\begin{cases}\frac{1}{m(B_{j})}\left(\int_{B_{i}\cap
B_{j}}(\max\\{v(t_{l},x)\cdot n_{ij},0\\})\
dm_{1}(x)\right)+\epsilon^{2}/(2\ell_{i,{\rm lon}}^{2}),&\quad i\neq
j,\mbox{$n_{ij}$ is longitudinal}\\\ \frac{1}{m(B_{j})}\left(\int_{B_{i}\cap
B_{j}}(\max\\{v(t_{l},x)\cdot n_{ij},0\\})\
dm_{1}(x)\right)+\epsilon^{2}/(2\ell_{i,{\rm lat}}^{2}),&\quad i\neq
j,\mbox{$n_{ij}$ is latitudinal}\\\ \sum_{j\neq
i}-\frac{m(B_{j})}{m(B_{i})}\mathbf{G}^{t_{l}}_{ij},&\quad\mbox{otherwise.}\end{cases}$
(21)

Because our box sizes vary with latitude, we denote the median box side length
(computed over longitudinal and latitudinal directions) by $\bar{\ell}$.
Analogously to the choice in Section 4.1 the small spatial diffusion parameter
$\epsilon$ is set to $\epsilon=\sqrt{0.1\bar{v}\bar{\ell}}$. The initial guess
for the temporal diffusion strength parameter $a$ is correspondingly modified
from the discussion in Section 4.2.1; we set
$a=\tau\sqrt{1.1\bar{v}\bar{\ell}}/L_{\rm max}$.

####  6.2.1 East Block (July-August 2003)

The first blocking event we consider occurred in Eastern Europe between 00:00
29/07/03 to 00:00 03/08/03 (UTC). To identify this block, we let
$M=\left[15^{\circ}\mathrm{E},60^{\circ}\mathrm{E}\right]\times\left[30^{\circ}\mathrm{N},75^{\circ}\mathrm{N}\right]$.
The temporal domain includes the full duration of the block with a three day
extension applied to either side of this time span; giving an interval ranging
from 00:00 26/07/03 to 00:00 06/08/03 (UTC) (inclusive) of length $\tau=11$
days. We use the full temporal resolution of 6 hours, and a spatial resolution
of $1^{\circ}$ in both latitude and longitude; providing a spatial grid of
$45\times 45$ boxes.

Following the instructions described in Section 4.2.1, we calculate
$\epsilon\approx 319.2460\ \mathrm{ms}^{-1/2}$, given a median box side length
$\bar{\ell}\approx 103618\ \mathrm{m}$ and median speed $\bar{v}\approx
9.8360\ \mathrm{ms^{-1}}$, and let $a=0.0032$; our heuristic yielded
$a=0.0023$. The leading ten eigenvalues for the inflated generator are shown
in Figure 5.

Figure 5: The leading ten eigenvalues $\Lambda_{k}$ (ordered according to
their real parts in decreasing order) of the inflated generator calculated for
the European atmospheric flow (eastern region) over the date range
26/07/03-06/08/03, plotted in the complex plane. Spatial and temporal
eigenvalues are distinguished using separate plot markers and colours.

We then visually inspect each of the real-valued spatial eigenvectors
individually to identify which of these generates a quasi-stationary object
that best resembles the east block. The deep red elliptical objects identified
from the ninth eigenvector $\mathbf{F}_{9}^{*}$ shown in Figure 6 (upper row)
were found to provide the best visualisation of this block. By comparison, in
results not shown, $\mathbf{F}_{2}^{*}$ identifies a combination of the east
block and another quasi-stationary object positioned further east; while
$\mathbf{F}_{4}^{*}$ and $\mathbf{F}_{5}^{*}$ each identify a cyclonic object
nested underneath the east block. The latter object is the cut-off low
expected to appear underneath a block [56]. We then apply the SEBA algorithm
to all real-valued spatial eigenvectors from $\mathbf{F}_{1}^{*}$ up to and
including $\mathbf{F}_{9}^{*}$ (namely
$\mathbf{F}_{1}^{*},\mathbf{F}_{2}^{*},\mathbf{F}_{4}^{*},\mathbf{F}_{5}^{*},\mathbf{F}_{9}^{*}$),
and then inspect the five SEBA vectors to identify which of these gives the
best visualisation of the east block. The third SEBA vector $\mathbf{S}_{3}$,
shown in Figure 6 (lower row) isolates this block.

Figure 6: Temporal evolution of the spatial eigenvector $\mathbf{F}_{9}^{*}$
and the SEBA vector $\mathbf{S}_{3}$ of the inflated generator for the
European atmospheric flow (eastern region) at select dates taken over the date
range 26/07/03-06/08/03. The time of day is 00:00 (UTC) for all of the above
dates.

Figure 6 (lower row) shows that the east block gradually begins to form around
July 26, with the block taking its expected shape just after July 29. The
block is maintained until August 3, by which time the block has reduced to a
smaller elliptical object. No trace of the block can be identified from
$\mathbf{F}_{9}^{*}$ or $\mathbf{S}_{3}$ beyond August 3.

We now compare $\mathbf{S}_{3}$ to the geophysical block diagnostic methods
listed in Section 6.1, as well as the wind streamline field. Figure 7 shows
time slices of $\mathbf{S}_{3}$ and all of the other diagnostic quantities
taken at the same time instances as those used in Figure 6.

Figure 7: Several block diagnostic tools for the European atmospheric flow
(eastern region) taken at specific dates over the date range
26/07/03-06/08/03. Included are (from top to bottom): the wind speed scalar
field and corresponding streamlines; blocked regions $\Gamma_{t}$ identified
using the $\mathrm{Z}_{500}$ gradient, the $\mathrm{Z}_{500}^{*}$ scalar
field, the $\overline{\mathrm{VPV}}^{*}$ scalar field, and the third SEBA
vector $\mathbf{S}_{3}$ as shown in Figure 6. The time of day is 00:00 (UTC)
for all of the above dates. Blocks are identified through medium red coloured
objects in the $\Gamma_{t}$ and $\mathrm{Z}_{500}^{*}$ fields; deep red
objects in the $\mathbf{S}_{3}$ field; medium blue objects in the
$\overline{\mathrm{VPV}}^{*}$ scalar field, and closed clockwise rotating
streamline loops in the wind field.

All of the other block identification methods are able to identify the east
block to some extent. This is evidenced by closed contours of clockwise
oriented wind streamlines [1], blocked regions $\Gamma_{t}$ obtained from the
$\mathrm{Z}_{500}$ gradient, and objects of strongly positive
$\mathrm{Z}_{500}^{*}$ and negative $\overline{\mathrm{VPV}}^{*}$. Regarding
false positives, the $\mathrm{Z}_{500}$ gradient method shows a block on July
26–27 that is not indicated by any of the other methods. After August 3, a
block cannot be identified using $\mathbf{S}_{3}$. In contrast, the results
from each of the other methods give the impression that after August 3 the
block moves Eastwards until August 6 where the block is roughly centred at
$55^{\circ}$ E. This Eastward movement is not consistent with the quasi-
stationary nature of blocking events [45, 46, 50] and we regard these temporal
extensions of the block arising from the non-generator techniques as false
positives.

We emphasise that the other blocking detection methods use instantaneous
geophysical proxies to diagnose a block at individual times, and there is no
explicit temporal connection between the diagnoses at distinct times. This
disconnection in time not only ignores the important quasi-stationary aspect
of blocks, but can make it difficult to identify separate blocking events and
to piece together block identification across time, including when the block
appears and disappears. In contrast, the inflated generator approach
explicitly connects different time instances and enforces quasi-stationarity
of the identified blocks, while the SEBA algorithm isolates individual blocks.
The inflated generator approach requires some prior knowledge of the lifespan
and spatial extent of a block to set up the spatial and temporal domain, but
as one can see from Figure 7, closed, clockwise rotating wind streamlines can
provide this rough initial information.

####  6.2.2 West Block (August 2003)

Next we apply the inflated generator method to a block that occurred further
west in Europe from 00:00 02/08/03 to 00:00 12/08/03 (UTC) (inclusive). This
significant blocking event was responsible for some of the worst extreme heat
recorded in Western Europe (particularly in France) during the summer of 2003
[5, 59]. For this longer-duration block, we alter the spatial domain $M$
slightly to
$\left[10^{\circ}\mathrm{W},40^{\circ}\mathrm{E}\right]\times\left[30^{\circ}\mathrm{N},75^{\circ}\mathrm{N}\right]$,
and define a time interval covering the full duration of the block as well as
five days immediately preceding and following this block, ranging from 00:00
28/07/03 to 00:00 17/08/03 (UTC) (inclusive), yielding $\tau=20$ days. We
maintain a temporal resolution of 6 hours and a spatial resolution of
$1^{\circ}$ in both latitude and longitude, providing a spatial grid of
$50\times 45$ boxes.

We calculate $\epsilon\approx 327.5440\ \mathrm{ms}^{-1/2}$ from a median box
length $\bar{\ell}\approx 103618\ \mathrm{m}$ and median speed $\bar{v}\approx
10.3539\ \mathrm{ms^{-1}}$, and let $a=0.0045$; our heuristic yielded
$0.0043$. The first eleven eigenvalues of the inflated generator spectrum are
in Figure 8.

Figure 8: The leading eleven eigenvalues $\Lambda_{k}$ (ordered according to
their real parts in decreasing order) of the inflated generator calculated for
the European atmospheric flow (western region) over the date range
28/07/03-17/08/03, plotted in the complex plane. Spatial and temporal
eigenvalues are distinguished using separate plot markers and colours.

Similarly to the east block experiment, we inspect each non-trivial, real-
valued spatial eigenvector until we find one that provides the best
visualisation of the west block. In this example we find that the eighth
eigenvector $\mathbf{F}_{8}^{*}$ (see Figure 9 (upper row)) best isolates the
west block as a separate quasi-stationary object covering the expected spatial
and temporal extent of the block. The block is indicated by the deep red
elliptical objects in the $\mathbf{F}_{8}^{*}$ fields, which form around
August 4 and are positioned over France between August 6-12, the period over
which France’s most extreme heat was recorded [5, 59]. The block is still
present on August 12 and begins to disappear after this date before almost
completely vanishing after August 14. To isolate this block from all other
objects in the $\mathbf{F}_{8}^{*}$ field, we apply SEBA to all real-valued
spatial eigenvectors from $\mathbf{F}_{1}^{*}$ up to and including
$\mathbf{F}_{8}^{*}$. The first SEBA vector $\mathbf{S}_{1}$ isolates the west
block and is plotted in Figure 9 (lower row).

Figure 9: Temporal evolution of the spatial eigenvector $\mathbf{F}_{8}^{*}$
and the SEBA vector $\mathbf{S}_{1}$ of the inflated generator for the
European atmospheric flow (western region) at select dates taken over the date
range 28/07/03-17/08/03. The time of day is 00:00 (UTC) for all of the above
dates.

We now compare $\mathbf{S}_{1}$ to the other block diagnostics in Figure 10.

Figure 10: Several block diagnostic tools for the European atmospheric flow
(western region) taken at select dates over the date range 28/07/03-17/08/03.
Included are (from top to bottom): the wind speed scalar field and
corresponding streamlines; blocked regions $\Gamma_{t}$ identified using the
$\mathrm{Z}_{500}$ gradient, the $\mathrm{Z}_{500}^{*}$ scalar field, the
$\overline{\mathrm{VPV}}^{*}$ scalar field, and the first SEBA vector
$\mathbf{S}_{1}$ as shown in Figure 9. The time of day is 00:00 (UTC) for all
of the above dates. Blocks are identified through medium red coloured objects
in the $\Gamma_{t}$ and $\mathrm{Z}_{500}^{*}$ fields; deep red objects in the
$\mathbf{S}_{1}$ field; medium blue objects in the
$\overline{\mathrm{VPV}}^{*}$ scalar field, and closed clockwise rotating
streamline loops in the wind field.

The west block is considerably harder to identify using these geophysical
proxies than the east block. The $\mathrm{Z}_{500}^{*}$ and
$\overline{\mathrm{VPV}}^{*}$ fields roughly identify a similar object to that
identified using $\mathbf{S}_{1}$, though this object is centred over the
United Kingdom and Ireland rather than France. The $\mathrm{Z}_{500}$ gradient
struggles to identify the west block, initially identifying it on August 3,
but only managing to identify random fragments of it from August 7 onwards.
There appears to be no consistent agreement between the three geophysical
block detection methods regarding the location, size and duration of the west
block, and none were able to provide a clear identification of the west block.

In contrast, the deep red elliptical object found in $\mathbf{S}_{1}$ strongly
resembles the west block, with the largest values of $\mathbf{S}_{1}$
pinpointing the presence of a roughly elliptical structure over France during
the expected period of the west block. The paler red periphery of this object
extends beyond the expected range of the block (particularly south west into
Northwest Africa). This extension of the west block was more strongly
identifiable in all real-valued spatial eigenvectors preceding
$\mathbf{F}_{8}^{*}$, suggesting that another blocking event was present over
Northwest Africa between July 28 and August 3. Only the $\mathrm{Z}_{500}$
gradient on July 28 suggests the presence of this other block.

##  7 Conclusion

We have developed a robust methodology for the identification of time-
parameterised regions in phase space that are (i) almost-invariant or
metastable at each time instant and (ii) vary slowly in space across time (are
quasi-stationary). Our approach also determined the timing of the onset and
disappearance of such quasi-stationary families of almost-invariant sets from
a single eigencomputation of an inflated generator. In addition to introducing
the concept of the inflated generator, we developed theory to underpin the
almost-invariance and quasi-stationarity properties of the eigenvectors of the
inflated generator, and proposed an efficient computational approach that
requires only an estimate of the underlying time-varying velocity field (in
particular, no trajectory integration is needed). We illustrated our new
approach by successfully identifying two atmospheric blocking events in the
European 2003 summer, and compared our results with common blocking diagnosis
techniques based on geophysical proxies. We found that when compared to a
ground truth, our inflated generator results were considerably more accurate
than all geophysical proxies. Advantages of our approach over existing
geophysical methods for atmospheric block diagnosis include (i) the direct use
of wind fields in assessing almost-invariance, rather than less direct
proxies, and (ii) the explicit requirement that the diagnosed family of
almost-invariant sets is quasi-stationary. We hope that our spacetime
methodology will aid a better understanding of the physics behind critical
blocking phenomena through its complete lifecycle description of these events.

##  8 Acknowledgements

AB is supported by an ARC Discovery Project and GF is partially supported by
an ARC Discovery project. The authors acknowledge the use of UNSW’s Katana
compute cluster for some of the computations.

##  9 Code

The code used to produce the results featured in this article can be obtained
from the Inflated-Generator repository on GitHub, available at:
https://github.com/gfroyland/Inflated-Generator

## References

  * [1] anticyclone.  https://glossary.ametsoc.org/wiki/Anticyclone, 2024.  American Meteorological Society Glossary of Meteorology.
  * [2] J. Atnip, G. Froyland, and P. Koltai.  An inflated dynamic Laplacian to track the emergence and disappearance of semi-material coherent sets.  https://arxiv.org/abs/2403.10360, 2024.
  * [3] E. A. Barnes, J. Slingo, and T. Woollings.  A methodology for the comparison of blocking climatologies across indices, models and climate scenarios.  Climate Dynamics, 38(11–12):2467–2481, 2012.
  * [4] D. Barriopedro, R. García-Herrera, and R. Trigo.  Application of blocking diagnosis methods to General Circulation Models. Part I: a novel detection scheme.  Climate Dynamics, 35(1):1373–1391, 2010.
  * [5] E. Black, M. Blackburn, G. Harrison, B. Hoskins, and J. Methven.  Factors contributing to the summer 2003 European heatwave.  Weather (Royal Meteorological Society), 59(8):217–223, 2004.
  * [6] T. Buehler, C. C. Raible, and T. F. Stocker.  The relationship of winter season North Atlantic blocking frequencies to extreme cold or dry spells in the ERA-40.  Tellus A: Dynamic Meterology and Oceanography, 63(2):212–222, 2011.
  * [7] W. Cai, X. Xu, X. Cheng, F. Wei, X. Qiu, and W. Zhu.  Impact of “blocking” structure in the troposphere on the wintertime persistent heavy air pollution in Northern China.  Science of the Total Environment, 741(1):140325, 2020.
  * [8] M. Carrera, R. Higgins, and V. Kousky.  Downstream Weather Impacts Associated with Atmospheric Blocking over the Northeast Pacific.  Journal of Climate, 17(24):4823–4839, 2004.
  * [9] M. Dellnitz and O. Junge.  On the Approximation of Complicated Dynamical Behavior.  SIAM Journal on Numerical Analysis, 36(2):491–515, 1999.
  * [10] P. Deuflhard and M. Weber.  Robust Perron cluster analysis in conformation dynamics.  Linear Algebra and its Applications, 398(1):161–184, 2005.
  * [11] R. M. Dole and N. D. Gordon.  Persistent Anomalies of the Extratropical Northern Hemisphere Wintertime Circulation: Geographical Distribution and Regional Persistence Characteristics.  Monthly Weather Review, 111(8):1567–1586, 1983.
  * [12] E. Dunn-Sigouin and S.-W. Son.  Northern Hemisphere blocking frequency and duration in the CMIP5 models.  Journal of Geophysical Research: Atmospheres, 118(3):1179–1188, 2013.
  * [13] N. Ehstand, R. V. Donner, C. López, and E. Hernández-García.  Characteristic signatures of Northern Hemisphere blocking events in a Lagrangian flow network representation of the atmospheric circulation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9):093128, 2021.
  * [14] D. Faranda, G. Masato, N. Moloney, Y. Sato, F. Daviaud, B. Dubrulle, and P. Yiou.  The switching between zonal and blocked mid-latitude atmospheric circulation: a dynamical system perspective.  Climate Dynamics, 47(5–6):1587–1599, 2015.
  * [15] G. Froyland.  Statistically optimal almost–invariant sets.  Physica D: Nonlinear Phenomena, 200(3–4):205–219, 2005.
  * [16] G. Froyland.  An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems.  Physica D: Nonlinear Phenomena, 250(1):1–19, 2013.
  * [17] G. Froyland and M. Dellnitz.  Detecting and locating near-optimal almost-invariant sets and cycles.  SIAM Journal on Scientific Computing, 24(6):1839–1863, 2003.
  * [18] G. Froyland, O. Junge, and P. Koltai.  Estimating long-term behavior of flows without trajectory integration: The infinitesimal generator approach.  SIAM Journal on Numerical Analysis, 51(1):223–247, 2013.
  * [19] G. Froyland and P. Koltai.  Detecting the birth and death of finite-time coherent sets.  Communications on Pure and Applied Mathematics, 76(12):3642–3684, 2023.
  * [20] G. Froyland, P. Koltai, and M. Stahn.  Computation and Optimal Perturbation of Finite-Time Coherent Sets for Aperiodic Flows Without Trajectory Integration.  SIAM Journal on Applied Dynamical Systems, 19(3):1659–1700, 2020.
  * [21] G. Froyland, S. Lloyd, and N. Santitissadeekorn.  Coherent sets for nonautonomous dynamical systems.  Physica D: Nonlinear Phenomena, 239(16):1527–1541, 2010.
  * [22] G. Froyland and K. Padberg.  Almost–invariant sets and invariant manifolds — Connecting probabilistic and geometric descriptions of coherent structures in flows.  Physica D: Nonlinear Phenomena, 238(16):1507–1523, 2009.
  * [23] G. Froyland, C. P. Rock, and K. Sakellariou.  Sparse eigenbasis approximation: Multiple feature extraction across spatiotemporal scales with application to coherent set identification.  Communications in Nonlinear Science and Numerical Simulation, 77(1):81–107, 2019.
  * [24] G. Froyland, N. Santitissadeekorn, and A. Monahan.  Transport in time-dependent dynamical systems: Finite–time coherent sets.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 20(4):043116, 2010.
  * [25] G. Froyland, R. M. Stuart, and E. van Sebille.  How well-connected is the surface of the global ocean?  Chaos: An Interdisciplinary Journal of Nonlinear Science, 24(3):033126, 2014.
  * [26] G. Haller and F. Beron-Vera.  Coherent Lagrangian vortices: the black holes of turbulence.  Journal of Fluid Mechanics, 731(1):R4, 2013.
  * [27] G. Haller and A. Poje.  Finite time transport in aperiodic flows.  Physica D: Nonlinear Phenomena, 119(3–4):352–380, 1998.
  * [28] G. Haller and G. Yuan.  Lagrangian coherent structures and mixing in two–-dimensional turbulence.  Physica D: Nonlinear Phenomena, 147(3–4):352–370, 2000.
  * [29] S. Hauser, F. Teubler, M. Riemer, P. Knippertz, and C. M. Grams.  Towards a holistic understanding of blocked regime dynamics through a combination of complementary diagnostic perspectives.  Weather and Climate Dynamics, 4(2):399–425, 2023.
  * [30] H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut.  ERA5 hourly data on pressure levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS).  https://doi.org/10.24381/cds.bd0915c6, 2023.
  * [31] M. Hirt, L. Schielicke, A. Müller, and P. Névir.  Statistics and dynamics of blockings with a point vortex model.  Tellus A: Dynamic Meteorology and Oceanography, 70(1):1458565, 2018.
  * [32] S. G. Johnson.  QuadGK.jl: Gauss–Kronrod integration in Julia: v2.9.4.  https://github.com/JuliaMath/QuadGK.jl, 2013.
  * [33] L.-A. Kautz, O. Martius, S. Pfahl, J. G. Pinto, A. M. Ramos, P. M. Sousa, and T. Woollings.  Atmospheric blocking and weather extremes over the Euro-Atlantic sector – a review.  Weather and Climate Dynamics, 3(1):305–336, 2022.
  * [34] M. Kittisopikul, T. E. Holy, and T. Aschan.  JuliaMath/Interpolations.jl: v0.14.7.  https://github.com/JuliaMath/Interpolations.jl, 2022.
  * [35] P. Koltai.  A stochastic approach for computing the domain of attraction without trajectory simulation.  American Institute of Mathematical Sciences Conference Publications, 2011(Special):854–863, 2011.
  * [36] H. Lejenäs and H. Økland.  Characteristics of northern hemisphere blocking as determined from a long time series of observational data.  Tellus A: Dynamic Meterology and Oceanography, 35(5):350–362, 1983.
  * [37] Q. Liu.  On the definition and persistence of blocking.  Tellus A: Dynamic Meterology and Oceanography, 46(3):286–298, 1994.
  * [38] T. Ma and E. M. Bollt.  Differential geometry perspective of shape coherence and curvature evolution by finite-time nonhyperbolic splitting.  SIAM Journal on Applied Dynamical Systems, 13(3):1106–1136, 2014.
  * [39] R. C. McOwen.  Partial Differential Equations: Methods and Applications.  Pearson, London, first edition, 1995.
  * [40] N. Nakamura and C. S. Y. Huang.  Atmospheric blocking as a traffic jam in the jet stream.  Science, 361(6397):42–47, 2018.
  * [41] T.-W. Park, C.-H. Ho, and S. Yang.  Relationship between the Arctic Oscillation and Cold Surges over East Asia.  Journal of Climate, 24(1):68–83, 2011.
  * [42] J. L. Pelly and B. J. Hoskins.  A New Perspective on Blocking.  Journal of the Atmospheric Sciences, 60(5):743–755, 2003.
  * [43] R. T. Pierrehumbert and H. Yang.  Global chaotic mixing on isentropic surfaces.  Journal of the Atmospheric Sciences, 50(15):2462–2480, 1993.
  * [44] M. C. Pinheiro, P. A. Ullrich, and R. Grotjahn.  Atmospheric blocking and intercomparison of objective detection methods: flow field characteristics.  Climate Dynamics, 53(7):4189–4216, 2019.
  * [45] D. F. Rex.  Blocking Action in the Middle Troposphere and its Effect upon Regional Climate Part I: An Aerological Study of Blocking Action.  Tellus A: Dynamic Meterology and Oceanography, 2(3):196–211, 1950.
  * [46] D. F. Rex.  Blocking Action in the Middle Troposphere and its Effect upon Regional Climate Part II: The Climatology of Blocking Action.  Tellus A: Dynamic Meterology and Oceanography, 2(4):237–324, 1950.
  * [47] V. Rom-Kedar, A. Leonard, and S. Wiggins.  An analytical study of transport, mixing and chaos in an unsteady vortical flow.  Journal of Fluid Mechanics, 214(1):347–394, 1990.
  * [48] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson.  Spectral analysis of nonlinear flows.  Journal of Fluid Mechanics, 641(1):115–127, 2009.
  * [49] R. Sausen, W. König, and F. Sielmann.  Analysis of blocking events from observations and ECHAM model simulations.  Tellus A: Dynamic Meterology and Oceanography, 47(4):421–438, 1995.
  * [50] S. C. Scherrer, M. Croci-Maspoli, C. Schwierz, and C. Appenzeller.  Two-dimensional indices of atmospheric blocking and their statistical relationship with winter climate patterns in the Euro-Atlantic region.  International Journal of Climatology, 26(2):233–249, 2006.
  * [51] S. Schubert and V. Lucarini.  Dynamical analysis of blocking events: spatial and temporal fluctuations of covariant Lyapunov vectors.  Quarterly Journal of the Royal Meteorological Society, 142(698):2143–2158, 2016.
  * [52] C. Schütte.  Conformational Dynamics: Modelling, Theory, Algorithm, and Application to Biomolecules.  Habilitation thesis, Zuse Institute, Freie Universität Berlin, July 1999.
  * [53] C. Schwierz, M. Croci-Maspoli, and H. C. Davies.  Perspicacious indicators of atmospheric blocking.  Geophysical Research Letters, 31(6):L06125, 2004.
  * [54] E. Ser-Giacomi, R. Vasile, I. Recuerda, E. Hernández-García, and C. López.  Dominant transport pathways in an atmospheric blocking event.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(8):087413, 2015.
  * [55] S. C. Shadden, F. Lekien, and J. E. Marsden.  Definition and properties of Lagrangian coherent structures from finite–time Lyapunov exponents in two–dimensional aperiodic flows.  Physica D: Nonlinear Phenomena, 212(3–4):271–304, 2005.
  * [56] G. J. Shutts.  The propagation of eddies in diffluent jetstreams: Eddy vorticity forcing of ‘blocking’ flow fields.  Quarterly Journal of the Royal Meteorological Society, 109(462):737–761, 1983.
  * [57] D. Steinfeld and S. Pfahl.  The role of latent heating in atmospheric blocking dynamics: a global climatology.  Climate Dynamics, 53(9–10):6159–6180, 2019.
  * [58] A. Tantet, F. R. van der Burgt, and H. A. Dijkstra.  An early warning indicator for atmospheric blocking events using transfer operators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(3):036406, 2015.
  * [59] C. Thomas, A. Voulgarakis, G. Lim, J. Haigh, and P. Nowack.  An unsupervised learning approach to identifying blocking events: the case of European summer.  Weather and Climate Dynamics, 2(3):581–608, 2021.
  * [60] S. Tibaldi and F. Molteni.  On the operational predictability of blocking.  Tellus A: Dynamic Meterology and Oceanography, 42(3):343–365, 1990.
  * [61] S. Tibaldi, P. Ruti, E. Tosi, and M. Maruca.  Operational predictability of winter blocking at ECMWF: an update.  Annales Geophysicae, 13(1):305–317, 1995.
  * [62] J. von Lindheim, A. Harikrishnan, T. Dörffel, R. Klein, P. Koltai, N. Mikula, A. Müller, P. Névir, G. Pacey, R. Polzin, and N. Vercauteren.  Definition, detection, and tracking of persistent structures in atmospheric flows.  Journal of Mathematical Physics, pages 1–64, 2021.
  * [63] S. Yun and C. Yoo.  The Effects of Spring and Winter Blocking on PM10 Concentration in Korea.  Atmosphere, 10(7):410, 2019.