# A Physics-Informed Auto-Learning Framework for Developing Stochastic
Conceptual Models for ENSO Diversity

Yinling Zhang  Department of Mathematics, University of Wisconsin-Madison,
Madison, USA  Nan Chen  Department of Mathematics, University of Wisconsin-
Madison, Madison, USA  Corresponding author: Nan Chen<EMAIL_ADDRESS>Jérôme Vialard  LOCEAN-IPSL, CNRS/IRD/MNHN/Sorbonne Université, Paris, France
Xianghui Fang  Department of Atmospheric and Oceanic Sciences and Institute of
Atmospheric Sciences, Fudan University, Shanghai, China

###### Summary

Understanding ENSO dynamics has tremendously improved over the past decades.
However, one aspect still poorly understood or represented in conceptual
models is the ENSO diversity in spatial pattern, peak intensity, and temporal
evolution. In this paper, a physics-informed auto-learning framework is
developed to derive ENSO stochastic conceptual models with varying degrees of
freedom. The framework is computationally efficient and easy to apply. Once
the state vector of the target model is set, causal inference is exploited to
build the right-hand side of the equations based on a mathematical function
library. Fundamentally different from standard nonlinear regression, the auto-
learning framework provides a parsimonious model by retaining only terms that
improve the dynamical consistency with observations. It can also identify
crucial latent variables and provide physical explanations. Exploiting a
realistic six-dimensional reference recharge oscillator-based ENSO model, a
hierarchy of three- to six-dimensional models is derived using the auto-
learning framework and is systematically validated by a unified set of
validation criteria assessing the dynamical and statistical features of the
ENSO diversity. It is shown that the minimum model characterizing ENSO
diversity is four-dimensional, with three interannual variables describing the
western Pacific thermocline depth, the eastern and central Pacific sea surface
temperatures (SSTs), and one intraseasonal variable for westerly wind events.
Without the intraseasonal variable, the resulting three-dimensional model
underestimates extreme events and is too regular. The limited number of weak
nonlinearities in the model are essential in reproducing the observed extreme
El Niños and nonlinear relationship between the eastern and western Pacific
SSTs.

##  1 Introduction

El Niño-Southern Oscillation (ENSO) dominates the planet’s interannual
variability with significant impacts on global weather and climate through
atmospheric teleconnections (Philander, 1983, Ropelewski and Halpert, 1987,
Klein et al., 1999, McPhaden et al., 2006, Dai and Wigley, 2000). In the
traditional viewpoint, ENSO consists of alternating periods of anomalously
warm El Niño conditions and cold La Niña conditions that peak in the
equatorial eastern Pacific (EP). However, in recent decades, many El Niño
events have been observed to peak in the central Pacific (CP) (Ashok et al.,
2007, Kao and Yu, 2009, Kim et al., 2012). The discovery of this new type of
ENSO events leads to the El Niño diversity concept (Capotondi et al., 2015),
which underlines that El Niño events vary continuously between the CP and EP
ideal types (Larkin and Harrison, 2005, Yu and Kao, 2007, Ashok et al., 2007,
Kao and Yu, 2009, Kug et al., 2009, Jin, 2022). The shift of the warming
center can cause significant differences in the air-sea coupling over the
equatorial Pacific, changing how ENSO affects the global climate (Chen and
Cane, 2008, Jin et al., 2008, Barnston et al., 2012, Hu et al., 2012, Zheng et
al., 2014, Fang et al., 2015, Sohn et al., 2016, Santoso et al., 2019,
Stuecker et al., 2013). In addition to these two major categories, individual
ENSO events further exhibit diverse characteristics in spatial pattern, peak
intensity, and temporal evolution. This is known as the ENSO complexity
(Timmermann et al., 2018).

Conceptual models of ENSO are useful due to their low dimensionality, which
makes them physically explainable, and mathematically tractable, then further
provide the test beds for examining different physical hypotheses and
interdependence between state variables (Jin, 1997a, b). Several low-order
conceptual models have been developed, including the recharge-discharge
oscillator (Jin, 1997a), the delayed oscillator (Suarez and Schopf, 1988,
Battisti and Hirst, 1989), the western-Pacific oscillator (Weisberg and Wang,
1997), and the advective-reflective oscillator (Picaut et al., 1997). Later, a
unified ENSO oscillator motivated by the dynamics and thermodynamics of Zebiak
and Cane’s coupled ocean-atmosphere model has been built (Wang, 2001). Some
other conceptual models have also been constructed to emphasize the nonlinear
dynamics of ENSO and explore the ENSO features in the decadal time scale
(Timmermann et al., 2003, Roberts et al., 2016, Timmermann and Jin, 2002).
These models were mainly proposed based on physical intuition and highlighted
one or two specific dynamical features of the ENSO as building blocks. They
primarily focused on EP El Niños, and led to many successes in understanding
and predicting their key characteristics.

Comparatively, few models have been built to describe CP events or capture
ENSO diversity and complexity. It was suggested in an early work (Ren and Jin,
2013) that EP and CP dynamics can both be treated as recharge-discharge
oscillators but with different strengths of the thermocline and the zonal
advective feedback. Along this direction, the conceptual model in (Geng et
al., 2020) generalized the original discharge-recharge oscillation for the EP
El Niño (Jin, 1997a) by adding a prognostic equation for the CP SST.
Leveraging the recharge-discharge oscillator concept, (Thual and Dewitte,
2023) proposed a low-dimensional model representing ENSO diversity by
replacing the fixed SST index with a zonally variable warm-pool edge index. On
the other hand, (Fang and Mu, 2018) extended the deterministic two-region
model to a three-region one, including the western Pacific (WP) heat content,
CP SST and zonal current, as well as EP SST. Recently, (Chen et al., 2022)
extended Fang and Mu model (Fang and Mu, 2018) by including seasonality and
two stochastic equations representing the effect of state-dependent synoptic
atmosphere forcing such as (Jin et al., 2007) and decadal fluctuations in the
strength of the Walker circulation. The decadal variability modulates the
amplitude of the zonal advection in the CP region. It leads to the system
alternating between EP-dominant and CP-dominant periods as observed (Capotondi
et al., 2015). This stochastic multiscale model reproduces many important
dynamical and statistical properties of the observed ENSO complexity,
including occurrence frequency (Chen et al., 2022).

While recent modeling efforts have advanced the preliminary understanding of
the diverse features of ENSO, fundamental issues remain in further
conceptualizing ENSO diversity and complexity. First, obtaining the minimal
model reproducing key observed features of ENSO is crucial to understanding
the most fundamental mechanisms contributing to ENSO diversity. Despite its
importance, such a model is still lacking. Second, instead of relying entirely
on human knowledge for model development, building an automatic machine-based
tool that systematically derives appropriate models satisfying specific
requirements is of practical significance. Third, by quantifying model error
and missing information in an existing model, it is vital to systematically
supplement an effective minimal component that facilitates improving the model
performance.

To address these gaps, a physics-informed auto-learning framework is developed
in this paper. This framework can automatically derive stochastic conceptual
models with different complexities that focus on capturing ENSO diversity.
Compared to traditional manual model building, the auto-learning approach
enables efficient exploration of model structures without being limited by
human biases. Specifically, the model structures are determined using rigorous
causal inference that discovers the underlying physics from a comprehensive
inference exploiting a balance between the observational data and prior
knowledge. The physics-informed auto-learning framework is applied to derive a
hierarchy of ENSO stochastic conceptual models with different state variables
and degrees of freedom focusing on ENSO diversity. These models will be
systematically explored and compared using a unified set of validation
criteria metrics that assess how essential ENSO properties, including
diversity, are reproduced. Several issues will be studied to understand
processes crucial for ENSO diversity. First, the framework will be utilized to
test model hierarchies to identify the minimal sufficient model that captures
the essential dynamical and statistical features of the ENSO complexity.
Second, the roles of linear and nonlinear dynamics in describing the ENSO
diversity (Chen and Majda, 2017, An and Jin, 2004, Timmermann, 2003) will be
explored. Such a study aims to discover the fundamental dynamics contributing
to the ENSO diversity in different regions. Third, the necessity of
representing interactions across timescales (Chen et al., 2022, Jin et al.,
2020), including intraseasonal westerly wind (WWB) events, will be analyzed.
Lastly, the framework enables identifying effective latent variables referring
to additional undetermined variables constructed to improve model performance
when the models with predetermined state variables show deficiencies. The time
series and resultant properties of these latent variables guide identifying
suitable missing physical processes to determine the minimal set of variables
needed to capture key ENSO features.

This paper is organized as follows. Section 2 describes the physics-informed
auto-learning framework for generating stochastic conceptual models of ENSO.
Section 3 describes the observational data sets, the reference model whose
synthetic data will be used in the auto-learning procedure, and the
definitions of various types of ENSO. Section 4 focuses on describing the
validation metrics. The hierarchy of models derived from the auto-learning
framework and their intercomparisons are presented in Section 5. Section 6
exploits the set of models in understanding the mechanisms and exploring the
minimum model of ENSO diversity. The paper is concluded in Section 7.

##  2 The Physics-Informed Auto-Learning Framework

###  2.1 Overview

The physics-informed auto-learning framework aims to extract the known
physical knowledge from existing models with the help of additional data to
develop new models automatically. It determines the model structure based on
causal inference and provides a physics-informed parsimonious model. An
overview of the framework is as follows.

  1. Step 1.

An appropriate existing model is needed to provide a long time series for the
auto-learning framework. This model is pre-determined by physical knowledge
and is calibrated by observations. The model simulation is expected to meet as
many desirable validation criteria as possible. In contrast to short
observational data, the long model simulation provides robust learning
outcomes.

  2. Step 2.

The quantities of interests, namely the state variables, are pre-determined to
derive a new model using the framework.

  3. Step 3.

A library of candidate functions is developed as a prerequisite for the auto-
learning framework. These candidate functions are potential terms appearing on
the right-hand side of one or a few equations of the target model. Physical
knowledge can be used to assist the development of such a candidate set.

  4. Step 4.

A causal inference technique is utilized to reveal the underlying relationship
between the time evolution of the state variables and candidate functions. It
determines the right-hand side of the model equations and leads to a
physically explainable parsimonious model.

  5. Step 5.

A simple maximum likelihood method is adopted to estimate model parameters and
characterize the model uncertainty utilizing appropriate stochastic
parameterizations.

It is worthwhile highlighting several significant features of this framework.
First, the framework differs fundamentally from nonlinear regression in
determining model structure. Typical nonlinear regression retains many
nonphysical terms with small coefficients to minimize error. Sparse regression
can eliminate some terms (Santosa and Symes, 1986, Tibshirani, 1996), but
lacks physical constraints. In contrast, the causal inference exploits the
causal relationship, which reflects the underlying physics, to automatically
determine the model structure that provides a physically explainable
parsimonious model (AlMomani et al., 2020, Fish et al., 2021, Kim et al.,
2017, AlMomani and Bollt, 2020, Elinger, 2021, Elinger and Rogers, 2021).

Second, the framework utilizes relationships and insights from the established
models as a starting point rather than fitting models directly to short
observational data. Candidate functional forms for the model structure are
constrained by incorporating physical dependencies learned from the existing
models. Unlike the pure data-driven method, the framework imposes known
dynamical couplings from the known knowledge to inform model derivation.
Physical insights from the established models guide the exploration of model
structures and dependencies.

Third, analytic formulae are available to determine both the model structure
via causal inference and the parameters using maximum likelihood (Chen, 2020).
In other words, the auto-learning framework developed here is wholly explicit
and highly efficient. It differs from many other machine learning methods that
contain black boxes and require expensive training. Therefore, the auto-
learning framework is easy to use and generalizable to many scientific topics.

The details of Steps 2-5 are provided below, where a schematic illustration is
shown in Figure 2.1. The established model used in Step 1, also known as the
reference model, will be described at the end of this section.

Figure 2.1: A schematic illustration of the auto-learning framework.

###  2.2 Determining model structure using causal inference

####  2.2.1 Determining the state variables

The state variables of the model to be derived from the auto-learning
framework are pre-determined. These variables are denoted by an
$n$-dimensional column vector $\mathbf{U}=(u_{1},\ldots,u_{n})^{\mathtt{T}}$,
where each $u_{i}$ is one state variable. To characterize ENSO diversity,
$T_{E}$ and $T_{C}$ are the two variables appearing in all the models to be
derived. Depending on the interests, additional variables are included in
different models, such as the zonal current in the central Pacific to
represent advection processes and the extra heat content variables to account
for more ocean adjustment timescales.

####  2.2.2 Developing a function library

After determining the state variables of the target model, a library
$\mathbf{f}$ consisting of in total $M$ possible candidate functions to
describe the right-hand side of the target model is developed,

$\mathbf{f}=\\{f_{1},\ldots,f_{m-1},f_{m},f_{m+1},\ldots,f_{M}\\}.$ (2.1)

Typically, a large number of candidate functions is included in the library to
allow for coverage over different possible dynamical features. Each $f_{m}$ is
given by a linear or nonlinear function containing a few components of
$\mathbf{U}$. In describing ENSO diversity, the library should at least
include $T_{E}$, $T_{C}$, and some of their linear and nonlinear combinations,
for example, $T_{E}^{3}$ and $T_{C}^{3}$. These functions may not appear in
the resulting model from the auto-learning framework. But the library can
provide as many potential candidate functions as possible that allow the
causal inference to determine the useful ones.

####  2.2.3 Determining model structure using causation inference

Next, causal inference is utilized to evaluate the relevance of the candidate
functions to the underlying dynamics and thus determine the model structure.
To this end, the so-called causation entropy
$C_{f_{m}\rightarrow\dot{u}_{i}\mid\left[\mathbf{f}\backslash{f}_{m}\right]}$
is computed to detect if the candidate function $f_{m}$ contributes to the
right-hand side of the equation for each $u_{i}$, namely ${\,\rm
d}u_{i}/{\,\rm d}t:=\dot{u}_{i}$. The causation entropy is given by (Elinger
and Rogers, 2021, AlMomani et al., 2020, Fish et al., 2021):

