# On the archetypal ‘flavours’, indices and teleconnections of ENSO revealed
by global sea surface temperatures

###### Abstract

El Niño-Southern Oscillation global (ENSO) imprint on sea surface temperature
comes in many guises. To identify its tropical fingerprints and impacts on the
rest of the climate system, we propose a global approach based on archetypal
analysis (AA), a pattern recognition method based on the identification of
extreme configurations in the dataset under investigation. Relying on
detrended sea surface temperature monthly anomalies over the 1982 to 2022
period, the technique recovers central and eastern Pacific ENSO types
identified by more traditional methods and allows one to hierarchically add
extra flavours and nuances to both persistent and transient phases of the
phenomenon. Archetypal patterns found compare favorably to phase
identification from K-means, fuzzy C-means and recently published network-
based machine-learning algorithms. The AA implementation is modified for the
identification of ENSO phases in sub-seasonal-to-seasonal prediction systems
and complements current alert systems in characterising the diversity of ENSO
and its teleconnections. Tropical and extra-tropical teleconnection composites
from various oceanic and atmospheric fields derived from the analysis are
shown to be robust and physically relevant. Extending AA to sub-surface ocean
fields improves the discrimination between phases when the characterisation of
ENSO based on sea surface temperature is uncertain. We show that AA on
detrended sea-level monthly anomalies provides a clearer expression of ENSO
types.

##  1 Significance statement

ENSO comes in many guises and is challenging to classify in a changing
climate. This issue is becoming critical as operational centers recognize the
shortcomings of simple two category ENSO systems. Archetypal analysis, a
pattern recognition technique, is applied to sea surface temperature and sea-
level anomalies to provide a robust characterisation of more flavours of ENSO.
The additional flavours are related to the state of the Pacific Ocean. The
inclusion of sea-level anomalies resolves issues in the classification of ENSO
when sea surface temperature alone is inconclusive. The ability of the method
to recover previously detected ENSO flavours and to extend its application to
the classification of nuances in the onset, main phase, and decay of the
phenomenon, is assessed.

##  2 Introduction

Operational centers have long realised that their climate driver alert systems
in general, and ENSO alert protocols in particular, may suffer from non-
stationary processes, driven, not only by anthropogenic climate change, but
also by long-period internal climate modes and their likely interactions
(Nguyen et al., 2022). The advice given by operational centers mainly relies
on simple indices derived from various oceanic and atmospheric anomaly fields
reaching and exceeding critical thresholds combined with a persistence
criterion. This advice may also be conditional on multivariate alert
mechanisms derived from additional atmospheric and oceanic diagnostics
reviewed by a panel of experts111For example,
https://www.cpc.ncep.noaa.gov/products/analysis_monitoring/enso_advisory/ensodisc.pdf
and http://www.bom.gov.au/climate/enso/outlook/#tabs=About-ENSO-and-the-
Outlook. For ENSO alert and prediction systems, operational centers typically
issue an all-encompassing type of advice for La Niña, La Niña watch, Neutral,
El Niño watch, and El Niño phases. A detailed inspection of ENSO global
imprint on sea surface temperature (SST) alone over the 1982-2022 period
reveals a diversity of expressions of ENSO and its teleconnections with the
rest of the climate system and their environmental impact, unlikely to be
adequately described by a single index. The non-stationarity of SST at long-
time scales also interferes with alert protocols rendering them prone to
mischaracterisation (L’Heureux et al., 2013; Timmermann et al., 2018;
Turkington et al., 2019; Oldenborgh et al., 2021).

Hereafter, we propose an approach based on archetypal analysis (AA) adding
more flavours222Here, we, the authors, understand flavour in the ENSO context
as a range of clearly differentiated (arche)typical patterns manifest in the
tropical SST anomalies with well defined centers of action or loci of absolute
maximum amplitudes, which expression may have a distinct dynamical origin. to
the reclassification of historical ENSO events (Hannachi and Trendafilov,
2017; Black et al., 2022) and readily applicable to the identification of
future ENSO phases in subseasonal-to-seasonal (S2S) predictions. The method
allows the unsupervised detection of Eastern and Central Pacific (EP and CP)
ENSO phases amongst others, and can easily be extended not only to capture the
diversity in the onset and decay of the phenomenon, but also to other climate
drivers such as the North Atlantic Oscillation (NAO), the Pacific–North
American pattern (PNA), the Pacific-South American pattern (PSA), the Southern
Annular Mode (SAM) or the Indian Ocean Dipole (IOD) for example when applied
to specific region of interests (Rogers and Loon, 1982; Lau et al., 1994; Saji
et al., 1999; Thompson and Wallace, 2000; Katz, 2002, for a review of Sir
Gilbert Walker contributions to the identification of large-scale climate
teleconnections such as the Southern Oscillation or the NAO). To further frame
the issue, a short and far from exhaustive history of ENSO follows with a
focus on the seminal works leading to the identification and naming of the
different flavours of ENSO.

Following Bjerknes (1966, 1969) early work on El Niño atmospheric
teleconnections linking it to the Sir Gilbert Walker ‘Southern Oscillation’,
Wyrtki (1975) proposed that El Niño episodes were due to anomalously high
southeast trades present in the central Pacific, over the two years preceding
El Niño, intensifying the subtropical gyre and associated South Equatorial
Current of the South Pacific, leading up to a build-up of water in the western
equatorial Pacific and increasing the sea level east-west slope across the
basin. Wyrtki postulated that the relaxation of these enhanced and persistent
trades would then trigger the opposite reaction whereby the water accumulated
in the West of the basin, over two years prior, would flow back eastwards
initiating El Niño, leading to a warming of the eastern equatorial SST and a
deepening of the equatorial thermocline. It is interesting to note that, at
the time, the La Niña phase of ENSO had not been coined neither by Bjerknes
nor Wyrtki in their seminal papers. However comparing and contrasting the
onset, persistence and decay of the 1957-1958, 1965-1966 and 1972-1973 El Niño
events, Wyrtki (1975) was led to conclude

> El Niño also appears in varying intensity and, like any other large-scale
> event in the ocean-atmosphere system, El Niño certainly does not have a
> single cause; no two El Niño events are quite alike

intimating at this early stage the diversity or flavours of El Niño. As far as
La Niña is concerned, the realisation that the condition corresponded to a
distinct phase of the ENSO phenomenon came some fifteen years later, although
Wyrtki may have alluded to it when mentioning the anomalously high southeast
trades present in the central Pacific as a precursor to El Niño, as the
enhancement of the trades would likely cool the Eastern equatorial Pacific
below its climatological values and lead to a further shoaling of the
thermocline there. According to the Oxford English Dictionary search tool,
Guillén et
al.333https://www.oed.com/view/Entry/252722?isAdvanced=true&result=2&rskey=RRdzwh&.
Interestingly, it has been rather challenging to find the peer-reviewed
publication where La Niña is first mentioned as the opposite phase of El Niño.
Izumo (2023) notes that the term La Niña may have emanated from George
Philander work through his Spanish speaking South-American wife. reported on
SST anomalies having a negative effects on primary production in Glantz and
Thompson (1981) and, it appears, mentioned for the first time the La Niña
moniker for the El Niño colder counterpart. The weather and climate community
had to wait for the aftermath of the intense 1997-1998 Niño, which was
followed by a rapid swing to a strong Niña condition, for directing its
attention to La Niña and its impacts (Glantz, 2002). The peculiar
characteristics and teleconnections of the 2002-2003 El Niño, when compared to
the conventional El Niños (Rasmusson and Carpenter, 1982, for an early
review), prompted Ashok et al. (2007) to undertake a comprehensive review of
the event and led to the definition of a distinct type of El Niño named El
Niño Modoki. The colder El Niño Modoki counterpart, La Niña Modoki, was not
mentioned in Ashok et al. (2007), but was introduced by Cai and Cowan (2009)
two years later when discussing the asymmetric impact on Australian rainfall
of what appears to be four distinct phases of ENSO. By 2009, four possibly
distinct ENSO conditions have been identified in the literature: El Niño, La
Niña, El Niño Modoki and La Niña Modoki.

Multiple approaches have been applied to the characterisation of the flavours
and naming of ENSO types. The specificity of the ENSO types discovered by each
of these methods is inherent to the methodological choices made, as it is with
our approach, but is not necessarily inherent to the phenomenon these methods
intend to capture. Hereafter, we will opt for the CP and EP naming convention
based on resemblance, not equivalence. We could have used Central/Eastern
Pacific (Kug et al., 2009; Takahashi et al., 2011; Sullivan et al., 2016),
Canonical/Modoki (Cane, 1983; Ashok et al., 2007; Cai and Cowan, 2009) or
Canonical/Dateline (Trenberth and Stepaniak, 2001) ENSO terminology
interchangeably.

With the reliance on both supervised and unsupervised machine learning (ML)
algorithms applied to weather and climate research, the number of new ENSO
classes is increasing, leading some to question if the inflationary nature of
the ML discovery of these classes corresponds to truly different climatic
processes (Christiansen, 2007; Johnson, 2013; Karnauskas, 2013; Kashinath et
al., 2021). Although there is still some debate whether or not the different
flavours of ENSO correspond to dynamically distinct phenomena (Ashok et al.,
2007; Kao and Yu, 2009; Karnauskas, 2013; Johnson, 2013; Capotondi et al.,
2015; Liguori and Di Lorenzo, 2019; Amaya, 2019; Capotondi et al., 2020; Feng
et al., 2020), various classification methods (Leloup et al., 2007; Takahashi
et al., 2011; Sullivan et al., 2016; Freund et al., 2019; Lemmon and
Karnauskas, 2019; Wiedermann et al., 2021; Black et al., 2022; Schlör et al.,
2023) allow to discriminate between the manifestation and evolution of either
a singular mode, plural modes or even a continuum of modes (Trenberth and
Stepaniak, 2001; Johnson, 2013; Williams and Patricola, 2018; Feng et al.,
2020; Schulte et al., 2021; Dieppois et al., 2021) of tropical variability
based on sea surface temperature anomalies (SSTAs).

In this work, we build on the AA results reported by Hannachi and Trendafilov
(2017), Black et al. (2022) and Chapman et al. (2022) for SSTA, who showed AA
unique ability in identifying extreme SST configurations. Section 3 introduces
all the datasets considered in this study and justifies some of the AA data
pre-processing steps. Section 4 summarises the AA method, together with
Section S1 of the supplementary material, provides an overview of AA in
layperson’s terms. It also compares and contrasts it to two commonly used
clustering algorithms: K-mean and fuzzy C-means. Section 5 presents AA results
applied to global SSTAs and discusses the resulting AA ENSO patterns and
indices, which are compared to more traditional methods of classification in
Section 6. Section 7 demonstrates that AA can be applied to S2S forecasts and
could complement current ENSO alert systems in characterising ENSO diversity.
Section 8 illustrates the reliability and relevance of the AA classification
of ENSO types by showing that extreme configurations detected in detrended
SSTAs (dSSTAs) correspond to robust and explainable teleconnection patterns of
various atmospheric and oceanic fields. We also present AA results performed
on detrended sea level anomalies showing that they provide a clearer
expression of ENSO types. ENSO phase transitions characterised by AA are
illustrated in Section 9. Section 10 concludes the paper by summarising our
findings and discusses potential applications of AA to possibly advance
weather and climate research.

We would like to conclude this introduction with an important caveat. Although
we will borrow from the ENSO nomenclature to label the AA discovered phases
with the CP and EP qualifiers, we do not claim that these phases correspond
exactly to those mentioned in the cited literature, but we do show that they
resemble these types as much as these types resemble one another. As mentioned
previously, we characterise the flavour or type in the ENSO context as a range
of clearly differentiated patterns manifest in the SST anomalies in the
Pacific with well defined centers of action or loci of absolute maximum
amplitudes and broad spatial extent.

##  3 Data

We apply AA to linearly detrended Optimum Interpolation Sea Surface
Temperature (OISST) v2.1 high resolution dataset (Reynolds et al., 2007)
global SST anomalies (SSTAs) over the 1982-2022 period. We mainly focus here
on monthly anomalies represented as the departure from monthly climatological
values defined as the time-mean of SST values for each month-of-year across
all years considered in the 1982-2022 interval with the local linear trend
computed over the same period removed at every grid point. Only complete
calendar years are considered.

Our justification for going global and not reducing the geographical domain in
the analysis is manifold. Firstly, limiting the domain to a given latitude-
longitude box is subjective, and may incur the risk of the box boundaries
cutting through important teleconnection pathways potentially useful to the
classification of ENSO events and their impacts.

Secondly and often without being explicitly stated, reducing the domain is
tantamount to infer a priori that the geophysical processes under
consideration within it are the most important drivers of interactions with
the rest of the climate system outside the initially chosen domain. As causal
links are particularly difficult to ascertain, these choices are at time
problematic. For high temporal resolution, from 1 to 1/5 days, Bach et al.
(2019) (Fig. 8) show that the local atmosphere-ocean predictability, based on
a spectral causality argument, is mainly driven by atmospheric processes apart
from a narrow tropical band, where the ocean drives the atmosphere. However,
for frequencies lower than 1/month, the regions where the ocean drives
predictability expand to the sub-tropics and higher latitudes. These findings
are based on local interactions between atmosphere and ocean only. The picture
is likely to change when remote teleconnections are taken into account. A
multivariate approach to assess non-local interactions between ocean and
atmosphere is challenging. Lagged correlation and linear response theoretical
arguments have been put forward to explain this relationship (Capotondi and
Sardeshmukh, 2017; Pegion et al., 2020; Capotondi and Ricciardulli, 2021).
Lagged-regression, -correlation or -covariance analyses are often used to
assess causality in weather and climate studies to investigate relationships
between variables. These methods often suffer from drawbacks where variables’
auto-correlation and their non-linear interactions have a confounding effect
as discussed recently by McGraw and Barnes (2018), for example. They advocate
for the implementation of dedicated causality-probing methods (Granger, 1969;
Barnett and Seth, 2014) as in Bach et al. (2019). However, these approaches
are based on time-mean (expectation) statistics and are difficult to compare
directly with AA results, which are not. As consequences of the monthly
temporal resolution used in this study and the potential causal and non-local
relationships between tropical, extra-tropical and high-latitude SSTs,
enlarging our domain of study can be justified. This choice can help to link
the ENSO tropical signature - its location, extent and intensity - to
concurrent extra-tropical conditions and to inform on its potential remote
impacts.

Thirdly, the choice of fixed geographical boxes, apart from being
subjective444The choice of climate indices, based on restricted geographical
domains and commonly used to benchmark climate driver strength and phases, may
suffer from 1) the non-stationarity of the underlying field under study and 2)
a lack of sample diversity when initially established as the history of ENSO
diversity exposed in Section 2 attests., corresponds to a crude way to capture
dynamical processes and is probably not suitable for a direct comparison with
ocean-atmosphere global climate model (OAGCM) simulations due to unresolved
model biases in the model representation of geophysical processes in general,
and in the spatiotemporal representation of ENSO in particular (Guilyardi et
al., 2020, for a recent exposé). As a case in point, Feng et al. (2020)
recently embarked in a thorough examination of the implication of these
choices for the detection and their power of discrimination when applied to
the ENSO diversity based on SSTAs. We refer the reader to it as an
illustration on how arbitrary and critical these choices can be in the
definition of various ENSO indices and associated patterns.

As we are interested in detecting the different flavours of ENSO, removing a
linear trend to SSTA prior to AA implementation helps to approximately
separate the anthropogenic greenhouse forcing and natural secular trends from
sub-seasonal-to-decadal internal climate mode variability, de facto assuming
that the internal variability of the climate system can be separated from
these processes operating at longer timescales and ignoring their potential
interaction. The assumption of separability is difficult to test given the
short high-quality global SST satellite records of 41 years at our disposal
and would be probably untenable for centennial to millennial timescales (Bunge
and Clarke, 2009; Levine et al., 2017; Han et al., 2020; Emile-Geay et al.,
2020; Jiang et al., 2021; Power et al., 2021, for ENSO variability over the
last millenium).

