# A Simple Multiscale Intermediate Coupled Stochastic Model for El Niño
Diversity and Complexity

Nan Chen  Department of Mathematics, University of Wisconsin-Madison, Madison,
WI, USA  Xianghui Fang  Department of Atmospheric and Oceanic Sciences and
Institute of Atmospheric Sciences, Fudan University, 220 Handan Rd., Shanghai
200433, China  Innovation Center of Ocean and Atmosphere System, Zhuhai Fudan
Innovation Research Institute, Zhuhai 518057, China  CMA-FDU Joint Laboratory
of Marine Meteorology, Shanghai 200438, China<EMAIL_ADDRESS>
###### Abstract

El Niño-Southern Oscillation (ENSO) is the most prominent interannual climate
variability in the tropics and exhibits diverse features in spatiotemporal
patterns. In this paper, a simple multiscale intermediate coupled stochastic
model is developed to capture the ENSO diversity and complexity. The model
starts with a deterministic and linear coupled interannual atmosphere, ocean
and sea surface temperature (SST) system. It can generate two distinct
dominant linear solutions that represent the eastern Pacific (EP) and the
central Pacific (CP) El Niños, respectively. In addition to adopting a
stochastic model for characterizing the intraseasonal wind bursts, another
simple stochastic process is developed to describe the decadal variation of
the background Walker circulation. The latter links the two dominant modes in
a simple nonlinear fashion and advances the modulation of the strength and
occurrence frequency of the EP and the CP events. Finally, a cubic nonlinear
damping is adopted to parameterize the relationship between subsurface
temperature and thermocline depth. The model succeeds in reproducing the
spatiotemporal dynamical evolution of different types of the ENSO events. It
also accurately recovers the strongly non-Gaussian probability density
function, the seasonal phase locking, the power spectrum and the temporal
autocorrelation function of the SST anomalies in all the three Niño regions
(3, 3.4 and 4) across the equatorial Pacific. Furthermore, both the composites
of the SST anomalies for various ENSO events and the strength-location
bivariate distribution of equatorial Pacific SST maxima for the El Niño events
from the model simulation highly resemble those from the observations.

##  1 Introduction

El Niño-Southern Oscillation (ENSO) is the most prominent interannual
variability in the tropics. It also affects the global climate, ecosystem and
socioeconomic development through the atmospheric teleconnections [1, 2].
Therefore, understanding and predicting ENSO is a central problem in
contemporary meteorology with large societal impacts. Bjerknes [3] first
suggested that ENSO is the product of tropical air-sea interaction. Since
then, great achievements have been made in its simulation and prediction
abilities [4, 5].

From the traditional point of view, El Niño is defined as the anomalous warm
sea surface temperature (SST) in the equatorial eastern Pacific (EP) region.
Zebiak and Cane [6] developed the first coupled ocean-atmosphere model of
intermediate complexity that successfully characterizes and predicts these EP
warming events. Several deterministic and linear conceptual models were also
proposed to explain the slow physics of ENSO. Among them, the delayed
oscillator [7] describes the delayed effects of oceanic wave reflection at the
ocean western boundary on the EP SST anomalies while the recharge-discharge
oscillator [8] combines SST dynamics and ocean adjustment dynamics into a
coupled basinwide recharge oscillator that relies on the non-equilibrium
between the zonal mean equatorial thermocline depth and wind stress.

With the continuously improved understanding of nature, the spatiotemporal
diversity and complexity of ENSO have been progressively highlighted [9, 10].
In particular, the observational data shows that the center of anomalous SST
is mainly located in the EP from 1980 to 2000, whereas it lies more towards
the central Pacific (CP) after 2000 [11, 12, 13]. See Panel (b) of Figure 1.
The emergence of these distinct ENSO events is called the El Niño diversity
[9]. It suggests the existence of at least two types of El Niños, which are
named as the EP and the CP ones when the peak of the SST anomaly locates in
the cold tongue and near the dateline region, respectively [14, 15, 11, 12,
16]. It is important to notice that the shift of the warming center can make
significant differences in the air-sea coupling over the equatorial Pacific,
which changes the way ENSO affects the global climate and brings serious
challenges to its prediction [17, 18, 19, 20, 21, 22, 23, 24]. In addition to
these two major categories, individual ENSO events further exhibit diverse
characteristics in spatial pattern, peak intensity, and temporal evolution,
which are known as the ENSO complexity [10]. Thus, developing effective
dynamical models that capture the ENSO complexity is of practical importance,
not only for improving the understanding of the formation mechanisms of ENSO
but advancing the prediction of different ENSO events and the associated
varying climatic impacts as well.

The main physical mechanisms of the EP and the CP El Niños are very different.
Due to the strong zonal asymmetry of the tropical Pacific air-sea system, the
thermocline, which is a thin layer that separates the upper warm water from
the cold water in the lower layer, is deep in the western Pacific (WP) while
shallow in the EP region. Such a structure is consistent with the easterly
trade wind. As a result, the SST in the EP is more susceptible to the oceanic
vertical processes, i.e., the thermocline feedback. On the other hand, the
background mean state suggests that the CP is the region with the largest
zonal SST gradient. Consequently, the anomalous zonal current can
significantly affect the local SST variations, which means the development of
the CP type of ENSO is largely influenced by the zonal advective feedback [16,
25, 26, 27, 28, 29, 30]. It is worthwhile to point out that, since the state-
of-art coupled general circulation models (CGCMs) in general have difficulties
in accurately describing the mean state of the tropical Pacific (e.g., the
unrealistic westward extension of the cold tongue), the simulations often
contain biases in reproducing the equatorial SST gradient and the relevant
zonal advective feedback [31, 30, 32, 33]. Such an issue is one of the main
reasons that result in great challenges for the CGCMs to correctly simulate
the El Niño diversity [34, 26, 9, 22, 30].

In addition to the interannual state variables, it is essential to further
take into consideration several other key variabilities belonging to different
temporal scales in the modeling procedure to realistically generate the ENSO
complexity. On one hand, the intraseasonal atmospheric variability, e.g., the
westerly wind burst (WWB) [35, 36, 37] and the Madden-Julian oscillation (MJO)
[38, 39], has been understood as one of the main sources that lead to the ENSO
irregularity and extreme events. Specifically, WWBs can influence the ENSO
development by stimulating eastward-propagating oceanic Kelvin waves,
generating surface zonal currents and weakening the evaporation. In fact, many
modeling works have been attempted to incorporate semistochastic
parameterization for the WWBs [40, 37, 41, 42, 43, 44] and suggest that the
coupled feedbacks between the interannual SST and the intraseasonal WWBs is
sufficient to transfer a damped system to a semi-regular self-sustained
oscillator. Likewise, in light of an intermediate coupled model (ICM), Lian et
al. [45] found that the WWBs are responsible for the existence of the
irregularity and intensity of El Niño, whose specific characteristics depend
on the timing of the WWBs relative to the phase of the recharge–discharge
cycle. On the other hand, ENSO is also modulated by the decadal variation of
the background mean state. Particularly, McPhaden et al. [46] and Xiang et al.
[47] revealed that the changes in the equatorial Pacific around 2000, i.e., a
La Niña-like background state with enhanced trade winds and a more tilted
thermocline, are in favor of the occurrence of more frequent CP El Niño
events. Power et al. [48] also emphasized the role of the decadal variability
in affecting the equatorial Pacific. In addition, by extending the recharge
oscillator into a three-region (i.e., WP, CP and EP) conceptual model that
contains a set of 6 stochastic ordinary differential equations and includes
both the thermocline and zonal advective feedbacks, Chen et al. [49]
demonstrated that the decadal variability plays a crucial role in modulating
the occurrence of the CP and EP El Niños.

The conceptual model in Chen et al. [49] succeeds in capturing many desirable
large-scale features of the ENSO complexity. It thus provides an essential
theoretical basis for the development of a more sophisticated dynamical model,
namely an ICM, that aims at realistically reproducing detailed spatiotemporal
patterns of the ENSO complexity. Different from the conceptual models, the
ICMs have unique advantages of incorporating more elaborate underlying physics
and spatially-extended dynamics into the model development that facilitate the
understanding and prediction of nature. The ICMs also serve as a bridge that
connects the low-order conceptual models and the more complicated CGCMs with a
relatively low computational cost.

To this end, a stochastic ICM for the ENSO complexity is developed in this
paper. The dynamical core of this new stochastic ICM is deterministic and
linear, which involves a coupled interannual atmosphere, ocean and SST system
that drives the basic ENSO dynamics in a simple fashion. Here, the latent
heating that is proportional to the SST is depleted from the ocean and forces
an atmospheric circulation. The resulting zonal wind stress in turn forces
ocean dynamics that provides a feedback to the SST through the thermocline
depth anomalies and the ocean zonal advection. The coupled linear and
deterministic interannual starting model can generate two distinct dominant
linear solutions that represent the EP and the CP El Niños, respectively,
which are essential for simulating the ENSO complexity. The interannual
components are then coupled with the intraseasonal and the decadal
variabilities, which are both described by suitable stochastic processes. The
former is the main contributor to the ENSO irregularity and extreme events
while the latter links the two dominant modes in a simple nonlinear fashion
and advances the modulation of the strength and occurrence frequency of the EP
and the CP events. Seasonal synchronization is further incorporated into the
model, facilitating the ENSO events to have a tendency to peak in boreal
winter. Finally, a cubic nonlinear damping is adopted to parameterize the
relationship between the subsurface temperature and thermocline depth. See
Panel (a) of Figure 1 for a schematic illustration of the model structure and
the key components. Note that, the originally pioneering Zebiak and Cane [6]
ICM was not designed for characterizing the ENSO complexity. A revised version
was recently developed that captures certain diversity features of the ENSO
[50]. Nevertheless, the new simple ICM to be developed in this paper is very
different from the Zebiak and Cane [6] model and its revised version. In fact,
the new model highlights the interactions between variabilities in different
time scales, where only a minimum nonlinearity is adopted to maintain the
model in a simple fashion. The model also exploits suitable stochastic
processes to effectively characterize the dynamical properties and accurately
reproduce the non-Gaussian statistics of the ENSO complexity in different Niño
regions across the equatorial Pacific. The latter is particularly crucial to
realistically simulate different ENSO events. It is a key prerequisite for the
unbiased statistical forecast of the ENSO complexity as well [51, 52].