$C_{f_{m}\rightarrow\dot{u}_{i}\mid\left[\mathbf{f}\backslash{f}_{m}\right]}=H(\dot{u}_{i}|\left[\mathbf{f}\backslash{f}_{m}\right])-H(\dot{u}_{i}|\mathbf{f}),$
(2.2)

where $\mathbf{f}\backslash{f}_{m}$ represents the set that contains all
functions in $\mathbf{f}$ except $f_{m}$. In other words,
$\mathbf{f}\backslash{f}_{m}$ contains $M-1$ candidate functions and is
defined as

$\mathbf{f}\backslash{f}_{m}=\\{f_{1},\ldots,f_{m-1},f_{m+1},\ldots,f_{M}\\}.$
(2.3)

The term $H(\cdot|\cdot)$ is the conditional entropy, which is related to
Shannon’s entropy $H(\cdot)$ and the joint entropy $H(\cdot,\cdot)$. For two
multi-dimensional random variables $\mathbf{X}$ and $\mathbf{Y}$ (with the
corresponding states being $\mathbf{x}$ and $\mathbf{y}$), they are defined as
(Cover, 1999):

$\begin{split}H(\mathbf{X})&=-\int_{x}p(\mathbf{x})\log(p(\mathbf{x})){\,\rm
d}\mathbf{x},\\\
H(\mathbf{Y}|\mathbf{X})&=-\int_{\mathbf{x}}\int_{\mathbf{y}}p(\mathbf{x},\mathbf{y})\log(p(\mathbf{y}|\mathbf{x})){\,\rm
d}\mathbf{y}{\,\rm d}\mathbf{x},\\\
H(\mathbf{X},\mathbf{Y})&=-\int_{\mathbf{x}}\int_{\mathbf{y}}p(\mathbf{x},\mathbf{y})\log(p(\mathbf{x},\mathbf{y})){\,\rm
d}\mathbf{y}{\,\rm d}\mathbf{x}.\end{split}$ (2.4)

On the right-hand side of (2.2), the difference between the two conditional
entropies indicates the information in $\dot{u}_{i}$ contributed by the
specific function $f_{m}$ given the contributions from all the other
functions. Thus, it tells if $f_{m}$ provides additional information to
$\dot{u}_{i}$ conditioned on the other potential terms in the dynamics. It is
worthwhile to highlight that the causation entropy in (2.2) is fundamentally
different from directly computing the correlation between $\dot{u}_{i}$ and
$f_{m}$, as the causation entropy also considers the influence of the other
library functions. If both $\dot{u}_{i}$ and $f_{m}$ are caused by a common
factor $f_{m^{\prime}}$ in the function library, then $\dot{u}_{i}$ and
$f_{m}$ can be highly correlated. Yet, in such a case, the causation entropy
$C_{f_{m}\rightarrow\dot{u}_{i}\mid\left[\mathbf{f}\backslash{f}_{m}\right]}$
will be zero as $f_{m}$ is not the causation of $\dot{u}_{i}$.

The causation entropy is computed from each of the candidate functions in
$\mathbf{f}$ to each $\dot{u}_{i}$. Thus, there are in total $NM$ causation
entropies, which can be written as a $N\times M$ matrix, called the causation
entropy matrix. Note that the dimension of $\mathbf{X}$ in (2.4) is $M$ when
it is applied to compute the difference between two terms on the right-hand
side of the causation entropy in (2.2). This implies that the direct
calculation of the entropies in (2.4) involves a high-dimensional numerical
integration, which is a well-known computationally challenging issue (Bellman,
1961). To circumvent the direct numerical integration, the entropy calculation
approximates all the joint and marginal distributions as Gaussians. In such a
way, the causation entropy can be approximated by

$\begin{split}C_{\mathbf{Z}\rightarrow\mathbf{X}|\mathbf{Y}}&=H(\mathbf{X}|\mathbf{Y})-H(\mathbf{X}|\mathbf{Y},\mathbf{Z})\\\
&=H(\mathbf{X},\mathbf{Y})-H(\mathbf{Y})-H(\mathbf{X},\mathbf{Y},\mathbf{Z})+H(\mathbf{Y},\mathbf{Z})\\\
&\approx\frac{1}{2}\ln(\operatorname{det}(\mathbf{R}_{\mathbf{X}\mathbf{Y}}))-\frac{1}{2}\ln(\operatorname{det}(\mathbf{R}_{\mathbf{Y}}))-\frac{1}{2}\ln(\operatorname{det}(\mathbf{R}_{\mathbf{X}\mathbf{Y}\mathbf{Z}}))+\frac{1}{2}\ln(\operatorname{det}(\mathbf{R}_{\mathbf{Y}\mathbf{Z}})),\end{split}$
(2.5)

where $\mathbf{R}_{\mathbf{XYZ}}$ denotes the covariance matrix of the state
variables $(\mathbf{X},\mathbf{Y},\mathbf{Z})$ and similar for other
covariances. The notations $\ln(\cdot)$ and $\operatorname{det}(\cdot)$ are
the logarithm of a number and determinant of a matrix, respectively. Note that
the second equality in (2.5) is derived by using the chain rule of conditional
entropy.

The simple and explicit expression in (2.5) based on the Gaussian
approximation can efficiently compute the causation entropy. It allows the
computation of the causation entropy with a moderately large dimension,
sufficient for deriving conceptual models. It is worth noting that the
Gaussian approximation may lead to certain errors in computing the causation
entropy if the true distribution is highly non-Gaussian. Nevertheless, the
primary goal is not to obtain the exact value of the causation entropy.
Instead, it suffices to detect if the causation entropy
$C_{f_{m}\rightarrow\dot{u}_{i}\mid\left[\mathbf{f}\backslash{f}_{m}\right]}$
is nonzero (or practically above a small threshold value). In most
applications, if a significant causal relationship is detected in the higher-
order moments, it is very likely in the Gaussian approximation. This allows us
to efficiently determine the sparse model structure, where the exact values of
the nonzero coefficients on the right-hand side of the model will be
calculated via a simple maximum likelihood estimation to be discussed in the
following, which will be elaborated in the following. Note that the Gaussian
approximation has been widely applied to compute various information
measurements and leads to reasonably accurate results (Majda and Chen, 2018,
Tippett et al., 2004, Kleeman, 2011, Branicki and Majda, 2012).

With the $N\times M$ causation entropy matrix in hand, the next step is
determining the model structure. This can be done by setting up a threshold
value of the causation entropy and retaining only those candidate functions
with the causation entropies exceeding the threshold. Causation entropy allows
the resulting model to contain only functions that significantly contribute to
the dynamics and facilitates a sparse model structure. Sparsity is crucial to
discovering the correct underlying physics and prevents overfitting (Ying,
2019, Brunton et al., 2016). It will also guarantee the robustness of the
model in response to perturbations and allow the model to apply to certain
extrapolation tests.

####  2.2.4 Parameter estimation

The final step is to estimate the parameters in the resulting model. The model
of $\mathbf{U}$ can be written in the vector form:

$\frac{{\,\rm d}\mathbf{U}}{{\,\rm
d}t}=\boldsymbol{\Phi}(\mathbf{U})+\boldsymbol{\sigma}\dot{\mathbf{W}}(t),$
(2.6)

where $\boldsymbol{\Phi}(\mathbf{U})$ is an $n$-dimensional column vector with
each entry representing the summation of the remaining nonlinear candidate
functions from the causal inference for the corresponding component of
$\mathbf{U}$. The matrix $\boldsymbol{\sigma}\in\mathbb{R}^{n\times d}$ is the
noise amplitude, and $\dot{\mathbf{W}}(t)\in\mathbb{R}^{d\times 1}$ is a white
noise. Denote by $\boldsymbol{\Theta}\in\mathbb{R}^{s}$ a column vector
containing all the parameters in the model with $s\gg n$ in a typical
situation. The first term on the right-hand side of (2.6) can be written as

$\boldsymbol{\Phi}(\mathbf{U})=\mathbf{M}\boldsymbol{\Theta}+\mathbf{Q},$
(2.7)

where $\mathbf{M}\in\mathbb{R}^{n\times s}$ is a matrix, each entry of which
is a function of $\mathbf{U}$ while $\mathbf{Q}\in\mathbb{R}^{n}$ is a column
vector that depends on $\mathbf{U}$. Only $\mathbf{M}$ is multiplied by the
parameters $\boldsymbol{\Theta}$. The parameters $\boldsymbol{\Theta}$ can be
easily determined using a maximum likelihood estimator. Meanwhile, the noise
coefficients $\boldsymbol{\sigma}$ are computed. See (Chen, 2020) for the
technical details. Notably, the entire parameter estimation can be solved via
closed analytic formulae, making the procedure efficient and accurate.

###  2.3 Auto-learning with latent variables

Determining the state variables is a prerequisite for model identification.
Yet, it is not always guaranteed that pre-determined state variables are
appropriate for the entire system to capture the target features. Therefore,
it is of practical importance to systematically recognize the sources of the
model deficiency and suggest appropriate additional state variables that
improve the model performance by augmenting the existing model with these
additional quantities.

These can be achieved by incorporating latent variables into the auto-learning
framework. Latent variables are additional processes whose observed time
series are not directly available. They may act as surrogates for the missing
physics due to the lack of a complete understanding of nature. They may also
effectively characterize the state variables at the unresolved scales or play
the role of additional stochastic parameterizations.

In the presence of latent variables, the auto-learning framework needs to be
slightly modified by introducing a joint learning procedure to derive the
equations for both groups of variables. Note that the time series is only
available for the pre-determined state variables (also called observed
variables in the following context). Therefore, an iterative method is used to
alternate between inferring the time series of the latent variables and
deriving the entire model structure. The model structure is determined using
the causal inference discussed in Section 2.2 when the inferred time series of
the latent variables is updated. The recovered time series of the latent
variables, based on the current structure of the latent processes, is obtained
by a Bayesian sampling technique (Chen, 2020). The iterative procedure
terminates when the model structure and parameters converge, which is ensured
by the information revealed from observed variables. See Figure 2.2 for the
illustration of incorporating this iterative process into the auto-learning
framework. The technical details of this iterative method can be found in
(Chen and Zhang, 2022).

One crucial task is to explain the possible physical meanings of the latent
variables once their time series is simulated or sampled from the augmented
system. A stochastic conceptual model with a latent variable representing the
variability in a fast time scale will be presented in Section 6.2.

Figure 2.2: Overview of the auto-learning framework with latent variables.

##  3 Reference Model and Observational data

###  3.1 Observational data sets

The following data sets are utilized in this study to provide some estimates
of the observed ENSO properties, including its diversity. For the assessment
of the model dynamical properties, monthly ocean temperatures, zonal current,
thermocline depth and daily zonal winds at 850 hPa are used. For this, the
oceanic variables are derived from the GODAS dataset (Behringer and Xue, 2004)
(https://www.ncdc.noaa.gov/oisst), with thermocline depth calculated from
potential temperature as the depth of the $20^{o}$C isotherm, while the
atmospheric variable is obtained from the NCEP-NCAR reanalysis (Kalnay et al.,
1996) (http://www.esrl.noaa.gov/psd/). The analysis period is from 1982 to
2019. This period is used because it is during the satellite era and its
observations are more reliable than those before. Anomalies are calculated by
removing the monthly mean climatology of the whole period. For the assessment
of the model performance on simulating different types of ENSO events, a
longer period for SST (i.e., 1950-2020) is also used, which is obtained from
the Extended Reconstructed Sea Surface Temperature version 5 (Huang et al.,
2017) (https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html). The Niño 4
($T_{C}$ in the model) and Niño 3 ($T_{E}$ in the model) indices are the
averages of SST anomalies over the regions 160oE-150oW, 5oS-5oN (CP) and
150o-90oW, 5oS-5oN (EP), respectively. The $h_{W}$ index in the model is the
mean thermocline depth anomaly over the WP region (120oE-180o, 5oS-5oN). The
$u$ index in the model is the mean mixed layer zonal current in the CP region.

###  3.2 The reference model

Instead of using short observational data, an appropriate physics-informed
reference model is adopted to generate a long time series that will be used as
the input data for the auto-learning framework. This reference model captures
many observed features and is calibrated using observational data. Therefore,
it can be regarded as an appropriate surrogate for observations. Various new
models will be derived using the auto-learning framework with simulated data
from this model in Section 5.

The reference model is given by a recently developed three-region multiscale
stochastic model (Chen et al., 2022). This model includes intraseasonal,
interannual, and decadal processes and captures many observed dynamical and
statistical features of ENSO diversity. The model is as follows,

$\displaystyle\frac{{\,\rm d}u}{{\,\rm
d}t}=-ru-\frac{a_{1}b_{0}\mu}{2}\left(T_{C}+T_{E}\right)+\beta_{u}\tau+\sigma_{u}\dot{W}_{u},$
(3.1a) $\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm
d}t}=-rh_{W}-\frac{a_{2}b_{0}\mu}{2}\left(T_{C}+T_{E}\right)+\beta_{h}\tau+\sigma_{h}\dot{W}_{h},$
(3.1b) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=\left(\frac{\gamma
b_{0}\mu}{2}-c_{1}\left(T_{C},t\right)\right)T_{C}+\frac{\gamma
b_{0}\mu}{2}T_{E}+\gamma
h_{W}+c_{I}Iu+C_{u}+\beta_{C}\tau+\sigma_{C}\dot{W}_{C},$ (3.1c)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=\gamma
h_{W}+\left(\frac{3\gamma b_{0}\mu}{2}-c_{2}(t)\right)T_{E}-\frac{\gamma
b_{0}\mu}{2}T_{C}+\beta_{E}\tau+\sigma_{E}\dot{W}_{E},$ (3.1d)
$\displaystyle\frac{{\,\rm d}\tau}{{\,\rm
d}t}=-d_{\tau}\tau+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\tau},$ (3.1e)
$\displaystyle\frac{{\,\rm d}I}{{\,\rm
d}t}=-\lambda(I-m)+\sigma_{I}(I)\dot{W}_{I}.$ (3.1f)