In the following, we justify the removal of the 41-year long linear trend from
SSTAs by referring to the recent publications of Wills et al. (2018), Wills et
al. (2020) and Jebri et al. (2020). Wills et al. (2018) propose supervised ML
pattern recognition methods to separate the anthropogenic forced response from
internal variability in both OAGCMs and observations using a combination of
low-pass filtering to separate long-timescale signal – the forced response –
from short-timescale noise – the internal variability – and a fingerprinting
technique to maximize the associated signal-to-noise variance. As
observational records, 100 years of Cowtan and Way (2014) infilled Hadley
Center Climate Research Unit version 4 surface temperature (HadCRUT4) is used.
Not only the global mean HadCRUT4 temperature, but also the basin mean SST
difference between eastern and western equatorial Pacific (Wills et al., 2020,
Figures 15a and 15b) for example, displays strong non-linear trends when
considered over the entire 1920 to 2020 period. However from 1980 onward,
these trends are approximately linear in Wills et al. (2020). In comparison,
Figures 1a and 1b display the first four empirical orthogonal function (EOF)
and principal component (PC) modes for both non-detrended and linearly
detrended SSTAs for OISST v2.1 (Reynolds et al., 2007). The first EOF/PC pair
in Figure 1a, corresponding to 14.1% of the variance, captures the global
warming pattern as indicated by the upward trending PC combined to the overall
positive EOF pattern, although the South-Eastern Pacific shows significant
cooling in a wedge-shaped domain, which has been linked to atmospherically
driven enhanced up-welling in response to both increasing greenhouse gas and
ozone depletion since the 1980s by Jebri et al. (2020) through dedicated
detection and attribution experiments. The first EOF/PC pair is followed by
three EOF/PC pairs highly ($>0.9$) and significantly (p-value $<10^{-6}$ at
0.05 significance level) correlated to the first three EOF/PC pairs of Figure
1b with a slight increase in variance for the individual modes of the
detrended case.

As we strive for simplicity in the data preparation, we posit that the forced
response to anthropogenic climate change and secular natural trends can be
approximated by a linear trend and removed from SSTAs over the 1982-2022
period. Furthermore, the choice of OISST v2.1 dataset over longer gridded SST
records such as NOAA Extended Reconstructed Sea Surface Temperature (ERSST)
Version 5555Boyin Huang, Peter W. Thorne, Viva F. Banzon, Tim Boyer, Gennady
Chepurin, Jay H. Lawrimore, Matthew J. Menne, Thomas M. Smith, Russell S.
Vose, and Huai-Min Zhang (2017): NOAA Extended Reconstructed Sea Surface
Temperature (ERSST), Version 5. NOAA National Centers for Environmental
Information. doi:10.7289/V5T72FNM [2023-07-07]. or COBE-SST Version 2
(Hirahara et al., 2014) is motivated by concerns that reconstructed SST
products prior to the satellite area may suffer from spatiotemporal sampling
inhomogeneities due to the disparity in the number of available observations
between the Southern and Northern hemispheres for example, potentially leading
to the mischaracterisation of ENSO phases or even the detection of spurious
ones, as illustrated and discussed in Sections S2 and S3 of the supplementary
material.

Figure 1: First four (rows) empirical orthogonal function (EOF, left column)
and principal component (PC, right column) modes of global a) non-detrended
and b) detrended SSTAs over the 1982 to 2022 period. The MEI index (black
line) is overlaid on the PCs. For each mode, the percentage of variance
explained is given on the bottom left corner of each subplot on the right and,
on the left. Reported in red, is the correlation coefficient computed at zero
lag over the 1982 to 2022 period between MEI and $PC_{2}$ and $PC_{1}$ for a)
and b), respectively.

The archetypes spatial and temporal imprints, and their linkages to extra-
tropical atmospheric circulation are revealed by compositing Japanese 55-year
Reanalysis (JRA-55) monthly averaged atmospheric fields (Kobayashi et al.,
2015) over overlapping time periods as described in Black et al. (2022). The
imprint of ENSO on global precipitation is further ascertained using Global
Precipitation Climatology Project (GPCP) Version 2.3 (Adler et al., 2018)
gridded data for the entire globe. Sea level anomalies (SLA) from AVISO+
Satellite altimetry data666The Ssalto/Duacs altimeter products were produced
and distributed by the Copernicus Marine and Environment Monitoring Service
(CMEMS) (http://www.marine.copernicus.eu) and Simple Ocean Data Assimilation
(SODA) Version
3.15.2777https://www2.atmos.umd.edu/~ocean/index_files/soda3.15.2_mn_download_b.htm
(Carton and Giese, 2008) ocean temperature composites will provide additional
evidence of the relevance of AA to the study of ENSO diversity for sub-surface
ocean fields. The Multivariate ENSO (Wolter and Timlin, 2011) and the Oceanic
ENSO Indices (MEI888https://psl.noaa.gov/enso/mei and
ONI999https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php)
will be compared to indices derived by AA.

##  4 Archetypal analysis

Archetypal analysis is an unsupervised pattern recognition method. Since its
development by Cutler and Breiman (1994), AA has only recently been adopted by
the weather and climate community to identify extreme geophysical
configurations and to investigate their impacts at the regional and global
scales (Steinschneider and Lall, 2015; Hannachi and Trendafilov, 2017;
Richardson et al., 2021; Risbey et al., 2021; Black et al., 2022; Chapman et
al., 2022).

The AA representation involves stochastic or probability matrices, the
matrices C and S defined hereafter. These matrices correspond to the
probability of expression of dSSTA records in archetypes (cluster definition),
themselves extreme geophysical configurations of the dataset as a whole, and
the probability of expression of archetypes in each and every dSSTA record
(cluster membership). This unique property is crucial and leads to a
probabilistic expression of the identified phases definition and membership.
It can be also applied to composites of any geophysical fields associated with
ENSO.

For a detailed review of the method and its various implementations, the
reader is referred to the following monographs, Hannachi (2021) and
Trendafilov and Gallo (2021), and specifically to Hannachi and Trendafilov
(2017) and Black et al. (2022) when applied to global SST. Hereafter, AA is
briefly introduced and compared to the commonly used principal component
analysis (PCA) in Table 1, PCA being also used hereafter as a dimension
reduction first step to speed-up computation as in Richardson et al. (2021),
Risbey et al. (2021), Black et al. (2022) and Chapman et al. (2022).

Singular value decomposition (SVD) |  Archetypal analysis (AA)
---|---
$\textbf{X}=\textbf{U}\Lambda\textbf{V}^{T}$ |  $\textbf{X}\approx\textbf{XCS}$
(Deterministic) solution of |  (Non-deterministic) solution of
$\underset{\textbf{U},\textbf{S},\textbf{V}}{\arg\min}\|\textbf{X}-\textbf{U}\Lambda\textbf{V}^{T}\|_{F}^{2}$ |  $\underset{\textbf{C},\textbf{S}}{\arg\min}\|\textbf{X}-\textbf{XCS}\|_{F}^{2}$
with (quadratic) constraints |  with (convexity) constraints
$\sum_{s}U_{sr}U_{st}=\delta_{rt}$ |  $\sum_{p}S_{pt}=\textbf{1}_{t},\,S_{pt}\geq 0,\,\forall p,t$
$\sum_{t}V_{rt}V_{st}=\delta_{rs}$ |  $\sum_{t}C_{tp}=\textbf{1}_{p},\,C_{tp}\geq 0,\,\forall p,t$
$\text{Trace}(\Lambda\Lambda^{T})=\text{Trace}(\textbf{XX}^{T})=\|\textbf{X}\textbf{X}^{T}\|_{F}^{2}$ |
Table 1: Compare and contrast SVD and AA matrix factorisation methods. The
symbol $\|\textbf{X}\|_{F}=\sqrt{\sum_{i}\sum_{j}|X_{ij}|^{2}}$ denotes the
Froebenius norm of the matrix X. The superscript $T$ denotes the matrix
transpose.

AA is not only a matrix factorisation method (Cichocki et al., 2009; Elad,
2010; Elden, 2019; Gan et al., 2020; Gillis, 2020; Hannachi, 2021; Trendafilov
and Gallo, 2021), but, like PCA, is also a dimension reduction procedure
(Izenman, 2008; Nguyen and Holmes, 2019), whereby a dataset, represented by a
data matrix X , $\textbf{X}=\textbf{X}_{s\times t}$ of space-time dimensions
($s,t$), can be reduced to or approximated by a product of factors such that
$\textbf{X}\approx\textbf{XCS}=\textbf{X}_{s\times t}\textbf{C}_{t\times
p}\textbf{S}_{p\times t}=\textbf{XC}_{s\times p}\textbf{S}_{p\times t}$, where
the summation over the repeated indices $t$ and $p$ is implied. The factors,
$\textbf{XC}_{s\times p}=\textbf{X}_{s\times t}\textbf{C}_{t\times p}$, are
called archetypes. They are located on the convex hull of the dataset X and
are convex combinations of data points, the matrix
$\textbf{C}=\textbf{C}_{t\times p}$ being a right-stochastic matrix with each
row summing to 1, $\sum_{t}C_{tp}=\textbf{1}_{p},\,C_{tp}\geq 0,\,\forall
p,t$. The number of archetypes, $p$, the cardinality of the archetypal
representation, is chosen a priori. The other factors,
$\textbf{S}=\textbf{S}_{p\times t}$, can be viewed as soft-clustering
affiliation weights such that points in X can be approximated by a convex
combination of archetypes, $\textbf{X}\approx\textbf{XC}_{s\times
p}\textbf{S}_{p\times t}$, with S being also a right-stochastic matrix with
$\sum_{p}S_{pt}=\textbf{1}_{t},\,S_{pt}\geq 0,\,\forall p,t$101010Here,
$\textbf{1}_{p}$ and $\textbf{1}_{t}$ stand for unit vectors of dimensions $p$
and $t$..

Both matrices, C and S, are found by minimizing the square distance between
the original dataset X and its reduced representation XCS,

$\underset{\textbf{C},\textbf{S}}{\arg\min}\|\textbf{X}-\textbf{XCS}\|_{F}^{2},$
(1)

where the distance function $\|\cdot\|_{F}$ stands here for the Froebenius
norm, the convexity constraints on C and S being enforced in the minimization
algorithm by adding penalty terms to Expression 1 (Cutler and Breiman, 1994;
Mørup and Hansen, 2012).

In their original paper, Cutler and Breiman (1994) proved through a reductio
ad absurdum111111https://en.wikipedia.org/wiki/Reductio_ad_absurdum. argument
that the archetypes, XC, can only be on the convex hull of the dataset under
consideration. In broad terms, the main reasons for it lay in the convexity
constraints such that the archetypes are a convex combination of data points
and that the data points themselves are approximated by a convex combination
of archetypes, and in the mathematical form of the cost function. This
remarkable property can be used therefore to identify extreme configurations
in geophysical observations (Hannachi and Trendafilov, 2017; Risbey et al.,
2021; Black et al., 2022; Chapman et al., 2022) and, conversely, to express
observations in term of these extreme configurations.

To illustrate the nature of AA, Figure 2 summarises its main characteristics
for cardinality 4 and a highly truncated 2-dimensional dataset, where the data
matrix $\textbf{X}_{2\times t}=[\lambda_{1}PC_{1},\lambda_{2}PC_{2}]$ is
chosen to be the first 2 EOF/PC pairs of dSSTAs, illustrated in Figure 1b,
rows 1 and 2. Its convex hull is represented by the black polygon
encapsulating all the data points in the set. The archetypes
$\textbf{XC}_{2\times 4}$, the coloured circles located on the hull,
correspond to extreme global configurations. The associated archetypal
patterns are reported in each quadrant and recover the familiar warm/cold and
eastern/central phases of ENSO. The green tetragon, whose vertices are the
archetypes $\textbf{XC}_{2\times 4}$, is an approximation of the dataset
convex hull. It encapsulates all the points of the representation
$\textbf{XCS}_{2\times t}$. These points are found by solving Expression 1 for
C and S by minimising the distances, illustrated by the thin black segments in
Figure 2, between the original data points in grey and their approximations in
colour. The coloured squares stand for the composites of X,
$\textbf{X}\tilde{\textbf{S}}^{T}$, constructed from the archetype mean
probability of expression given by the S as in Equation 3. In the Section S1
of the supplementary material, we show that this example can be extended to 3
dimensions (3D) with datasets such as the ‘Stanford
Bunny’121212https://en.wikipedia.org/wiki/Stanford_bunny, a computer graphics
3D test model, and the first 3 EOF/PC pairs of dSSTAs, for selected
cardinalities. In the rest of the paper, we use 492 PCs of dSSTAs and so a
492-dimensional dataset.

Figure 2: AA results for cardinality 4 based only on dSSTAs first 2 PCs
multiplied by their respective eigenvalues, $\textbf{X}_{2\times
t}=[\lambda_{1}PC_{1},\lambda_{2}PC_{2}]$. The 2D convex hull is the black
polygon. Points in the pointset are represented by coloured dots if they lay
within the convex hull approximation - the green tetragon - which 4 vertices
correspond to the 4 archetypes, $\textbf{XC}_{2\times 4}$ (large coloured
circles). Individual data points located outside the AA convex hull
approximation are coloured grey, the data points within the AA approximation
of the convex hull are coloured with a mixture of pure colours corresponding
to the 4 archetypes. The respective archetypal maps corresponding to
$\textbf{XC}_{2\times 4}$ are reported in each quadrant. The coloured squares
mark the composites of X, $\textbf{X}\tilde{\textbf{S}}^{T}$ (Equation 3),
constructed from the archetype mean probability of expression given by the S.
The arrows link the archetypes $\textbf{XC}_{2\times 4}$ to their respective
composites $\textbf{X}\tilde{\textbf{S}}^{T}$.

Hereafter, AA solutions will be compared to results from two commonly used
pattern recognition or clustering methods: K-means (KM) and fuzzy C-means
(FCM) (Steinhaus, 1957; Forgy, 1965; Lloyd, 1982). As in Mørup and Hansen
(2012), Table 2 summarises the main characteristics of KM and FCM factors with
those of PCA and AA. The main difference between AA and both clustering
methods, KM and its extension FCM, is that the KM and FCM cluster centres do
not generally correspond to extreme configurations, as they are not
necessarily located on the convex hull of the dataset.

Table 2: Relation between the AA model and unsupervised methods such as SVD/PCA, fuzzy C-means, and K-means (Steinhaus, 1957; Forgy, 1965; Lloyd, 1982) for the factorisation of X adapted from Mørup and Hansen (2012), where the symbol $\|\textbf{v}\|_{1}=\sum_{i}|v_{i}|$ denotes the Taxicab or Manhattan norm $l_{1}$ and $\|\cdot\|$ the Euclidean distance or $l_{2}$ norm. The $s\times t$ elements of the data matrix X are denoted by $X_{ij},\forall i=1,\ldots,s\,(\text{spatial dimension}),j=1,\ldots,t\,(\text{time dimension})$, whereas the index $p$ corresponds to the clustering cardinality, the exponent $m\in(1,+\infty)$ stands for the C-means ‘fuzzyfier’ (Dunn, 1973; Bezdek et al., 1984), $\textbf{x}_{t}=X_{\cdot t}$ and $\textbf{c}_{c}$ are the data records and the cluster centers in vector form, with $c=1,\ldots,p$. The symbols $\mathbb{R}$ and $\mathbb{B}$ respectively stand for the sets of real numbers and booleans, consisting of two elements, false-0 and true-1, used for categorical classification. SVD/PCA | AA/PCH | Fuzzy C-means131313https://en.wikipedia.org/wiki/Fuzzy_clustering | K-means
---|---|---|---
$\textbf{C}\in\mathbb{R}$ | $\|\textbf{c}_{p}\|_{1}=1,\textbf{C}\geq 0$ | $C^{m}_{tc}=\frac{S_{ct}^{m}}{\sum_{t}S_{ct}^{m}}$ | $\|\textbf{c}_{p}\|_{1}=1,\textbf{C}\geq 0$
$\textbf{S}\in\mathbb{R}$ | $\|\textbf{s}_{t}\|_{1}=1,\textbf{S}\geq 0$ | $S^{m}_{ct}=\left\\{\sum_{k}^{p}\left(\frac{\left\|\textbf{x}_{t}-\textbf{c}_{c}\right\|}{\left\|\textbf{x}_{t}-\textbf{c}_{k}\right\|}\right)^{\frac{2}{m-1}}\right\\}^{-1}$ | $\|\textbf{s}_{t}\|_{1}=1,\textbf{S}\in\mathbb{B}$

Time-mean composites of any geophysical field can be straightforwardly
constructed from the AA factors, C and S, as weighted or conditional means,
where the matrix elements $C_{tp}$ and $S_{pt}$ can be viewed as 1) the
probability of a data record at time $t$ to contribute to the definition of
archetype $p$ for C or 2) the probability of archetype $p$ to be expressed in
a data record at time $t$ for S, respectively. Let’s recall the AA
approximation,