The rest of the paper is organized as follows. Section 2 presents the details
of the simple stochastic ICM, including the deterministic and linear
interannual components, the stochastic intraseasonal parameterization, the
stochastic decadal process, the seasonal synchronization and the nonlinearly
coupled multiscale system. Section 3 contains the observational datasets and
the definitions of different types of ENSO events. The model simulations are
presented in Section 4 and are compared with the observations. In addition to
showing the spatiotemporal patterns of different ENSO events, the skill of
reproducing several key statistics in different Niño regions and the composite
analysis are also highlighted in this section. Finally, Section 5 contains the
conclusions and discussion.

##  2 The Simple Stochastic ICM

###  2.1 The starting deterministic and linear interannual model

The starting interannual model is a deterministic and linear coupled
atmosphere-ocean-SST system:

Atmosphere:

$\begin{split}&-yv-\partial_{x}\theta=0\\\ &yu-\partial_{y}\theta=0\\\
&-(\partial_{x}u+\partial_{y}v)=E_{q}/(1-\overline{Q})\end{split}$ (1)

Ocean:

$\begin{split}&\partial_{t}U-c_{1}YV+c_{1}\partial_{x}H=c_{1}{\tau_{x}}\\\
&YU+\partial_{Y}H=0\\\
&\partial_{t}H+c_{1}(\partial_{x}U+\partial_{Y}V)=0\end{split}$ (2)

SST:

$\begin{split}&\partial_{t}T=-c_{1}\zeta
E_{q}+c_{1}\eta_{1}H+c_{1}\eta_{2}U.\end{split}$ (3)

The coupled system (1)–(3) consists of a non-dissipative Matsuno–Gill type
atmosphere model [53, 54], a simple shallow-water ocean model [55] and a SST
budget equation [56]. Here, the state variables $u$ and $v$ are the zonal and
meridional wind speeds, $\theta$ is the potential temperature, $U$ and $V$ are
the zonal and meridional ocean currents, $H$ is the thermocline depth and $T$
is the SST. All of them are anomalies. For the coordinate variables, $t$ is
the interannual time coordinate and $x$ is the zonal coordinate, while $y$ and
$Y$ are the meridional coordinates for the atmosphere and ocean components,
respectively. The reason to adopt two distinct meridional axes is that the
deformation radii of atmosphere and ocean are different. In these equations,
$E_{q}=\alpha_{q}T$ is the latent heat with $\overline{Q}$ a constant
representing the background vertical moisture gradient [57], $\tau_{x}=\gamma
u$ is the wind stress, $\zeta$ is the latent heating exchange coefficient,
$\eta_{1}$ and $\eta_{2}$ are the strengths of the thermocline and zonal
advective feedback, respectively. Here, $\eta_{1}$ is stronger in the EP due
to the shallower thermocline there while $\eta_{2}$ is stronger in the CP
because of the larger zonal gradient of the background SST in that region. The
constant $c_{1}$ is related to the ratio between the ocean and atmosphere
phase speeds. The atmosphere extends over the entire equatorial belt $0\leq
x\leq L_{A}$ with periodic boundary conditions, namely
$u(0,y,t)=u(L_{A},y,t)$, and similar for other atmospheric variables. The
Pacific ocean extends over $0\leq x\leq L_{O}$ with reflection boundary
conditions $\int_{-\infty}^{\infty}U(0,Y,t)\mathrm{d}Y=0$ and $U(L_{O},Y,t)=0$
[58, 56]. The detailed parameter values are listed in Appendix.

The above model retains a few essential ingredients that couple the
interannual atmosphere, ocean and SST components and drive the basic ENSO
dynamics in a simple fashion. Specifically, the latent heating $E_{q}$ that is
proportional to the SST $T$ is depleted from the ocean and forces an
atmospheric circulation. The resulting zonal wind stress $\tau_{x}$ in turn
forces ocean dynamics that provides a feedback to the SST through the
thermocline depth anomalies $H$ and the zonal current $U$. See the dashed box
in Panel (a) of Figure 1 that depicts the interannual components.

To facilitate the computation of the model solution, a meridional projection
and truncation is applied to the coupled system, which is known to have
meridional bases in the form of parabolic cylinder functions [59, 43]. For the
purpose of developing a simple ICM, only the leading basis function is kept
for both the atmosphere and the ocean, which are denoted by $\phi_{0}(y)$ and
$\psi_{0}(Y)$, respectively. Both $\phi_{0}(y)$ and $\psi_{0}(Y)$ have
Gaussian profiles and are centered at the equator but the meridional span of
$\phi_{0}(y)$ is larger than that of $\psi_{0}(Y)$. The details of these basis
functions are included in Appendix. The meridional truncations trigger
atmosphere Kelvin, Rossby waves $K_{A},R_{A}$, and ocean Kelvin, Rossby waves
$K_{O},R_{O}$. Once the system (1)–(3) is projected to the leading meridional
bases, the dependence of $y$ and $Y$ are eliminated. The resulting system is
only a function of $t$ and $x$. It reads:

Atmosphere:

$\begin{split}&\partial_{x}K_{A}=-\chi_{A}E_{q}(2-2\bar{Q})^{-1}\\\
&-\partial_{x}R_{A}/3=-\chi_{A}E_{q}(3-3\bar{Q})^{-1}\\\
(B.C.)\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
&K_{A}(0,t)=K_{A}(L_{A},t)\\\ (B.C.)\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ &R_{A}(0,t)=R_{A}(L_{A},t)\end{split}$ (4)

Ocean:

$\begin{split}&\partial_{t}K_{O}+c_{1}\partial_{x}K_{O}=\chi_{O}c_{1}\tau_{x}/2\\\
&\partial_{t}R_{O}-(c_{1}/3)\partial_{x}R_{O}=-\chi_{O}c_{1}\tau_{x}/3\\\
(B.C.)\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
&K_{O}(0,t)=r_{W}R_{O}(0,t)\\\ (B.C.)\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ &R_{O}(L_{O},t)=r_{E}K_{O}(L_{O},t)\end{split}$ (5)

SST:

$\begin{split}&\partial_{t}T=-c_{1}\zeta
E_{q}+c_{1}\eta_{1}(K_{O}+R_{O})+c_{1}\eta_{2}(K_{O}-R_{O}),\end{split}$ (6)

where $r_{W}$ and $r_{E}$ are the reflection coefficients associated with the
ocean reflection boundary conditions (B.C.) while $\chi_{A}$ and $\chi_{O}$
are the meridional projection coefficients. Once these waves are solved, the
physical variables can be reconstructed,

$\begin{split}&u=(K_{A}-R_{A})\phi_{0}+(R_{A}/\sqrt{2})\phi_{2}\\\
&\theta=-(K_{A}+R_{A})\phi_{0}-(R_{A}/\sqrt{2})\phi_{2}\\\
&U=(K_{O}-R_{O})\psi_{0}+(R_{O}/\sqrt{2})\psi_{2}\\\
&H=(K_{O}+R_{O})\psi_{0}+(R_{O}/\sqrt{2})\psi_{2}\end{split}$ (7)

where $\phi_{2}$ and $\psi_{2}$ are the third meridional bases of atmosphere
and ocean, respectively.

After applying a spatial discretization in the $x$ direction, the coupled
system (4)–(6) is solved numerically via an upwind finite difference scheme.
Since the coupled system is linear and deterministic, its final solution,
after the numerical discretization, can be written as a superposition of a set
of non-interacting linear modes (the so-called linear solutions). Each linear
solution is associated with one eigenmode of the system. In the numerical
discretization here, the entire equatorial band is divided into $N_{A}$
equidistance grids and correspondingly there are $N_{O}$ grid points in the
Pacific ocean. In the simulations of this paper, $N_{A}=128$ and $N_{O}=56$
are utilized. In other words, the distance between every two grid points is
$312.5$km, as the entire equatorial band and the span of the Pacific ocean are
$40,000$km and $17,500$km, respectively. With an appropriate choice of the
physical parameters (see Appendix), all the eigenvalues have negative real
parts, indicating the decaying nature of these linear solutions. It is
important to highlight that although the eigenvalues of the strongly decaying
and fast oscillating small-scale modes may vary by changing the resolution in
the spatial discretization, the leading two eigenmodes with the slowest
decaying rate remain almost unchanged as long as the spatial discretization is
not too coarse. The leading two eigenmodes appear as a pair and the associated
eigenvalues are complex conjugate to each other, where the associated
oscillation frequency lies on the interannual time scale. Due to the slowest
decaying rate, the full solution of the coupled system (4)–(6) is dominated by
these eigenmodes.

Figure 2 shows the spatiotemporal evolution of the leading two eigenmodes. For
the purpose of illustration, the decaying rate is manually set to be zero in
demonstrating the spatiotemporal pattern of such linear solutions in this
figure. By varying the strength of the zonal advective feedback coefficient
$\eta_{2}$, the dominant eigenmodes can have distinct behavior (see Appendix
for the detailed parameter values). Specifically, if the role of the zonal
advection is weakened, then the leading linear solution exhibits
spatiotemporal patterns with the EP El Niño being dominated. See Panel (a). In
such a situation, the thermocline feedback is the main mechanism for the SST
development. In addition, the convergence center of the atmospheric wind lies
in the eastern Pacific. On the other hand, if the zonal advective feedback
becomes stronger, then the CP El Niño pattern becomes dominant.
Correspondingly, the ocean zonal current leads to the warming in the CP region
and the associated convergence center of the atmospheric wind shifts westward
as well. See Panel (b). It is worthwhile to highlight that the occurrence
frequency of the CP events (every 2.5 years) is higher than that of the EP
events (every 4.5 years) in these dominant linear solutions, which is
consistent with the observations. Note that these two linear solutions are the
necessary conditions and mechanisms for the model to capture the ENSO
complexity.

###  2.2 Simple stochastic models for the intraseasonal and decadal
variabilities

As was seen in Section 2.1, the coupled system (4)–(6) can generate basic
linear solutions that exhibit regular patterns of the EP and CP El Niño events
in different situations. However, the irregularity and complexity of ENSO
require extra mechanisms beyond the deterministic and linear dynamics. In
particular, the ENSO variability is often triggered or inhibited by a broad
range of random atmospheric disturbances in the tropics, such as the WWBs [35,
36, 37], the easterly wind bursts (EWBs) [60], as well as the convective
envelope of the MJO [38]. On the other hand, it has been shown that the EP and
CP events were alternatively prevalent every 10 to 20 years over the past
century [61, 62]. For example, the EP events were dominant in the 1980s while
the CP El Niños were more frequently identified after 2000 [63, 64]. These
findings imply that the decadal variability plays a crucial role in driving
the transitions between the CP- and EP-dominant regimes. Thus, it is also
essential to incorporate the decadal effect into the coupled ENSO model to
link the separate linear solutions.