Here, the interannual components, i.e., (3.1a)–(3.1d), depict the main
dynamics for ENSO. The intraseasonal variable $\tau$ in (3.1e) represents the
amplitude of the random wind bursts (Thual et al., 2016) while the decadal
variability $I$ in (3.1f) represents the strength of the background Walker
circulation. In the model, $T_{C}$ and $T_{E}$ are the SSTs in the CP and EP,
$u$ is the zonal ocean current in the CP, and $h_{W}$ is the thermocline depth
in the WP. As was discussed in (Chen et al., 2022), $I$ also stands for the
zonal SST difference between the WP and CP, which affects the interannual
component by modulating the efficiency of the zonal advection through $Iu$ in
(3.1c).

In this model, stochasticity plays a crucial role in coupling variables at
different time scales and parameterizing the unresolved features in the model.
First, intraseasonal wind bursts $\tau$ are modeled by a stochastic
differential equation with state-dependent $\sigma_{\tau}(T_{C})$. The wind
burst $\tau$ is then coupled to the processes of the interannual part serving
as an external forcing. In addition, four Gaussian random noises
$\sigma_{u}\dot{W}_{u}$, $\sigma_{h}\dot{W}_{h}$, $\sigma_{C}\dot{W}_{C}$ and
$\sigma_{E}\dot{W}_{E}$ are further added to the processes, which effectively
parameterize the untracked factors, such as the subtropical atmospheric
forcing at the Pacific Ocean. More broadly, the noise increases variability
for statistics to match observations (Palmer et al., 2009). Second, since the
details of the background Walker circulation consist of uncertainties and
randomness (Chen et al., 2015), a simple but effective stochastic process is
used to describe the temporal evolution of the decadal variability $I$ (Yang
et al., 2021). The multiplicative noise in the process of $I$ is aimed at
guaranteeing its positivity because the long-term average of the background
Walker circulation is non-negative. Besides, the effects of seasonality are
added to both the wind activity and the collective damping to depict the
seasonal phase-locking characteristics realistically, which manifests as the
tendency of ENSO events to peak during boreal winter (Tziperman et al., 1997,
Stein et al., 2014, Fang and Zheng, 2021).

The model parameters are listed in Table 1.

$[h]$ | $150\mathrm{~{}m}$ | $[T]$ | $7.5^{\circ}\mathrm{C}$
---|---|---|---
$[u]$ | $1.5\mathrm{~{}m}\mathrm{~{}s}^{-1}$ | $[t]$ | 2 months
$[\tau]$ | $5\mathrm{~{}m}\mathrm{~{}s}^{-1}$ | $d_{\tau}$ | 2 (1 month-1)
$\gamma$ | 0.45 | $r$ | 0.15
$\alpha_{1}$ | 0.0375 | $\alpha_{2}$ | 0.075
$b_{0}$ | 2.5 | $\mu$ | 0.5
$\sigma$ | $0.12I$ | $\lambda$ | 2/60 (5 years-1)
$p(I)$ | 0.25 in $I\in(0,4)$ | $c_{U}$ | 0.018
$\Phi(x)$ | $\int_{b}^{x}(y-m)p(y)dy$ | $\sigma_{I}(I)$ | $\sqrt{\frac{2}{p(I)}[-\lambda\Phi(I)]}$
$\beta_{E}$ | $0.1239\times(2-0.2I)$ | $\beta_{u}$ | $-0.2\beta_{E}$
$\beta_{h}$ | $-0.4\beta_{E}$ | $\beta_{C}$ | $0.8\beta_{E}$
$\sigma_{u}$ | 0.0310 | $\sigma_{h}$ | 0.0155
$\sigma_{C}$ | 0.0310 | $\sigma_{E}$ | 0.0232
$\sigma_{\tau}\left(T_{C},t\right)$ | $0.9\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(\frac{2\pi}{6}t\right)\right]$ |  |
$c_{1}\left(T_{C},t\right)$ | $\left[15.6\left(T_{C}+\frac{0.75}{7.5}\right)^{2}+0.57\right]\left[1+0.4\sin\left(\frac{2\pi}{6}t\right)\right]$ |  |
$c_{2}(t)$ | $0.9\left[1+0.4\sin\left(\frac{2\pi}{6}t+\frac{2\pi}{6}\right)+0.2\sin\left(\frac{2\pi}{3}t+\frac{2\pi}{6}\right)\right]$ |  |
Table 1: Summary of the non-dimensional units and reference model parameters.
The square brackets around the $h$, $u$, $\tau$, $T$, and $t$ mean the unit of
each variable. Key parameters include the seasonally modulated damping rates
$c_{1}$ and $c_{2}$, where $c_{1}$ has a $T_{C}$ nonlinearity justified by a
heat budget analysis to capture asymmetric CP statistics. The wind burst noise
$\sigma_{\tau}$ also depends on $T_{C}$ to generate irregularity and non-
Gaussian SST distributions which is known as a multiplicative noise forcing
and represents the effect of interannual SST anomalies on wind bursts (Jin et
al., 2007). $\sigma_{I}$ defines stochastic decadal variability for $I$, and
the parameter $m$ is the mean of $I$ which can be computed from its PDF.

###  3.3 Definitions of different types of the ENSO events

The definitions of different El Niño and La Niña events are based on the
average SST anomalies during boreal winter, i.e., from December to February
(DJF). Following the definition in (Kug et al., 2009), EP events are defined
as situations where the EP is warmer than the CP with EP SST anomalies above
0.5 oC. Furthermore, based on the definitions in (Wang et al., 2019), an
extreme El Niño event corresponds to a maximum EP SST anomaly above 2.5oC
between April and the following March. Similarly, CP events are defined as
situations where the CP is warmer than the EP, and CP SST is above 0.5 oC. A
La Niña event is also defined as either the CP or EP SST anomaly is below -0.5
oC. Finally, when El Niño or La Niña occurs for two or more consecutive years,
it is defined as a multi-year El Niño or La Niña, respectively.

##  4 Validation Metrics for Assessing the Model Performance

The models derived from the auto-learning framework will be validated by a
systematic set of metrics with observational data. These metrics measure the
model performance in characterizing the essential features of ENSO diversity
from both the dynamical and the statistical perspectives. The validation
metrics are partitioned into three categories, which are summarized in Figure
4.1. They will be discussed in detail below.

Figure 4.1: A summary of validation metrics for assessing the model
performance.

It is worth mentioning that the auto-learning framework does not explicitly
exploit any of these metrics in deriving the stochastic conceptual models.
Therefore, these metrics are objective for validating the performance of the
resulting models.

###  4.1 A metric for assessing the fundamental dynamical properties

#### (1a) Nonlinear dependence between state variables

One of the main goals of developing a physical model is to reveal the
underlying dynamics. This can, for instance, be achieved by quantifying the
nonlinear inter-dependency between different state variables. Standard model
evaluation metrics such as the lagged correlation are linear and may not
properly reveal such nonlinear relationships.

Instead, a Bayesian inference approach is utilized. This idea is to test the
ability of the model to reconstruct a subset of its state variables
$\mathbf{u}(t)$ given the remaining subset $\mathbf{v}(t)$, using a Bayesian
approach. Further denote the observational time series of these two sets of
variables by $\mathbf{u}^{obs}(t)$ and $\mathbf{v}^{obs}(t)$. Loosely
speaking, the Bayesian inference treats the time series of
$\mathbf{v}^{obs}(t)$ as the constraint for the model to infer the time series
$\mathbf{u}(t)$ and is then compared with $\mathbf{u}^{obs}(t)$. Technically,
this can be achieved by computing the conditional distribution
$p(\mathbf{u}(t)|\mathbf{v}^{obs}(t))$ at each time $t$ using the Bayes
formula,

$p(\mathbf{u}(t)|\mathbf{v}^{obs}(t))\sim p(\mathbf{u}(t))\cdot
p(\mathbf{v}^{obs}(t)|\mathbf{u}(t)),$ (4.1)

where $p(\mathbf{u}(t))$ is obtained from the model simulation, while
$p(\mathbf{v}^{obs}(t)|\mathbf{u}(t))$ is computed based on both the model
output and the available observation $\mathbf{v}^{obs}(t)$. In this work, a
nonlinear Kalman smoother (Emerick and Reynolds, 2013, Evensen, 2003) is
utilized to recover $p(\mathbf{u}(t)|\mathbf{v}^{obs}(t))$. Then the mean
values of $p(\mathbf{u}(t)|\mathbf{v}^{obs}(t))$ at different time instants
provide the most likely time evolution of the model, given the observational
time series. This allows point-wise comparisons between model simulations and
observations. If recovering $\mathbf{u}(t)$ conditioned on
$\mathbf{v}^{obs}(t)$ closely matches $\mathbf{u}^{obs}(t)$, it indicates that
the model accurately captures the nonlinear dependence between these two sets
of state variables. The performance of a model in capturing the fundamental
dynamical properties is then assessed by analyzing whether individual
simulation events match observations. In addition to comparing individual
events, the correlation coefficient between $\mathbf{u}(t)$ and
$\mathbf{u}^{obs}(t)$ can be adopted as a simple metric to quantify such
dynamical consistency.

###  4.2 Metrics for fundamental properties of the ENSO

ENSO exhibits fundamental properties that conceptual models should aim to
reproduce, including its amplitude, asymmetry, seasonal phase locking, and
characteristic timescales. This category of validation metrics focuses on
examining model skill in reproducing these fundamental ENSO properties.

To examine the performance of the models in reproducing fundamental ENSO
properties, a time series of 20000 years from each model is generated.
Omitting the first 1000 years as possibly the burning period, the remaining
time series of 19000 years is divided into 271 non-overlapping segments, each
having a length of 70 years as the observational time series. These different
segments allow us to estimate the 95% confidence interval when computing these
statistics.

#### (2a) Seasonal synchronization

The seasonal variation of the SST variance is a standard statistic reflecting
the seasonal synchronization of a model. This statistic is computed in both
the EP and CP regions.

#### (2b) Probability density functions (PDFs)

ENSO is asymmetrical. Its PDFs are positively skewed in the EP and negatively
skewed in the CP. These asymmetry properties are caused by ENSO extreme
events. These extreme events, known as the super El Niño, have become more
frequent in recent decades and have significantly impacted the global climate
(Thual et al., 2019, Chen et al., 2017, Hameed et al., 2018). The super El
Niño usually occurs in the EP region while most CP events are moderate with
almost no extreme CP El Niño. Therefore, the SST PDFs of these two regions are
crucial for examining the performance of the model in capturing ENSO extreme
events.

#### (2c) Autocorrelation functions (ACFs)

The ACF measures the overall memory of a chaotic system and describes the
averaged convergence rate of the statistics toward the climatology (Gubner,
2006). It is a valuable quantity to characterize the transient behavior and
directly affects the model predictability. If the ACF of a model is strongly
biased, then the model will fail to describe and predict any ENSO events, even
if it has perfectly non-Gaussian PDFs.

Note that despite being widely used for many other studies of time series
analysis, the power spectrum is not included as a validation metric in this
paper since the information delivered by the spectrum overlaps with that of
the temporal correlations (Priestley, 1981, Kay, 1988, Percival and Walden,
1993).

###  4.3 Metrics for ENSO diversity

Besides the fundamental properties of the ENSO, two additional quantitative
metrics are designed here to evaluate model fidelity in reproducing the
observed spatial pattern diversity of ENSO categories and proportions.

#### (3a) Bivariate diagram of location and strength of ENSO events

The first criterion in this set is a bivariate diagram that depicts the
location and the strength of the ENSO events, which was first proposed in
(Capotondi et al., 2015). This is an essential criterion to characterize the
spatiotemporal pattern of the ENSO events.

This metric requires the spatiotemporal reconstruction of the SST data. Denote
by $x$ and $t$ the longitude and the time, $\mbox{SST}(x,t)$ can be obtained
from a multivariate linear regression based on the two time series of the CP
SST $T_{C}$ and the EP SST $T_{E}$ from the model:

$\mbox{SST}(x,t)=r_{C}(x)T_{C}(t)+r_{E}(x)T_{E}(t),$ (4.2)

where $r_{C}$ and $r_{E}$ are the two regression coefficients between the
observational time series at the longitude $x$ and the $T_{C}$ and $T_{E}$
time series, respectively.

#### (3b) Occurrence frequency of different types of events

The second criterion is the occurrence frequency of different ENSO events.
This is a refined metric that counts the number of events in multiple sub-
categories, including all the El Niño events, the EP El Niños, the CP El
Niños, the extreme El Niños, the multi-year El Niños, all the La Niña events,
and the multi-year La Niñas (Chen et al., 2022) according to the definitions
mentioned in Section 3.3.

##  5 A Hierarchy of Stochastic Conceptual Models for ENSO Diversity

###  5.1 Overview of the models

The auto-learning framework is applied to derive five models with different
state variables and degrees of freedom, summarized in Table 2. The detailed
expressions of these models are included in the Appendix. Depending on the
prescribed state variables in the model, the corresponding time series
generated from the reference model is utilized for model identification in the
auto-learning framework.