$\textbf{X}\approx\textbf{XCS},$

with convexity constraints $\sum_{p}S_{pt}=\textbf{1}_{t}$, $S_{pt}\geq 0$ and
$\sum_{t}C_{tp}=\textbf{1}_{p}$, $C_{tp}\geq 0,\,\forall p,t$. Comparing
composites of the original dataset X as in Black et al. (2022) Figures 6 to 8,
XC and $\textbf{X}\bar{\textbf{S}}^{T}$, separately built on C and S factors,
is instructive as these composites may inform on the robustness of the AA
factorisation; the more similar they are, the stronger the discrimination or
separation between archetypes. Typically, for any given archetype $p$, the
non-zero elements of the sparser matrix $C_{tp}$ at record $t$ should be
reflected in the increased value of the element $S_{pt}$ at record $t$ for the
denser matrix S when compared to the time-mean
$\bar{S}_{p}=\frac{1}{N}\sum_{t}S_{pt},\,\forall p$, where $N$ is the total
number of records. However, the reverse assertion does not necessarily hold
true as the sparsity of C, resulting from the optimization, may well favour a
single data snapshot as archetype even if the corresponding condition is
repeated at multiple occasions across the time records considered.

For any spatiotemporal field, $\textbf{F}=\textbf{F}_{s\times t}$141414The
spatial dimensions of F and X do not need to be identical. We have kept the
same symbol $s$ for $\textbf{F}=\textbf{F}_{s\times t}$ to simplify the
notation., AA (time-)mean composite for each archetype $p$ can be constructed
in similar fashion as weighted means

$\textbf{F}^{C}_{s\times p}=\textbf{FC}=\sum_{t}F_{st}C_{tp}$ (2)

for C-composites, or as

$\textbf{F}^{S}_{s\times
p}=\textbf{F}\bar{\textbf{S}}^{T}=\frac{\sum_{t}F_{st}S^{T}_{pt}}{\sum_{t}S_{pt}}$
(3)

for S-composites. To generate seasonal composites, the summation over time in
Equations 2 and 3 is only performed on the appropriate month-of-year records,
hereafter December-January-February (DJF), March-April-May (MAM), June-July-
August (JJA) and September-October-November (SON).

The corresponding ‘in-sample’ AA composite uncertainty or spread is expressed
as
follows151515https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights..
With the AA composite weighted mean $\textbf{F}^{S}_{s\times p}$ constructed
as above, $\sum_{t}S_{pt}=W^{S}_{p}$ and the weights
$w_{pt}=S_{pt}/W^{S}_{p}$, the covariance matrix ${\textbf{C}^{S}_{F,s\times
p}}$ reads

$\textbf{C}^{S}_{F,s\times
p}=\frac{\sum_{t}w_{pt}}{\left(\sum_{t}w_{pt}\right)^{2}-\sum_{t}w_{pt}^{2}}\sum_{t}w_{pt}\left(F_{st}-F^{S}_{sp}\right)^{T}\left(F_{st}-F^{S}_{sp}\right).$
(4)

A similar, but simpler expression, can be derived for the covariance matrix
${\textbf{C}^{C}_{F,s\times p}}$, corresponding to C-composites,

$\begin{split}\textbf{C}^{C}_{F,s\times
p}&=\frac{\sum_{t}C_{tp}}{\left(\sum_{t}C_{tp}\right)^{2}-\sum_{t}C_{tp}^{2}}\sum_{t}C_{tp}\left(F_{st}-F^{C}_{sp}\right)^{T}\left(F_{st}-F^{C}_{sp}\right),\\\
&=\frac{1}{1-\sum_{t}C_{tp}^{2}}\sum_{t}C_{tp}\left(F_{st}-F^{C}_{sp}\right)^{T}\left(F_{st}-F^{C}_{sp}\right),\end{split}$
(5)

as the weights reduce to $C_{tp}$ as $W^{C}_{p}=\sum_{t}C_{tp}=1,\forall p$.
Of note is that similar expressions for the KM and FCM cluster centers and
associated composites can be constructed.

The introduced ‘in-sample’ measure of significance based on weighted
covariance, Equations 5 and 4, indicates whether or not the archetypal
patterns discovered are repeatable or reoccur over the period considered in a
conditional time-mean sense by comparing pattern signal-to-noise (or mean-to-
spread) ratios. We note that the ‘in-sample’ signal-to-noise ratios are a very
conservative measure of reliability, especially for composites based of the
affiliation weights, the stochastic matrix S. For ‘out-of-sample’ AA composite
uncertainty or spread determination, a less stringent measure of reliability
or significance, the problem rapidly becomes computationally intractable for
(spatial or reduced) dimensions larger than 9161616http://www.qhull.org
(Barber et al., 1996), as one is required to randomly sample the convex hull
of the original or reduced dataset. An even less stringent approach as in
Chapman et al. (2022), would be to employ a brute-force Monte-Carlo approach,
1) to generate an ensemble of random synthetic stochastic matrices to
replicate the features of both the C and S, 2) to compute the composites on
these, and 3) to declare each point of the composites ‘significant’ if not
falling within a certain percentile range constructed from the synthetic
ensemble. As archetypes correspond by construction to extreme configurations,
the latter method mentioned is too permissive. The composites generated from
random stochastic matrices are not necessarily located on the convex hull of
the dataset and cannot therefore be considered ‘legitimate’ samples as
archetypes are located on the convex hull by construction.

##  5 Archetypal flavours and indices of ENSO

AA is applied to linearly detrended monthly mean SSTAs (dSSTA) over the
1982-2022 period. To increase computational speed, we implement the reduced
space AA (RSAA) methodology reviewed in Black et al. (2022), whereby a PCA
dimension reduction step is performed on the dataset prior to apply AA, to
factor out of the cost function – Expression 1, the spatial dependency given
by the orthonormal empirical functions U in Table 1. AA is done on all
$12\times 41=492$ PCs multiplied by their corresponding eigenvalues,
$\Lambda\textbf{V}^{T}$, where 100% of variance of dSSTA is retained for the
ensuing AA step. The factors
$(\textbf{XC})_{RSAA}=\textbf{U}(\Lambda\textbf{V}^{T})\textbf{C}$ and
$\textbf{S}_{RSAA}$ are now the archetypal patterns and affiliation
probabilities. It can be shown that $(\textbf{XC})_{RSAA}$ and
$\textbf{S}_{RSAA}$ are almost identical to the AA solutions of the non-
reduced problem (Black et al., 2022, Figure 3). To simplify the notation, the
subscript $RSAA$ is omitted from this point onwards. An identical dimension
reduction strategy is applied to both KM and FCM clustering algorithms. For a
given cardinality or cluster number $p$ and in all three cases, AA, KM and
FCM, multiple trials are necessary to make sure the solutions are close to
optimal. Typically, 1000 trials are performed with different random
initialisations and the solutions correspond to the lowest cost, Expression 1,
across all trials.

Figure 3: Reduced space (RS) detrended SSTA cluster centers for K-means, fuzzy
C-means, AA C-composites XC and S-composites $\textbf{X}\textbf{\~{S}}^{T}$
(columns). The cluster centers for each method (rows) are referenced to the AA
results for a cardinality of 4 (last column) and assigned using the euclidean
distance. Each pattern has been normalised by its maximum absolute value for
increased contrast across clusters, reported on the Eurasian continent in °C.
The ‘fuzzyfier’ exponent $m$ was set to 1.05 for C-means. The labels $A_{p}$
in each row indicate the archetype ranks based on the time mean of the AA
stochastic matrix $\textbf{S}_{p\times t}$,
$\overline{S}_{p}=\frac{1}{N}\sum_{t}S_{pt}$ with $p=1,\ldots,4$, where $N$ is
the total number of records.

As in Black et al. (2022), the dSSTA data matrix, $\textbf{X}_{s\times t}$, is
weighted spatially by W with elements $W_{ij}=\sqrt{cos(\theta_{i})}$,
$\theta_{i}$ being the latitude of the grid cell $i=1,\ldots,s$ for all times
$j=1,\ldots,t$ prior applying PCA, AA, KM and FCM to $\textbf{XW}_{s\times
t}=\textbf{X}_{s\times t}\odot\textbf{W}_{s\times t}$171717Element-wise or
Hadamard matrix product, $(A\odot B)_{ij}=A_{ij}B_{ij}$. rather than
$\textbf{X}_{s\times t}$ in order to compensate for unequal grid spacing.
Therefore the tropical data points are carrying more weight than mid- and
high-latitude ones. As not to render our notation too cumbersome, we omit the
weights, $W_{ij}$. However, when KM, FCM and AA composites are computed,
unweighted snapshots based on the dSSTA cluster affiliation sequences are
used.

The choice of the domain driving the analysis is motivated by minimizing the
number of arbitrary choices with the explicit goal of capturing both the
tropical character and remote teleconnection of the ENSO phenomenon at zero
monthly lag. In Section 2 of the supplementary material, we have repeated the
analysis for Pacific Ocean regions over a wide range of latitude bands ranging
from $\pm 10^{\circ}$ to $\pm 90^{\circ}$ centered on the Equator with a fixed
longitudinal sector from $120^{\circ}$E to $280^{\circ}$E and compared
archetypes, XC, for cardinality 8 using pattern correlation between
overlapping regions of the restricted and global domains in Tables S1 and S2.
From the latitudinal band $\pm 50^{\circ}$ and higher, individual regional
archetypes can be uniquely matched to the global ones, with high correlation
values. Alternatively, when comparing archetypes from one latitudinal band to
the next in increments of $\pm 5^{\circ}$, starting with the Pacific region
$\pm 10^{\circ}$ latitudinal band as in Tables S3 and S4, individual regional
archetypes can only be uniquely matched from the latitudinal band $\pm
55^{\circ}$ and higher, an indication that to discriminate between extreme
dSSTAs conditions, the ‘latitudinal reach’ has to be increased.

Both matching strategies show that, to discriminate between archetypes across
bands, it is necessary to select Pacific regions which latitudinal extent
larger than about $\pm 50^{\circ}$ centered on the Equator. As we strive for
parsimony, but also would like to capture the global teleconnections (at lag
zero) across oceanic basins, we have therefore opted for a global domain.

Figure 3 compares KM, FCM and AA results for a cardinality of 4. At first
glance, the 4 cluster centers and archetypes (rows), labelled
$A_{i},\,i=1,\ldots,4$, share similarities across methods (columns). The
patterns in first two rows resemble La Niña and El Niño Modoki – or Central
Pacific (CP) – ENSO conditions with large anomalous SSTs located in the
equatorial central Pacific, whereas the third and fourth rows display cold and
warm anomalies in the Eastern Pacific, reminiscent of the canonical – or
Eastern Pacific (EP) – ENSO patterns for La Niña and El Niño, respectively.
The AA CP-Niño pattern detected for cardinality 4, $A_{2}$ in Figure 3, shares
commonalities with the so-called Pacific Meridional Mode (PMM) as defined by
Chiang and Vimont (2004). A clearer discrimination between PMM- and CP-Niño-
like conditions will only be revealed for higher cardinalities. A closer
inspection of these patterns for KM and FCM reveals both Central Pacific cold
and warm patterns are almost anti-symmetric; a property shared neither by the
AA C\- nor S-composites, Equations 2 and 3 where F has been replaced by X.

Figure 4: dSSTA cluster affiliation sequences for K-means, fuzzy C-means and
AA. The coloured patches denote K-means ‘categorical’ states corresponding to
CP Niña, CP Niño, EP Niña and EP Niño intervals. Thick (thin) coloured lines
stand for the AA (C-means) ‘soft’ affiliation probabilities
$S_{pt}\,(S^{m}_{pt})\in[0,1]$ with $p=1,\ldots,4$ corresponding to the
spatial patterns depicted in rows of Fig. 3.

The asymmetry in amplitudes, spatial patterns, timing and teleconnection
characteristics between the cold and warm phases of ENSO has been widely
reported (Hoerling et al., 1997; Sardeshmukh et al., 2000; Hannachi et al.,
2003, for early references on the subject) and is reflected by the affiliation
probabilities reported in Figure 4. However, the anti-symmetric nature of the
KM and FCM clusters is an artefact of cardinality 4 as is evident for higher
cluster cardinalities reported in Figure S7.

KM corresponds to a categorical clustering method where each data record is
associated with a single pattern or cluster center with affiliation
$\|\textbf{s}_{t}\|_{1}=1,\textbf{S}\in\mathbb{B}$ as in Table 2. Figure 5
compares KM to the soft-clustering methods, FCM and AA, whereby the ENSO phase
categorical or ‘hard’ affiliations are given by the states indexed by
$i_{max}(t)=\underset{i={1,\ldots,p}}{\arg\max}(S_{it})$ and
$i^{m}_{max}(t)=\underset{i={1,\ldots,p}}{\arg\max}(S^{m}_{it})$ corresponding
to the maximum values of the AA and FCM affiliation time-series $S_{it}$ and
$S^{m}_{it}$, respectively labeled by $i_{max}(t)$ and $i_{max}^{m}(t)$, and
contributing with highest probability to the snapshot $X_{st}$ for each time
record $t$181818A tie, although always possible, is very unlikely given that
the algorithms used implement double-precision arithmetic., indicating the
‘winning’ archetype or FCM cluster - the archetype or FCM cluster
corresponding to the highest probability of expression given by $i_{max}(t)$
or $i^{m}_{max}(t)$.

Figure 5: ENSO phase ‘hard’ affiliation sequences for AA, fuzzy C-means, and
K-means. The coloured patches denote ‘categorical’ states corresponding to CP
Niña, CP Niño, EP Niña and EP Niño intervals. For AA and fuzzy C-means
clustering, the phases labeled by $i_{max}(t)$ and $i_{max}^{m}(t)$ correspond
to highest affiliation probabilities, $S_{it}$ and $S^{m}_{it}$, are colour-
coded.

Overall, the ‘hard’ ENSO state affiliation mostly agrees across methods for
cardinality 4, especially over EP Niño (dark red) and CP Nina (light blue)
intervals in Figure 5. Of interest is the 2001-2007 interval, where both the
KM and its soft-clustering extension FCM disagree with the AA categorical
representation. However, increasing the cardinality improves the concordance
between KM, FCM and AA as shown for the patterns displayed in Figure S7, where
not only the issue of KM and FCM phase anti-symmetry raised previously is
resolved, but also a better categorical cluster agreement across methods is
achieved, as shown in Figures 6 and 7.

Due to the archetypal ‘nestedness’191919Archetypal ‘nestedness’, although not
expected apriori in the original formulation of Cutler and Breiman (1994),
describes the remarkable property that archetypes can be matched across
cardinality or number of clusters. originally reported by Black et al. (2022)
for dSSTAs, the 4 new distinct archetypes in Figure S7 for cardinality 8, not
only adds more diversity to the ENSO states already identified in Figure 4,
but also reveals the complexity of the ENSO phenomenon and its sequential
evolution over the 1982-2022 period. The 2001-2007 interval, where substantial
discordance between KM, FCM and AA categorical affiliation was observed, shows
now a better agreement across methods.

Figure 6: dSSTA cluster affiliation sequences for K-means, fuzzy C-means and
AA for cardinality 8. The coloured patches denote K-means ‘categorical’ states
corresponding to CP Niño, $A_{2}$, CP Niña, $A_{4}$, EP Niño, $A_{6}$, EP
Niña, and $A_{8}$ intervals. Thick (thin) coloured lines stand for the AA
(C-means) ‘soft’ affiliation probabilities $S_{pt}\in[0,1]$ with
$p=1,\ldots,8$ corresponding to the spatial patterns depicted in Fig. S7 rows.
Figure 7: ENSO phase ‘hard’ affiliation sequences for AA, fuzzy C-means, and
K-means for cardinality 8. The coloured patches denote ‘categorical’ states
corresponding to CP Niño, $A_{2}$, CP Niña, $A_{4}$, EP Niño, $A_{6}$, EP
Niña, and $A_{8}$ intervals. For AA and fuzzy C-means clustering, the phases
corresponding to highest affiliation probabilities are colour-coded.