To this end, two stochastic processes are developed and coupled to the
starting interannual model (4)–(6). These two stochastic processes
characterize the intraseasonal random wind bursts and the decadal variability,
respectively. The former is a natural component that depicts the random
atmospheric disturbances. The latter describes the decadal variation of the
background Walker circulation. It is related to the climate change scenario
and plays an important role in modulating the strength and the occurrence
frequency of the EP and the CP events [63, 49].

First, with the stochastic wind bursts, the wind stress $\tau_{x}$ now
contains two components $\tau_{x}=\gamma(u+u_{p})$, where $u$ remains the same
as the atmospheric circulation in (2) while $u_{p}$ is the contribution from
the stochastic wind bursts, which is assumed to have the following structure,

$u_{p}(x,y,t)=a_{p}(t)s_{p}(x)\phi_{0}(y),$

where $\phi_{0}(y)$ is again the leading meridional basis while $s_{p}(x)$ is
a fixed spatial structure that is localized in the western Pacific, due to the
fact that most of the observed wind bursts are active in the western Pacific.
The time series $a_{p}(t)$ describes the wind burst amplitude and is governed
by a simple one-dimensional real-valued stochastic process [65]

$\frac{\mathrm{d}a_{p}}{\mathrm{d}t}=-d_{p}a_{p}+{\sigma_{p}(T_{C})}\dot{W}_{p},$
(8)

where $d_{p}$ is the damping term that is chosen such that the decorrelation
time of the wind is about 1 month. In (8), $\dot{W}_{p}$ is a white noise
source while $\sigma_{p}(T_{C})$ is its strength. When $a_{p}$ is positive and
negative, it represents the WWB and EWB, respectively. It is important to
highlight that the noise strength ${\sigma_{p}(T_{C})}$ is state-dependent
(the so-called multiplicative noise), as a function of the interannual SST
from (6) averaged over the western-central Pacific, namely Niño 4 region. The
reason of choosing such a state-dependent noise coefficient is that wind burst
activity is usually more active with warmer SST in the western-central Pacific
due to the strengthening or eastward extension of the warm pool [36, 38],
which is modeled here in a simple parameterized fashion. Note that the
enhanced SST only increases the amplitude of the wind bursts while the
individual wind burst event generated from the stochastic process in (8) does
not have a preference towards westerly or easterly. This allows an equal
chance to create both the WWB and the EWB as individual events, which is
consistent with the observations. Due to the state-dependent noise
coefficient, the modeling procedure here indicates that the intraseasonal wind
bursts not only affect the interannual variability but also are modulated by
the latter as well.

Next, the decadal variability is driven by another simple stochastic process,

$\frac{\mathrm{d}I}{\mathrm{d}t}=-\lambda(I-m)+\sigma_{I}(I)\dot{W}_{I},$ (9)

where the damping $\lambda$ is set to be $5$ years-1 representing the decadal
time scale. Similar to (8), $\sigma_{I}(I)$ and $\dot{W}_{I}$ here are the
state-dependent noise strength and the white noise source. The reason to adopt
a state-dependent noise coefficient, which is a function of $I$ itself, is to
allow the distribution of $I$ to be non-Gaussian. In particular, the trade
wind in lower level of the Walker circulation in the decadal time scale is
easterly, which means the sign of $I$ should stay the same throughout time and
thus the distribution of $I$ is obviously non-Gaussian. This feature can be
easily incorporated into the process of $I$ with the state-dependent noise
coefficient [66, 67]. Based on the limited observational data and the theory
of inferring the least unbiased maximum entropy solution for a distribution, a
uniform distribution between $[0,1]$ is adopted for $I$ in this work. Here, a
larger $I$ corresponds to a stronger easterly trade wind. The details of the
maximum entropy solution and the way to determine $\sigma_{I}(I)$ are included
in Appendix. Note that the decadal variability $I$ also stands for the zonal
SST gradient between the WP and CP regions that directly determines the
strength of the zonal advective feedback, which is the main interaction
between decadal and interannual variabilities in the coupled system. In fact,
in Kang et al. [68], a Walker circulation strength index is defined as the sea
level pressure difference over the CP/EP region ($160^{o}$W-$80^{o}$W,
$5^{o}$S-$5^{o}$N) and over the Indian ocean/WP region ($80^{o}$E-$160^{o}$E,
$5^{o}$S-$5^{o}$N). The monthly zonal SST gradient between the WP and CP
region is highly correlated with this Walker circulation strength index
(correlation coefficient of $\sim 0.85$), suggesting significant air–sea
interaction over the equatorial Pacific. Since the latter is more directly
related to the zonal advective feedback strength over the CP region, the
decadal variable mainly illustrates such a feature.

###  2.3 Seasonal synchronization

Seasonal phase locking is one of the remarkable features of ENSO, which
manifests in the tendency of ENSO to peak during boreal winter and is mainly
related to the pronounced seasonal cycle of the mean state [69, 70]. The
seasonal synchronization is incorporated into the multiscale coupled model
developed above through two simple parameterizations.

First, it has been shown that the climatological SST in the central-eastern
Pacific cools in boreal fall and warms in spring as a result of the seasonal
motion of the Intertropical Convergence Zone (ITCZ), which also modulates the
strength of the upwelling and horizontal advection processes to influence the
evolution of the SST anomalies [71]. Since the cool (warm) SSTs tend to
coincide with decreased (increased) convective activity and upper cloud cover,
a season-dependent damping term, which represents the cloud radiative
feedback, is included for describing such a seasonal variation [44]. More
specifically, two sinusoidal functions are utilized for parameterizing the
otherwise constant $\alpha_{q}$, which appears as $E_{q}=\alpha_{q}T$ in (6).
One sinusoidal function has a period of one year that naturally describes the
seasonal cycle. The other sinusoidal function has a period of half a year that
represents a semiannual contribution to the seasonally modulated variance, as
was suggested by Stein et al. [70].

Second, the increased wind burst activity in the western Pacific during the
boreal winter as a direct response to the increased atmospheric intraseasonal
variability, such as the MJO, is another main contributor of the seasonal
synchronization [38, 72]. Therefore, a sinusoidal function with period of one
year is utilized for parameterizing the seasonal variation of the wind burst
strength coefficient $\sigma_{p}$ in (8).

###  2.4 The nonlinearly coupled multiscale system

The coupled model developed so far is a linear model, despite the state-
dependent noise. However, the linear nature of the model is insufficient in
characterizing some of the key observed dynamical and statistical features of
the ENSO complexity.

From the dynamical point of view, there are at least two major nonlinearities
that are expected to be added to the starting linear model. First, the decadal
variability determines the strength of the zonal advective feedback.
Therefore, it is natural to treat the modulation of the decadal variability on
the ENSO dynamics to be nonlinear, where the decadal variability plays the
role of a multiplicative factor of the zonal advection coefficient. In other
words, the decadal variability $I$ is incorporated into the SST budget
equation (6) and appear in front of the zonal advection coefficient,

$\partial_{t}T=-c_{1}\zeta
E_{q}+c_{1}\eta_{1}(K_{O}+R_{O})+c_{1}I\eta_{2}(K_{O}-R_{O})+c_{1}\eta_{2}c_{2},$
(10)

such that a quadratic nonlinearity is introduced from $I(K_{O}-R_{O})$ (recall
from (7) that $K_{O}$ and $R_{O}$ are the linear combination of $H$ and $U$).
This nonlinearity represents the mechanism that the strengthening of the
Walker circulation in the decadal time scale will trigger more CP events. It
is crucial in simulating the correct occurrence frequencies of both the CP and
the EP El Niños. One additional small constant $c_{2}$ is further added to
(10), which guarantees all the variables have climatology with zero mean since
otherwise the nonlinearity can cause a slight shift of the mean state.

Another nonlinearity incorporated here is the damping coefficient in the SST
equation. Recall that $E_{q}=\alpha_{q}T$ and therefore
$-c_{1}\zeta\alpha_{q}$ is the damping coefficient. Here $\alpha_{q}$ is
parameterized by a quadratic nonlinear function of the CP SST and the spatial
structure of such a nonlinear function is concentrated in the CP area. In
addition, the symmetric axis of this quadratic function has a negative value,
which means a stronger damping is imposed when the CP SST is positive. This
effectively gives a cubic damping in the CP region. The reason to introduce
this nonlinearity is twofold. On one hand, the relationship between the
subsurface temperature and the thermocline depth is more complicated in the CP
region [73] while only a simple shallow water model is utilized here. Thus,
such a nonlinear damping is introduced to parameterize the additional
relationship beyond the capability of the shallow water model. On the other
hand, it is justified from a simple statistical analysis of the observational
data that a linear relationship between the damping and SST anomaly in the CP
region is broken while a cubic nonlinearity fits the data in a more accurate
fashion [49]. The nonzero symmetric axis in parameterizing $\alpha_{q}$ is
also crucial for recovering the correct non-Gaussian statistics of the SST in
the CP region.

##  3 Observational Data Sets and the Definitions of Different Types of the
ENSO Events

###  3.1 Data

The monthly SST data is taken from the GODAS dataset [74]. Anomalies are
calculated by removing the monthly mean climatology of the entire period. The
Niño 4, Niño 3.4, and Niño 3 indices are the average of SST anomalies over the
zonal regions $160^{o}$E-$150^{o}$W, $170^{o}$W-$120^{o}$W and
$150^{o}$W-$90^{o}$W, respectively, together with a meridional average over
$5^{o}$S-$5^{o}$N.

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

The definitions of different El Niño and La Niña events for studying the ENSO
complexity follow those in Kug et al. [16], which are based on the average SST
anomalies during the boreal winter (December–January–February; DJF). When the
EP is warmer than the CP and is greater than $0.5^{o}$C, it is classified as
the EP El Niño. Among this, based on the definitions used by Wang et al. [75],
an extreme El Niño event corresponds to the situation that the maximum of EP
SST anomaly from April to the next March is larger than $2.5^{o}$C. When the
CP is warmer than the EP and is larger than $0.5^{o}$C, the event is then
defined as a CP El Niño. Finally, when either the CP and EP SST anomaly is
cooler than $-0.5^{o}$C, it is defined as a La Niña event.