Model Name | Acronym | Dimension | State Variables
---|---|---|---
Interannual Intraseasonal Decadal Model with Advection | $I_{a}I_{s}DMA$ | 6 | $u,h_{W},T_{C},T_{E},\tau,I$
Interannual Intraseasonal Model with Advection | $I_{a}I_{s}MA$ | 5 | $u,h_{W},T_{C},T_{E},\tau$
Interannual Intraseasonal Decadal Model | $I_{a}I_{s}DM$ | 5 | $h_{W},T_{C},T_{E},\tau,I$
Interannual Intraseasonal Model | $I_{a}I_{s}M$ | 4 | $h_{W},T_{C},T_{E},\tau$
Interannual Model | $I_{a}M$ | 3 | $h_{W},T_{C},T_{E}$
Table 2: Summary of the five models derived from the auto-learning framework

  1. 1.

Interannual Intraseasonal Decadal Model with Advection ($I_{a}I_{s}DMA$): This
6-dimensional model involves three time scales, consistent with the multiscale
system in (Chen et al., 2022).

  2. 2.

Interannual Intraseasonal Model with Advection ($I_{a}I_{s}MA$): This
5-dimensional model excludes the decadal variable $I$ and only contains
interactions between the interannual and intraseasonal variables. This model
will be used to examine if the decadal modulation significantly impacts the
dynamics and statistics of ENSO diversity.

  3. 3.

Interannual Intraseasonal Decadal Model ($I_{a}I_{s}DM$): This 5-dimensional
model involves three time scales but excludes the variable describing the
ocean zonal current. This current variable was introduced to resolve
advection, which plays a strong role in the central Pacific (Vialard et al.,
2001). However, (Jin and An, 1999) pointed out that advection is implicitly
included in the standard Recharge Oscillator framework, suggesting that the
$u$ variable may not be necessary.

  4. 4.

Interannual Intraseasonal Model ($I_{a}I_{s}M$): This 4-dimensional model
excludes the decadal variability and the ocean zonal current. It aims to use a
minimum number of state variables to characterize the interactions between
interannual and intraseasonal variabilities.

  5. 5.

Interannual Model ($I_{a}M$): This 3-dimensional model further excludes the
intraseasonal wind variable. Therefore, it reduces to a single time scale
model describing the interannual variabilities. Note that this model still
includes stochastic forcing. However, the intraseasonal variability, which is
characterized by a state-dependent red noise is omitted in this model. This
model will be compared with the $I_{a}I_{s}M$ to reveal the crucial role of
intraseasonal variability in contributing to ENSO diversity (Chen et al.,
2015, Eisenman et al., 2005, Jin et al., 2007). Despite the lack of skill to
fully characterize ENSO diversity, this $I_{a}M$ will serve as the basis for
the auto-learning framework to include latent variables for model improvement.

The entire library of candidate functions includes all the linear and
quadratic functions of the six state variables, the cubic function of each
variable, and a constant term representing the forcing:

$\begin{gathered}u,\ h_{W},\ T_{C},\ T_{E},\ \tau,\ I,\ u^{2},\ h_{W}^{2},\
T_{C}^{2},\ T_{E}^{2},\ \tau^{2},\ I^{2},\\\ uh_{W},\ uT_{C},\ uT_{E},\
u\tau,\ uI,\ h_{W}T_{C},\ h_{W}T_{E},\ h_{W}\tau,\ h_{W}I,\\\ T_{C}T_{E},\
T_{C}\tau,\ T_{C}I,\ T_{E}\tau,\ T_{E}I,\ \tau I,\ u^{3},\ h_{W}^{3},\
T_{C}^{3},\ T_{E}^{3},\ \tau^{3},\ I^{3},\ 1.\end{gathered}$ (5.1)

For models with only a subset of the six possible variables, the corresponding
candidate functions involving the missing state variables are removed from the
library.

It is worth mentioning that the derivation of the 6-dimensional interannual
intraseasonal decadal model with advection ($I_{a}I_{s}DMA$) is also
considered a proof-of-concept twin experiment to validate the auto-learning
framework. As seen in the following (e.g., Figure 6.1), the identified model
has almost an identical structure to the reference model. The simulations from
these two models also resemble each other. This evidence justifies the auto-
learning framework.

###  5.2 Intercomparison between different models

To provide an overview of the model performance in characterizing different
features of ENSO diversity, let us start with an intercomparison between these
models.

####  5.2.1 Criterion category 1: assessment of the fundamental dynamical
properties

Figure 5.1 shows the results using the Bayesian inference to access the
fundamental dynamical properties. The correlation between the recovered mean
time series and observations is used to quantify the performance of different
models. Panel (a) and (b) illustrate the results of two experiments. In the
first experiment, the SST variables $T_{C}$, $T_{E}$, and decadal variability
$I$ are observed to recover the time series of the remaining variables in each
model. Note that $I$ is the slowest variable and is always treated as
observations. The results show that all models, except the 3-dimensional
interannual model, provide comparably high correlations in recovering $h_{W}$,
which are around 0.7. In contrast, the 3-dimensional interannual model only
has a correlation around 0.5. One interesting finding in Panel (a) is that
recovery of $u$ seems quite challenging because of the weak coupling
relationship between the SST and current $u$ in the models. This also
indicates that $u$ may only have a weak direct contribution to modulating the
SST variables. It is consistent with the fact that $u$ is heavily driven by
the atmospheric wind, and its role has primarily been included by the wind.
Such a finding suggests that the process of $u$ may not be entirely necessary
for building the minimum model of ENSO diversity.

Panel (b) shows the results of the second experiment. Here, the two SST
variables $T_{C}$ and $T_{E}$ are unobserved, while the time series of all the
remaining variables in each model are utilized as observations. All first four
models (from the 6-dimensional interannual intraseasonal decadal model with
advection to the 4-dimensional interannual intraseasonal model) perform
comparably well in recovering the two SST variables, where the correlation is
around 0.75 for $T_{C}$ and 0.82 for $T_{E}$, respectively. These values are
similar to those of the reference model in the inference. In contrast, the
correlations using the 3-dimensional interannual model reduce to 0.53 for
$T_{C}$ and 0.73 for $T_{E}$. These numbers are still above the standard
threshold for skillful results, with a correlation of 0.5. However, compared
with other models, the 3-dimensional interannual is less accurate in
characterizing the interdependence between different state variables.

Two dashed boxes on the left side in Figure 5.1 show the direct results from
Bayesian inference, including the ensemble mean time series (blue), the
confidence interval, and the observations (red) of the reference model as an
example. The reference model has a good performance in characterizing the
nonlinear relationship between the state variables. Specifically, the
recovered $T_{E}$ time series shows significant consistency with the
observation in depicting the events around 1982-1983 and 1997-1998, and the
2015-2016 super El Niños.

Figure 5.1: Comparison between different models using Criterion (1a) –
assessing the fundamental dynamical properties via Bayesian inference. The
correlations between the recovered mean time series and observations are used
to quantify the performance of these models in the right box chart. Panel (a):
correlations of the results for recovery of $u$, $h_{W}$ and $\tau$ by
observing $T_{C}$, $T_{E}$ and $I$. Panel (b): correlations of the results for
recovery of $T_{C}$ and $T_{E}$ by observing other variables. On the left
side, there is an example showing the recovered time series of the reference
model in two experiments. The red lines represent the observation, the solid
blue lines represent the ensemble mean simulation, and the blue shading areas
are the confidence intervals of the ensemble simulations. The correlation $r$
between the truth (red line) and the recovered ensemble mean time series
(solid blue line) for each state variable is also listed.

####  5.2.2 Criterion category 2: quantification of fundamental ENSO
properties

The first row in Figure 5.2 shows the seasonal dependence of the $T_{C}$ and
$T_{E}$ amplitude (Criterion 2a). Only the results of three models are shown
to avoid redundancy – the first four models (from the 6-dimensional
interannual intraseasonal decadal model with advection to the 4-dimensional
interannual intraseasonal model) perform comparably similarly. All the models
have the skill to reproduce seasonal synchronization.

Next, the comparison of the non-Gaussian PDFs (Criterion 2b) is included in
the second row of Figure 5.2. All the models have similar profiles regarding
these PDFs, especially the sign of the skewness for $T_{C}$ and $T_{E}$. The
only noticeable difference among different models is that the PDF of $T_{C}$
associated with the 3-dimensional interannual model is slightly biased
compared with other models. Yet, it is worth remarking that the higher-order
moments that characterize non-Gaussian features may not be easily seen from
the standard PDF plots. Therefore, Figure 5.3 shows the box plots for the
variance, skewness, and kurtosis from the model based on the 271 simulations.
Despite the correct signs of the skewness being correctly recovered by all
these models, the 3-dimensional interannual model has significant biases in
reproducing the skewness and kurtosis values. For $T_{C}$, both the skewness
and kurtosis are stronger than the observations. In particular, the kurtosis
exceeds $3$, which means the probability tail is fatter than a Gaussian
distribution that will cause more occurrence of extreme CP La Niña events in
such a model. In contrast, both the skewness and kurtosis are weaker than
observations for $T_{E}$. This corresponds to the underestimation of the
extreme EP El Niño events. These findings indicate that the state-dependent
red noise is crucial in reproducing these strong non-Gaussian features. The
simple additive white noise is insufficient to capture the observed features
in nature.

Finally, the third row in Figure 5.2 compares different models in recovering
the ACFs (Criterion 2c). Similar to the variances and PDFs, only the
3-dimensional interannual model demonstrates different features compared with
other models. The ACFs of the 3-dimensional interannual model have stronger
oscillation patterns for both $T_{C}$ and $T_{E}$. Such a feature is related
to the strong peak in the power spectrums, as the noise in the temporal
direction is weaker than in the other models.

Figure 5.2: Comparison between different models using Criteria 2a–2c including
the seasonal variation of the variance, PDFs, and ACFs of $T_{C}$ and $T_{E}$.
Figure 5.3: Additional supporting information of Criterion 2b. The second,
third, and fourth moments, i.e., the variance, skewness, and kurtosis, of
$T_{C}$ and $T_{E}$ for each model. These statistics can be computed using the
PDFs in Figure 5.2. In each box plot, the central mark on each box denotes the
median, while the bottom and top edges of the box signify the 5th and 95th
percentiles, respectively. The horizontal dashed line in each panel shows the
true value using observations.

####  5.2.3 Criterion category 3: evaluation of ENSO diversity

Figure 5.4 shows the bivariate diagram of the location and strength of El Niño
events (Criterion 3a). Again, the first four derived models have similar
performance. Thus, only the 4-dimensional interannual intraseasonal model and
3-dimensional interannual model are shown. A continuous spectrum of events
across the equatorial Pacific is obtained for each model. In general, the CP
events always exhibit weak to moderate amplitudes, while amounts of the strong
events can be found in the EP region, which is consistent with observations
(Capotondi et al., 2015). Among different models, the 3-dimensional
interannual model is the only one that underestimates the number of extreme
events in EP. This is related to the inaccuracy in the skewness and kurtosis
shown in Figure 5.3.

Figure 5.4: Comparison between different models using Criterion 3a: bivariate
diagram of location and strength of El Niño events. Each panel shows the
distribution of equatorial Pacific SSTa maxima for El Niño events of the
standard simulations. For each of the El Niño events, the winter-mean SST
anomalies are averaged over the equatorial zone, and then the Pacific zonal
maximum is located. For the presentation purpose, the blue dots show the El
Niño events based on the time series of only the first 4000 years of
simulation in each panel.

Figure 5.5 presents the occurrence frequency of different types of events
(Criterion 3b). All the models, except the 3-dimensional interannual model,
can capture the observed occurrence frequency of all the seven types of ENSO
events. This guarantees the skill of these models in simulating ENSO diversity
and its complex features. Note that the 3-dimensional interannual model
nevertheless reproduces the occurrence frequency of most ENSO events. The only
two exceptions are the underestimation of extreme events, which has already
been discussed in Figures 5.3 and 5.4, and the overestimation of the number of
multi-year El Niño events due to the longer memory shown in the third row of
Figure 5.2.

Figure 5.5: Comparison between different models using Criterion 3b: occurrence
frequency of different types of events. In the box plot, the central mark on
each box denotes the median, while the bottom and top edges of the box signify
the 5th and 95th percentiles, respectively. The horizontal dashed line in each
panel shows the true value using observations.

##  6 Mechanisms and the Minimum Model for ENSO Diversity

Figure 6.1 illustrates the main structure of each model determined by the
auto-learning framework. Note that the same number of columns is utilized for
all the models to keep a unified format in the figure. The columns with dark
blue shadings indicate the functions not built in the library for the
corresponding model. Additional additive noise is added to the equations of
$u$, $h_{W}$, $T_{E}$, and $T_{C}$, while multiplicative noise is imposed on
the $\tau$ and $I$. The details of the models can be found in the Appendix.

Figure 6.1: A summary of the model structures for the reference model (Panel
(a)) and the five derived models from the auto-learning framework (Panels
(b)–(f)). The yellow grids indicate the terms that appear in the model, while
the blue ones mean the functions do not exist on the right-hand side of the
derived model. The columns with dark blue shadings indicate the functions not
built in the library for the corresponding model. The pink boxes mark the
terms that are present in the derived models but not in the reference model.

###  6.1 The minimum model

The 4-dimensional $I_{s}I_{a}M$ is the minimum model that succeeds in
capturing the target features of ENSO diversity characterized by the
validation criteria in Figures 5.1–5.5. A further reduction of the model
dimension to $3$ by eliminating the intraseasonal variability will
significantly decrease the model’s capacity to reproduce many observed
features.

The model, according to Panel (e) of Figure 6.1 or (8.4), also has a
relatively simple structure. It is almost a linear model containing the
interactions between the three interannual variabilities. The only nonlinear
terms are the cubic damping in the process of $T_{C}$, which have been shown
to be crucial in reproducing the non-Gaussian features in $T_{C}$ (Chen et
al., 2022). Similarly, the multiplicative noise in the wind burst process
contributes to the extreme events and non-Gaussian features in $T_{E}$ (Thual
et al., 2016, Levine and Jin, 2017).