Figures S15a and S15b, updated from Black et al. (2022), illustrate the impact
of cardinality on the affiliation probabilities expressed by the matrix S and
the patterns ‘nestedness’, demonstrating that no spurious blending occurs
across dSSTA states or clusters when extreme conditions are ‘well-expressed’,
as can be seen when a single archetype corresponding to the 3 major EP Niño
intervals of 1982-1983, 1997-1998 and 2015-2016, matches across 2-8
cardinalities as the dark red coloured intervals indicate. For a cardinality
of 2 (the supplementary material, top row in Figure S15a, the other records
when not corresponding to pure CP La Niña intervals in light blue, being
convex combinations of archetypes, have to be expressed as ‘blended’ patterns
by construction.

Taking advantage of the convexity constraint on S, the extent of how ‘well-
expressed’ or not an archetype is, can be measured by a discrimination score,
introduced for AA by Risbey et al. (2021),

$\Delta^{p}_{i_{max}}(t)=1-\frac{1}{p-1}\left\\{\frac{1}{\underset{i=1,\ldots,p}{\max}(S_{it})}-1\right\\}=1-\frac{1}{p-1}\left\\{\frac{1}{S_{i_{max}t}}-1\right\\},$
(6)

where $\Delta^{p}_{i_{max}}(t)\in[0,1]$ with
$S_{i_{max}t}=\underset{i=1,\ldots,p}{\max}(S_{it})$ is obtained from maximum
values of $\textbf{S}_{p\times t}$ corresponding to the archetype labelled by
$i_{max}=i_{max}(t)$ contributing with highest probability to the snapshot
$\textbf{X}_{s\times t}$ for each time record $t$202020A similar expression
can be devised for FCM affiliation probabilities $S^{m}_{pt}$ in Table 2..
Here, the symbol $i_{max}$ indicates that the discrimination score not only
corresponds to a value ranging from 0 to 1, but also to the label of the
‘winning’ archetype. A low discrimination score, $\Delta^{p}_{i_{max}}(t)$
close to zero, corresponds to a mixed or blended state at time $t$ with
probabilities being approximately the same and equal to $1/p$ for all
archetypes for a given cardinality $p$, whereas a discrimination score close
to 1 indicates a pure archetype or ‘well-expressed’ state at time $t$.

Figures 8a and 8b compare AA affiliation sequences, the AA ENSO indices, for
cardinalities $p=4$ and $p=8$. The discrimination scores,
$\Delta^{p}_{i_{max}}(t)\in[0,1]$ for all monthly records $t$, are represented
by archetypal colour-coded dots of the archetype $i_{max}(t)$. Increasing the
cardinality from 4 to 8 leads to overall higher discrimination scores where
they were initially low for cardinality 4 over the periods 1983-1997,
1998-2015 and 2016-2022, and introduces 4 new and distinct patterns
corresponding archetypes $A_{2}$, $A_{4}$, $A_{6}$ and $A_{8}$ in Figures S7,
for overall time-mean, and S8, for seasonal composites, respectively. As such,
increasing the cardinality tends to improve the discrimination between states
and reduces the number of ‘blended’ states. When checked against Figure 9, the
CP and EP ENSO patterns previously identified for cardinality 4 remain.
However, the onset and decay of the 3 major EP Niño intervals of 1982-1983,
1997-1998 and 2015-2016, are now different with $A_{4}$ $\rightarrow$ EP Niño
$\rightarrow$ EP Niña, $A_{4}$ $\rightarrow$ EP Niño $\rightarrow$ CP Niña,
and CP Niño $\rightarrow$ EP Niño $\rightarrow$ EP Niña, respectively. Apart
from the 1998-2000 period, the long lasting CP Niña intervals detected for
cardinality 4 over 1988-1989, 2007-2008, 2010-2012, and 2021-2022, are now
further split by the introduction of an EP Niña-like and EP Niño-like
patterns, $A_{2}$ and $A_{6}$. The interval commonly referred in the
literature (Capotondi et al., 2015; Wiedermann et al., 2021) as the (Niño)
Modoki period over 2001-2007 is now mainly captured by CP Niño-like archetypes
$A_{4}$ and $A_{8}$.

Figure 8: Comparison of AA affiliation sequences for cardinalities of a) $p=4$
and b) $p=8$. The discrimination scores (Risbey et al. (2021))
$\Delta^{p}_{i_{max}}(t)\in$ [0, 1], $\,\forall i=1,\ldots,p$ for all monthly
records $t$ are represented by archetypal colour-coded dots of the winning
archetype. A higher discrimination score indicates less blended state. Both
MEI (black) and ONI (grey) time series were scaled to fit within the
$S_{4/8\times 492}$ and $\Delta_{4/8}$ bounds.

The dSSTA composites based on time averages of the S-affiliation matrix over
the 1982-2022 period for a cardinality of 4, Figure 9, show that archetypal
patterns across seasons resemble the annual S\- and C-composites displayed in
Figure 3, but vary in amplitude from the tropics to higher latitudes. Of note
are the strong seasonal differences expressed by archetypes corresponding CP-
Niño ($A_{2}$) and EP-Niña ($A_{3}$), possibly indicative of more transitory –
but still extreme – global SST conditions. For CP-Niño, the positive SST
anomalies approximately centered around the date-line in JJA, migrate
westwards in SON, DJF and MAM, and are accompanied by a strengthening,
relative to the Tropics that is, of a PMM-like pattern in the North Pacific
(Chiang and Vimont, 2004; Amaya, 2019; Amaya et al., 2020; Richter et al.,
2022). This is reminiscent of the North Pacific SST conditions of 2013-2015
preceding the 2015-2016 EP Niño event, as Figures 8a and 8b attest, and also
reported by Tseng et al. (2017), for example. The 2013-2015 CP-Niño (PMM-like)
conditions persist across AA cardinalities ranging from 3 to 8 as can be seen
on Figures S15a (rows 2 to 8) and S15b (row 1). The EP-Niña ($A_{3}$) phases
are not as persistent for cardinality 4. An extra EP-Niña like pattern,
$A_{2}$, is now expressed for cardinality greater than 4, with the original
EP-Niña previously detected for a cardinality 4, now solely well expressed
over 2016-2017 interval (dark blue line) in Figure 8b.

Figure 9: AA seasonal lag-0 composites (columns) and signal-to-noise ratios
($\textbf{+}>0.15,\bullet>0.3$, $(F^{2}_{S}/C^{S}_{F})^{1/2}$) based on time
averages of the S-affiliation matrix over the 1982-2022 period for a
cardinality of 4. As each individual subplot is scaled to increase contrast
across archetypal patterns (rows), the boldfaced numbers on the Eurasian
continent label the maximum absolute SSTA for each archetypal pattern
$\textbf{A}_{i},i=1,...,4$ and seasons (DJF, MAM, JJA and SON).

In summary, the comparison between AA results across cardinalities shows that
the ENSO phenomenon takes multiple guises. AA, applied here to global dSSTAs,
is able to distinguish between both recurrent patterns such as those observed
over the well-known, but short, EP-Niño 1982-1983, 1997-1998 and 2015-2016
intervals, and a series of intermittent ones between these.

Recent work by Schlör et al. (2023) derived from fuzzy unsupervised clustering
also classifies the AA-detected EP-Niños as Extreme El Niños and our results,
based on the detection of extreme configurations in global SSTA, corroborate
these findings. Over the period 1982-2022 under study, the EP-Niño phases are
typically followed by Niña-like conditions of various types and interspersed
by longer periods of somewhat muted CP-Niño conditions. The archetypal
sequences over these ‘inter’–EP-Niño intervals are all different as can be
seen on Figures 7 and 8b. This suggests that ENSO likely corresponds to a
‘continuum’ of types as advocated by Johnson (2013) using a self-organizing-
map (SOM, Kohonen (1990)) or by Dieppois et al. (2021) when analysing changes
in the distributions of ENSO phases’ location and intensity. The question of
statistical distinguishability of these patterns raised in Johnson (2013) is a
difficult one to assess given the short interval of high-quality satellite
data with global coverage available. Three EP-Niño events over 41 years can
barely be considered as a large enough sample size for statistical inferences.
However, we will show in Section 8 that archetypes, expressing extreme global
dSSTA conditions, correspond to distinct teleconnection composite patterns in
atmospheric, surface and sub-surface fields. From this point onward,
archetypal patterns and geophysical field composites pertaining to a
cardinality of 8, will be given in the supplementary material.

An attentive reader would notice the unusual feature of AA that the event
categories do not include a ‘neutral’ condition. The ‘neutral’ concept stems
mainly from the adoption of tercile-type (negative/neutral/positive)
classification based on univariate indices used by operational centers to
qualify the current state of the climate. These have the advantage of
simplicity in the messaging, but the over-reliance on these can be fraught and
controversial as recently illustrated by the difficulty of identifying the
2023 El Niño event.

The ‘neutral’ ENSO condition can easily be implemented with AA by inspecting
the discrimination score $\Delta^{p}_{i_{max}}(t)$, Equation 6, implemented
for AA in Risbey et al. (2021) and Chapman et al. (2022), where archetypes
dominance and persistence are illustrated. This score is a measure how
representative archetypes are in any given data record $t$, ranging from 0 –
for low discrimination or a ‘neutral’ ENSO state, the record corresponding
then to a blend of archetypes with the equal probability of expression ($1/p$
for p archetypes) – to 1 – for high discrimination, the record being expressed
by a single archetype $i_{max}$ or a ‘pure’ ENSO state, with (maximum)
probability of expression equal to 1. As with any neutral condition
identification, thresholds selection based on expert consensus will have to be
tested and implemented.

Together with a persistence criterion, neutral ENSO states could be defined as
period of low archetype discrimination below a certain threshold over a
certain length of time. However, we note that the power discrimination across
archetypes from cardinalities ranging from 2 to 20 increases monotonically
with cardinality as illustrated in Figure S9 of the supplementary material,
making what one considers ‘neutral’ somewhat moot.

##  6 Traditional approaches and other clustering methods

Several studies (Capotondi et al., 2015, 2020; Wiedermann et al., 2021; Schlör
et al., 2023, for recent examples) have compared different approaches in the
characterisation and identification of ENSO phases based on SSTAs. These
approaches can be broadly categorised as methods based on 1) principal
component or multiple correspondence analysis, as in Chiang and Vimont (2004),
Takahashi et al. (2011), and Feng et al. (2020), 2) composites based on region
of influence or box-type indices, as in Ashok et al. (2007), Kug et al.
(2009), Karnauskas (2013), Cai and Cowan (2009), and Sullivan et al. (2016),
and 3) pattern recognition, as in Johnson (2013), Black et al. (2022),
Wiedermann et al. (2021) and Schlör et al. (2023).

Table 3: Adapted from Schlör et al. (2023), classification of El Niño events over the 1982-2022 DJF interval. The approaches of Kug et al. (2009), Takahashi et al. (2011), Schlör et al. (2023), and AA are compared. Similar to Schlör et al. (2023), AA categorize El Niño events based on the archetype corresponding to the highest affiliation probabilities for cardinalities 4 and 8, with archetypes $A_{1,4,8}$ standing for CP Niño and $A_{5}$ for EP Niño, as in Figures S7 and S8 for cardinality 8. The EP Niño archetype, $A_{5}$ for cardinality 8, corresponds to the Extreme El Niños type, EP*, identified by Schlör et al. (2023). The last column indicates when AA8 aligns with the consensus characterisation of El Niño events across the 3 other methods listed. Niño | Kug et al. | Takahashi et al. | Schlör et al. | AA4 | AA8 | Consensus
---|---|---|---|---|---|---
Events | 2009 | 2011 | 2023 |  |  |
1982–1983 | EP | EP | EP* | EP | EP | Yes
1986–1987 | EP | EP | EP | CP/EP | CP | No
1987–1988 | CP | CP | EP/CP | EP | CP | Yes
1991–1992 | EP | CP | EP | CP/EP | CP/EP | Yes
1994–1995 | CP | CP | EP | CP | CP | Yes
1997–1998 | EP | EP | EP* | EP | EP | Yes
2002–2003 | CP | CP | CP/EP | CP/EP | CP | Yes
2004–2005 | CP | CP | CP | CP | CP | Yes
2006–2007 | CP | CP | EP | EP | CP | Yes
2009–2010 | CP | CP | CP/EP | CP/EP | CP | Yes
2014–2015 | CP | - | CP | CP | CP | Yes
2015–2016 | EP | CP | EP* | EP | EP | Yes
2018–2019 | CP | CP | EP/CP | CP/EP | CP | Yes
2019-2020 | CP | CP | CP | CP | CP | Yes

Table 3 compares the classification of El Niño conditions over the satellite
era 1982 to 2020 for the region of influence, Niño 3 – Niño 4, of Kug et al.
(2009), the original EOF/PC based approach of Takahashi et al. (2011), the
unsupervised pattern recognition methods of Schlör et al. (2023) and AA for
cardinalities 4 and 8. Both unsupervised pattern recognition methods have been
‘rendered’ categorical for comparison by only listing the pattern or archetype
corresponding to the highest and second highest affiliation probabilities for
Schlör et al. (2023) and the ‘winning’ archetype based on the discrimination
scores $\Delta^{4}_{i_{max}}$ and $\Delta^{8}_{i_{max}}$, i.e. the archetype
contributing to any given record with the highest probability of expression.
The primary aim of using global SSTAs is not only to categorise the ENSO
phenomenon in the Tropics, but also via its impacts through teleconnections.
Although, the choice of going global with AA may incur the risk of ‘blending’
different processes and be potentially prone to confounding effects, it leads
to similar types when compared to other categorization methods and associated
indices, which are all based on SSTAs in the tropical and sub-tropical
Pacific. These indices are derived from separated, but at time overlapping,
domains and linearly combined to form EP- and CP-like ENSO indices. It could
be argued that such combinations are not necessarily justifiable from a
dynamical standpoint as they mix values, coincident in time, from spatially
separated regions. Climate and weather research practitioners also rely on
prior knowledge gathered from historical records for the choice and extent of
the domains used to construct these metrics; this corresponds to a form of
supervised learning. AA does not suffer from these drawbacks as the resulting
patterns are close to the observed (global) dSSTA snapshots without any non-
local manipulation needed to discriminate between CP and EP conditions or a
priori choices of domains. It is therefore expected that noticeable
differences may emerge when comparing AA categorical ENSO phase affiliation to
other methods. The non-local manipulations will also impact teleconnection
patterns and associated dynamical inferences derived typically from
correlation maps, normalised expectations in statistical terms, where the
contributing records of both the ENSO indices and geophysical fields under
investigation are equally weighted. This is not the case with AA composites,
where teleconnection patterns of any fields built on AA indices (S-matrix),
correspond to weighted averages, or conditional expectations, of field records
multiplied by the probabilities of expression of the associated dSSTA
archetypes as per Equation 3.

Interestingly, the majority of the comparative studies mentioned earlier only
focus on differentiating between El Niño conditions and rarely attempt to
classify La Niña flavours, with Cai and Cowan (2009) and Black et al. (2022)
being notable exceptions. Results discussed in Cai and Cowan (2009) highlight
the importance of CP-Niña or Modoki conditions affecting Australian rainfall.
Black et al. (2022) work reveals the ENSO teleconnection global imprint on
surface and atmospheric fields, unequivocally showing in particular that CP-
Niña and EP-Niño conditions strongly influence the circulation at mid- and
high-latitude as identified by large amplitude composites when compared to
opposite CP-Niño and EP-Niña phases. Using a reduced geographical domain,
Chapman et al. (2022) recently illustrate the value of AA methodology in
potentially disentangling the impacts of local and remote large-scale drivers
when applied to marine heatwaves (MHWs) and cold spells (MCSs) in the
Australasian region by comparing AA to the conventional geographical point-
based methods described in Oliver et al. (2021). They show that AA
unambiguously identifies concurrent CP- and EP-Niña conditions associated with
MHWs in South-Eastern Indian Ocean of the Western Australian coast and in
South Pacific close to New Zealand, respectively. As with Niña flavours, the
different Niño types have varied teleconnection impacts. These will be further
examined in Section 8.