##  4 Model Simulation Results

The numerical solution of the model is calculated utilizing the forward Euler
time integration scheme with a time step of $0.5$ days for the interannual
variabilities. The Euler-Maruyama scheme is adopted to compute the stochastic
processes of the decadal variability and the intraseasonal wind bursts, with a
numerical integration time step being $0.5$ days and $0.05$ days,
respectively. The monthly averaged model output for the interannual and
decadal variabilities is utilized in presenting the dynamical and statistical
results. This has almost no difference with the direct model solution but is
adopted mainly for the purpose of being consistent with the monthly averaged
observational data. On the other hand, the monthly average is not applied to
the wind burst data.

###  4.1 Model simulation of the ENSO complexity

Figure 3 shows a 50-year model simulation. With the random and nonlinear
components in the model, the resulting atmosphere-ocean-SST fields exhibit
irregular spatiotemporal patterns, mimicking the observations shown in Panel
(b) of Figure 1. To begin with, the model simulation succeeds in reproducing
both the realistic CP (e.g., years 154, 157, 160) and EP (e.g., years 163,
175, 192) events as well as some mixed events (e.g., year 188). In fact, the
two separate linear solutions presented in Figure 2 are now linked by the
decadal variability, which directly modulates the strength of the zonal
advective feedback. In other words, the decadal variability advances the model
to have a preference towards either the EP or the CP mode at each time
instant, although the details of each single ENSO event are still largely
affected by other stochastic and nonlinear effects. Next, the spatiotemporal
fields of the interannual atmosphere wind $u$, ocean current $U$, thermocline
depth $H$ and the SST $T$ in Panels (a)–(d) as well as the wind bursts
strength $a_{p}$ in Panel (f) reveal distinct formation mechanisms for the CP
and EP El Niños. The EP El Niño, especially the extreme EP El Niño, is
triggered by the random wind bursts and the thermocline depth plays an
important role in the event development. In contrast, the zonal advection is
the dominant contributor to the CP events. In addition to the response of the
CP and EP events to the zonal advective and thermocline feedbacks, the
convergence center of the interannual atmosphere wind locates in the CP and EP
regions when these two types of the events occur, respectively. These causal
relationships are consistent with observations and the previous findings [12,
16, 47, 21, 76]. Furthermore, as in the observations, the strength of the CP
El Niños is overall weaker than that of the EP ones [21]. Particularly,
extreme El Niño events are only observed in the EP region, due to the
anomalously strong wind bursts. It is also noticed that the probability of
generating CP events increases as the decadal variability becomes stronger.
This is again consistent with the observations, for example, the CP events
becoming more frequent as the strengthening of the Walker circulation in the
21st century [46, 47]. Nevertheless, regardless of the strength of the decadal
variability, the model always allows both the CP and the EP events to be
triggered with a certain chance. Finally, the Niño indices shown in Panel (f)
mimic the reality, where the Niño 4 index has a slightly larger value than
Niño 3 at the CP El Niño phases while the Niño 3 index becomes much bigger
than Niño 4 during the occurrence of the extreme EP events.

Figures 4–5 show a simulation of the SST field for 200 consecutive years,
accompanying with the associated wind bursts and the decadal variability. To
summarize the findings in these figures, Table 1 lists examples of different
ENSO events belonging to 9 refined categories in such a long model simulation
and are compared with observations. The results indicate the ability of the
model in reproducing the realistic ENSO complexity. First, the model is able
to simulate various EP El Niño events with different strengths. In addition to
the moderate EP El Niños, the extreme El Niño events, which appear as a result
of the strong WWBs generated from the intraseasonal model, are also reproduced
by the model. It is worthwhile to highlight that the so-called delayed super
El Niño, as was observed in 2014-2015 [60, 77, 78, 79], are realistically
simulated by the model, for example during model years 905-906. In fact, the
model succeeds in recovering the associated peculiar westerly-easterly-
westerly wind burst structure that is the key mechanism to trigger such an El
Niño event. Here, the initial WWB tends to trigger a strong El Niño but the
subsequent EWB kills the event and postpones it until the next year when
another series of strong WWBs occur. Next, the model generates many realistic
CP El Niño events. In particular, both single-year (e.g., years 764 and 799)
and multi-year (e.g., years 760-761 and 899-900) CP El Niño events can be
reproduced from the model. The latter mimics the observed CP episodes, for
example during 2018-2020. In addition to those events that clearly belong to
either the EP or the CP categories, the model also creates some mixed CP-EP
events (e.g., years 785 and 939), which are similar to the observed ones in
early 1990s (e.g., the one occurred in year 1992). Finally, the La Niña events
from the model usually follow the El Niño ones as the discharge phase. Some La
Niña events have cold SST in the CP region while others have cold centers
locating around the EP area. The model is also able to simulate multi-year La
Niña events, namely, a La Niña transiting to another La Niña, such as the one
spans over years 774-775 and 902-904, mimicking the observed events 1999-2000
and 1984-1986, respectively.

###  4.2 Comparison of the statistics between model simulations and
observations

In addition to the dynamical properties, the model statistics is another
important measurement for assessing its skill in reproducing the realistic
ENSO features. Since the focus is on the ENSO complexity, it is essential to
study various statistics that represent unique aspects of the ENSO
characteristics in different Niño regions across the equatorial Pacific. To
this end, four statistical quantities with respect to the SST anomalies of the
model simulation are compared with those of the observations in Niño 4, Niño
3.4 and Niño 3 regions, respectively. They are 1) the probability density
function (PDF), 2) the seasonal variance, 3) the power spectrum and 4) the
autocorrelation function (ACF). Here, the statistics of the observations are
computed based on the observed SST between 1951 and 2020, which contains 70
years. On the other hand, a long simulation with in total 3500 years is
utilized for computing the model statistics. The total simulation is divided
into 50 non-overlapping subperiods, each having the same length as the
observation. The statistics is then computed for each of these 50 subperiods,
the difference among which reflects the uncertainty in computing these
statistics.

Panel (a) of Figure 6 shows that the model almost perfectly recovers the
strong non-Gaussian statistics of the SST anomalies in all the three Niño
regions. In particular, the observed Niño 3 SST has a positive skewness and a
one-sided fat tail that is associated with the occurrence of the super El
Niños in the EP region. Due to the state-dependent noise in the wind burst
process, the model is able to trigger the extreme events and thus accurately
reproduce such a highly non-Gaussian distribution. On the other hand, the
skewness of the observed Niño 4 SST is negative, and the kurtosis is less than
the standard Gaussian value 3, indicating the suppression of extreme El Niño
events in the CP region. With the help of the cubic and non-centered damping
in the CP area, the model succeeds in capturing such a skewed and light-tailed
distribution. Similarly, the model is able to reproduce the PDF of the Niño
3.4 SST, which demonstrates a slight positive skewness. It is worthwhile to
highlight that, despite the skill in recovering many dynamical features of
ENSO, the CGCMs and many other intermediate models often have difficulties in
simultaneously capturing the highly non-Gaussian PDFs in all the three Niño
regions, which are nevertheless one of the most important and necessary
conditions for simulating the realistic ENSO complexity. Panel (b) of Figure 6
reveals the capability of the model in recovering the observed seasonal
synchronizations of ENSO, which is represented by the monthly variance of the
SST in different Niño regions. The observed ENSO events have their onset in
boreal spring, develop in summer, and peak in the subsequent winter. These
features are overall well captured by the model, especially given the fact
that the model only exploits simple sinusoidal functions for parameterizing
the seasonal effects.

Next, Panel (a) of Figure 7 shows the power spectrums of the SST. It can be
seen that the major signal of the power associated with the Niño 4 SST is
between 2 and 4 years. The power decreases rapidly outside this window. In
contrast, the signal of the Niño 3 (and Niño 3.4) SST has a wider range, that
is, the power remains significant between the range of 2 and 7 years. These
features are well captured by the model simulations, which are the key
requirements for the model to reproduce the observed irregular oscillations.
In addition, as is shown in Panel (b), the model is able to create very
similar ACFs in different Niño regions as the observations. It indicates that
the model has realistic decaying rate and memory, which are consistent with
nature in the equatorial Pacific and are important prerequisites for the
forecast of the ENSO complexity.

Finally, Figure 8 shows the scatter plot of the equatorial Pacific boreal
winter (DJF) mean SST maxima along the equator for the El Niño events from the
model simulation, which is compared with the observations. Each point
displayed here is a function of the maximum SST and its corresponding
longitude. Panel (a) shows the distribution of the locations of these El Niño
events. The distribution exhibits two major centers. One is near the dateline
and another is in the cold tongue region. These two peaks correspond to the CP
and EP El Niño events, respectively. It should be noticed that the model
distribution here is in general more consistent with the observations than
most of the state-of-art CGCMs, which often exhibit a common bias for a
farther west position [9]. Another desirable feature to highlight is that,
despite the bimodality, there remains a relatively large probabilities of the
event occurrence in the region spanning from the dateline to $120^{o}$W (240).
This reveals that the ENSO diversity is not simply composed by the events that
belong to two separate categories. Instead, there are many mixed EP-CP events.
In fact, according to Panel (b), there are several observed El Niño events
(red dots) locating in this region, indicating that the distribution should be
in the form of a continuum rather than two disjoint sets [80]. This, however,
seems not the case in many CGCM results, as was pointed out by Capontondi et
al. [9]. Next, in terms of the strength, the events with CP SST anomaly peaks
are overall weaker than the corresponding EP ones. While the strongest events
always locate in the eastern Pacific, the EP events can exhibit a wide range
of amplitudes. These characteristics are physically reasonable since there is
more potential for the warming amplitude in the EP due to its climatological
SST being much less than the radiative-convective equilibrium temperature of
about 30oC while the potential for the warming amplitude in the CP is
relatively smaller [81].

###  4.3 Composite analysis