Figure 6.2 shows the time series and the spatiotemporal pattern of the
4-dimensional $I_{s}I_{a}M$. The time series in Panel (a) illustrates a clear
relationship between the wind $\tau$ and the EP events in the model, for
example, around $t=230$, $236$, and $254$. The coupled relationship between
the three interannual variabilities in the model is also quite similar to
those in the observations. The spatiotemporal patterns of different types of
ENSO events are presented in Panel (b), including EP El Niño ($t=202$), CP El
Niño ($t=138$), La Niña ($t=203$), extreme El Niño ($t=147$), multi-year El
Niño ($t=206$ to $t=207$) and multi-year La Niña ($t=215$ to $t=216$).

Figure 6.2: Simulation of the 4-dimensional $I_{s}I_{a}M$. Panel (a): time
series of the model (blue), which is compared with observations (red) in a
qualitative way. Panel (b): Hovmoller diagram of the spatiotemporal pattern of
the model simulation. Several different types of ENSO events are marked by
dashed boxes, where the orange ones stand for El Niño events and the blue ones
represent the La Niño events.

###  6.2 Improving the model performance by incorporating latent variables

It has been shown that the 3-dimensional interannual model fails to reproduce
several crucial features of ENSO diversity. This is because the model lacks
certain critical state variables. In general, perfectly pre-determining the
most appropriate set of state variables is not always an easy task. Therefore,
it is beneficial if the auto-learning framework can identify the missing
components in the target model and include additional suitable processes to
compensate for the missing information automatically. As the physical meaning
of the additional processes is unknown at this stage, the associated variables
are treated as latent variables. Once the time series of these latent
variables are recovered, physical justifications can be provided to explain
the resulting model.

To improve the performance of the 3-dimensional interannual model, one latent
variable is added to the system. The auto-learning framework then involves a
joint learning procedure to derive the equations for the three pre-determined
state variables $(h_{W},T_{C},T_{E})$ and the latent variable from the
synthetic data generated from the reference model. Notably, the time series is
only available for the three pre-determined state variables, as the meaning of
the latent variable remains unclear at this stage. Therefore, the learning
procedure here fundamentally differs from deriving the 4-dimensional
$I_{a}I_{s}M$, in which all the four state variables have observations. An
iterative method is applied to alternate between inferring the time series of
the latent variables and deriving the model structure (Chen and Zhang, 2022).
The output of this method is a new model involving existing state variables
and the latent variable after iteration, as shown in Figure 2.2. For
simplicity, only additive noise is assumed in the latent process.

To figure out if the latent variable mimics the behavior of a certain
unobserved state variable, a comparison between the latent variable and this
unobserved state variable in a path-wise way is necessary. Here, synthetic
data generated from the reference model and real observations are taken as
truth in the following two experiments, respectively. Based on the newly
derived 4-dimensional model from the auto-learning framework and true
trajectories of existing state variables, the time series of the latent
variable could be sampled via the conditional sampling method. Figure 6.3
shows the observed time series of the three interannual variables
$(h_{W},T_{C},T_{E})$ and one sampled time series of the latent variable using
the ensemble Kalman smoother data assimilation approach (see also Section
4.1). Panel (a) shows a sampled trajectory based on the time series of
$(h_{W},T_{C},T_{E})$ generated from the reference model, while Panel (b)
utilizes the observational data set as the input. A significant finding is
that the sampled trajectory of the latent variable has a very similar pattern
with the atmospheric wind in both cases. Therefore, the auto-learning
framework succeeds in identifying the dominant component missing in the
3-dimensional interannual model. The latent variable thus provides a natural
surrogate of the atmospheric wind in this new model. For this reason, the
latent variable is denoted by $\tau$. Note that the amplitude of the sampled
trajectories differs from that of the actual observed wind burst $\tau$ (see
the last row of Figure 6.3). This is not surprising as $\tau$ multiplied by
the wind stress coefficient is the actual forcing to the SST processes; for
example, the forcing is given by $L_{34}{\tau}$ in (8.7c). Neither $L_{34}$
nor ${\tau}$ is uniquely determined but the actual wind forcing $L_{34}{\tau}$
can be learned accurately. If the wind stress coefficient is known, then the
latent variable is naturally the wind burst.

Figure 6.3: Comparison of observed time series for interannual variables
$(h_{W},T_{C},T_{E})$ and a sampled latent variable trajectory derived from
the auto-learning framework. Panel (a): the trajectories (pink) of four
variables $(h_{W},T_{C},T_{E},\tau)$ from the 6-dimensional reference model
and a sampled trajectory (blue) of the latent variable $\tau$ based on the new
4-dimensional model from the learning framework when the time series of
$(h_{W},T_{C},T_{E})$ from the reference model is given. Panel (b): Similar to
Panel (a) except the time series from the reference model is replaced by the
observational data (red).

Figure 6.4 displays the structure of the resulting new 4-dimensional model
(Panel (b)) and compares it with the 6-dimensional reference model (Panel
(a)). The latent variable is denoted by $\tau$ in this figure since it has
been identified as a surrogate for the atmospheric wind. The resulting new
4-dimensional model provides a very similar structure to the reference model
in terms of the dynamics $(h_{W},T_{C},T_{E})$, which are dominated by linear
terms but contain a cubic damping for the process of $T_{C}$. It is worth
highlighting that $T_{C}\tau$ appears on the right-hand side of the processes
for the three interannual variables. Since only an additive noise is used in
the latent process, this $T_{C}\tau$ term provides an equivalent
multiplicative noise, indicating that the effective strength of the wind
depends on the SST.

Figure 6.4: Comparison of the model structures of the 6-dimensional reference
model (Panel (a)) and the resulting new 4-dimensional model containing one
latent variable $\tau$ (Panel (b)). Panel (c) shows the causation entropy of
the new model. The dark blue shadings show the terms not included in the pre-
determined function library.

Figure 6.5 includes the dynamical and statistical features of this new
4-dimensional model. These features are mostly missed by the 3-dimensional
interannual model (see Figures 5.1–5.5) but are overall accurately recovered
by the new 4-dimensional model. Panel (a) shows the results of two experiments
utilizing Bayesian Inference. In the first experiment, the observed time
series of $T_{C}$, $T_{E}$ are used to estimate the state of $h_{W}$ and the
latent variable (the surrogate for the wind burst $\tau$ in observational
data), where the amplitude of the recovered wind time series has been adjusted
to match the observations. In the second one, the observed time series of
$h_{W}$, $\tau$ are used to recover the state of $T_{C}$ and $T_{E}$. Among
the two experiments, the recovered time series of $h_{W}$ is more accurate
than that of the 3-dimensional model, which indicates that the new
4-dimensional model captures a better oscillation mechanism between SST and
thermocline with the contribution of the latent variable. Next, Panel (b)
presents the bivariate plot of the location and strength of El Niño events.
The new 4-dimensional model captures the feature that the average strength of
CP events is lower than those of EP events. Besides, the model can generate a
sufficient number of extreme events. This is also supported by the box plot of
the occurrence frequency of different ENSO events in Panel (c). The occurrence
frequencies of different ENSO events from observational data (the red
diamonds) lie in the corresponding box. Similarly, the model performs well in
reproducing long-term statistics, especially the PDFs, which can be seen in
Panel (d). Recall that a noticeable bias was found in the PDF of $T_{C}$ from
the 3-dimensional model in Panel (c) of Figure 5.2. With the addition of the
latent variable, the new 4-dimensional model can accurately capture the non-
Gaussian features. Further validation is shown in Panel (e), which contains
box plots of the variance, skewness, and kurtosis. It illustrates similar
profiles among these aspects with observations (black dashed line). In
conclusion, the latent variable serving as the wind process provides essential
feedback to the interannual variables. It allows the new model to reproduce
the observed features of ENSO diversity.

Figure 6.5: Validation tests for the new 4-dimensional model derived from the
auto-learning framework. Panel (a): assessing the fundamental dynamical
properties. Note that the wind time series has been plotted as a monthly
average for illustration purposes. Panel (b): the bivariate diagram of the
locations and strengths of different ENSO events. Panel (c): the occurrence
frequency of different ENSO events. In the box plot, the central mark on each
box denotes the median, and the bottom and top edges present the 5th and 95th
percentiles, respectively. Panel (d): long-time statistics including PDFs and
ACFs of $T_{C}$ and $T_{E}$. Panel (e): additional supporting information of
PDFs in Panel (d). It presents the statistical features of SST, containing the
variance, skewness, and kurtosis of $T_{C}$ and $T_{E}$.

##  7 Further Discussions and Conclusion

###  7.1 Summary

In this paper, a physics-informed auto-learning framework is built to derive a
hierarchy of stochastic conceptual models with different state variables and
degrees of freedom for ENSO diversity. These models are intercompared based on
validation criteria. Except for the 3-dimensional interannual model, all other
models have comparably skillful results in reproducing many observed features
of ENSO diversity. Therefore, the minimum model for ENSO diversity should
contain at least three interannual variables describing the thermocline depth
in the western Pacific, SST in the EP and the CP, and one intraseasonal
variable for wind bursts. Without the latter, the resulting 3-dimensional
interannual model suffers from underestimating extreme events and biased PDFs.
In addition to the dimensions of the system, the minimum complexity of the
model is also discussed. It is also shown that certain nonlinearities are
crucial in reproducing the observed non-Gaussian features and realistic
amplitudes of different ENSO events. Finally, the auto-learning framework can
also recognize the sources of model error and systematically incorporate
physically explainable latent variables into the model derivation. These
latent variables significantly improve the performance of the resulting
models. Therefore, the auto-learning framework offers a powerful tool to
derive models when only limited prior knowledge is available.

###  7.2 The multiscale nature in modeling ENSO diversity

It has already been shown in the previous section that the 3-dimensional
interannual model fails to characterize several crucial features of ENSO
diversity. In contrast, the new 4-dimensional model that includes the
intraseasonal atmospheric wind variations significantly improves the resulting
dynamics and statistics. This confirms the necessity of having intraseasonal
variability in modeling ENSO diversity.

On the other hand, the role of decadal variability in contributing to the ENSO
diversity can be understood by comparing the 4-dimensional $I_{a}I_{s}M$ with
the 5-dimensional $I_{a}I_{s}DM$ and the 6-dimensional $I_{a}I_{s}DMA$. It
turns out that the results based on the validation metrics do not have a
significant change when the decadal variable is excluded from the models.
However, caution is needed to explain such a result. All the results are based
on 70-year segments of data as same as the length of observations, which could
be too short for the random modulation of the decadal variability to play a
significant role in changing the model statistics. In other words, its role in
modifying the statistics may not be distinguishable from certain random noise.
It is worth highlighting that decadal variability is essential (Dieppois et
al., 2021) in situations with regime switching, such as the climate change
scenario (Yu and Kao, 2007, Yu and Kim, 2013). This was the guideline for
designing the multiscale model in (Chen et al., 2022), where the decadal
variability determines the CP- or EP-favored regimes. If the mean state of the
decadal variability is changed, then it will affect the occurrence frequency
of the CP and EP events (Chen et al., 2022), which cannot be easily reproduced
by simple random noises. The model performance for such a sensitivity
experiment is not included in the model validation criteria adopted here.
These criteria are primarily designed to examine the skill of the models in
reproducing the current climate.

###  7.3 Applicability to CMIP and paleo databases

The auto-learning framework utilizes information theory to identify model
structures requiring long-time datasets. While real-world observations are
limited, model simulations and paleo reconstructions provide ample data for
training conceptual models. Coupled Model Intercomparison Project (CMIP) model
outputs contain multi-century projections enabling robust structure extraction
(Eyring et al., 2016). Paleo proxy records also have long climate time series
that could inform minimal model development for different climate epochs
(Brown et al., 2008). Therefore, applying the auto-learning approach with CMIP
and paleo databases is an intriguing future direction. The framework may
extract underlying dynamics from complex model projections. Paleo-informed
models could also give insights into ENSO evolution under past climate regimes
(Ford et al., 2012). The feasibility of extracting structures from such
alternative datasets merits in-depth exploration.

###  7.4 Relation to recharge oscillator models

The recharge oscillator paradigm established in (Jin, 1997a) represented a
foundation in ENSO conceptual modeling. Recent studies have sought to extend
the recharge oscillator framework by incorporating additional processes and
dimensions. The 6-dimensional conceptual model used here to generate synthetic
data follows this extended recharge oscillator direction. This work
demonstrates that the minimum requirements for a skillful ENSO diversity model
should include three interannual states and intraseasonal wind bursts. The
results from this work provide helpful guidance on potential pathways to
improve recharge oscillator models through appropriate dimensionality
expansion. One guideline for advancing recharge oscillator models includes
incorporating diversity in ENSO spatial patterns (Geng and Jin, 2022, Chen and
Fang, 2023). Thus, the conceptual modeling approach furnishes valuable tools
and insights to guide recharge oscillator model development toward more
comprehensive ENSO simulation.

###  7.5 Potential future directions

In addition to guiding the development of more sophisticated models directly
capturing the spatial patterns of ENSO diversity, the set of stochastic
conceptual models developed in this work has several potential applications.
These models can be used for multi-model data assimilation and forecast for
the large-scale features of ENSO diversity. In particular, the stochastic
nature of these models allows statistical tools to analyze the predictability
and sensitivity analysis (Fang and Chen, 2023, Chen, 2023).

## Acknowledgments

The research of N.C. is partially funded by ARO W911NF-23-1-0118. Y.Z. is
partially supported by this grant. The research of X.F. is supported by
Guangdong Major Project of Basic and Applied Basic Research (Grant No.
2020B0301030004), the Ministry of Science and Technology of the People’s
Republic of China (Grant No. 2020YFA0608802) and the National Natural Science
Foundation of China (Grant Nos. 42192564 and 41805045). The research of J.V.
was supported by the French Agence Nationale pour la Recherche (ANR) ARiSE
grant (ANR-18-CE01-0012). The authors thank Dr. Fei-Fei Jin for useful
discussions.