##  7 AA on-line updates

AA can lend itself to the characterisation of forecasted ENSO conditions and
corresponding teleconnections in S2S or multi-year prediction systems by first
training AA on observed or reanalysed historical global SST patterns for a
given cardinality $p$ as follows. For $M$ historical time records, AA
solutions for the dataset X, here dSSTAs, are first computed,

$X_{st}\approx\sum_{r=1}^{M}\sum_{i=1}^{p}X_{sr}C_{ri}S_{it}=\sum_{j=1}^{p}XC_{sj}S_{jt},\forall
t\in 1,\ldots,M$

where $\textbf{XC}=\sum_{r=1}^{M}X_{sr}C_{rj}$ are the previously found
archetypal patterns. The on-line affiliation probability updates,
$\hat{S}_{it^{\prime}},\,\forall t^{\prime}\in M+1,\ldots,M+n$ for $n$ new
forecasted records of $Y_{st^{\prime}}$, are then approximated by the
solutions of

$\underset{\hat{\textbf{S}}}{\arg\min}\|Y_{st^{\prime}}-\sum_{i=1}^{p}XC_{si}\hat{S}_{it^{\prime}}\|_{F},$
(7)

where $p$ stands for the archetype cardinality considered a priori. The
constraints, $\sum_{j=1}^{p}\hat{S}_{jt^{\prime}}=\textbf{1}_{t^{\prime}}$ and
$\hat{\textbf{S}}\geq 0$, still apply. By modifying Mørup and Hansen (2012)
algorithm, Expression 1 reduces to the simpler least squares problem with
convexity constraints of finding the minimum for Expression 7 as a function of
$\hat{\textbf{S}}$ for given archetypes computed over the training, here XC,
and the forecasted records $Y_{st^{\prime}}$.

$\hat{\textbf{S}}_{4\times 12}$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
$\hat{S}_{1,\cdot}$ | 0.13 | 0.26 | 0.31 | 0.36 | 0.45 | 0.48 | 0.40 | 0.30 | 0.21 | 0.06 | 0.00 | 0.07
$\hat{S}_{2,\cdot}$ | 0.62 | 0.56 | 0.56 | 0.41 | 0.42 | 0.36 | 0.38 | 0.19 | 0.39 | 0.69 | 0.68 | 0.68
$\hat{S}_{3,\cdot}$ | 0.25 | 0.19 | 0.08 | 0.22 | 0.13 | 0.09 | 0.14 | 0.47 | 0.35 | 0.21 | 0.32 | 0.25
$\hat{S}_{4,\cdot}$ | 0.00 | 0.00 | 0.05 | 0.00 | 0.00 | 0.07 | 0.08 | 0.03 | 0.06 | 0.04 | 0.00 | 0.00
$\textbf{S}_{4\times 12}$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
$S_{1,\cdot}$ | 0.13 | 0.26 | 0.31 | 0.36 | 0.45 | 0.48 | 0.40 | 0.29 | 0.19 | 0.04 | 0.00 | 0.08
$S_{2,\cdot}$ | 0.65 | 0.57 | 0.57 | 0.43 | 0.44 | 0.39 | 0.40 | 0.22 | 0.44 | 0.78 | 0.75 | 0.77
$S_{3,\cdot}$ | 0.22 | 0.17 | 0.07 | 0.21 | 0.11 | 0.06 | 0.12 | 0.45 | 0.31 | 0.13 | 0.25 | 0.16
$S_{4,\cdot}$ | 0.00 | 0.00 | 0.05 | 0.00 | 0.00 | 0.07 | 0.08 | 0.03 | 0.06 | 0.04 | 0.00 | 0.00
$|\Delta\textbf{S}_{4\times 12}|$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
$|S_{1,\cdot}-\hat{S}_{1,\cdot}|$ | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.01 | 0.02 | 0.00 | 0.00
$|S_{2,\cdot}-\hat{S}_{2,\cdot}|$ | 0.03 | 0.02 | 0.01 | 0.02 | 0.02 | 0.03 | 0.03 | 0.03 | 0.05 | 0.09 | 0.07 | 0.09
$|S_{3,\cdot}-\hat{S}_{3,\cdot}|$ | 0.03 | 0.02 | 0.01 | 0.01 | 0.02 | 0.03 | 0.03 | 0.02 | 0.04 | 0.08 | 0.07 | 0.09
$|S_{4,\cdot}-\hat{S}_{4,\cdot}|$ | 0.00 | 0.02 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.00 | 0.00
Table 4: $|\Delta\textbf{S}_{4\times 12}|$ absolute differences for the 2021
monthly records between the ‘truth’ based on 1982-2021 AA,
$\textbf{S}_{4\times 12}$, and the ‘approximation’ update based on 1982-2020
AA, $\hat{\textbf{S}}_{4\times 12}$.

Table 4 compares absolute differences, $|\Delta\textbf{S}_{4\times 12}|$, for
the 2021 monthly records between the ‘truth’ based on 1982-2021 AA,
$\textbf{S}_{4\times 12}$, and the ‘approximated’ update based on 1982-2020
AA, $\hat{\textbf{S}}_{4\times 12}$. The 12 month updates for 2021, based on
1982-2020 archetypes, were computed with a single initialization trial and
application of the algorithm. The 2021 monthly records set aside for testing,
$Y_{st^{\prime}}$, are not used for training. The largest absolute mismatches
of the affiliation matrix S elements over the forecast period are about 0.09
and are mainly seen on affiliation time-series for CP- and EP-Niña
conditions212121As in Black et al. (2022) for AA based on dSSTAs over the
1982-2020 period, archetypes $A_{1}$ and $A_{2}$ correspond to CP-Niño and CP-
Niña, respectively., which are difficult to classify for cardinality 4 as
illustrated in Section 5. When applied to S2S forecasts, the AA on-line update
is likely to degrade as model representation and observations may diverge when
the forecast lead-time increases. The AA on-line update procedure presented
here is a proof-of-concept only. It requires thorough benchmarking for S2S
applications. Such an exercise would not only include characterising S2S
forecasts with AA, but also would test the ability of any S2S OAGCM to
represent ENSO, and is beyond the scope of this work. As an aside, we note
that AA could be added to the suite of diagnostics commonly used to check the
OAGCM representation fidelity of various climate modes across geophysical
realms, an approach illustrated recently by Chapman et al. (2022) for
2500-year-long control run of the GFDL Climate Model 2.1 (CM2.1) Coupled Model
Intercomparison Project (CMIP) phase 3 – class OAGCM, used in older versions
of the Australian Community Climate and Earth System Simulator (ACCESS).

##  8 Teleconnections

An indirect way to test the reliability and relevance of the AA classification
of ENSO types, is to ascertain that the extreme configurations detected in
dSSTAs correspond to robust and explainable teleconnection patterns of various
atmospheric and oceanic fields. In this section, we explore the archetypes
spatial and temporal imprints, and their linkages to extra-tropical
atmospheric circulation and rainfall. Composites of depth-integrated
quantities such as SLA and sub-surface ocean temperature anomaly profiles in
the Tropics will provide additional evidences of the realistic AA
representation of ENSO teleconnections. Only results pertaining to the AA
cardinality of 4, will be shown, with signal-to-noise ratios for individual
fields and results for cardinality 8 given in the supplementary material.

Figure 10: AA seasonal lag-0 composites (columns) of GPCP precipitation
centered undeciles (shaded) overlaid with anomalies of velocity potential
differences between isobaric pressure levels of 150 and 850 hPa (contours)
based on time averages of the S-affiliation matrix over the 1982-2022 period
for a cardinality of 4. Each individual subplot is scaled to increase contrast
across archetypal patterns (rows) and seasons (column). The scaling factors
are given at the bottom of each panel on the right for rainfall undeciles and
on the left for velocity potential anomalies. Green (brown) shadings
correspond to enhanced (reduced) precipitation above (below) the median
rainfall in white. Dashed (full) contours correspond to enhanced (reduced)
convection as diagnosed by $\Delta\chi_{150-850}$. Signal-to-noise ratios
$(F^{2}_{S}/C^{S}_{F})^{1/2}$ for both fields are reported in Figures S13 and
S14 of the supplementary material.

To begin with, ENSO conditions lead to enhanced or reduced atmospheric
convective activity over the tropical and associated convergence zones, not
only modulating rainfall over these regions, but also interacting with other
modes of variability such as the Madden-Julian Oscillation (MJO) (Madden and
Julian, 1972), for example. Figure 10 displays the lag-0 seasonal imprints of
ENSO archetypes on monthly GPCP precipitation centered undeciles222222Centered
undeciles or elftiles label the rank order of values of a distribution
partitioned in 11 frequency classes of equal width measured with respect to
‘median’ class 6, hence ranging from $-5$ to $+5$. (shaded) overlaid with
anomalies of JRA-55 velocity potential difference (contours) between isobaric
pressure levels of 150 and 850 hectoPascal (hPa), $\Delta\chi_{150-850}$,
based on time averages of the S-affiliation matrix over the 1982-2022 period
for a cardinality of 4. As reported by Adames and Wallace (2014) Table 2, a
MJO index based on the first two PCs of $\Delta\chi_{150-850}$ explains 57% of
its total variance is not only an excellent proxy for probing the propagation
of tropical Kelvin and Rossby-waves and their impact on extra-tropical
circulation, but also marks regions of enhanced and reduced convective
activity. Based on these findings, we are using the velocity potential
difference $\Delta\chi_{150-850}$ anomalies to diagnose the imprint of ENSO
phases on vertical atmospheric convection. Overall, both the precipitation and
the velocity potential composites agree well with the expected behaviour
associated with SST anomalies driven by ENSO. EP- and CP-Niño phases.
Associated SST signatures are characterised by enhanced precipitation (in
blue) in the Eastern and Central Pacific, respectively, collocated with
enhanced upward convection marked by dashed contours on Figure 10. On the
western side of the Pacific Basin, drier conditions (in red) prevail and the
convective activity reduces as shown by full contours in the
$\Delta\chi_{150-850}$ composites. The situation is reversed for EP- and CP-
Niña phases. Marked seasonal differences are also observed in both fields, but
especially in the $\Delta\chi_{150-850}$ composites, where the enhanced and
reduced convective center of action locations change with latitudes and
longitudes. Clear extra-tropical signatures along the Pacific Rim impacting
the North- and South-American, the Australian and Far-East Asian continents,
extending at times to the Indian sub-continent, the African and the European
side of the Eurasian continents, are apparent for CP-Niño and EP-Niña
composites. The corresponding composites from cardinality 8 are reproduced in
the supplementary material, Figure S9.

Figure 11: AA seasonal lag-0 composites (columns) of detrended sea-level
anomalies and signal-to-noise ratios ($\textbf{+}>0.15,\bullet>0.3$,
$(F^{2}_{S}/C^{S}_{F})^{1/2}$) based on time averages of the S-affiliation
matrix over the 1993-2021 period for a cardinality of 4. Each individual
subplot is scaled to increase contrast across archetypal patterns (rows). The
maximum absolute value of the detrended sea-level anomaly composites in a $\pm
20^{\circ}$ tropical band is reported in meter on the Eurasian continent for
each panel.

In the tropics, detrended sea-level anomalies (dSLAs) provide a first-order
approximation for (detrended) steric height. As shown by Gasparin and Roemmich
(2017) and Izumo et al. (2019), the dSLA signal is mainly driven by sub-
surface anomalies in the thermohaline expansion of seawater integrated from
the surface to below the thermocline down to about a depth of 300 meters. The
seasonal dSLA S-composites displayed in Figure 11 discriminate between CP and
EP ENSO sub-surface expressions, not only along the Equator, but also in the
inter-tropical convergence and South-Pacific convergence zones (ITCZ and
SPCZ), and even further afield, in the tropical Indian Ocean basin. The
largest anomalies are observed for the CP-Niña and EP-Niño conditions peaking
in DJF, with values twice as high as their CP-Niño and EP-Niña counterparts,
with their maximum values now expressed in MAM and JJA, respectively. The
behaviour is confirmed by sub-surface SODA profiles in potential temperature,
$\theta$, averaged 1) in latitudes over the $\pm 5\degree$ tropical band,
Figure 12, and 2) in longitudes over the Niño 3.4, 190$\degree$\- to
240$\degree$-East, bounds, Figure 13. The signal-to-noise ratios between the
composite means and standard deviations,
$(\theta^{2}_{S}/C^{S}_{\theta})^{1/2}$, are 2 or 3 time larger than the GPCP
precipitation composite ratios displayed on Figures S13 and S18, for example,
and provide an indication of the strength of the influence of ENSO cycles (as
determined by the AA metric), especially between 50- to 250-meter depth. For
EP-Niño, the strong anomalies in ocean temperature reach all the way to the
surface and have the largest signal-to-noise ($>$0.6) across all archetypes in
the Niño 3.4 longitude band displayed in Figure 12.

Figure 12: Detrended anomalies of temperature AA lag-0 composites longitudinal
profiles (left) from SODA v3.15.2 reanalysis averaged over the $-5^{\circ}$ to
$5^{\circ}$ latitude band (dashed lines) and corresponding detrended SSTA
archetypes (right) for a cardinality of 4. The shadings on the left correspond
to time-mean temperature composites and the contours to the (time) mean to
standard deviation ratio $(F^{2}_{S}/C^{S}_{F})^{1/2}$ in 0.2 increment and
the green values in the lower right corners are the maximum absolute
temperature anomalies in ∘Celsius. Each individual subplot is scaled by its
maximum amplitude to increase contrast across archetypal patterns (rows).
Figure 13: Detrended anomalies of temperature AA lag-0 composites latitudinal
profiles (left) from SODA v3.15.2 reanalysis averaged over the 190$\degree$ to
240$\degree$ (Niño 3.4 longitudes) longitude band (dashed lines) and
corresponding detrended SSTA archetypes (right) for a cardinality of 4. The
shadings on the left correspond to time-mean temperature composites and the
contours to the (time) mean to standard deviation ratio
$(F^{2}_{S}/C^{S}_{F})^{1/2}$ in 0.2 increment and the green values in the
lower right corners are the maximum absolute temperature anomalies in
∘Celsius. Each individual subplot is scaled by its maximum amplitude to
increase contrast across archetypal patterns (rows).

The larger signal-to-noise results, obtained for both sea-level height and
potential temperature composites, indicate that AA based on ocean fields is
likely to be less prone to atmospheric weather noise and, therefore, possibly
more suited to the detection and classification of ENSO conditions.

To test this hypothesis, we performed AA on detrended SLAs directly. Figure 14
displays AA spatial pattern and time series results using detrended monthly
SLA anomalies over 1982-2021 for cardinality 4. The first two columns on the
left show archetypes and dSLA composites, XC and
$\textbf{X}\overline{\textbf{S}}^{T}$, with
$\overline{S}_{pt}=S_{pt}/\sum^{t}_{i=1}S_{pi}$. The archetypal patterns are
ranked, from top-to-bottom, according to the time mean of S,
$\frac{1}{N}\sum^{t}_{i=1}S_{pi},\forall p$, the mean probability of
expression of archetype $p$ over all $N$ records. The last two columns on the
right show C and S matrix time series, respectively, with MEI and ONI time
series overlaid on both. The value in bold, third column top row, corresponds
to the maximum C-weights across all archetypes. A cursory comparison with
Figure 11 indicates that the seasonal dSLA composites at 0-lag built on dSSTA
AA affiliation weights resemble the dSLA archetypes. The time-mean affiliation
weights, $\frac{1}{N}\sum^{t}_{i=1}S^{SLA}_{pi}$ for dSLA and
$\frac{1}{N}\sum^{t}_{i=1}S^{SST}_{pi}$ for dSSTA, are nevertheless different,
so archetype patterns do not follow the same top-to-bottom row order in
Figures 11 and 14.