To provide a more quantitative assessment of the model performance on
simulating each type of the ENSO events, Table 2 summarizes the occurrence
frequency of different El Niño and La Niña events (as were defined in Section
3.2) per 70 years. For the El Niño events, although the occurrence frequency
of both the EP and the CP El Niño events from the model (18.2 and 14.7) is
higher than that from the observations (14 and 10), the gap in counting both
types of the El Niños between the model and the observations is just around
one standard deviation of different model simulation segments, which is
nevertheless within a relatively reasonable range. Such a difference mainly
comes from overestimating the number of the multi-year El Niño events in the
model. Except for this overestimation issue, other statistics from the model
simulation are all very similar to those from the observations. First, the
ratio between the numbers of the EP and the CP El Niño events from the model
simulation (55% v.s. 45%) is almost the same as that in the observations (58%
v.s. 42%), which indicates the skill of the model in capturing the overall El
Niño diversity. Second, a total of four extreme El Niño events have occurred
since 1951, namely 1972–1973, 1982–1983, 1997–1998 and 2015–2016, while a
comparable number of 3.1$\pm$2.3 events is found in the model simulations.
Third, the occurrence frequency for the La Niña events (26.1$\pm$4.6) from the
model is very close to that in the observations (24). In particular, the model
and the observations share approximately the same numbers of both the single-
year and multi-year La Niña events. It is worthwhile to remark that, as the
classification of El Niño events in observations is subject to the limited
sample size and suffers from uncertainties associated with varying datasets
[82, 83], the perfect agreement of the occurrence frequency with observations
should not be a strict metric on evaluating the model performance. Therefore,
it can be concluded that the model is overall skillful in reproducing
reasonably accurate numbers of different types of the ENSO events.

Next, Figure 9 exhibits the composites of the DJF mean SST anomalies on the
equatorial Pacific for the EP El Niño, the CP El Niño and the La Niña with
respect to the spatial distribution. Here, all the ENSO events (i.e., during
the total 3500 model years) are used to compare with the observations, instead
of separating it into the smaller segments. It is seen that the composites
from the model simulation are almost identical with those from the
observations in terms of both the spatial patterns and the amplitudes.
Specifically, the warming center locates in the EP and CP regions for the EP
and CP El Niños, respectively, although the simulated EP events are closer to
the coastline of the South America. The cooling center of La Niña is located
between the warming centers of the EP and CP El Niño, which is also in
accordance with the observations. Next, the model succeeds in recovering the
ENSO asymmetry. In other words, the amplitude of the EP El Niño is overall
stronger than those of the CP El Niño and the La Niña [84]. It should be
noticed that accurately reproducing the spatial distributions of the ENSO
events and the ENSO asymmetry is still one of the main challenges for the
state-of-art CGCMS [33].

###  4.4 Sensitivity analysis

What remains is to study the role of each key process in the coupled model,
which facilitates the understanding of the model dynamics. To this end,
several sensitivity tests are carried out in the following.

Let us begin with investigating the role of the decadal variability. In the
standard run, $I$ is driven by a simple stochastic process (9) and its value
of the decadal variability time series varies between $0$ and $1$. In the
sensitivity tests, the model simulations with a fixed $I$ of either $I\equiv
0$ or $I\equiv 1$ are studied. In each test, the strength of the wind bursts
is slightly tuned by multiplying a constant such that the variance of the SST
remains the same as the standard run, which allows a fair comparison for the
occurrence frequency in different scenarios. Note that the decadal variability
is related to the climate change and climate projection. Therefore, the
response of the interannual variability due to the variation of the decadal
variability is of great interest.

First, the decadal variability is set to be zero (i.e., $I\equiv 0$), which
corresponds to the situation with a weakened background Walker circulation as
the period between 1980 and 2000 but towards the more extreme case. In such a
scenario, the model simulation leads to an increase of the El Niño events
(from 32.9 to 42.1 per 70 years) and a decrease of the La Niña ones (from 26.1
to 19.4). More specifically, among the El Niño events, 65.8% events are the EP
type while only 34.2% events remain as the CP El Niño. This means the scenario
with an weakened Walker circulation is more favorable for the EP than the CP
El Niño events, as the zonal advective feedback in this situation is reduced
and the role played by the thermocline feedback becomes dominant. Note that
the CP El Niño events can still be generated because of the stochastic noise.
Similar to the overall occurrence frequency, more multi-year El Niño (from
10.5 to 15.9) and less multi-year La Niña (from 6.9 to 4.6) events are found
in such a case. In addition, about 7.6 extreme El Niño events are produced per
70 years. This is nearly twice as many as that in the standard run and is in
accordance with the observations, where 3 out of the total 4 extreme El Niño
events occurred before 2000. Notably, such a result is consistent with the
climate projection that an increased frequency of extreme El Niño events will
appear due to the greenhouse warming, since a projected surface warming over
the EP is faster than that in the surrounding ocean waters [85]. On the other
hand, if the decadal variability is set to be one (i.e., $I\equiv 1$), then
the model mimics the situation when the Walker circulation and zonal
thermocline slope are relatively strong, similar to the period after 2000. In
such a scenario, more CP El Niño events (from 10 to 20.6) are found in the
model simulation, as a natural consequence of the strengthened zonal advective
feedback. Similarly, less multi-year El Niño (from 10.5 to 9.6) and multi-year
La Niña (from 6.9 to 5.2) events are generated. In addition, there remains
only 1.3 extreme El Niño events per 70 years since the overall occurrence of
the EP events becomes lower. These findings further indicate that the
difference between the positive and negative phases of ENSO is weakened. As a
result, the PDF of the EP SST becomes more towards nearly Gaussian, which is
fundamentally different from the observed non-Gaussian PDF with a fat tail.

The next analysis study is about the effect of the multiplicative noise
$\sigma_{p}(T_{C})$ in the stochastic wind burst process. The multiplicative
noise is one of the main contributors to the asymmetry of the EP type of El
Niño. If an additive noise (i.e., setting $\sigma_{p}$ as a constant) is
adopted instead, then the PDF of the simulated EP SST becomes more symmetric
and Gaussian. This is very different from the observations, where the
amplitude of the extreme El Niño is typically larger than that of the
strongest La Niña.

Finally, the nonlinear damping in the CP region is crucial to the ENSO
dynamics and statistics. According to the observations, the asymmetry with
respect to the SST PDF in the CP is reversed compared with that in the EP.
That is, the amplitude of the negative phase of the CP SST is in general
stronger than that of the positive one, which leads to a negative skewness of
the CP SST PDF. Such a negative skewness is accurately recovered with the help
of the nonlinear function of $\alpha_{q}$ with the non-zero symmetric axis. On
the other hand, the nonlinear damping in the CP region also plays an important
role in suppressing the amplitude of the strong CP events, since the damping
becomes larger as the amplitude of the SST anomaly. As a result, the system is
more in favor of the small and moderate SST anomalies and a reduced kurtosis
appears for the CP SST PDF.

##  5 Conclusions and Discussion

In this paper, a simple multiscale stochastic ICM is developed to capture the
ENSO diversity and complexity. The model highlights the interconnections
between intraseasonal, interannual, and decadal variabilities. It also
exploits suitable stochastic processes to facilitate the realistic simulation
of the ENSO. The model succeeds in reproducing the spatiotemporal dynamical
evolution of different types of the ENSO events. It also accurately recovers
the strongly non-Gaussian probability density function, the seasonal phase
locking, the power spectrum and the temporal autocorrelation function of the
SST anomalies in all the three Niño regions (3, 3.4 and 4) across the
equatorial Pacific. Furthermore, both the composites of the SST anomalies for
various ENSO events and the strength-location bivariate distribution of
equatorial Pacific SST maxima for the El Niño events from the model simulation
highly resemble those from the observations. These desirable features of the
model are particularly important for the purpose of realistically simulating
different ENSO events. They are also the key prerequisites for the unbiased
statistical forecast of the ENSO complexity.

It is worthwhile to point out that the stochastic ICM developed here share
many common features as the conceptual model in Chen et al. [49]. Both models
include three time scales, where the decadal variability modulates the
solution that alternates between the EP- and the CP-dominant regimes while the
intrasesonal wind bursts triggers most of the irregularities and extreme
events. The underlying principles of incorporating stochastic wind bursts and
the nonlinearity into both models also appear in a similar fashion.
Nevertheless, the ICM emphasizes more sophisticated physics and includes many
additional dynamical properties. It also involves spatially-extended
structures, which allow a better understanding and potentially an improved
forecast of the spatiotemporal patterns.

A few important topics are remained as future work. First, the intraseasonal
model adopted here is a one-dimensional simple stochastic process, where the
spatial structure is prescribed and fixed. A more realistic intraseasonal
model can be a spatially-extended (stochastic) model for the wind bursts and
the MJO, for example one of the models in Adames et al. [86], Wang et al. [87]
and Thual et al. [88]. A dynamical intraseasonal model allows the wind bursts
and the MJO to have realistic spatial propagation mechanisms, spanning from
the Indian ocean to the WP. Such a coupled model will also be helpful in
understanding of the coupling between the MJO and ENSO. Second, the model
developed here has a symmetric meridional structure, due to the fact that only
the leading meridional basis function is utilized in the meridionally
truncated system. Yet, both the wind bursts and ENSO have certain meridionally
asymmetric features. Therefore, incorporating additional meridional basis
functions into the model is a natural extension of the current system. Third,
the ICM developed here is applicable to the forecast of different ENSO events.
In particular, the ICM can be combined with the conceptual model in Chen et
al. [49] (and possibly the coupled MJO-ENSO model to be developed as well) for
the multi-model data assimilation and forecast, which advances the
understanding of the role of each model and each component in improving the
forecast of the ENSO complexity.