##  8 Appendix: The Details of the Stochastic Conceptual Models and Parameters
Resulting from the Auto-Learning Framework

###  8.1 The interannual intraseasonal decadal model with advection

The interannual intraseasonal decadal model with advection ($I_{a}I_{s}DMA$)
is given as follows:

$\displaystyle\frac{{\,\rm d}u}{{\,\rm d}t}=$ $\displaystyle
L_{11}u+L_{12}\tau+L_{13}+\sigma_{u}\dot{W}_{u},$ (8.1a)
$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}T_{E}+L_{23}\tau+L_{24}\tau I+L_{25}+\sigma_{h}\dot{W}_{h},$
(8.1b) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}(t)T_{C}+L_{33}T_{E}+L_{34}\tau+L_{35}(t)T_{C}^{2}+L_{36}uI+L_{37}\tau
I$ $\displaystyle+L_{38}(t)T_{C}^{3}+L_{39}+\sigma_{C}\dot{W}_{C},$ (8.1c)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{41}h_{W}+L_{42}T_{C}+L_{43}(t)T_{E}+L_{44}\tau+L_{45}\tau
I+L_{46}+\sigma_{E}\dot{W}_{E},$ (8.1d) $\displaystyle\frac{{\,\rm
d}\tau}{{\,\rm d}t}=$ $\displaystyle
L_{51}\tau+L_{52}+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\mathrm{\tau}},$
(8.1e) $\displaystyle\frac{{\,\rm d}I}{{\,\rm d}t}=$ $\displaystyle
L_{61}I+L_{62}+\sigma_{I}(I)\dot{W}_{I},$ (8.1f)

where the parameters are listed in Table 3.

$L_{11}$ | $-0.1400$ | $L_{31}$ | $0.4540$ | $L_{41}$ | $0.4493$
---|---|---|---|---|---
$L_{12}$ | $-0.0428$ | $L_{32}$ | $-0.4505-0.2989s_{1}$ | $L_{42}$ | $-0.2830$
$L_{13}$ | $0.0000$ | $L_{33}$ | $0.2850$ | $L_{43}$ | $-0.0558-0.3618s_{2}-0.1788s_{3}$
$L_{21}$ | $-0.1663$ | $L_{34}$ | $0.1983$ | $L_{44}$ | $0.2470$
$L_{22}$ | $-0.0694$ | $L_{35}$ | $-3.1142-1.2090s_{1}$ | $L_{45}$ | $-0.0245$
$L_{23}$ | $-0.1007$ | $L_{36}$ | $0.1218$ | $L_{46}$ | $0.0001$
$L_{24}$ | $0.0100$ | $L_{37}$ | $-0.0196$ | $L_{51}$ | $-1.9942$
$L_{25}$ | $0.0001$ | $L_{38}$ | $-15.6559-6.2024s_{1}$ | $L_{52}$ | $0.0045$
$\sigma_{u}$ | $0.0310$ | $L_{39}$ | $0.0177$ | $L_{61}$ | $-0.0323$
$\sigma_{h}$ | $0.0155$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$ | $L_{62}$ | $0.0639$
$\sigma_{C}$ | $0.0310$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$
$\sigma_{E}$ | $0.0232$ | $\sigma_{\tau}$ | $0.8999\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(2\pi/6t\right)\right]$

Table 3: Parameters of the interannual intraseasonal decadal model with
advection ($I_{a}I_{s}DMA$).

Note that the functions of multiplicative noises are not included in the
library since they could not be identified from the framework. The
multiplicative noise of $\tau$ and $I$ are assumed to be known and they are
the same as the reference model.

###  8.2 The interannual intraseasonal model with advection

The interannual intraseasonal model with advection ($I_{a}I_{s}MA$) is given
as follows:

$\displaystyle\frac{{\,\rm d}u}{{\,\rm d}t}=$ $\displaystyle
L_{11}u+L_{12}\tau+L_{13}+\sigma_{u}\dot{W}_{u},$ (8.2a)
$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}T_{C}+L_{23}T_{E}+L_{24}\tau+L_{25}+\sigma_{h}\dot{W}_{h},$
(8.2b) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{31}u+L_{32}h_{W}+L_{33}(t)T_{C}+L_{34}T_{E}+L_{35}\tau+L_{36}(t)T_{C}^{2}+L_{37}(t)T_{C}^{3}$
$\displaystyle+L_{38}+\sigma_{C}\dot{W}_{C},$ (8.2c)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{41}h_{W}+L_{42}T_{C}+L_{43}(t)T_{E}+L_{44}\tau+L_{45}+\sigma_{E}\dot{W}_{E},$
(8.2d) $\displaystyle\frac{{\,\rm d}\tau}{{\,\rm d}t}=$ $\displaystyle
L_{51}\tau+L_{52}+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\tau},$ (8.2e)

where the parameters are listed in Table 4.

$L_{11}$ | $-0.1400$ | $L_{31}$ | $0.2183$ | $L_{41}$ | $0.4443$
---|---|---|---|---|---
$L_{12}$ | $-0.0428$ | $L_{32}$ | $0.4222$ | $L_{42}$ | $-0.2741$
$L_{13}$ | $0.0000$ | $L_{33}$ | $-0.4061-0.2989s_{1}$ | $L_{43}$ | $-0.0553-0.3597s_{2}-0.1791s_{3}$
$L_{21}$ | $-0.1481$ | $L_{34}$ | $0.2455$ | $L_{44}$ | $0.1976$
$L_{22}$ | $-0.0501$ | $L_{35}$ | $0.1568$ | $L_{45}$ | $0.0001$
$L_{23}$ | $-0.0456$ | $L_{36}$ | $-2.8721-1.1140s_{1}$ | $L_{51}$ | $-1.9942$
$L_{24}$ | $-0.0793$ | $L_{37}$ | $-14.5992-5.6033s_{1}$ | $L_{52}$ | $0.0045$
$L_{25}$ | $0.0000$ | $L_{38}$ | $0.0166$ |  |
$\sigma_{u}$ | $0.0310$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$ |  |
$\sigma_{h}$ | $0.0155$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$ |  |
$\sigma_{C}$ | $0.0310$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$ |  |
$\sigma_{E}$ | $0.0232$ | $\sigma_{\tau}$ | $0.8999\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(2\pi t/6\right)\right]$
Table 4: Parameters of the interannual intraseasonal model with advection
($I_{a}I_{s}MA$).

###  8.3 The interannual intraseasonal decadal model

The interannual intraseasonal decadal model ($I_{a}I_{s}DM$) is given as
follows:

$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{11}h_{W}+L_{12}T_{E}+L_{13}\tau+L_{14}\tau I+L_{15}+\sigma_{h}\dot{W}_{h},$
(8.3a) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}(t)T_{C}+L_{23}T_{E}+L_{24}\tau+L_{25}(t)T_{C}^{2}+L_{26}T_{C}I+L_{27}T_{E}I+L_{28}\tau
I$ $\displaystyle+L_{29}(t)T_{C}^{3}+L_{2,10}+\sigma_{C}\dot{W}_{C},$ (8.3b)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}T_{C}+L_{33}(t)T_{E}+L_{34}\tau+L_{35}\tau
I+L_{36}+\sigma_{E}\dot{W}_{E},$ (8.3c) $\displaystyle\frac{{\,\rm
d}\tau}{{\,\rm d}t}=$ $\displaystyle
L_{41}\tau+L_{42}+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\tau},$ (8.3d)
$\displaystyle\frac{{\,\rm d}I}{{\,\rm d}t}=$ $\displaystyle
L_{51}I+L_{52}+\sigma_{I}(I)\dot{W}_{I},$ (8.3e)

where the parameters are listed in Table 5.

$L_{11}$ | $-0.1663$ | $L_{27}$ | $-0.0563$ | $L_{41}$ | $-1.9942$
---|---|---|---|---|---
$L_{12}$ | $-0.0694$ | $L_{28}$ | $-0.0203$ | $L_{42}$ | $0.0045$
$L_{13}$ | $-0.1007$ | $L_{29}$ | $-15.3784-6.2253s_{1}$ | $L_{51}$ | $-0.0323$
$L_{14}$ | $0.0100$ | $L_{2,10}$ | $0.0166$ | $L_{52}$ | $0.0639$
$L_{15}$ | $0.0001$ | $L_{31}$ | $0.4493$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$
$L_{21}$ | $0.5030$ | $L_{32}$ | $-0.2830$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$
$L_{22}$ | $-0.4785-0.2993s_{1}$ | $L_{33}$ | $-0.0558-0.3518s_{2}-0.1788s_{3}$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$
$L_{23}$ | $0.3142$ | $L_{34}$ | $0.2470$ | $\sigma_{h}$ | $0.0155$
$L_{24}$ | $0.1981$ | $L_{35}$ | $-0.0245$ | $\sigma_{C}$ | $0.0310$
$L_{25}$ | $-2.9602-1.1545s_{1}$ | $L_{36}$ | $0.0001$ | $\sigma_{E}$ | $0.0232$
$L_{26}$ | $0.0555$ | $\sigma_{\tau}$ | $0.8999\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(2\pi t/6\right)\right]$
Table 5: Parameters of the interannual intraseasonal decadal model
($I_{a}I_{s}DM$).

###  8.4 The interannual intraseasonal model

The interannual intraseasonal model ($I_{a}I_{s}M$) is given by

$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{11}h_{W}+L_{12}T_{C}+L_{13}T_{E}+L_{14}\tau+L_{15}+\sigma_{h}\dot{W}_{h},$
(8.4a) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}(t)T_{C}+L_{23}T_{E}+L_{24}\tau+L_{25}(t)T_{C}^{2}+L_{26}(t)T_{C}^{3}$
$\displaystyle+L_{27}+\sigma_{C}\dot{W}_{C},$ (8.4b)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}T_{C}+L_{33}(t)T_{E}+L_{34}\tau+L_{35}+\sigma_{E}\dot{W}_{E},$
(8.4c) $\displaystyle\frac{{\,\rm d}\tau}{{\,\rm d}t}=$ $\displaystyle
L_{41}\tau+L_{42}+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\tau},$ (8.4d)

where the parameters are listed in Table 6.

$L_{11}$ | $-0.1481$ | $L_{25}$ | $-2.7519-1.0865s_{1}$ | $L_{41}$ | $-1.9942$
---|---|---|---|---|---
$L_{12}$ | $-0.0501$ | $L_{26}$ | $-14.2956-5.7679s_{1}$ | $L_{42}$ | $0.0045$
$L_{13}$ | $-0.0456$ | $L_{27}$ | $0.0158$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$
$L_{14}$ | $-0.0793$ | $L_{31}$ | $0.4446$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$
$L_{15}$ | $0.0000$ | $L_{32}$ | $-0.2741$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$
$L_{21}$ | $0.4861$ | $L_{33}$ | $-0.0553-0.3597s_{2}-0.1791s_{3}$ | $\sigma_{h}$ | $0.0155$
$L_{22}$ | $-0.3443-0.3014s_{1}$ | $L_{34}$ | $0.1976$ | $\sigma_{C}$ | $0.0310$
$L_{23}$ | $0.1813$ | $L_{35}$ | $0.0000$ | $\sigma_{E}$ | $0.0232$
$L_{24}$ | $0.1558$ | $\sigma_{\tau}$ | $0.8999\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(2\pi t/6\right)\right]$
Table 6: Parameters of the interannual intraseasonal model ($I_{a}I_{s}M$).

###  8.5 The interannual model

The interannual model ($I_{a}M$) is given by

$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{11}h_{W}+L_{12}T_{C}+L_{13}T_{E}+L_{14}T_{C}^{2}+L_{15}+\sigma_{h}\dot{W}_{h},$
(8.5a) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}(t)T_{C}+L_{23}T_{E}+L_{24}\mathrm{h}_{W}T_{C}+L_{25}T_{C}T_{E}+L_{26}+\sigma_{C}\dot{W}_{C},$
(8.5b) $\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}(t)T_{E}+L_{33}T_{C}^{2}+L_{34}(t)T_{E}^{2}+L_{35}(t)T_{E}^{3}+L_{36}+\sigma_{E}\dot{W}_{E},$
(8.5c)

where the parameters are listed in Table 7.

$L_{11}$ | $-0.0678$ | $L_{24}$ | $0.6729$ | $L_{36}$ | $-0.0078$
---|---|---|---|---|---
$L_{12}$ | $-0.1927$ | $L_{25}$ | $-1.1864$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$
$L_{13}$ | $-0.0593$ | $L_{26}$ | $0.0070$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$
$L_{14}$ | $-0.6729$ | $L_{31}$ | $0.2719$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$
$L_{15}$ | $0.0049$ | $L_{32}$ | $0.0523-0.2549s_{2}-0.1792s_{3}$ | $\sigma_{h}$ | $0.0156$
$L_{21}$ | $0.3146$ | $L_{33}$ | $1.1864$ | $\sigma_{C}$ | $0.0311$
$L_{22}$ | $-0.2553-0.1799s_{1}$ | $L_{34}$ | $-0.0134+0.0284s_{2}+0.0381s_{3}$ | $\sigma_{E}$ | $0.0236$
$L_{23}$ | $0.2340$ | $L_{35}$ | $-0.2611-0.3158s_{2}+0.1207s_{3}$ |  |
Table 7: Parameters of the interannual model ($I_{a}M$).

###  8.6 The linear model

The linear 6-dimensional model is given as follows:

$\displaystyle\frac{{\,\rm d}u}{{\,\rm d}t}=$ $\displaystyle
L_{11}u+L_{12}\tau+L_{13}+\sigma_{u}\dot{W}_{u},$ (8.6a)
$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}T_{E}+L_{23}\tau+L_{24}\tau I+L_{25}+\sigma_{h}\dot{W}_{h},$
(8.6b) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}(t)T_{C}+L_{33}T_{E}+L_{34}\tau+L_{35}uI+L_{36}\tau I$
$\displaystyle+L_{37}+\sigma_{C}\dot{W}_{C},$ (8.6c)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{41}h_{W}+L_{42}T_{C}+L_{43}(t)T_{E}+L_{44}\tau+L_{45}\tau
I+L_{46}+\sigma_{E}\dot{W}_{E},$ (8.6d) $\displaystyle\frac{{\,\rm
d}\tau}{{\,\rm d}t}=$ $\displaystyle
L_{51}\tau+L_{52}+\sigma_{\tau}\left(T_{C}\right)\dot{W}_{\mathrm{\tau}},$
(8.6e) $\displaystyle\frac{{\,\rm d}I}{{\,\rm d}t}=$ $\displaystyle
L_{61}I+L_{62}+\sigma_{I}(I)\dot{W}_{I},$ (8.6f)

where the parameters are listed in Table 8.

$L_{11}$ | $-0.1400$ | $L_{31}$ | $0.4461$ | $L_{41}$ | $0.4493$
---|---|---|---|---|---
$L_{12}$ | $-0.0428$ | $L_{32}$ | $-0.5564-0.2799s_{1}$ | $L_{42}$ | $-0.2830$
$L_{13}$ | $0.0000$ | $L_{33}$ | $0.1990$ | $L_{43}$ | $-0.0558-0.3618s_{2}-0.1788s_{3}$
$L_{21}$ | $-0.1663$ | $L_{34}$ | $0.1766$ | $L_{44}$ | $0.2470$
$L_{22}$ | $-0.0694$ | $L_{35}$ | $0.1000$ | $L_{45}$ | $-0.0245$
$L_{23}$ | $-0.1007$ | $L_{36}$ | $-0.0162$ | $L_{46}$ | $0.0001$
$L_{24}$ | $0.0100$ | $L_{37}$ | $-0.0007$ | $L_{51}$ | $-1.9942$
$L_{25}$ | $0.0001$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$ | $L_{52}$ | $0.0045$
$\sigma_{u}$ | $0.0310$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$ | $L_{61}$ | $-0.0323$
$\sigma_{h}$ | $0.0155$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$ | $L_{62}$ | $0.0639$
$\sigma_{C}$ | $0.0310$ | $\sigma_{\tau}$ | 0.8999$\left[\tanh\left(4.5T_{C}\right)+1\right]\left[1+0.25\cos\left(2\pi/6t\right)\right]$
$\sigma_{E}$ | $0.0232$ |  |  |  |
Table 8: Parameters of the linear 6-dimensional model.

###  8.7 The new 4-dimensional model with the latent variable

The new 4-dimensional model with the latent variable is given by

$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}=$ $\displaystyle
L_{11}h_{W}+L_{12}T_{C}+L_{13}\hat{\tau}+L_{14}T_{C}\hat{\tau}+L_{15}+\sigma_{h}\dot{W}_{h},$
(8.7a) $\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}=$ $\displaystyle
L_{21}h_{W}+L_{22}(t)T_{C}+L_{23}\hat{\tau}+L_{24}(t)T_{C}^{2}+L_{25}T_{C}\hat{\tau}+L_{26}(t)T_{C}^{3}$
$\displaystyle+L_{27}+\sigma_{C}\dot{W}_{C},$ (8.7b)
$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}=$ $\displaystyle
L_{31}h_{W}+L_{32}T_{C}+L_{33}(t)T_{E}+L_{34}\hat{\tau}+L_{35}T_{C}\hat{\tau}+L_{36}+\sigma_{E}\dot{W}_{E},$
(8.7c) $\displaystyle\frac{{\,\rm d}\hat{\tau}}{{\,\rm d}t}=$ $\displaystyle
L_{41}\hat{\tau}+L_{42}+\sigma_{\hat{\tau}}\dot{W}_{\hat{\tau}},$ (8.7d)

where the $\hat{\tau}$ denotes the latent variable, the parameters are listed
in Table 9.

$L_{11}$ | $-0.1003$ | $L_{25}$ | $0.1834$ | $L_{41}$ | $-1.5815$
---|---|---|---|---|---
$L_{12}$ | $-0.0819$ | $L_{26}$ | $-10.1500-1.8734s_{1}$ | $L_{42}$ | $-0.1297$
$L_{13}$ | $-0.0289$ | $L_{27}$ | $0.0122$ | $s_{1}$ | $\sin\left(2\pi t/6\right)$
$L_{14}$ | $-0.1069$ | $L_{31}$ | $0.3222$ | $s_{2}$ | $\sin(2\pi t/6+2\pi/6)$
$L_{15}$ | $0.0012$ | $L_{32}$ | $-0.1863$ | $s_{3}$ | $\sin(2\pi t/3+2\pi/6)$
$L_{21}$ | $0.3470$ | $L_{33}$ | $-0.2110-0.1402s_{2}-0.0768s_{3}$ | $\sigma_{h}$ | $0.0155$
$L_{22}$ | $-0.2554-0.2584s_{1}$ | $L_{34}$ | $0.0782$ | $\sigma_{C}$ | $0.0310$
$L_{23}$ | $0.0589$ | $L_{35}$ | $0.2580$ | $\sigma_{E}$ | $0.0232$
$L_{24}$ | $-2.1976-0.5077s_{1}$ | $L_{36}$ | $-0.0017$ | $\sigma_{\hat{\tau}}$ | $2.2034$
Table 9: Parameters of the new 4-dimensional model with one latent variable.