Figure 14: AA spatial pattern and time series results using detrended monthly
SLA anomalies over 1982-2021 for cardinality 4. The two left columns show
archetypes and dSLA composites, XC and $\textbf{X}\overline{\textbf{S}}^{T}$,
with $\overline{S}_{p\times t}=S_{p\times t}/\sum^{t}_{i=1}S_{pi}$,
respectively. The two right columns show C and S matrix time series,
respectively, with MEI (black) and ONI (grey) time series overlaid on both.
The values in the upper left corner of each panel, in the first two columns,
correspond to the maximum dSLA anomalies range scale expressed in meters,
whereas the value in bold in third column first row is the maximum C-weights
across all archetypes. Both MEI (black) and ONI (grey) time series were scaled
to fit within the C and S weights bounds.

Figures 15a and 15b try to match the affiliation sequences for both dSSTAs and
dSLAs AA based on the ‘nestedness’ information displayed in Figures S15 and
S16. However, a perfect match between dSSTA and dSLA archetypes is not to be
expected. The main differences are found in the archetype persistence across
time records expressed by the discrimination scores, $\Delta^{8}_{i_{max}}$.
In any given phases, the scores are more elevated and persist over longer
periods for dSLA archetypes when compared to those of dSSTA. The dSLA score
persistence indicates 1) that sub-surface quantities are less exposed to the
vagaries of the weather and associated drivers over both the Tropics and Sub-
Tropics, and 2) that the ENSO phenomenon operates over longer timescales than
evident from surface-based representations. There are clear periods where the
ENSO state is poorly expressed in dSSTA, but well expressed in dSLA in Figures
15a and 15b. For example, the period 2001–2004 where there is a mix of
archetypes expressed in dSSTA but a well expressed preference for CP Nino in
dSLA. More generally, when the MEI and ONI indices have small amplitude
variation, dSSTA ENSO types are poorly expressed, but there is still a well
expressed dSLA type.

Although all flavours of El Niño may have similar early subsurface origins at
long lead times as claimed by Ramesh and Murtugudde (2013), the extreme
configurations revealed by AA indices in both S time-series on Figure 15 and
in sub-surface composites at 0-lag on Figures S11, S12 and S13, not only point
to significantly different expressions across ENSO types, but also to fast –
over 6 to 12 months – transitions for both the initiation and termination of
events across all flavours. These extreme configurations are interspersed by
longer lasting CP Niño and Niña intervals, especially evident in the sub-
surface.

Many different extra-tropical modes have been hypothesized as potential ENSO
precursors, both in the northern and southern hemispheres such as in Yu et al.
(2010), Wang et al. (2012) and Zhang et al. (2014). A statistical relationship
and the underlying mechanisms are very difficult to separate from each other
and from ENSO itself, as all modes are correlated and not independent one from
another, as nicely illustrated by Pegion et al. (2020) for example. Pegion et
al. (2020) and Capotondi and Ricciardulli (2021) work supports the global
approach adopted in this study. By extending our analysis domain to global
dSSTAs, we attempt to follow their lead, not only including the extra-tropics,
but also the Indian and Atlantic Ocean basins, where extreme dSSTA
configurations do have an imprint. Amongst other influencing factors, the so-
called ‘capacitor effect’ of the Indian and Atlantic Oceans discussed in
Okumura and Deser (2010), An and Kim (2018) and An et al. (2020), where the
ENSO driven anomalies imparted to these basins via the atmospheric bridge are
imprinted and eventually feeding back to the Pacific, may play a role.

Figure 15: Comparison of AA affiliation sequences for a) dSSTAs and b) dSLAs.
The discrimination scores (Risbey et al. (2021)) $\Delta^{p}_{i_{max}}(t)\in$
[0, 1], $\,\forall i=1,\ldots,8$ for all monthly records $t$ are represented
by archetypal colour-coded dots of the winning archetype. Here, the color
codes for dSLA archetypes are tentatively matched to those of dSSTA
archetypes. The corresponding archetypal patterns are given in Figures S7,
columns 3 and 4 for dSSTAs, and S13, columns 1 and 2 for dSLAs. Both MEI
(black) and ONI (grey) time series were scaled to fit within the $S_{8\times
492}$ and $\Delta_{8}$ bounds.

The results presented in Figures 15a and 15b not only reaffirm and extend the
earlier observations of Wyrtki (1975) that no ENSO events are quite the same
across geophysical realms, but also show for a given cardinality that dSLA
archetypes provide a clearer representation of ENSO types and that they are
less prone to mischaracterisation than dSSTA archetypes. For both dSSTA and
dSLA archetypes, increasing the AA cardinality adds variety to the phenomenon
and sharpens the discrimination across ENSO types.

##  9 ENSO phase transitions

Although hampered by the short sample size of the well-observed SST records
over the satellite area, the phase evolution of ENSO can be studied by
counting transitions from one well-expressed archetype to the next. Tables 5
and 6 tally 491 monthly transitions, including transitions to oneself,
$\pi_{ij}=\frac{1}{N-1}\sum_{t=1}^{N-1}\\{\underset{i=1,\ldots,4/8}{\arg\max}\,S_{it}\longrightarrow\underset{j=1,\ldots,4/8}{\arg\max}\,S_{jt+1}\\}$,
across the categorical AA representation of dSSTAs for cardinalities 4 and 8.
$N=492$ stands for the total number of months in the 41-year long record and
$\pi_{ij}$ the time-mean transition probability from archetypes $i$ to $j$.
The diagonal entries in both tables correspond to the number of self-
transitions. The individual phase mean persistence can be easily computed as
the diagonal entries divided by the number of separate intervals of a given
phase. For a cardinality of 8 for example, the EP-Niño intervals, labeled
$A_{5}$ in Table 6, persist in average for about $33/4\approx 7$ to 8 months,
where 4 corresponds to the number of separate EP-Niño intervals as shown in
Figure 8a.

| $\textbf{A}_{1}$ | $\textbf{A}_{2}$ | $\textbf{A}_{3}$ | $\textbf{A}_{4}$ | Phase
---|---|---|---|---|---
$\textbf{A}_{1}$ | 129 | 8 | 7 | 1 | CP Niña
$\textbf{A}_{2}$ | 3 | 143 | 8 | 8 | CP Niño
$\textbf{A}_{3}$ | 12 | 8 | 95 | 3 | EP Niña
$\textbf{A}_{4}$ | 2 | 3 | 7 | 54 | EP Niño
Table 5: Number of the transitions from winning $i$ (row) to winning $j$ (column) archetypes, $(N-1)\pi_{ij}=\sum_{t=1}^{N-1}\\{\underset{i=1,\ldots,4}{\arg\max}\,S_{it}\longrightarrow\underset{j=1,\ldots,4}{\arg\max}\,S_{jt+1}\\}$, where $N=492$ stands for the total number of months in the 41-year long record and $\pi_{ij}$ the transition probability from archetype $i$ to $j$. The table diagonal elements correspond here to the numbers of self transitions omitted in Figure 16. The winning archetypes are the archetypes contributing to any given record with the highest probability of expression, $\underset{i=1,\ldots,4}{\arg\max}\,S_{it},\,\forall t$. | $\textbf{A}_{1}$ | $\textbf{A}_{2}$ | $\textbf{A}_{3}$ | $\textbf{A}_{4}$ | $\textbf{A}_{5}$ | $\textbf{A}_{6}$ | $\textbf{A}_{7}$ | $\textbf{A}_{8}$ | Phase
---|---|---|---|---|---|---|---|---|---
$\textbf{A}_{1}$ | 89 | 4 | 1 | 2 | 3 | 2 | 1 | 3 | CP Niño
$\textbf{A}_{2}$ | 6 | 83 | 3 | 6 | 0 | 0 | 1 | 0 | –
$\textbf{A}_{3}$ | 0 | 5 | 62 | 2 | 0 | 4 | 2 | 1 | CP Niña
$\textbf{A}_{4}$ | 6 | 5 | 0 | 58 | 2 | 0 | 2 | 0 | –
$\textbf{A}_{5}$ | 1 | 0 | 2 | 0 | 33 | 1 | 1 | 0 | EP Niño
$\textbf{A}_{6}$ | 3 | 1 | 3 | 0 | 0 | 33 | 0 | 1 | –
$\textbf{A}_{7}$ | 0 | 1 | 3 | 4 | 0 | 0 | 20 | 2 | EP Niña
$\textbf{A}_{8}$ | 0 | 1 | 2 | 0 | 0 | 1 | 3 | 22 | –
Table 6: Number of the transitions from winning $i$ (row) to winning $j$
(column) archetypes,
$(N-1)\pi_{ij}=\sum_{t=1}^{N-1}\\{\underset{i=1,\ldots,8}{\arg\max}\,S_{it}\longrightarrow\underset{j=1,\ldots,8}{\arg\max}\,S_{jt+1}\\}$,
where $N$ stands for the number of months in the 41-year long record. The
table diagonal elements correspond here to the numbers of self transitions
omitted in Figure S14. Figure 16: The directed graph representation of time-
mean states transition probabilities from winning to winning archetypes,
$\pi_{ij}=\frac{1}{N-1}\sum_{t=1}^{N-1}\\{\underset{i=1,\ldots,4}{\arg\max}\,S_{it}\longrightarrow\underset{j=1,\ldots,4}{\arg\max}\,S_{jt+1}\\}$,
where $N$ is the total number of monthly records. The arrows thicknesses are
proportional to the transition probabilities. Self transitions, $\pi_{ii}$
with $i=1,\ldots,4$, have been omitted. The colour-coded graph nodes
correspond to the ENSO phase colours in Fig. 4 and 5.

In Figures 16 and S14, the time-mean probability of transitions between ENSO
phases is represented by directed graphs illustrating the transitions from
‘winning’ to ‘winning’ archetypes,
$\pi_{ij}=\frac{1}{N-1}\sum_{t=1}^{N-1}\\{\underset{i=1,\ldots,4}{\arg\max}\,S_{it}\longrightarrow\underset{j=1,\ldots,4}{\arg\max}\,S_{jt+1}\\}$,
where $N$ is the total number of monthly records. The arrows thicknesses are
proportional to the transition probabilities. Self transitions, $\pi_{ii}$
with $i=1,\ldots,4$, have been omitted given the persistence of ENSO
conditions. The colour-coded graph nodes correspond to the ENSO phase colours
in Figures 4 and 5. Such a representation suggests that AA could be used as a
dimension reduction strategy, whereby the evolution of the ENSO phenomenon
(and its teleconnection impacts) is reduced to and represented by the S-matrix
time-series alone, our AA ENSO indices. Based on these indices, machine
learning methods such long short-term memory (LSTM) deep learning methods
(Hochreiter and Schmidhuber, 1997; Chollet, 2021) or Hidden Markov Models
(HMM) (Baum and Petrie, 1966; Rabiner, 1989) could be deployed to forecast
ENSO, keeping in mind that the power232323Here, ‘power’ is understood as the
predictive accuracy of the statistical algorithms employed as in Breiman
(2001). of any statistical forecasting method is contingent on the number of
high-quality historical records available for training. In principle, the
approach could be tested on long OAGCM simulations for models with an adequate
ENSO representation.

##  10 Summary and conclusion

Since the discovery of the Southern Oscillation by Sir Gilbert Walker, our
understanding of the ENSO has advanced in leaps and bounds. As illustrated in
the introduction, every unusual event departing from the accepted historical
paradigm at the time, triggered a surge of activity to further research, and
lead to an extension of ocean-atmosphere observing systems in order to better
capture the intricacies of the phenomenon given its global environmental and
societal impacts. These successive changes in the understanding of ENSO and,
not only associated with, but also driving, improvements in spatiotemporal
observational coverage in both the atmosphere and oceans, all contributed to
the current state of ENSO research. To date, ENSO events are commonly
categorized by anomalous patterns of SST and their centers of action (regions
of absolute anomaly maxima) in the tropical Pacific Ocean; these typically
correspond to the central and eastern Pacific El Niño and La Niña conditions.

Built on recently published work by Hannachi and Trendafilov (2017), Black et
al. (2022) and Chapman et al. (2022), we show that the approach based on
archetypal analysis in Sections 4 and 5 and is able to recover the CP- and EP-
like phases of ENSO using detrended global SST anomalies. The favorable
comparison of AA to the categorical and soft clustering methods, K-means and
fuzzy-C-means, underlines the robustness of the method.

The choice of going global for the analysis domain is intentional as both
tropical and extra-tropical information allows to capture the essence of the
phenomenon and, importantly, its remote impact. We posit in Section 7 that AA,
when applied to S2S forecasts, may complement and possibly sharpen the advice
given by operational centers by further characterising ENSO teleconnections.
Given the ‘nestedness’ properties as a function of AA cardinality, extra
flavours can be hierarchically introduced to study onset and decay of
different phases, without changing the archetypal patterns uncovered for lower
cardinalities as shown for dSSTAs and dSLAs in the supplementary material S15
and S16; the higher the cardinality, the more discriminatory AA becomes, as
illustrated in Figure S17 for dSSTAs.

The archetype affiliation probabilities, the S-matrix, associated with each
SST records, allow forecasting centers to set explicit criteria for onset,
persistence, and decay of events. They correspond to the archetypal flavours
and ENSO indices. We consider the probabilistic nature of the AA factorisation
as a strong point. It allows the identification of teleconnections based on
conditional expectations, in both mean and variance, by deferentially weighing
the records of any geophysical fields of interest as a function of how well-
expressed dSSTA global configurations are.

The non-stationarity in event sequences can be captured by changes in event
probabilities of expression throughout the period under study as illustrated
by Figures 8b, 15a and 15b, showing that no two identical sequences of ENSO
events are observed over the 1982 to 2022 period, even after removing a linear
trend at every point approximately representing anthropogenic forcing. The
classification of events via AA is reasonably robust to changes in the
selection of time periods and geographical domains. The identification of new
archetypes, where and when the dynamics or the underlying observing systems
have changed, allows for data-driven regime shift detection and
identification, either real or spurious, as AA captures to extreme
configurations. AA performed on ocean sub-surface fields, such as dSLAs,
provides a clearer expression of ENSO types, especially when the
categorisation based on surface fields lacks discrimination.

Finally, when considered as a dimension reduction strategy in Section 9, AA
results, the affiliation probability matrix S in particular, can be used as
input to statistical prediction models applied to time-series forecasting to
speed up deep learning algorithms.

###### Acknowledgements.