## References

  * [1] Ropelewski, C. F. & Halpert, M. S.  Global and regional scale precipitation patterns associated with the El Niño/Southern Oscillation.  _Monthly Weather Review_ 115, 1606–1626 (1987).
  * [2] McPhaden, M. J., Zebiak, S. E. & Glantz, M. H.  ENSO as an integrating concept in earth science.  _Science_ 314, 1740–1745 (2006).
  * [3] Bjerknes, J.  Atmospheric teleconnections from the equatorial pacific.  _Monthly Weather Review_ 97, 163–172 (1969).
  * [4] Latif, M. _et al._ A review of the predictability and prediction of enso.  _Journal of Geophysical Research: Oceans_ 103, 14375–14393 (1998).
  * [5] Neelin, J. D. _et al._ ENSO theory.  _Journal of Geophysical Research: Oceans_ 103, 14261–14290 (1998).
  * [6] Zebiak, S. E. & Cane, M. A.  A model El Niño–Southern Oscillation.  _Monthly Weather Review_ 115, 2262–2278 (1987).
  * [7] Schopf, P. S. & Suarez, M. J.  Vacillations in a coupled ocean–atmosphere model.  _Journal of Atmospheric Sciences_ 45, 549–566 (1988).
  * [8] Jin, F.-F.  An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model.  _Journal of the Atmospheric Sciences_ 54, 811–829 (1997).
  * [9] Capotondi, A. _et al._ Understanding ENSO diversity.  _Bulletin of the American Meteorological Society_ 96, 921–938 (2015).
  * [10] Timmermann, A. _et al._ El Niño-Southern Oscillation complexity.  _Nature_ 559, 535–545 (2018).
  * [11] Ashok, K., Behera, S. K., Rao, S. A., Weng, H. & Yamagata, T.  El Niño Modoki and its possible teleconnection.  _Journal of Geophysical Research: Oceans_ 112 (2007).
  * [12] Kao, H.-Y. & Yu, J.-Y.  Contrasting eastern-Pacific and central-Pacific types of ENSO.  _Journal of Climate_ 22, 615–632 (2009).
  * [13] Kim, J.-S., Kim, K.-Y. & Yeh, S.-W.  Statistical evidence for the natural variation of the central pacific el niño.  _Journal of Geophysical Research: Oceans_ 117 (2012).
  * [14] Larkin, N. K. & Harrison, D.  Global seasonal temperature and precipitation anomalies during El Niño autumn and winter.  _Geophysical Research Letters_ 32 (2005).
  * [15] Yu, J.-Y. & Kao, H.-Y.  Decadal changes of ENSO persistence barrier in SST and ocean heat content indices: 1958–2001.  _Journal of Geophysical Research: Atmospheres_ 112 (2007).
  * [16] Kug, J.-S., Jin, F.-F. & An, S.-I.  Two types of El Niño events: cold tongue El Niño and warm pool El Niño.  _Journal of Climate_ 22, 1499–1515 (2009).
  * [17] Chen, D. & Cane, M. A.  El Niño prediction and predictability.  _Journal of Computational Physics_ 227, 3625–3640 (2008).
  * [18] Jin, E. K. _et al._ Current status of ENSO prediction skill in coupled ocean–atmosphere models.  _Climate Dynamics_ 31, 647–664 (2008).
  * [19] Barnston, A. G., Tippett, M. K., L’Heureux, M. L., Li, S. & DeWitt, D. G.  Skill of real-time seasonal ENSO model predictions during 2002–11: Is our capability increasing?  _Bulletin of the American Meteorological Society_ 93, 631–651 (2012).
  * [20] Hu, Z.-Z. _et al._ An analysis of warm pool and cold tongue El Niños: Air–sea coupling processes, global influences, and recent trends.  _Climate Dynamics_ 38, 2017–2035 (2012).
  * [21] Zheng, F., Fang, X., Yu, J.-Y. & Zhu, J.  Asymmetry of the Bjerknes positive feedback between the two types of El Niño.  _Geophysical Research Letters_ 41, 7651–7657 (2014).
  * [22] Fang, X., Zheng, F. & Zhu, J.  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, 4357–4369 (2015).
  * [23] Sohn, S.-J., Tam, C.-Y. & Jeong, H.-I.  How do the strength and type of ENSO affect SST predictability in coupled models.  _Scientific Reports_ 6, 1–8 (2016).
  * [24] Santoso, A. _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, 403–420 (2019).
  * [25] Kug, J.-S., Choi, J., An, S.-I., Jin, F.-F. & Wittenberg, A. T.  Warm pool and cold tongue El Niño events as simulated by the GFDL 2.1 coupled GCM.  _Journal of Climate_ 23, 1226–1239 (2010).
  * [26] Ham, Y.-G. & Kug, J.-S.  How well do current climate models simulate two types of El Niño?  _Climate Dynamics_ 39, 383–398 (2012).
  * [27] Kug, J.-S., Ham, Y.-G., Lee, J.-Y. & Jin, F.-F.  Improved simulation of two types of El Niño in CMIP5 models.  _Environmental Research Letters_ 7, 034002 (2012).
  * [28] Chen, N. & Majda, A. J.  Simple dynamical models capturing the key features of the central pacific el niño.  _Proceedings of the National Academy of Sciences_ 113, 11732–11737 (2016).
  * [29] Fang, X. & Mu, M.  A three-region conceptual model for central Pacific El Niño including zonal advective feedback.  _Journal of Climate_ 31, 4965–4979 (2018).
  * [30] Fang, X. & Zheng, F.  Simulating eastern-and central-pacific type ENSO using a simple coupled model.  _Advances in Atmospheric Sciences_ 35, 671–681 (2018).
  * [31] Lin, J.-L.  The double-itcz problem in ipcc ar4 coupled gcms: Ocean–atmosphere feedback analysis.  _Journal of Climate_ 20, 4497–4525 (2007).
  * [32] Xie, R. & Jin, F.-F.  Two leading enso modes and el niño types in the zebiak–cane model.  _Journal of Climate_ 31, 1943–1962 (2018).
  * [33] Planton, Y. Y. _et al._ Evaluating climate models with the clivar 2020 ENSO metrics package.  _Bulletin of the American Meteorological Society_ 102, E193–E217 (2021).
  * [34] Wittenberg, A. T.  Are historical records sufficient to constrain enso simulations?  _Geophysical Research Letters_ 36 (2009).
  * [35] Harrison, D. & Vecchi, G. A.  Westerly wind events in the tropical pacific, 1986–95.  _Journal of Climate_ 10, 3131–3156 (1997).
  * [36] Vecchi, G. A. & Harrison, D.  Tropical pacific sea surface temperature anomalies, el niño, and equatorial westerly wind events.  _Journal of Climate_ 13, 1814–1830 (2000).
  * [37] Tziperman, E. & Yu, L.  Quantifying the dependence of westerly wind bursts on the large-scale tropical Pacific SST.  _Journal of Climate_ 20, 2760–2768 (2007).
  * [38] Hendon, H. H., Wheeler, M. C. & Zhang, C.  Seasonal dependence of the MJO–ENSO relationship.  _Journal of Climate_ 20, 531–543 (2007).
  * [39] Puy, M., Vialard, J., Lengaigne, M. & Guilyardi, E.  Modulation of equatorial pacific westerly/easterly wind events by the madden–julian oscillation and convectively-coupled rossby waves.  _Climate Dynamics_ 46, 2155–2178 (2016).
  * [40] Gebbie, G., Eisenman, I., Wittenberg, A. & Tziperman, E.  Modulation of westerly wind bursts by sea surface temperature: A semistochastic feedback for enso.  _Journal of the Atmospheric Sciences_ 64, 3281–3295 (2007).
  * [41] Gebbie, G. & Tziperman, E.  Predictability of sst-modulated westerly wind bursts.  _Journal of climate_ 22, 3894–3909 (2009).
  * [42] Levine, A., Jin, F.-F. & McPhaden, M. J.  Extreme noise–extreme El Niño: How state-dependent noise forcing creates El Niño–La Niña asymmetry.  _Journal of Climate_ 29, 5483–5499 (2016).
  * [43] Thual, S., Majda, A. J., Chen, N. & Stechmann, S. N.  Simple stochastic model for El Niño with westerly wind bursts.  _Proceedings of the National Academy of Sciences_ 113, 10245–10250 (2016).
  * [44] Thual, S., Majda, A. J. & Chen, N.  Seasonal synchronization of a simple stochastic dynamical model capturing El Niño diversity.  _Journal of Climate_ 30, 10047–10066 (2017).
  * [45] Lian, T., Chen, D., Tang, Y. & Jin, B.  A theoretical investigation of the tropical Indo-Pacific tripole mode.  _Science China Earth Sciences_ 57, 174–188 (2014).
  * [46] McPhaden, M. J., Lee, T. & McClurg, D.  El Niño and its relationship to changing background conditions in the tropical Pacific ocean.  _Geophysical Research Letters_ 38 (2011).
  * [47] Xiang, B., Wang, B. & Li, T.  A new paradigm for the predominance of standing central Pacific warming after the late 1990s.  _Climate Dynamics_ 41, 327–340 (2013).
  * [48] Power, S. _et al._ Decadal climate variability in the tropical Pacific: Characteristics, causes, predictability, and prospects.  _Science_ 374, eaay9165 (2021).
  * [49] Chen, N., Fang, X. & Yu, J.-Y.  A multiscale model for el niño complexity.  _npj Climate and Atmospheric Science_ 5, 1–13 (2022).
  * [50] Geng, L. & Jin, F.-F.  Enso diversity simulated in a revised cane-zebiak model.  _Frontiers in Earth Science_ 10, 899323 (2022).
  * [51] Majda, A. J. & Chen, N.  Model error, information barriers, state estimation and prediction in complex multiscale systems.  _Entropy_ 20, 644 (2018).
  * [52] Fang, X. & Chen, N.  Quantifying the predictability of enso complexity using a statistically accurate multiscale stochastic model and information theory.  _arXiv preprint arXiv:2203.02657_ (2022).
  * [53] Matsuno, T.  Quasi-geostrophic motions in the equatorial area.  _Journal of the Meteorological Society of Japan. Ser. II_ 44, 25–43 (1966).
  * [54] Gill, A. E.  Some simple solutions for heat-induced tropical circulation.  _Quarterly Journal of the Royal Meteorological Society_ 106, 447–462 (1980).
  * [55] Vallis, G. K.  Geophysical fluid dynamics: whence, whither and why?  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 472, 20160140 (2016).
  * [56] Jin, F.-F.  An equatorial ocean recharge paradigm for enso. part ii: A stripped-down coupled model.  _Journal of the Atmospheric Sciences_ 54, 830–847 (1997).
  * [57] Majda, A. J. & Stechmann, S. N.  The skeleton of tropical intraseasonal oscillations.  _Proceedings of the National Academy of Sciences_ 106, 8417–8422 (2009).
  * [58] Cane, M. A., Sarachik, E., Moore, D. & Sarachik, E.  The response of a linear baroclinic equatorial ocean to periodic forcing.  _Journal of Marine Research_ 39, 651–693 (1981).
  * [59] Majda, A. J.  _Introduction to PDEs and Waves for the Atmosphere and Ocean_ , vol. 9 (American Mathematical Soc., 2003).
  * [60] Hu, S. & Fedorov, A. V.  Exceptionally strong easterly wind burst stalling El Niño of 2014.  _Proceedings of the National Academy of Sciences_ 113, 2005–2010 (2016).
  * [61] Yu, J.-Y. & Kim, S. T.  Identifying the types of major El Niño events since 1870.  _International Journal of Climatology_ 33, 2105–2112 (2013).
  * [62] Dieppois, B. _et al._ ENSO diversity shows robust decadal variations that must be captured for accurate future projections.  _Communications Earth & Environment_ 2, 1–13 (2021).
  * [63] Chen, D. _et al._ Strong influence of westerly wind bursts on El Niño diversity.  _Nature Geoscience_ 8, 339–345 (2015).
  * [64] Freund, M. B. _et al._ Higher frequency of central Pacific El Niño events in recent decades relative to past centuries.  _Nature Geoscience_ 12, 450–455 (2019).
  * [65] Gardiner, C.  _Stochastic methods_ , vol. 4 (Springer Berlin, 2009).
  * [66] Averina, T. & Artemiev, S.  Numerical solution of systems of stochastic differential equations.  _Russian Journal of Numerical Analysis and Mathematical Modelling_ 3, 267–286 (1988).
  * [67] Yang, Q., Majda, A. J. & Chen, N.  ENSO diversity in a tropical stochastic skeleton model for the MJO, El Niño, and dynamic Walker circulation.  _Journal of Climate_ 1–56 (2021).
  * [68] Kang, S. M. _et al._ Walker circulation response to extratropical radiative forcing.  _Science Advances_ 6, eabd3021 (2020).
  * [69] Tziperman, E., Zebiak, S. E. & Cane, M. A.  Mechanisms of seasonal–ENSO interaction.  _Journal of the Atmospheric Sciences_ 54, 61–71 (1997).
  * [70] Stein, K., Timmermann, A., Schneider, N., Jin, F.-F. & Stuecker, M. F.  ENSO seasonal synchronization theory.  _Journal of Climate_ 27, 5285–5310 (2014).
  * [71] Mitchell, T. P. & Wallace, J. M.  The annual cycle in equatorial convection and sea surface temperature.  _Journal of Climate_ 5, 1140–1156 (1992).
  * [72] Seiki, A. & Takayabu, Y. N.  Westerly wind bursts and their relationship with intraseasonal variations and ENSO. Part I: Statistics.  _Monthly Weather Review_ 135, 3325–3345 (2007).
  * [73] Zhao, S., Jin, F.-F., Long, X. & Cane, M. A.  On the breakdown of enso’s relationship with thermocline depth in the central-equatorial pacific.  _Geophysical Research Letters_ 48, e2020GL092335 (2021).
  * [74] Behringer, D. & Xue, Y.  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).
  * [75] Wang, B. _et al._ 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, 22512–22517 (2019).
  * [76] Chen, N., Majda, A. J. & Thual, S.  Observations and mechanisms of a simple stochastic dynamical model capturing el niño diversity.  _Journal of Climate_ 31, 449–471 (2018).
  * [77] Capotondi, A., Sardeshmukh, P. D. & Ricciardulli, L.  The nature of the stochastic wind forcing of ENSO.  _Journal of Climate_ 31, 8081–8099 (2018).
  * [78] Thual, S., Majda, A. J. & Chen, N.  Statistical occurrence and mechanisms of the 2014–2016 delayed super El Niño captured by a simple dynamical model.  _Climate Dynamics_ 52, 2351–2366 (2019).
  * [79] Xie, R. & Fang, X.  The unusual 2014–2016 El Niño events: Dynamics, prediction and enlightenments.  _Science China Earth Sciences_ 63, 626–633 (2020).
  * [80] Johnson, N. C.  How many enso flavors can we distinguish?  _Journal of Climate_ 26, 4816–4827 (2013).
  * [81] Jin, F.-F., An, S.-I., Timmermann, A. & Zhao, J.  Strong El Niño events and nonlinear dynamical heating.  _Geophysical Research Letters_ 30, 20–1 (2003).
  * [82] Wiedermann, M., Radebach, A., Donges, J. F., Kurths, J. & Donner, R. V.  A climate network-based index to discriminate different types of el niño and la niña.  _Geophysical Research Letters_ 43, 7176–7185 (2016).
  * [83] Capotondi, A., Wittenberg, A. T., Kug, J.-S., Takahashi, K. & McPhaden, M. J.  Enso diversity.  _El Niño Southern Oscillation in a changing climate_ 65–86 (2020).
  * [84] Hayashi, M., Jin, F.-F. & Stuecker, M. F.  Dynamics for El Niño-La Niña asymmetry constrain equatorial-Pacific warming pattern.  _Nature communications_ 11, 1–10 (2020).
  * [85] Cai, W. _et al._ Increasing frequency of extreme El Niño events due to greenhouse warming.  _Nature Climate Change_ 4, 111–116 (2014).
  * [86] Adames, Á. F. & Kim, D.  The MJO as a dispersive, convectively coupled moisture wave: Theory and observations.  _Journal of the Atmospheric Sciences_ 73, 913–941 (2016).
  * [87] Wang, B., Liu, F. & Chen, G.  A trio-interaction theory for Madden–Julian oscillation.  _Geoscience Letters_ 3, 1–16 (2016).
  * [88] Thual, S., Majda, A. J. & Chen, N.  A tropical stochastic skeleton model for the MJO, El Niño, and dynamic Walker circulation: A simplified GCM.  _Journal of Climate_ 31, 9261–9282 (2018).
  * [89] Kapur, J. N. & Kesavan, H. K.  Entropy optimization principles and their applications.  In _Entropy and energy dissipation in water resources_ , 3–20 (Springer, 1992).
  * [90] Majda, A. J. & Wang, X.  _Nonlinear dynamics and statistical theories for basic geophysical flows_ (Cambridge University Press, 2006).
  * [91] Branicki, M., Chen, N. & Majda, A. J.  Non-Gaussian test models for prediction and state estimation with model errors.  _Chinese Annals of Mathematics, Series B_ 34, 29–64 (2013).