## References

  * Philander (1983) S George H Philander.  El Niño Southern Oscillation phenomena.  _Nature_ , 302(5906):295–301, 1983.
  * Ropelewski and Halpert (1987) Chester F Ropelewski and Michael S Halpert.  Global and regional scale precipitation patterns associated with the El Niño/Southern Oscillation.  _Monthly Weather Review_ , 115(8):1606–1626, 1987.
  * Klein et al. (1999) Stephen A Klein, Brian J Soden, and Ngar-Cheung Lau.  Remote sea surface temperature variations during ENSO: Evidence for a tropical atmospheric bridge.  _Journal of Climate_ , 12(4):917–932, 1999.
  * McPhaden et al. (2006) Michael J McPhaden, Stephen E Zebiak, and Michael H Glantz.  ENSO as an integrating concept in earth science.  _Science_ , 314(5806):1740–1745, 2006.
  * Dai and Wigley (2000) Aiguo Dai and TML Wigley.  Global patterns of ENSO-induced precipitation.  _Geophysical Research Letters_ , 27(9):1283–1286, 2000.
  * Ashok et al. (2007) Karumuri Ashok, Swadhin K Behera, Suryachandra A Rao, Hengyi Weng, and Toshio Yamagata.  El Niño Modoki and its possible teleconnection.  _Journal of Geophysical Research: Oceans_ , 112(C11), 2007.
  * Kao and Yu (2009) Hsun-Ying Kao and Jin-Yi Yu.  Contrasting eastern-Pacific and central-Pacific types of ENSO.  _Journal of Climate_ , 22(3):615–632, 2009.
  * Kim et al. (2012) Jin-Soo Kim, Kwang-Yul Kim, and Sang-Wook Yeh.  Statistical evidence for the natural variation of the central Pacific El Niño.  _Journal of Geophysical Research: Oceans_ , 117(C6), 2012.
  * Capotondi et al. (2015) Antonietta Capotondi, Andrew T Wittenberg, Matthew Newman, Emanuele Di Lorenzo, Jin-Yi Yu, Pascale Braconnot, Julia Cole, Boris Dewitte, Benjamin Giese, Eric Guilyardi, et al.  Understanding ENSO diversity.  _Bulletin of the American Meteorological Society_ , 96(6):921–938, 2015.
  * Larkin and Harrison (2005) Narasimhan K Larkin and DE Harrison.  Global seasonal temperature and precipitation anomalies during El Niño autumn and winter.  _Geophysical Research Letters_ , 32(16), 2005.
  * Yu and Kao (2007) Jin-Yi Yu and Hsun-Ying Kao.  Decadal changes of ENSO persistence barrier in SST and ocean heat content indices: 1958–2001.  _Journal of Geophysical Research: Atmospheres_ , 112(D13), 2007.
  * Kug et al. (2009) Jong-Seong Kug, Fei-Fei Jin, and Soon-Il An.  Two types of El Niño events: cold tongue El Niño and warm pool El Niño.  _Journal of Climate_ , 22(6):1499–1515, 2009.
  * Jin (2022) Fei-Fei Jin.  Toward understanding El Niño southern-oscillation’s spatiotemporal pattern diversity.  _Frontiers in Earth Science_ , 10:899139, 2022.
  * Chen and Cane (2008) Dake Chen and Mark A Cane.  El Niño prediction and predictability.  _Journal of Computational Physics_ , 227(7):3625–3640, 2008.
  * Jin et al. (2008) Emilia K Jin, James L Kinter, Bin Wang, C-K Park, I-S Kang, BP Kirtman, J-S Kug, A Kumar, J-J Luo, J Schemm, et al.  Current status of ENSO prediction skill in coupled ocean–atmosphere models.  _Climate Dynamics_ , 31(6):647–664, 2008.
  * Barnston et al. (2012) Anthony G Barnston, Michael K Tippett, Michelle L L’Heureux, Shuhua Li, and David G DeWitt.  Skill of real-time seasonal ENSO model predictions during 2002–11: Is our capability increasing?  _Bulletin of the American Meteorological Society_ , 93(5):631–651, 2012.
  * Hu et al. (2012) Zeng-Zhen Hu, Arun Kumar, Bhaskar Jha, Wanqiu Wang, Bohua Huang, and Boyin Huang.  An analysis of warm pool and cold tongue El Niños: Air–sea coupling processes, global influences, and recent trends.  _Climate Dynamics_ , 38(9):2017–2035, 2012.
  * Zheng et al. (2014) Fei Zheng, Xiang-Hui Fang, Jin-Yi Yu, and Jiang Zhu.  Asymmetry of the Bjerknes positive feedback between the two types of El Niño.  _Geophysical Research Letters_ , 41(21):7651–7657, 2014.
  * Fang et al. (2015) Xiang-Hui Fang, Fei Zheng, and Jiang Zhu.  The cloud-radiative effect when simulating strength asymmetry in two types of El Niño events using CMIP5 models.  _Journal of Geophysical Research: Oceans_ , 120(6):4357–4369, 2015.
  * Sohn et al. (2016) Soo-Jin Sohn, Chi-Yung Tam, and Hye-In Jeong.  How do the strength and type of ENSO affect SST predictability in coupled models.  _Scientific Reports_ , 6(1):1–8, 2016.
  * Santoso et al. (2019) Agus Santoso, Harry Hendon, Andrew Watkins, Scott Power, Dietmar Dommenget, Matthew H England, Leela Frankcombe, Neil J Holbrook, Ryan Holmes, Pandora Hope, et al.  Dynamics and predictability of El Niño–Southern Oscillation: an Australian perspective on progress and challenges.  _Bulletin of the American Meteorological Society_ , 100(3):403–420, 2019.
  * Stuecker et al. (2013) Malte F Stuecker, Axel Timmermann, Fei-Fei Jin, Shayne McGregor, and Hong-Li Ren.  A combination mode of the annual cycle and the El Niño/southern oscillation.  _Nature Geoscience_ , 6(7):540–544, 2013.
  * Timmermann et al. (2018) Axel Timmermann, Soon-Il An, Jong-Seong Kug, Fei-Fei Jin, Wenju Cai, Antonietta Capotondi, Kim M Cobb, Matthieu Lengaigne, Michael J McPhaden, Malte F Stuecker, et al.  El Niño-Southern Oscillation complexity.  _Nature_ , 559(7715):535–545, 2018.
  * Jin (1997a) Fei-Fei Jin.  An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model.  _Journal of the Atmospheric Sciences_ , 54(7):811–829, 1997a.
  * Jin (1997b) Fei-Fei Jin.  An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-down coupled model.  _Journal of the Atmospheric Sciences_ , 54(7):830–847, 1997b.
  * Suarez and Schopf (1988) Max J Suarez and Paul S Schopf.  A delayed action oscillator for ENSO.  _Journal of Atmospheric Sciences_ , 45(21):3283–3287, 1988.
  * Battisti and Hirst (1989) David S Battisti and Anthony C Hirst.  Interannual variability in a tropical atmosphere–ocean model: Influence of the basic state, ocean geometry and nonlinearity.  _Journal of the Atmospheric Sciences_ , 46(12):1687–1712, 1989.
  * Weisberg and Wang (1997) Robert H Weisberg and Chunzai Wang.  A western Pacific oscillator paradigm for the El Niño-Southern Oscillation.  _Geophysical Research Letters_ , 24(7):779–782, 1997.
  * Picaut et al. (1997) Joël Picaut, François Masia, and Y Du Penhoat.  An advective-reflective conceptual model for the oscillatory nature of the ENSO.  _Science_ , 277(5326):663–666, 1997.
  * Wang (2001) Chunzai Wang.  A unified oscillator model for the el niño–southern oscillation.  _Journal of Climate_ , 14(1):98–115, 2001.
  * Timmermann et al. (2003) Axel Timmermann, Fei-Fei Jin, and Jan Abshagen.  A nonlinear theory for El Niño bursting.  _Journal of the atmospheric sciences_ , 60(1):152–165, 2003.
  * Roberts et al. (2016) Andrew Roberts, John Guckenheimer, Esther Widiasih, Axel Timmermann, and Christopher KRT Jones.  Mixed-mode oscillations of Ell Nino–Southern oscillation.  _Journal of the Atmospheric Sciences_ , 73(4):1755–1766, 2016.
  * Timmermann and Jin (2002) Axel Timmermann and Fei-Fei Jin.  A nonlinear mechanism for decadal El Niño amplitude changes.  _Geophysical Research Letters_ , 29(1):3–1, 2002.
  * Ren and Jin (2013) Hong-Li Ren and Fei-Fei Jin.  Recharge oscillator mechanisms in two types of ENSO.  _Journal of Climate_ , 26(17):6506–6523, 2013.
  * Geng et al. (2020) Tao Geng, Wenju Cai, and Lixin Wu.  Two types of ENSO varying in tandem facilitated by nonlinear atmospheric convection.  _Geophysical Research Letters_ , 47(15):e2020GL088784, 2020.
  * Thual and Dewitte (2023) Sulian Thual and Boris Dewitte.  Enso complexity controlled by zonal shifts in the walker circulation.  _Nature Geoscience_ , 16(4):328–332, 2023.
  * Fang and Mu (2018) Xiang-Hui Fang and Mu Mu.  A three-region conceptual model for central Pacific El Niño including zonal advective feedback.  _Journal of Climate_ , 31(13):4965–4979, 2018.
  * Chen et al. (2022) Nan Chen, Xianghui Fang, and Jin-Yi Yu.  A multiscale model for El Niño complexity.  _npj Climate and Atmospheric Science_ , 5(1):1–13, 2022.
  * Jin et al. (2007) Fei-Fei Jin, L Lin, A Timmermann, and J Zhao.  Ensemble-mean dynamics of the ENSO recharge oscillator under state-dependent stochastic forcing.  _Geophysical Research Letters_ , 34(3), 2007.
  * Chen and Majda (2017) Nan Chen and Andrew J Majda.  Simple stochastic dynamical models capturing the statistical diversity of El Niño Southern Oscillation.  _Proceedings of the National Academy of Sciences_ , 114(7):1468–1473, 2017.
  * An and Jin (2004) Soon-Il An and Fei-Fei Jin.  Nonlinearity and asymmetry of ENSO.  _Journal of Climate_ , 17(12):2399–2412, 2004.
  * Timmermann (2003) Axel Timmermann.  Decadal enso amplitude modulations: A nonlinear paradigm.  _Global and Planetary Change_ , 37(1-2):135–156, 2003.
  * Jin et al. (2020) Fei-Fei Jin, Han-Ching Chen, Sen Zhao, Michiya Hayashi, Christina Karamperidou, Malte F Stuecker, Ruihuang Xie, and Licheng Geng.  Simple enso models.  _El Niño Southern Oscillation in a changing climate_ , pages 119–151, 2020.
  * Santosa and Symes (1986) Fadil Santosa and William W Symes.  Linear inversion of band-limited reflection seismograms.  _SIAM journal on scientific and statistical computing_ , 7(4):1307–1330, 1986.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 58(1):267–288, 1996.
  * AlMomani et al. (2020) Abd AlRahman R AlMomani, Jie Sun, and Erik Bollt.  How entropic regression beats the outliers problem in nonlinear system identification.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(1), 2020.
  * Fish et al. (2021) Jeremie Fish, Alexander DeWitt, Abd AlRahman R AlMomani, Paul J Laurienti, and Erik Bollt.  Entropic regression with neurologically motivated applications.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 31(11), 2021.
  * Kim et al. (2017) Pileun Kim, Jonathan Rogers, Jie Sun, and Erik Bollt.  Causation entropy identifies sparsity structure for parameter estimation of dynamic systems.  _Journal of Computational and Nonlinear Dynamics_ , 12(1):011008, 2017.
  * AlMomani and Bollt (2020) Abd AlRahman AlMomani and Erik Bollt.  Erfit: Entropic regression fit MATLAB package, for data-driven system identification of underlying dynamic equations.  _arXiv preprint arXiv:2010.02411_ , 2020.
  * Elinger (2021) Jared Elinger.  _Information Theoretic Causality Measures For Parameter Estimation and System Identification._ PhD thesis, Georgia Institute of Technology, Atlanta, GA, USA, 2021.
  * Elinger and Rogers (2021) Jared Elinger and Jonathan Rogers.  Causation entropy method for covariate selection in dynamic models.  In _2021 American Control Conference (ACC)_ , pages 2842–2847. IEEE, 2021.
  * Chen (2020) Nan Chen.  Learning nonlinear turbulent dynamics from partial observations via analytically solvable conditional statistics.  _Journal of Computational Physics_ , page 109635, 2020.
  * Cover (1999) Thomas M Cover.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Bellman (1961) Richard Ernest Bellman.  _Dynamic programming treatment of the traveling salesman problem._ RAND Corporation, 1961.
  * Majda and Chen (2018) Andrew J Majda and Nan Chen.  Model error, information barriers, state estimation and prediction in complex multiscale systems.  _Entropy_ , 20(9):644, 2018.
  * Tippett et al. (2004) Michael K Tippett, Richard Kleeman, and Youmin Tang.  Measuring the potential utility of seasonal climate predictions.  _Geophysical research letters_ , 31(22), 2004.
  * Kleeman (2011) Richard Kleeman.  Information theory and dynamical system predictability.  _Entropy_ , 13(3):612–649, 2011.
  * Branicki and Majda (2012) Michal Branicki and Andrew J Majda.  Quantifying uncertainty for predictions with model error in non-Gaussian systems with intermittency.  _Nonlinearity_ , 25(9):2543, 2012.
  * Ying (2019) Xue Ying.  An overview of overfitting and its solutions.  In _Journal of physics: Conference series_ , volume 1168, page 022022. IOP Publishing, 2019.
  * Brunton et al. (2016) Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the national academy of sciences_ , 113(15):3932–3937, 2016.
  * Chen and Zhang (2022) Nan Chen and Yinling Zhang.  A causality-based learning approach for discovering the underlying dynamics of complex systems from partial observations with stochastic parameterization.  _arXiv preprint arXiv:2208.09104_ , 2022.
  * Behringer and Xue (2004) DW Behringer and Yan Xue.  Evaluation of the global ocean data assimilation system at NCEP: The Pacific Ocean.  In _Proc. Eighth Symp. on Integrated Observing and Assimilation Systems for Atmosphere, Oceans, and Land Surface_ , 2004.
  * Kalnay et al. (1996) Eugenia Kalnay, Masao Kanamitsu, Robert Kistler, William Collins, Dennis Deaven, Lev Gandin, Mark Iredell, Suranjana Saha, Glenn White, John Woollen, et al.  The NCEP/NCAR 40-year reanalysis project.  _Bulletin of the American Meteorological Society_ , 77(3):437–472, 1996.
  * Huang et al. (2017) Boyin 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.  Extended reconstructed sea surface temperature, version 5 (ersstv5): upgrades, validations, and intercomparisons.  _Journal of Climate_ , 30(20):8179–8205, 2017.
  * Thual et al. (2016) Sulian Thual, Andrew J Majda, Nan Chen, and Samuel N Stechmann.  Simple stochastic model for El Niño with westerly wind bursts.  _Proceedings of the National Academy of Sciences_ , 113(37):10245–10250, 2016.
  * Palmer et al. (2009) Tim N Palmer, Roberto Buizza, F Doblas-Reyes, Thomas Jung, Martin Leutbecher, Glenn J Shutts, Martin Steinheimer, and Antje Weisheimer.  Stochastic parametrization and model uncertainty.  _ECMWF Technical Memorandum_ , 598, 2009.
  * Chen et al. (2015) Dake Chen, Tao Lian, Congbin Fu, Mark A Cane, Youmin Tang, Raghu Murtugudde, Xunshu Song, Qiaoyan Wu, and Lei Zhou.  Strong influence of westerly wind bursts on El Niño diversity.  _Nature Geoscience_ , 8(5):339–345, 2015.
  * Yang et al. (2021) Qiu Yang, Andrew J Majda, and Nan Chen.  ENSO diversity in a tropical stochastic skeleton model for the MJO, El Niño, and dynamic Walker circulation.  _Journal of Climate_ , pages 1–56, 2021.
  * Tziperman et al. (1997) Eli Tziperman, Stephen E Zebiak, and Mark A Cane.  Mechanisms of seasonal–ENSO interaction.  _Journal of the Atmospheric Sciences_ , 54(1):61–71, 1997.
  * Stein et al. (2014) Karl Stein, Axel Timmermann, Niklas Schneider, Fei-Fei Jin, and Malte F Stuecker.  ENSO seasonal synchronization theory.  _Journal of Climate_ , 27(14):5285–5310, 2014.
  * Fang and Zheng (2021) Xiang-Hui Fang and Fei Zheng.  Effect of the air–sea coupled system change on the ENSO evolution from boreal spring.  _Climate Dynamics_ , pages 1–12, 2021.
  * Wang et al. (2019) Bin Wang, Xiao Luo, Young-Min Yang, Weiyi Sun, Mark A Cane, Wenju Cai, Sang-Wook Yeh, and Jian Liu.  Historical change of El Niño properties sheds light on future changes of extreme El Niño.  _Proceedings of the National Academy of Sciences_ , 116(45):22512–22517, 2019.
  * Emerick and Reynolds (2013) Alexandre A Emerick and Albert C Reynolds.  Ensemble smoother with multiple data assimilation.  _Computers & Geosciences_, 55:3–15, 2013.
  * Evensen (2003) Geir Evensen.  The ensemble kalman filter: Theoretical formulation and practical implementation.  _Ocean dynamics_ , 53:343–367, 2003.
  * Thual et al. (2019) Sulian Thual, Andrew J Majda, and Nan Chen.  Statistical occurrence and mechanisms of the 2014–2016 delayed super El Niño captured by a simple dynamical model.  _Climate Dynamics_ , 52(3-4):2351–2366, 2019.
  * Chen et al. (2017) Lin Chen, Tim Li, Bin Wang, and Lu Wang.  Formation mechanism for 2015/16 super El Niño.  _Scientific reports_ , 7(1):1–10, 2017.
  * Hameed et al. (2018) Saji N Hameed, Dachao Jin, and Vishnu Thilakan.  A model for super El Niños.  _Nature Communications_ , 9(1):2528, 2018.
  * Gubner (2006) John A Gubner.  _Probability and random processes for electrical and computer engineers_.  Cambridge University Press, 2006.
  * Priestley (1981) Maurice Bertram Priestley.  Spectral analysis and time series.  _(No Title)_ , 1981.
  * Kay (1988) Steven M Kay.  _Modern spectral estimation_.  Pearson Education India, 1988.
  * Percival and Walden (1993) Donald B Percival and Andrew T Walden.  _Spectral analysis for physical applications_.  cambridge university press, 1993.
  * Vialard et al. (2001) Jérôme Vialard, Christophe Menkes, Jean-Philippe Boulanger, Pascale Delecluse, Eric Guilyardi, Michael J McPhaden, and Gurvan Madec.  A model study of oceanic mechanisms affecting equatorial pacific sea surface temperature during the 1997–98 el niño.  _Journal of Physical Oceanography_ , 31(7):1649–1675, 2001.
  * Jin and An (1999) Fei-Fei Jin and Soon-Il An.  Thermocline and zonal advective feedbacks within the equatorial ocean recharge oscillator model for enso.  _Geophysical research letters_ , 26(19):2989–2992, 1999.
  * Eisenman et al. (2005) Ian Eisenman, Lisan Yu, and Eli Tziperman.  Westerly wind bursts: ENSO’s tail rather than the dog?  _Journal of Climate_ , 18(24):5224–5238, 2005.
  * Levine and Jin (2017) Aaron FZ Levine and Fei Fei Jin.  A simple approach to quantifying the noise–enso interaction. part i: Deducing the state-dependency of the windstress forcing using monthly mean data.  _Climate Dynamics_ , 48(1-2):1–18, 2017.
  * Dieppois et al. (2021) Bastien Dieppois, Antonietta Capotondi, Benjamin Pohl, Kwok Pan Chun, Paul-Arthur Monerie, and Jonathan Eden.  ENSO diversity shows robust decadal variations that must be captured for accurate future projections.  _Communications Earth & Environment_, 2(1):1–13, 2021.
  * Yu and Kim (2013) Jin-Yi Yu and Seon Tae Kim.  Identifying the types of major El Niño events since 1870.  _International Journal of Climatology_ , 33(8):2105–2112, 2013.
  * Eyring et al. (2016) Veronika Eyring, Sandrine Bony, Gerald A Meehl, Catherine A Senior, Bjorn Stevens, Ronald J Stouffer, and Karl E Taylor.  Overview of the coupled model intercomparison project phase 6 (cmip6) experimental design and organization.  _Geoscientific Model Development_ , 9(5):1937–1958, 2016.
  * Brown et al. (2008) Josephine Brown, Alexander W Tudhope, Matthew Collins, and Helen V McGregor.  Mid-holocene enso: Issues in quantitative model-proxy data comparisons.  _Paleoceanography_ , 23(3), 2008.
  * Ford et al. (2012) Heather L Ford, A Christina Ravelo, and Steven Hovan.  A deep eastern equatorial pacific thermocline during the early pliocene warm period.  _Earth and Planetary Science Letters_ , 355:152–161, 2012.
  * Geng and Jin (2022) Licheng Geng and Fei-Fei Jin.  ENSO diversity simulated in a revised Cane-Zebiak model.  _Frontiers in Earth Science_ , 10:899323, 2022.
  * Chen and Fang (2023) Nan Chen and Xianghui Fang.  A simple multiscale intermediate coupled stochastic model for el niño diversity and complexity.  _Journal of Advances in Modeling Earth Systems_ , 15(4):e2022MS003469, 2023.
  * Fang and Chen (2023) Xianghui Fang and Nan Chen.  Quantifying the predictability of enso complexity using a statistically accurate multiscale stochastic model and information theory.  _Journal of Climate_ , 36(8):2681–2702, 2023.
  * Chen (2023) Nan Chen.  _Stochastic Methods for Modeling and Predicting Complex Dynamical Systems: Uncertainty Quantification, State Estimation, and Reduced-Order Models_.  Springer Nature, 2023.