This work is supported the Australian Climate Service, a partnership of the
Australian Bureau of Meteorology, Geoscience Australia, Commonwealth
Scientific and Industrial Research Organisation (CSIRO) and Australian Bureau
of Statistics (https://www.acs.gov.au). The collaboration between Abdelwaheb
Hannachi and Didier Monselesan was fostered by the International
Meteorological Institute at the University of Stockholm through its visiting
scientist program and Didier Monselesan would like to thank Henriette Valsö at
Stockholm University – Department of Meteorology and Leesa MacKay at the CSIRO
– Environment, for facilitating international travel to Sweden. Sophie
Cravatte, Bastien Pagli and Takashi Izumo contributions to the paper were
supported in part by the Franco-Australian Hubert Curien Program (FASIC).
FASIC is implemented jointly by the Ministry of Europe and Foreign Affairs
(MEAE) and the Ministry of Higher Education and Research (MESR) in France, and
by research partners in Australia. Christopher Chapman was supported by the
CSIRO Environment Climate Atmosphere Oceans Interaction Program. All coauthors
have contributed to the work since its inception. No artificial intelligence
(AI) chat-bots were deployed to either conceptualise, code or generate content
for this work. Sea surface temperature (SST) data are from the Optimum
Interpolation Sea Surface Temperature (OISST) v2.1 high resolution dataset
provided by the NOAA/OAR/ESRL Physical Sciences Laboratory. This data is
available at https://psl.noaa.gov. The multi-variate ENSO index (MEI) was
retrieved from https://psl.noaa.gov/enso/mei. The Oceanic Niño Index (ONI)
from the Climate Prediction Center (CPC) is provided by NOAA/PSL downloaded
from http://www.cpc.ncep.noaa.gov/data/indices/oni.ascii.txt. Global
atmospheric reanalysis data used to relate extreme events to large-scale
climate modes comes from the Japanese 55-year Reanalysis (JRA-55) project
carried out by the Japan Meteorological Agency (JMA). JRA-55 data is available
at https://jra.kishou.go.jp. Global precipitation composites use Global
Precipitation Climatology Project (GPCP) Climate Data Record (CDR), Version
2.3 (Monthly), from the National Centers for Environmental Information,
doi:10.7289/V56971M6, whereas Australian precipitation composites are
generated from the Australian Gridded Climate Data (AGCD), the Australian
Bureau of Meteorology’s official dataset for climate analyses covering the
variables of rainfall, temperature (maximum and minimum) as well as vapour
pressure at daily and monthly timescales accessible at
https://geonetwork.nci.org.au/geonetwork/srv/eng/catalog.search#/metadata/f7421_3667_8543_6351.
The Ssalto/Duacs altimeter products were produced and distributed by the
Copernicus Marine and Environment Monitoring Service (CMEMS), accessible at
http://www.marine.copernicus.eu. The Simple Ocean Data Assimilation (SODA)
dataset was retrieved from https://climatedataguide.ucar.edu/climate-
data/soda-simple-ocean-data-assimilation on 2023-05-23.

## References

  * Adames and Wallace (2014) Adames, A., and J. Wallace, 2014: Three-dimensional structure and evolution of the MJO and its relation to the mean flow. J. Atmos. Sci., 71 (6), 2007–2026.
  * Adler et al. (2018) Adler, R., and Coauthors, 2018: The Global Precipitation Climatology Project (GPCP) Monthly Analysis (New Version 2.3) and a Review of 2017 Global Precipitation. Atmosphere, 9 (4), 138, 10.3390/atmos9040138.
  * Amaya (2019) Amaya, D. J., 2019: The Pacific Meridional Mode and ENSO: A Review. Current Climate Change Reports, 5 (4), 296–307, 10.1007/s40641-019-00142-x.
  * Amaya et al. (2020) Amaya, D. J., A. J. Miller, S.-P. Xie, and Y. Kosaka, 2020: Physical drivers of the summer 2019 north pacific marine heatwave. 11 (1), 1903, 10.1038/s41467-020-15820-w, URL https://www.nature.com/articles/s41467-020-15820-w, number: 1 Publisher: Nature Publishing Group.
  * An and Kim (2018) An, S.-I., and J.-W. Kim, 2018: ENSO Transition Asymmetry: Internal and External Causes and Intermodel Diversity. Geophysical Research Letters, 45 (10), 5095–5104, 10.1029/2018GL078476, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2018GL078476, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2018GL078476.
  * An et al. (2020) An, S.-I., E. Tziperman, Y. M. Okumura, and T. Li, 2020: ENSO Irregularity and Asymmetry. El Niño Southern Oscillation in a Changing Climate, American Geophysical Union (AGU), 153–172, 10.1002/9781119548164.ch7, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/9781119548164.ch7.
  * Ashok et al. (2007) Ashok, K., S. K. Behera, S. A. Rao, H. Weng, and T. Yamagata, 2007: El Niño Modoki and its possible teleconnection. Journal of Geophysical Research: Oceans, 112 (C11), 10.1029/2006JC003798.
  * Bach et al. (2019) Bach, E., S. Motesharrei, E. Kalnay, and A. Ruiz-Barradas, 2019: Local Atmosphere–Ocean Predictability: Dynamical Origins, Lead Times, and Seasonality. Journal of Climate, 32 (21), 7507–7519, 10.1175/JCLI-D-18-0817.1, URL https://journals.ametsoc.org/view/journals/clim/32/21/jcli-d-18-0817.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Barber et al. (1996) Barber, C. B., D. P. Dobkin, and H. Huhdanpaa, 1996: The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software, 22 (4), 469–483, 10.1145/235815.235821, URL https://dl.acm.org/doi/10.1145/235815.235821.
  * Barnett and Seth (2014) Barnett, L., and A. K. Seth, 2014: The MVGC multivariate Granger causality toolbox: A new approach to Granger-causal inference. Journal of Neuroscience Methods, 223, 50–68, 10.1016/j.jneumeth.2013.10.018, URL https://www.sciencedirect.com/science/article/pii/S0165027013003701.
  * Baum and Petrie (1966) Baum, L. E., and T. Petrie, 1966: Statistical Inference for Probabilistic Functions of Finite State Markov Chains. The Annals of Mathematical Statistics, 37 (6), 1554–1563, 10.1214/aoms/1177699147, URL https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-37/issue-6/Statistical-Inference-for-Probabilistic-Functions-of-Finite-State-Markov-Chains/10.1214/aoms/1177699147.full, publisher: Institute of Mathematical Statistics.
  * Bezdek et al. (1984) Bezdek, J. C., R. Ehrlich, and W. Full, 1984: FCM: The fuzzy c-means clustering algorithm. Computers & Geosciences, 10 (2-3), 191–203, 10.1016/0098-3004(84)90020-7, URL https://linkinghub.elsevier.com/retrieve/pii/0098300484900207.
  * Bjerknes (1966) Bjerknes, J., 1966: A possible response of the atmospheric hadley circulation to equatorial anomalies of ocean temperature. 18 (4), 820–829, 10.1111/j.2153-3490.1966.tb00303.x, URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1966.tb00303.x, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1966.tb00303.x.
  * Bjerknes (1969) Bjerknes, J., 1969: Atmospheric teleconnections from the equatorial Pacific. Mon. Weather Rev., 97 (3), 163–172.
  * Black et al. (2022) Black, A. S., and Coauthors, 2022: Archetypal Analysis of Geophysical Data Illustrated by Sea Surface Temperature. Artificial Intelligence for the Earth Systems, 1 (3), 10.1175/AIES-D-21-0007.1.
  * Breiman (2001) Breiman, L., 2001: Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author). Statistical Science, 16 (3), 199–231, 10.1214/ss/1009213726, URL https://projecteuclid.org/journals/statistical-science/volume-16/issue-3/Statistical-Modeling–The-Two-Cultures-with-comments-and-a/10.1214/ss/1009213726.full, publisher: Institute of Mathematical Statistics.
  * Bunge and Clarke (2009) Bunge, L., and A. J. Clarke, 2009: A verified estimation of the el niño index niño-3.4 since 1877. 22 (14), 3979–3992, 10.1175/2009JCLI2724.1, URL https://journals.ametsoc.org/view/journals/clim/22/14/2009jcli2724.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Cai and Cowan (2009) Cai, W., and T. Cowan, 2009: La Niña Modoki impacts Australia autumn rainfall variability. Geophysical Research Letters, 36 (12), 10.1029/2009GL037885, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2009GL037885.
  * Cane (1983) Cane, M. A., 1983: Oceanographic Events During El Niño. Science, 222 (4629), 1189–1195, 10.1126/science.222.4629.1189, URL https://www.science.org/doi/10.1126/science.222.4629.1189, publisher: American Association for the Advancement of Science.
  * Capotondi and Ricciardulli (2021) Capotondi, A., and L. Ricciardulli, 2021: The influence of pacific winds on ENSO diversity. 11 (1), 18 672, 10.1038/s41598-021-97963-4, URL https://www.nature.com/articles/s41598-021-97963-4, number: 1 Publisher: Nature Publishing Group.
  * Capotondi and Sardeshmukh (2017) Capotondi, A., and P. D. Sardeshmukh, 2017: Is El Niño really changing? Geophysical Research Letters, 44 (16), 8548–8556, https://doi.org/10.1002/2017GL074515, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL074515.
  * Capotondi et al. (2020) Capotondi, A., A. T. Wittenberg, J.-S. Kug, K. Takahashi, and M. J. McPhaden, 2020: ENSO Diversity. El Niño Southern Oscillation in a Changing Climate, American Geophysical Union (AGU), 65–86, 10.1002/9781119548164.ch4, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/9781119548164.ch4.
  * Capotondi et al. (2015) Capotondi, A., and Coauthors, 2015: Understanding ENSO Diversity. Bulletin of the American Meteorological Society, 96 (6), 921–938, 10.1175/BAMS-D-13-00117.1, URL https://journals.ametsoc.org/view/journals/bams/96/6/bams-d-13-00117.1.xml, publisher: American Meteorological Society Section: Bulletin of the American Meteorological Society.
  * Carton and Giese (2008) Carton, J. A., and B. S. Giese, 2008: A Reanalysis of Ocean Climate Using Simple Ocean Data Assimilation (SODA). Monthly Weather Review, 136 (8), 2999–3017, 10.1175/2007MWR1978.1.
  * Chapman et al. (2022) Chapman, C. C., D. P. Monselesan, J. S. Risbey, M. Feng, and B. M. Sloyan, 2022: A large-scale view of marine heatwaves revealed by archetype analysis. Nature Communications, 13 (1), 7843, 10.1038/s41467-022-35493-x, URL https://www.nature.com/articles/s41467-022-35493-x, number: 1 Publisher: Nature Publishing Group.
  * Chiang and Vimont (2004) Chiang, J. C. H., and D. J. Vimont, 2004: Analogous Pacific and Atlantic Meridional Modes of Tropical Atmosphere–Ocean Variability. Journal of Climate, 17 (21), 4143–4158, 10.1175/JCLI4953.1, URL https://journals.ametsoc.org/view/journals/clim/17/21/jcli4953.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Chollet (2021) Chollet, F., 2021: Deep learning with Python. Second edition ed., Manning Publications, Shelter Island, oCLC: on1289290141.
  * Christiansen (2007) Christiansen, B., 2007: Atmospheric Circulation Regimes: Can Cluster Analysis Provide the Number? Journal of Climate, 20 (10), 2229–2250, 10.1175/JCLI4107.1.
  * Cichocki et al. (2009) Cichocki, A., R. Zdunek, A. H. Phan, and S.-i. Amari, 2009: Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation. Wiley Publishing.
  * Cowtan and Way (2014) Cowtan, K., and R. Way, 2014: Coverage bias in the HadCRUT4 temperature series and its impact on recent temperature trends. Q. J. R. Meteor. Soc., 140 (683), 1935–1944.
  * Cutler and Breiman (1994) Cutler, A., and L. Breiman, 1994: Archetypal analysis. Technometrics, 36 (4), 338–347.
  * Dieppois et al. (2021) Dieppois, B., A. Capotondi, B. Pohl, K. P. Chun, P.-A. Monerie, and J. Eden, 2021: ENSO diversity shows robust decadal variations that must be captured for accurate future projections. Communications Earth & Environment, 2 (1), 1–13, 10.1038/s43247-021-00285-6, URL https://www.nature.com/articles/s43247-021-00285-6.
  * Dunn (1973) Dunn, J. C., 1973: A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters. 3 (3), 32–57, 10.1080/01969727308546046, URL https://doi.org/10.1080/01969727308546046, publisher: Taylor & Francis _eprint: https://doi.org/10.1080/01969727308546046.
  * Elad (2010) Elad, M., 2010: Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. 1st ed., Springer Publishing Company, Incorporated.
  * Elden (2019) Elden, L., 2019: Matrix Methods in Data Mining and Pattern Recognition. 2nd ed., SIAM-Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.
  * Emile-Geay et al. (2020) Emile-Geay, J., K. M. Cobb, J. E. Cole, M. Elliot, and F. Zhu, 2020: Past ENSO variability. El Niño Southern Oscillation in a Changing Climate, American Geophysical Union (AGU), 87–118, 10.1002/9781119548164.ch5, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119548164.ch5, section: 5 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119548164.ch5.
  * Feng et al. (2020) Feng, Y., X. Chen, and K.-K. Tung, 2020: ENSO diversity and the recent appearance of Central Pacific ENSO. Climate Dynamics, 54 (1), 413–433, 10.1007/s00382-019-05005-7, URL https://doi.org/10.1007/s00382-019-05005-7.
  * Forgy (1965) Forgy, E. W., 1965: Cluster analysis of multivariate data : efficiency versus interpretability of classifications. Biometrics, 21, 768–769.
  * Freund et al. (2019) Freund, M. B., B. J. Henley, D. J. Karoly, H. V. McGregor, N. J. Abram, and D. Dommenget, 2019: Higher frequency of Central Pacific El Niño events in recent decades relative to past centuries. Nature Geoscience, 12 (6), 450–455, 10.1038/s41561-019-0353-3.
  * Gan et al. (2020) Gan, G., C. Ma, and J. Wu, 2020: Data Clustering: Theory, Algorithms and Applications. 2nd ed., SIAM Press, URL https://my.siam.org/Store/Product/viewproduct/?ProductId=32863116.
  * Gasparin and Roemmich (2017) Gasparin, F., and D. Roemmich, 2017: The seasonal march of the equatorial Pacific upper-ocean and its El Niño variability. Progress in Oceanography, 156, 1–16, 10.1016/j.pocean.2017.05.010, URL https://www.sciencedirect.com/science/article/pii/S0079661116301628.
  * Gillis (2020) Gillis, N., 2020: Nonnegative Matrix Factorization. Society for Industrial and Applied Mathematics, Philadelphia, PA, 10.1137/1.9781611976410, URL https://epubs.siam.org/doi/abs/10.1137/1.9781611976410, https://epubs.siam.org/doi/pdf/10.1137/1.9781611976410.
  * Glantz and Thompson (1981) Glantz, M., and J. Thompson, Eds., 1981: Resource Management and Environmental Uncertainty: Lessons From Coastal Upwelling Fisheries. John Wiley & Sons, New York.
  * Glantz (2002) Glantz, M. H., Ed., 2002: La Niña and its impacts: facts and speculation. United Nations University Press, Tokyo ; New York.
  * Granger (1969) Granger, C. W. J., 1969: Investigating Causal Relations by Econometric Models and Cross-spectral Methods. Econometrica, 37 (3), 424–438, 10.2307/1912791, URL https://www.jstor.org/stable/1912791, publisher: [Wiley, Econometric Society].
  * Guilyardi et al. (2020) Guilyardi, E., A. Capotondi, M. Lengaigne, S. Thual, and A. T. Wittenberg, 2020: ENSO Modeling. El Niño Southern Oscillation in a Changing Climate, American Geophysical Union (AGU), 199–226, 10.1002/9781119548164.ch9, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119548164.ch9, section: 9 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119548164.ch9.
  * Han et al. (2020) Han, T., K. Yu, H. Yan, W. Jiang, H. Yan, and S. Tao, 2020: Coral d180-based reconstruction of el niño-southern oscillation from the northern south china sea since 1851 AD. 550, 159–168, 10.1016/j.quaint.2020.04.032, URL https://www.sciencedirect.com/science/article/pii/S104061822030197X.
  * Hannachi (2021) Hannachi, A., 2021: Further Topics. Patterns Identification and Data Mining in Weather and Climate, A. Hannachi, Ed., Springer Atmospheric Sciences, Springer International Publishing, Cham, 367–413, 10.1007/978-3-030-67073-3_16, URL https://doi.org/10.1007/978-3-030-67073-3˙16.
  * Hannachi et al. (2003) Hannachi, A., D. Stephenson, and K. Sperber, 2003: Probability-based methods for quantifying nonlinearity in the ENSO. 20 (2), 241–256, 10.1007/s00382-002-0263-7, URL https://doi.org/10.1007/s00382-002-0263-7.
  * Hannachi and Trendafilov (2017) Hannachi, A., and N. Trendafilov, 2017: Archetypal analysis: mining weather and climate extremes. J. Clim., 30 (17), 6927–6944.
  * Hirahara et al. (2014) Hirahara, S., M. Ishii, and Y. Fukuda, 2014: Centennial-scale sea surface temperature analysis and its uncertainty. J. Clim., 27 (1), 57–75.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S., and J. Schmidhuber, 1997: Long Short-Term Memory. Neural Computation, 9 (8), 1735–1780, 10.1162/neco.1997.9.8.1735, URL https://doi.org/10.1162/neco.1997.9.8.1735.
  * Hoerling et al. (1997) Hoerling, M. P., A. Kumar, and M. Zhong, 1997: El niño, la niña, and the nonlinearity of their teleconnections. 10 (8), 1769–1786, 10.1175/1520-0442(1997)010¡1769:ENOLNA¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/clim/10/8/1520-0442˙1997˙010˙1769˙enolna˙2.0.co˙2.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Izenman (2008) Izenman, A. J., 2008: Linear Dimensionality Reduction. Modern Multivariate Statistical Techniques: Regression, Classification, and Manifold Learning, A. J. Izenman, Ed., Springer Texts in Statistics, Springer, New York, NY, 195–236, 10.1007/978-0-387-78189-1_7, URL https://doi.org/10.1007/978-0-387-78189-1˙7.
  * Izumo (2023) Izumo, T., 2023: Private communication.
  * Izumo et al. (2019) Izumo, T., M. Lengaigne, J. Vialard, I. Suresh, and Y. Planton, 2019: On the physical interpretation of the lead relation between Warm Water Volume and the El Niño Southern Oscillation. Climate Dynamics, 52 (5), 2923–2942, 10.1007/s00382-018-4313-1, URL https://doi.org/10.1007/s00382-018-4313-1.
  * Jebri et al. (2020) Jebri, B., M. Khodri, V. Echevin, G. Gastineau, S. Thiria, J. Vialard, and N. Lebas, 2020: Contributions of Internal Variability and External Forcing to the Recent Trends in the Southeastern Pacific and Peru–Chile Upwelling System. Journal of Climate, 33 (24), 10 555–10 578, 10.1175/JCLI-D-19-0304.1, URL https://journals.ametsoc.org/view/journals/clim/33/24/jcliD190304.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Jiang et al. (2021) Jiang, L., K. Yu, S. Tao, S. Wang, T. Han, and W. Jiang, 2021: ENSO variability during the medieval climate anomaly as recorded by porites corals from the northern south china sea. 36 (4), e2020PA004 173, 10.1029/2020PA004173, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2020PA004173, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2020PA004173.
  * Johnson (2013) Johnson, N. C., 2013: How Many ENSO Flavors Can We Distinguish? Journal of Climate, 26 (13), 4816–4827, 10.1175/JCLI-D-12-00649.1.
  * Kao and Yu (2009) Kao, H.-Y., and J.-Y. Yu, 2009: Contrasting Eastern-Pacific and Central-Pacific Types of ENSO. Journal of Climate, 22 (3), 615–632, 10.1175/2008JCLI2309.1.
  * Karnauskas (2013) Karnauskas, K. B., 2013: Can we distinguish canonical El Niño from Modoki? Geophysical Research Letters, 40 (19), 5246–5251, 10.1002/grl.51007, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/grl.51007, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/grl.51007.
  * Kashinath et al. (2021) Kashinath, K., and Coauthors, 2021: Physics-informed machine learning: case studies for weather and climate modelling. 379 (2194), 20200 093, 10.1098/rsta.2020.0093, URL https://royalsocietypublishing.org/doi/full/10.1098/rsta.2020.0093, publisher: Royal Society.
  * Katz (2002) Katz, R. W., 2002: Sir Gilbert Walker and a Connection between El Niño and Statistics. Statistical Science, 17 (1), 97–112, 10.1214/ss/1023799000, URL https://projecteuclid.org/journals/statistical-science/volume-17/issue-1/Sir-Gilbert-Walker-and-a-Connection-between-El-Ni%c3%b1o-and/10.1214/ss/1023799000.full, publisher: Institute of Mathematical Statistics.
  * Kobayashi et al. (2015) Kobayashi, S., and Coauthors, 2015: The JRA-55 reanalysis: General specifications and basic characteristics. J. Met. Soc. Japan, 93 (1), 5–48.
  * Kohonen (1990) Kohonen, T., 1990: The self-organizing map. Proceedings of the IEEE, 78 (9), 1464–1480, 10.1109/5.58325, URL https://ieeexplore.ieee.org/abstract/document/58325, conference Name: Proceedings of the IEEE.
  * Kug et al. (2009) Kug, J.-S., F.-F. Jin, and S.-I. An, 2009: Two types of el niño events: Cold tongue el niño and warm pool el niño. 22 (6), 1499–1515, 10.1175/2008JCLI2624.1, URL https://journals.ametsoc.org/view/journals/clim/22/6/2008jcli2624.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Lau et al. (1994) Lau, K., P. Sheu, and I. Kang, 1994: Multiscale low-frequency circulation modes in the global atmosphere. J. Atmos. Sci., 51 (9), 1169–1193.
  * Leloup et al. (2007) Leloup, J. A., Z. Lachkar, J.-P. Boulanger, and S. Thiria, 2007: Detecting decadal changes in ENSO using neural networks. Climate Dynamics, 28 (2), 147–162, 10.1007/s00382-006-0173-1, URL https://doi.org/10.1007/s00382-006-0173-1.
  * Lemmon and Karnauskas (2019) Lemmon, D. E., and K. B. Karnauskas, 2019: A metric for quantifying el niño pattern diversity with implications for ENSO–mean state interaction. 52 (12), 7511–7523, 10.1007/s00382-018-4194-3, URL https://doi.org/10.1007/s00382-018-4194-3.
  * Levine et al. (2017) Levine, A. F. Z., M. J. McPhaden, and D. M. W. Frierson, 2017: The impact of the AMO on multidecadal ENSO variability. 44 (8), 3877–3886, 10.1002/2017GL072524, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/2017GL072524, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL072524.
  * Liguori and Di Lorenzo (2019) Liguori, G., and E. Di Lorenzo, 2019: Separating the North and South Pacific Meridional Modes Contributions to ENSO and Tropical Decadal Variability. Geophysical Research Letters, 46 (2), 906–915, 10.1029/2018GL080320.
  * Lloyd (1982) Lloyd, S., 1982: Least squares quantization in PCM. IEEE Transactions on Information Theory, 28 (2), 129–137, 10.1109/TIT.1982.1056489, conference Name: IEEE Transactions on Information Theory.
  * L’Heureux et al. (2013) L’Heureux, M. L., D. C. Collins, and Z.-Z. Hu, 2013: Linear trends in sea surface temperature of the tropical Pacific Ocean and implications for the El Niño-Southern Oscillation. Climate Dynamics, 40 (5), 1223–1236, 10.1007/s00382-012-1331-2, URL https://doi.org/10.1007/s00382-012-1331-2.
  * Madden and Julian (1972) Madden, R. A., and P. R. Julian, 1972: Description of Global-Scale Circulation Cells in the Tropics with a 40–50 Day Period. Journal of the Atmospheric Sciences, 29 (6), 1109–1123, 10.1175/1520-0469(1972)029¡1109:DOGSCC¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/atsc/29/6/1520-0469˙1972˙029˙1109˙dogscc˙2˙0˙co˙2.xml.
  * McGraw and Barnes (2018) McGraw, M. C., and E. A. Barnes, 2018: Memory Matters: A Case for Granger Causality in Climate Variability Studies. Journal of Climate, 31 (8), 3289–3300, 10.1175/JCLI-D-17-0334.1, URL https://journals.ametsoc.org/view/journals/clim/31/8/jcli-d-17-0334.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Mørup and Hansen (2012) Mørup, M., and L. Hansen, 2012: Archetypal analysis for machine learning and data mining. Neurocomputing, 80, 54–63.
  * Nguyen et al. (2022) Nguyen, H., C. Lucas, M. Wheeler, and A. Watkins, 2022: Summary of a workshop on enso-iod alert systems for a warming world held 16-17 august 2022. Tech. rep., Bureau of Meteorology, Australia. URL http://www.bom.gov.au/research/publications/researchreports/BRR-072.pdf, Bureau Research Report - 72.
  * Nguyen and Holmes (2019) Nguyen, L. H., and S. Holmes, 2019: Ten quick tips for effective dimensionality reduction. PLOS Computational Biology, 15 (6), e1006 907, 10.1371/journal.pcbi.1006907, URL https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006907.
  * Okumura and Deser (2010) Okumura, Y. M., and C. Deser, 2010: Asymmetry in the Duration of El Niño and La Niña. Journal of Climate, 23 (21), 5826–5843, 10.1175/2010JCLI3592.1, URL https://journals.ametsoc.org/view/journals/clim/23/21/2010jcli3592.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Oldenborgh et al. (2021) Oldenborgh, G. J. v., H. Hendon, T. Stockdale, M. L’Heureux, E. C. d. Perez, R. Singh, and M. v. Aalst, 2021: Defining El Niño indices in a warming climate. Environmental Research Letters, 16 (4), 044 003, 10.1088/1748-9326/abe9ed, URL https://doi.org/10.1088/1748-9326/abe9ed, publisher: IOP Publishing.
  * Oliver et al. (2021) Oliver, E. C., J. A. Benthuysen, S. Darmaraki, M. G. Donat, A. J. Hobday, N. J. Holbrook, R. W. Schlegel, and A. Sen Gupta, 2021: Marine heatwaves. 13 (1), 313–342, 10.1146/annurev-marine-032720-095144, URL https://www.annualreviews.org/doi/10.1146/annurev-marine-032720-095144.
  * Pegion et al. (2020) Pegion, K., C. M. Selman, S. Larson, J. C. Furtado, and E. J. Becker, 2020: The impact of the extratropics on ENSO diversity and predictability. Climate Dynamics, 54 (9), 4469–4484, 10.1007/s00382-020-05232-3, URL https://doi.org/10.1007/s00382-020-05232-3.
  * Power et al. (2021) Power, S., and Coauthors, 2021: Decadal climate variability in the tropical pacific: Characteristics, causes, predictability, and prospects. 374 (6563), eaay9165, 10.1126/science.aay9165, URL https://www.science.org/doi/10.1126/science.aay9165, publisher: American Association for the Advancement of Science.
  * Rabiner (1989) Rabiner, L., 1989: A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 77 (2), 257–286, 10.1109/5.18626, conference Name: Proceedings of the IEEE.
  * Ramesh and Murtugudde (2013) Ramesh, N., and R. Murtugudde, 2013: All flavours of El Niño have similar early subsurface origins. Nature Climate Change, 3 (1), 42–46, 10.1038/nclimate1600.
  * Rasmusson and Carpenter (1982) Rasmusson, E. M., and T. H. Carpenter, 1982: Variations in tropical sea surface temperature and surface wind fields associated with the southern oscillation/el niño. Mon. Weather Rev., 110, 354–384.
  * Reynolds et al. (2007) Reynolds, R. W., T. M. Smith, C. Liu, D. B. Chelton, K. S. Casey, and M. G. Schlax, 2007: Daily High-Resolution-Blended Analyses for Sea Surface Temperature. Journal of Climate, 20 (22), 5473–5496, 10.1175/2007JCLI1824.1, URL https://journals.ametsoc.org/view/journals/clim/20/22/2007jcli1824.1.xml.
  * Richardson et al. (2021) Richardson, D., A. Black, D. Monselesan, T. Moore, J. Risbey, D. Squire, and C. Tozer, 2021: Identifying periods of forecast model confidence for improved subseasonal prediction of precipitation. J. Hydromet., 22 (2), 371–385.
  * Richter et al. (2022) Richter, I., M. F. Stuecker, N. Takahashi, and N. Schneider, 2022: Disentangling the north pacific meridional mode from tropical pacific variability. 5 (1), 1–9, 10.1038/s41612-022-00317-8, URL https://www.nature.com/articles/s41612-022-00317-8, number: 1 Publisher: Nature Publishing Group.
  * Risbey et al. (2021) Risbey, J., D. Monselesan, A. Black, T. Moore, D. Richardson, D. Squire, and C. Tozer, 2021: The identification of long-lived Southern Hemisphere flow events using archetypes and principal components. Mon. Weather Rev., 149 (6), 1987–2010, https://doi.org/10.1175/MWR-D-20-0314.1.
  * Rogers and Loon (1982) Rogers, J. C., and H. v. Loon, 1982: Spatial variability of sea level pressure and 500 mb height anomalies over the southern hemisphere. 110 (10), 1375–1392, 10.1175/1520-0493(1982)110¡1375:SVOSLP¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/mwre/110/10/1520-0493˙1982˙110˙1375˙svoslp˙2˙0˙co˙2.xml, publisher: American Meteorological Society Section: Monthly Weather Review.
  * Saji et al. (1999) Saji, N., B. Goswami, P. Vinayachandran, and T. Yamagata, 1999: A dipole mode in the tropical Indian Ocean. Nature, 401 (6751), 360–363.
  * Sardeshmukh et al. (2000) Sardeshmukh, P. D., G. P. Compo, and C. Penland, 2000: Changes of probability associated with el niño. 13 (24), 4268–4286, 10.1175/1520-0442(2000)013¡4268:COPAWE¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/clim/13/24/1520-0442˙2000˙013˙4268˙copawe˙2.0.co˙2.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Schlör et al. (2023) Schlör, J., F. Strnad, A. Capotondi, and B. Goswami, 2023: A multi-modal representation of El Niño Southern Oscillation Diversity. arXiv, URL http://arxiv.org/abs/2307.11552, arXiv:2307.11552 [physics].
  * Schulte et al. (2021) Schulte, J., F. Policelli, and B. Zaitchik, 2021: A Continuum Approach to Understanding Changes in the ENSO–Indian Monsoon Relationship. Journal of Climate, 34 (4), 1549–1561, 10.1175/JCLI-D-20-0027.1, URL https://journals.ametsoc.org/view/journals/clim/34/4/JCLI-D-20-0027.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Steinhaus (1957) Steinhaus, H., 1957: Sur la division des corps matériels en parties. Bulletin de l’Académie Polonaise des Sciences, Classe 3, 4, 801–804.
  * Steinschneider and Lall (2015) Steinschneider, S., and U. Lall, 2015: Daily precipitation and tropical moisture exports across the Eastern United States: An application of archetypal analysis to identify spatiotemporal structure. J. Clim., 28 (21), 8585–8602.
  * Sullivan et al. (2016) Sullivan, A., J.-J. Luo, A. C. Hirst, D. Bi, W. Cai, and J. He, 2016: Robust contribution of decadal anomalies to the frequency of central-Pacific El Niño. Scientific Reports, 6 (1), 38 540, 10.1038/srep38540.
  * Takahashi et al. (2011) Takahashi, K., A. Montecinos, K. Goubanova, and B. Dewitte, 2011: ENSO regimes: Reinterpreting the canonical and Modoki El Niño. Geophysical Research Letters, 38 (10), 10.1029/2011GL047364, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2011GL047364, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2011GL047364.
  * Thompson and Wallace (2000) Thompson, D., and J. Wallace, 2000: Annular modes in the extratropical circulation. part i: Month-to-month variability. J. Clim., 13 (5), 1000–1016.
  * Timmermann et al. (2018) Timmermann, A., and Coauthors, 2018: El Niño–Southern Oscillation complexity. Nature, 559 (7715), 535.
  * Trenberth and Stepaniak (2001) Trenberth, K. E., and D. P. Stepaniak, 2001: Indices of El Niño Evolution. Journal of Climate, 14 (8), 1697–1701, 10.1175/1520-0442(2001)014¡1697:LIOENO¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/clim/14/8/1520-0442˙2001˙014˙1697˙lioeno˙2.0.co˙2.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Trendafilov and Gallo (2021) Trendafilov, N., and M. Gallo, 2021: Data analysis on simplexes. Multivariate Data Analysis on Matrix Manifolds: (With Manopt), N. Trendafilov, and M. Gallo, Eds., Springer Series in the Data Sciences, Springer International Publishing, Cham, 373–402, 10.1007/978-3-030-76974-1_11.
  * Tseng et al. (2017) Tseng, Y.-H., R. Ding, and X.-m. Huang, 2017: The warm blob in the northeast pacific—the bridge leading to the 2015/16 el niño. 12 (5), 054 019, 10.1088/1748-9326/aa67c3, URL https://dx.doi.org/10.1088/1748-9326/aa67c3, publisher: IOP Publishing.
  * Turkington et al. (2019) Turkington, T., B. Timbal, and R. Rahmat, 2019: The impact of global warming on sea surface temperature based El Niño–Southern Oscillation monitoring indices. International Journal of Climatology, 39 (2), 1092–1103, 10.1002/joc.5864, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/joc.5864, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/joc.5864.
  * Wang et al. (2012) Wang, S.-Y., M. L’Heureux, and H.-H. Chia, 2012: ENSO prediction one year in advance using western North Pacific sea surface temperatures. Geophysical Research Letters, 39 (5), 10.1029/2012GL050909, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2012GL050909, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2012GL050909.