## Acknowledgements

The research of N.C. is partially funded by the Office of VCRGE at UW-Madison.
The research of X.F. is supported by 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 No. 41805045).

## Appendix-A: Variables and Parameters

The definitions and units of the model variables are listed in Table 3. The
parameter values are summarized in Table 4.

As was stated in Section 2.1, different parabolic cylinder functions in the
ocean and atmosphere were used in the coupled model. Their profiles are shown
in Figure 10. The atmospheric parabolic cylinder functions read
$\phi_{0}(y)=(\pi)^{-1/4}\exp(-y^{2}/2)$,
$\phi_{2}(y)=(4\pi)^{-1/4}(2y^{2}-1)\exp(-y^{2}/2)$. The ocean parabolic
cylinder functions read $\psi_{m}(Y)$, which have the same profiles as the
atmospheric ones but depend on the oceanic meridional axis $Y$.

In order to couple the ocean and atmosphere, projection coefficients are
introduced, which read
$\chi_{A}=\int_{-\infty}^{+\infty}{\phi_{0}(y)\phi_{0}(y/\sqrt{c})\mathrm{d}y}$
and
$\chi_{O}=\int_{-\infty}^{+\infty}{\psi_{0}(Y)\psi_{0}(Y/\sqrt{c}Y)\mathrm{d}Y}$.
The atmosphere uses a truncation of the Kelvin and first Rossby atmospheric
equatorial waves of amplitude $K_{A}$ and $R_{A}$. The ocean uses a truncation
of zonal wind stress forcing to $\psi_{0}$, $\tau_{x}=\tau_{x}\psi_{0}$. This
is known to excite only the Kelvin and first Rossby oceanic waves, of
amplitude $K_{O}$ and $R_{O}$. Similarly, for the SST model a truncation
$\psi_{0}$, $T=T\psi_{0}$ is utilized. Then, the deterministic and linear part
of the ENSO model truncated meridionally yields (4)–(6).

## Appendix-B: Parameters for the Two Different Linear Solutions

The model starts with a deterministic and linear coupled interannual
atmosphere, ocean and SST system. Before the two stochastic processes on the
other two time scales are further incorporated, it is crucial to confirm that
the linear model is able to generate the basic solutions of the two types of
ENSO events under different parameter settings. In addition, the potential
difference should be as small as possible for the physical interpretation.

In this model, the behavior of the leading mode is mainly determined by the
relative amplitude between the zonal advective feedback and the thermocline
feedback, which is consistent with the observational analyses. For this
purpose, the only change of the setting is the strength of the zonal advective
feedback, i.e., $\eta_{2}$. Specifically, for the linear solution
corresponding to the CP ENSO regime,
$\eta_{2}(x)=\max(0,4\exp(-(x-L_{O}/(7/3))^{2}/0.05)\times 0.9)$ is utilized
(Panel (b) of Figure 2), while for that corresponding to the EP ENSO regime,
$\eta_{2}(x)=\max(0,4\exp(-(x-L_{O}/(7/3))^{2}/0.05)\times 0.9)\times 0.3$ is
adopted (Panel (a) of Figure 2), which is $30\%$ of the one in the CP ENSO
regime.

## Appendix-C: Stochastic Process with Multiplicative Noise for the Decadal
Variability

The decadal variability influences the occurrence frequency of the two types
of El Niño and thus the ENSO complexity. A stochastic model is introduced for
the decadal variability, which depicts the strength of the background Walker
circulation and affects the related zonal advective feedback. In the decadal
model (9), a state-dependent (i.e., multiplicative) noise coefficient
$\sigma_{I}(I)$ is adopted that allows $I$ to be non-negative, which comes
from the fact that the trade wind in the lower level of the Walker circulation
in the decadal time scale is easterly. Here, as only limited data for the
decadal variability is available, a uniform distribution function of $I$ is
adopted in the model. This is based on the fact that the uniform distribution
is the maximum entropy solution for a function in the finite interval without
additional information [89, 90, 91]. The parameter $m$ is the mean of $I$,
which can be inferred directly from the data. The damping parameter $\lambda$
can be determined by taking the inverse of the decorrelation time, which is
defined as

$\tau=\lim_{T\to\infty}\int_{0}^{T}ACF(t)\mathrm{d}t\qquad\mbox{with}\qquad
ACF(t)=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{I(t+t^{\prime})I(t^{\prime})}{var(I)}\mathrm{d}t^{\prime}.$

In practice, a sufficiently large $T$ is taken as a numerical approximation.
Finally, the multiplicative noise coefficient $\sigma_{I}(I)$ is determined in
the following way [66]

$\sigma^{2}_{I}(I)=\frac{-2\lambda}{p(I)}\int_{-\infty}^{I}\left(s-\frac{m}{\lambda}\right)p(s)\mathrm{d}s.$

Coarse category | Refined category | Model (yrs 750-950) | Observations (yrs 1980-2020)
---|---|---|---
EP El Niño |  |  |
| Moderate EP El Niño | 792, 897 | 1987
| Super El Niño | 752, 862 | 1998
| Delayed super El Niño | 918-919, 905-906 | 2014-2015
CP El Niño |  |  |
| Single-year CP El Niño | 764, 799 | 2005
| Multi-year CP El Niño: | 760-761, 899-900 | 2018-2020
| CP La Niña | 770, 848 | 1989
Mixed events |  |  |
| Mixed EP-CP El Niño | 785, 939 | 1992
La Niña |  |  |
| Single-year La Niña | 798, 802 | 2006
| Multi-year La Niña | 774-775, 902-904 | 1999-2000
Table 1: Examples of different ENSO events in the model simulation (from year 750 to year 950; showing in Figures 4–5) and observations (from year 1980 to 2020). Here, two examples from the model simulation and one example from observations are listed for each type of the ENSO events, respectively.  | El Niño |  | La Niña
---|---|---|---
| EP | CP | Extreme | Multi-year |  | Total | Multi-year
Obs | 14 | 10 | 4 | 5 |  | 24 | 8
Model | $18.2\pm 3.5$ | $14.7\pm 3.5$ | $3.1\pm 2.3$ | $10.5\pm 1.6$ |  | $26.1\pm 4.6$ | $6.9\pm 2.6$
Table 2: Occurrence frequency of different ENSO events per 70 years. The observations are based on the period of 1951-2020, which contains 70 years. For the model simulation, the mean value plus and minus one standard deviation computed from the 50 segments is shown for each case. Note that the counted number of the EP and the CP El Niños contains both single-year and multi-year events as well as the extreme events. Therefore, the total number of the El Niño events is simply the summation of the numbers in the first two columns (e.g., 24 in observations).  Variable | Unit | Unit Value
---|---|---
$x$ zonal axis | $[y]/\delta$ | 15000km
$y$ atmospheric meridional axis | $\sqrt{c_{A}/\beta}$ | 1500km
$Y$ oceanic meridional axis | $\sqrt{c_{O}/\beta}$ | 330km
$t$ interannual time axis | $[t]$ | 34days
$u$ zonal wind speed | $\delta c_{A}$ | 5ms-1
$v$ meridional wind speed | $\delta[u]$ | 0.5ms-1
$\theta$ potential temperature | 15$\delta$ | 1.5K
$E_{q}$ latent heating | $[\theta]/[t]$ | 0.45K.day-1
$U$ zonal current speed | $c_{O}\delta_{O}$ | 0.25ms-1
$V$ meridional current speed | $\delta\sqrt{c}[U]$ | 0.56cms-1
$H$ thermocline depth | $H_{O}\delta_{O}$ | 20.8m
$T$ sea surface temperature | $[\theta]$ | 1.5K
$a_{p}$ wind burst amplitude | $[u]$ | 5ms-1
Table 3: Model variables, definitions and units. Parameter | Value
---|---
$\epsilon$ Froude number | 0.4
$\delta$ long-wave scaling | 0.1
$\delta_{O}$ arbitrary constant | 0.1
$c_{A}$ atmospheric phase speed | 50ms${}^{-}1$
$c_{O}$ oceanic phase speed | 2.5ms${}^{-}1$
$c$ ratio of oceanic/atmospheric phase speed | 0.05
$c_{1}$ modified ratio of phase speed | 0.15
$\beta$ beta-plane parameter | 2.28 $10^{-11}$m-1s-1
$g^{\prime}$ reduced gravity | 0.03ms-2
$H_{O}$ mean thermocline depth | 50m
$\rho_{O}$ ocean density | 1000 kg.m-3
$\chi_{A}$ atmospheric meridional projection coefficient | 0.31
$\chi_{O}$ oceanic meridional projection coefficient | 1.38
$L_{A}$ equatorial belt length | 8/3
$L_{O}$ equatorial Pacific length | 1.2
$\bar{Q}$ mean vertical moisture gradient | 0.9
$\bar{T}$ mean SST | 16.6 (which is 25oC)
$\alpha_{q}$ latent heating factor | $\alpha_{q}=q_{c}q_{e}\exp(q_{e}\bar{T})/\tau_{q}\times\beta_{1}(T)\times\beta_{2}(t)$
$\beta_{1}(T)$ state dependent component in $\alpha_{q}$ | $\beta_{1}(T)=1.8-\eta_{2}/3+(0.2+|T_{C}+0.4|\times\eta_{2})^{2}/5$
$\beta_{2}(T)$ seasonal dependent component in $\alpha_{q}$ | $\beta_{2}(T)=1+0.5\sin(2\pi(t-1/12))+0.1\sin(2\pi t)\eta_{2}$
| $\qquad\qquad-0.0625\sin(4\pi(t-3/12))\eta_{1}$
$q_{c}$ latent heating multiplier coefficient | 7
$q_{e}$ latent heating exponential coefficient | 0.093
$\tau_{q}$ latent heating adjustment rate | 15
$\gamma$ wind stress coefficient | 6.53
$r_{W}$ western boundary reflection coefficient | 0.5
$r_{E}$ eastern boundary reflection coefficient | 1
$\zeta$ latent heating exchange coefficient | 8.7
$c_{2}$ mean correction coefficient | 0.1
$\eta$ profile of thermocline feedback | $\eta(x)=1.3+(1.1\times\tanh(7.5(x-L_{O}/3)))$
$\eta_{2}$ profile of zonal advective feedback | $\eta_{2}(x)=\max(0,4\exp(-(x-L_{O}/(7/3))^{2}/0.1)\times 0.9)$
$d_{p}$ wind burst damping | 1.12 (which is 1mon-1)
$s_{p}$ wind burst zonal structure | $s_{p}(x)=\exp(-45(x-L_{O}/4)^{2})$
$\sigma_{p}(T_{C})$ wind burst noise coefficient | $\sigma_{p}(T_{C})=1.6(\tanh(T_{C})+1)(1+0.6\cos(2\pi t))$
| $(1-0.75I)$
$\lambda$ damping of decadal variability | 0.0186 (which is 5year-1)
$m$ mean of $I$ | 0.5
Table 4: Model parameter values. Figure 1: Panel (a): A schematic illustration
of the multiscale model developed here. Panel (b): The observational SST
anomaly from 1980 to 2020 (unit: oC). It is based on the GODAS dataset [74]
and is computed by averaging over 5oS to 5oN followed by removing the monthly
mean climatology of the entire period. Figure 2: Linear solutions of the
coupled system (4)–(6) reconstructed utilizing the leading two eigenmodes,
which have the slowest decaying rate. These two modes appear as a pair of
complex conjugate and therefore the reconstructed spatiotemporal pattern is
real-valued. Panel (a) shows the solution by multiplying a small number to the
ocean zonal advective feedback $\eta_{2}$ coefficient to lower its role and
thus gives a EP El Niño dominant mode. Panel (b) shows the solution of the
system with a stronger zonal advective feedback and leads to a CP El Niño
dominant mode. In both panels, the four columns present the hovmoller diagrams
of the interannual atmosphere wind $u$ (unit: m/s), the ocean current $U$
(unit: m/s), the thermocline depth $H$ (unit: m) and the SST $T$ (unit: oC).
The detailed parameter values corresponding to the results here are listed in
Appendix. Note that, for the purpose of illustration, the decaying rate is
manually set to be zero in demonstrating the spatiotemporal pattern of such
linear solutions here. Figure 3: A 50-year model simulation of different
fields. Panels (a)–(d): Hovmoller diagrams of the interannual atmosphere wind
$u$ (unit: m/s), the ocean current $U$ (unit: m/s), the thermocline depth $H$
(unit: m) and the SST $T$ (unit: oC). The longitude ranges from $120^{o}$E
(120) to $80^{o}$W (280). Panel (e): time series of the decadal variability
$I$. Panel (f): time series of the intraseasonal random wind bursts $a_{p}$
(unit: m/s). Panel (g): Niño 4, Niño 3.4 and Niño 3 indices.  Figure 4: Model
simulation: Hovmoller diagram of SST $T$ (unit: oC), time series of the
decadal variability $I$ and time series of the intraseasonal random wind
bursts $a_{p}$ (unit: m/s) from year 750 to year 800 and from year 800 to year
850. Figure 5: Similar to Figure 4 but for the periods from year 850 to year
900 and from year 900 to year 950.  Figure 6: Comparison of the PDFs (Panel
(a)) and the seasonal variances (Panel (b)) of SST anomalies between the model
simulation and the observations in different Niño regions. The observations
are based on the period of 1951-2020, which contains 70 years.
Correspondingly, the model simulation has the length in total 3500 years and
is divided into 50 equally long periods, each of which is 70 years. The blue
shading area is the one standard deviation of the statistics computed from
these 50 non-overlapped periods and the blue solid curve is the average value.
The observational statistics is shown in red solid curve. Figure 7: Similar to
Figure 6 but for the spectrums (Panel (a)) and the autocorrelation functions
(ACFs; Panel (b)). Figure 8: Distribution of equatorial Pacific SST maxima for
the El Niño events from the model simulation of 3500 years (blue) and the
observations (red). For each of the qualified El Niño events, the winter-mean
SST anomalies are averaged over the equatorial zone (from 5oS to 5oN), and
then the Pacific zonal maximum is located. (a) Distribution of peak SST
anomaly longitudes. (b) Scatter plot of the peak SST anomaly value v.s. the
longitude at which it occurs. The blue (red) dots are for the model results
(observations). (c) Distribution of peak SST anomaly values. Figure 9:
Composited winter (DJF) mean SST anomalies (lines) and the corresponding error
bars (i.e., shaded by one standard deviation) over the equatorial Pacific for
EP El Niño (black), CP El Niño (red) and La Niña (blue) events. (a) and (b)
are for model and observation, respectively. Figure 10: Panels (a)–(d):
Spatial structure functions of the meridional bases $\phi_{0}(y)$,
$\phi_{2}(y)$, $\psi_{0}(Y)$ and $\psi_{2}(Y)$. Panels (e)–(g): Spatial
structure functions of $s_{p}(x)$, $\eta_{1}(x)$ and $\eta_{2}(x)$.