# Statistical Response of ENSO Complexity to Initial Value and Model Parameter
Perturbations

###### Abstract

Studying the response of a climate system to perturbations has practical
significance. Standard methods in computing the trajectory-wise deviation
caused by perturbations may suffer from the chaotic nature that makes the
model error dominate the true response after a short lead time. Statistical
response, which computes the return described by the statistics, provides a
systematic way of reaching robust outcomes with an appropriate quantification
of the uncertainty and extreme events. In this paper, information theory is
applied to compute the statistical response and find the most sensitive
perturbation direction of different El Niño-Southern Oscillation (ENSO) events
to initial value and model parameter perturbations. Depending on the initial
phase and the time horizon, different state variables contribute to the most
sensitive perturbation direction. While initial perturbations in sea surface
temperature (SST) and thermocline depth usually lead to the most significant
response of SST at short- and long-range, respectively, initial adjustment of
the zonal advection can be crucial to trigger strong statistical responses at
medium-range around 5 to 7 months, especially at the transient phases between
El Niño and La Niña. It is also shown that the response in the variance
triggered by external random forcing perturbations, such as the wind bursts,
often dominates the mean response, making the resulting most sensitive
direction very different from the trajectory-wise methods. Finally, despite
the strong non-Gaussian climatology distributions, using Gaussian
approximations in the information theory is efficient and accurate for
computing the statistical response, allowing the method to be applied to
sophisticated operational systems.

The purpose of this work is to better understand how the El Niño-Southern
Oscillation (ENSO) responds to changes in its initial state and internal
dynamics or external forcings. A statistical quantification of this response
allows for the comprehension of the triggering conditions and the effect of
climate change in the occurrence frequency and strength of each type of ENSO
event. Such a study also allows to detect the most dangerous perturbation
directions, which has practical significance in guiding anthropogenic
activities. The approach used to study the response in this work is through
the framework of information theory, which allows for an unbiased and robust
assessment of the statistical response that is not affected by the turbulent
dynamics of the system.

##  1 Introduction

El Niño-Southern Oscillation (ENSO) is the dominant interannual variability
over the equatorial central to eastern Pacific Ocean (Philander, 1983;
Ropelewski and Halpert, 1987; Klein et al., 1999; McPhaden et al., 2006),
which is characterized by its irregular and quasi-periodic anomalies in
atmospheric wind and sea surface temperatures (SST). While the immediate
effects are observed in the tropics and subtropics, ENSO has a significant
impact on the climate, ecosystems, economies, and societies around the globe
via atmospheric pathways (Ropelewski and Halpert, 1987; Klein et al., 1999),
making it a global climate phenomenon leading to tropical cyclones, floods,
and droughts. Under the classical viewpoint, ENSO is regarded as a
climatological phenomenon with oscillatory behavior between two roughly mirror
phases based on its features during its mature phase in the eastern Pacific
(EP) (Jin, 1997). In the warming phase of the SST, it is known as El Niño,
while in the cooling phase, it is known as La Niña. In recent decades, many El
Niño events have been observed to occur in the central Pacific (CP) area.
Therefore, the El Niño phenomenon is further divided into the EP and CP types
(Ashok et al., 2007; Yu and Kao, 2007; Kao and Yu, 2009), where the most
significant SST anomaly (SSTa) is located near the coast of South America and
the dateline, respectively (Yu and Kao, 2007; Kao and Yu, 2009). This is known
as the ENSO diversity (Capotondi et al., 2015). Although ENSO was initially
viewed as an essentially symmetric phenomenon with oscillatory and quasi-
periodic behavior, it demonstrates significant asymmetric spatial patterns and
irregularities, as well as diverse characteristics in its temporal evolution
and peak intensities. The former usually leads to different ENSO categories,
including EP El Niño, CP El Niño, mixed EP-CP events, and La Niña. The latter
consists of standard single-year moderate ENSO events, multi-year events (Yu
and Fang, 2018), extreme El Niños (Chen et al., 2015; Levine et al., 2016;
Capotondi et al., 2018; Sun and Yu, 2009) and delayed El Niños (Hu and
Fedorov, 2016, 2017). These spatiotemporal irregularities are called ENSO
complexity (Timmermann et al., 2018; Hayashi and Watanabe, 2017; Boucharel et
al., 2021).

Due to the strong connections with other climate variabilities, it is of
practical importance to compute the response of each ENSO event to various
perturbations of the climate system. On the one hand, predicting the
corresponding spatiotemporal patterns of an ENSO event due to the perturbed
initial state helps understand its precursors and analyze the triggering
conditions of the event. Such a study is vital in advancing an improved
understanding of ENSO predictability. It also facilitates discovering the
formation mechanisms of extreme El Niños and thus provides possible guidance
to prevent or cope with the associated natural hazards. On the other hand, a
perturbation of either the internal dynamics or the external forcing will also
lead to a change in the resulting ENSO spatiotemporal patterns. Such a
perturbation can be caused by potential climate change and may result in
regime switching. The perturbed system may bring about a different occurrence
frequency for each type of ENSO event. It may also increase the strength of
extreme events and enhance the probability of generating multi-year events.
With an appropriate climate model, a natural way to study these two types of
perturbation problems is to compute the model response to the perturbed
initial conditions and model parameters, respectively. Particularly, detecting
the most dangerous perturbation directions, namely the fastest-growing
perturbations, has practical significance in guiding anthropogenic activities.

A hierarchy of approaches has been developed to study the response to these
perturbations. One of the simplest methods for studying the initial
perturbation is to approximate the nonlinear governing equation by the tangent
linear model and then apply a linear singular vector to find the fastest-
growing perturbation (Lorenz, 1965; Samelson and Tziperman, 2001). To take
into account crucial nonlinear features in computing the model response,
methods based on the nonlinear singular vector have been built (Mu, 2000; Mu
and Wang, 2001). The approaches were later extended to the development of
conditional nonlinear optimal perturbation for studying the response of the
ENSO and other climate phenomena (Mu and Duan, 2003; Duan and Mu, 2018). In
addition, many studies exploited intermediate or global circulation models to
carry out numerical experiments that compare the model realizations under
various perturbations (Cai et al., 2020; Mayer et al., 2016; Toniazzo et al.,
2008; Callahan et al., 2021). Most existing methods aim to study the
trajectory-wise difference caused by perturbations. These methods provide
valuable insights for short-term behavior and lead to many successful results.
However, since nature is turbulent, applying trajectory-wise methods for
studying the response may not always be an optimal approach. One of the
fundamental issues in many complex turbulent systems is model error, which is
inevitable in practice. Due to the turbulent nature, model error can easily be
amplified after a very short term (Chen, 2023). Consequently, model error
dominates the intrinsic dynamics in driving the time evolution of the model
trajectory. In other words, a large portion of the computed response is
attributed to the model error. Furthermore, as stochastic parameterizations
have been incorporated into many climate systems (Palmer et al., 2009; Berner
et al., 2017; Franzke et al., 2015), each single model trajectory contains
randomness, which raises difficulties in using standard trajectory-wise
approaches to evaluate the response.

Statistical response, which computes the deviation of the model statistics
instead of trajectories, provides an alternative way to study the model
response due to the perturbation of the initial state or parameters (Majda and
Qi, 2018, 2019). One significant advantage of the statistical response is that
although each model trajectory is chaotic, the time evolution of the
statistics is deterministic and more predictable (Gardiner et al., 1985; Chen,
2023). These statistics are also robust to the random noises in the underlying
system. As a result, a small model error will not substantially impact the
response of the statistics, which is fundamentally different from its
trajectory-wise counterpart. Notably, the statistical response not only
describes the shift of the mean state but also measures the increase or
decrease of the level of uncertainty. The latter can be roughly reflected in
the variance or, more precisely, characterized by the probability density
function (PDF) of the state variables (Billingsley, 2017). Understanding the
corresponding change of the uncertainty as a response to perturbations is
essential in studying the predictability in ensemble prediction and advancing
the probabilistic forecast of extreme events (Fang and Chen, 2023). Although
several statistical methods have been utilized to assess prediction skill,
intrinsic predictability, and model error (DelSole, 2004; DelSole and Tippett,
2007; Kleeman, 2011; Majda et al., 2005), analyzing the statistical response
of ENSO complexity has not been systematically studied.

In this paper, a mathematical framework for computing the statistical response
of a complex system is developed, where information theory is utilized to
measure the strength of the response. The method is then applied to study the
response of different ENSO events, namely the ENSO complexity, to the
perturbations of initial conditions and model parameters. It is also used to
find the most sensitive perturbation direction for each ENSO event. The focus
is on highlighting the advantage of the statistical response over the
trajectory-wise approaches, especially for computing the response in
characterizing the uncertainty and extreme events. A recently developed
conceptual multiscale stochastic model is utilized to describe the ENSO
complexity (Chen et al., 2022). The model can reproduce many crucial observed
dynamical and statistical features, including the non-Gaussian climatology
statistics. It has also been the building block for developing an intermediate
coupled stochastic dynamical model for the ENSO complexity (Chen and Fang,
2023). Therefore, the conceptual model is an appropriate choice as a first
path for exploring the statistical response of different ENSO events. In
addition to the direct numerical algorithm based on the definition, several
approximate schemes are derived to facilitate the practical calculations of
the statistical response. These methods are applicable to more sophisticated
operational models. The conceptual model will serve as a testbed to validate
these computationally efficient methods.

The rest of the paper is organized as follows. The general framework of the
statistical response and the associated computationally efficient approximate
schemes are introduced in Section 2. The observational data set and the
conceptual multiscale stochastic model for the ENSO complexity are described
in Section 3. The results of the statistical response to initial value and
model parameters are presented in Section 4. The paper is concluded in Section
5.

##  2 Methods

###  2.1 Standard trajectory-wise response

Let us start reviewing the standard trajectory-wise methods for computing the
model response. Denote by $\mathbf{x}_{t}$ the state variable $\mathbf{x}$ at
time $t$ from the original model. Denote by $\mathbf{x}^{\delta}_{t}$ the
corresponding variable once a perturbation $\bm{\delta}$ is imposed on either
the initial conditions or the model parameters. Note that, although the
superscript $\delta$ in $\mathbf{x}^{\delta}_{t}$ is written in the scalar
form for notational simplicity, the perturbation $\bm{\delta}$ is usually a
vector. The response of the system at time $t$ to such a perturbation can be
defined as the distance between these two vectors (Samelson and Tziperman,
2001), namely

$d(\mathbf{x}^{\delta}_{t},\mathbf{x}_{t})=\|\mathbf{x}_{t}-\mathbf{x}^{\delta}_{t}\|,$
(1)

where a standard Euclidean norm is taken to obtain a scalar value of the
distance.

Due to the chaotic nature of the system, one fundamental issue is that model
error can dominate the intrinsic dynamics after a short term. Therefore, it
remains unclear in many situations if the computed distance is attributed to
the actual response or the model error.

###  2.2 Statistical response and the associated information measurement

Different from using the trajectories as the quantity for measuring the model
response, the statistical response focuses on the change in the model
statistics as a response to the perturbations. To this end, denote by
$p(\mathbf{x}_{t})$ the PDF of $\mathbf{x}_{t}$ for the original unperturbed
system and $p^{\delta}(\mathbf{x}_{t})$ the corresponding PDF after the
perturbation. The moments, such as the mean and the variance, can be easily
obtained from the PDF. Unlike the case with two vectors where the point-wise
difference as in (1) can be naturally used as the distance, such a direct
discrepancy between the two PDFs will significantly underestimate the role of
extreme events corresponding to the tail probability. Information theory
provides a practical way to characterize the distance between the two PDFs via
the following relative entropy (Majda and Gershgorin, 2010; Majda et al.,
2005; Kleeman, 2011),

$\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))=\int_{\mathbf{x}_{t}}p^{\delta}(\mathbf{x}_{t})\log\left(\frac{p^{\delta}(\mathbf{x}_{t})}{p(\mathbf{x}_{t})}\right)\mathrm{d}\mathbf{x}_{t},$
(2)

which is also known as Kullback-Leibler divergence or information divergence
(Kullback and Leibler, 1951; Kullback, 1987, 1959). The ratio between the two
PDFs inside the logarithm function quantifies the gap in the tail probability,
resulting in an unbiased way of characterizing the statistical difference. It
allows the relative entropy to be widely utilized to quantify model error,
predictability, and prediction skill (Majda and Gershgorin, 2010, 2011; Majda
and Branicki, 2012; Branicki and Majda, 2012, 2014; Kleeman, 2011, 2002;
DelSole, 2004; Branicki et al., 2013; Branstator and Teng, 2010). Despite the
lack of symmetry, the relative entropy has two attractive features. First,
$\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))\geq 0$ with
equality if and only if $p^{\delta}(\mathbf{x}_{t})=p(\mathbf{x}_{t})$.
Second, $\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))$ is
invariant under general nonlinear changes of variables. These provide an
attractive framework for assessing the discrepancy between the two statistical
quantities. A larger value of
$\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))$ means the
statistical response to the perturbation is more significant.

Since the relative entropy $\mathcal{P}$ in (2) is unbounded, it is
practically useful to introduce a rescaled version defined as

$\mathcal{E}=1-\exp(-\mathcal{P}),$ (3)

It rescales the original relative entropy $\mathcal{P}$ to the interval
$[0,1)$. The rescaled relative entropy $\mathcal{E}$ takes the value of $0$ if
and only if $\mathcal{P}=0$. It approaches $1$ when $\mathcal{P}$ becomes
infinity. The rescaled relative entropy $\mathcal{E}$ remains a monotonically
increasing function in terms of the difference between the two PDFs. The
rescaled relative entropy $\mathcal{E}$ will be used in all the numerical
results shown in this work.

Given a perturbation, the strength of the corresponding statistical response
is computed from the relative entropy in (2). However, unless the
distributions have desirable features, numerical integration is needed to
calculate the relative entropy, which is a computationally challenging issue.
Furthermore, in the situation of seeking the most dangerous perturbation
directions, an exhaustive search of the entire state space of $\mathbf{x}_{t}$
is needed based on the direct definition of the relative entropy in (2). This
becomes computationally prohibitive when $\mathbf{x}_{t}$ is high dimensional
(Robert et al., 2010; Kuo and Sloan, 2005). The following subsection aims to
provide alternative ways to accelerate the calculations.

###  2.3 Practical numerical approaches for computing the statistical response

####  2.3.1 Gaussian approximation

One practical setup for utilizing the framework of information theory in many
applications arises when both the distributions are Gaussian so that
$p^{\delta}(\mathbf{x}_{t})\sim\mathcal{N}(\bar{\mathbf{x}}_{t}^{\delta},\mathbf{R}_{t}^{\delta})$
and $p(\mathbf{x}_{t})\sim\mathcal{N}(\bar{\mathbf{x}}_{t},\mathbf{R}_{t})$.
In the Gaussian framework,
$\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))$ has the following
explicit formula (Majda and Gershgorin, 2010; Majda and Wang, 2006)

$\begin{split}\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))=&\left[\frac{1}{2}(\bar{\mathbf{x}}^{\delta}_{t}-\bar{\mathbf{x}}_{t})^{\mathtt{T}}(\mathbf{R}_{t})^{-1}(\bar{\mathbf{x}}^{\delta}_{t}-\bar{\mathbf{x}}_{t})\right]\\\
&\qquad+\left[-\frac{1}{2}\log\det(\mathbf{R}^{\delta}_{t}\mathbf{R}_{t}^{-1})+\frac{1}{2}(\mbox{tr}(\mathbf{R}^{\delta}_{t}\mathbf{R}_{t}^{-1})-\mbox{Dim}(\mathbf{x}_{t}))\right],\end{split}$
(4)

where $\mbox{Dim}(\mathbf{x}_{t})$ is the dimension of $\mathbf{x}_{t}$. The
first term in brackets in (4) is called ‘signal’, reflecting the information
gain in the mean but weighted by the inverse of the model variance,
$\mathbf{R}_{t}$, whereas the second term in brackets, called ‘dispersion’,
involves only the covariance ratio,
$\mathbf{R}^{\delta}_{t}\mathbf{R}^{-1}_{t}$. The signal and dispersion terms
are individually invariant under any (linear) change of variables which maps
Gaussian distributions to Gaussians.

For non-Gaussian PDFs, a Gaussian fit using the mean and covariance can always
be adopted to build the approximate Gaussian distributions. Then the explicit
formula in (4) is used to find the approximate statistical response. It is
worth highlighting two things. First, the Gaussian approximation in (4) is
very different from using a linear approximation of the original dynamics,
such as the linear tangent model. The full nonlinear model is still utilized
to obtain the non-Gaussian PDF as the first step. Only the Gaussian statistics
of the non-Gaussian distribution are used in the explicit formula (4).
Therefore, the statistical response still reflects the nonlinear features of
the underlying dynamics. Second, although the Gaussian approximation may lead
to errors in approximating the PDF itself, it may become a valuable surrogate
for finding the most sensitive perturbation direction, corresponding to the
strongest statistical response at a given forecast lead time. Therefore, one
task below compares the statistical response computed from the explicit
formula with a Gaussian approximation with the exact value. The conclusions
based on the conceptual model tests can provide valuable guidelines for more
sophisticated models.

####  2.3.2 Leading-order approximation via Fisher information

Recall that $\bm{\delta}$ is the perturbation vector, which can be the
perturbation of a subset of the state variables for the initial values or a
few selected parameters. It is assumed that the possible range for the
perturbation is within physical meanings and that it is further standardized
as for $\bm{\delta}=\mathbf{0}$ to correspond to the unperturbed
climatological system. Since the perturbation is usually small, the perturbed
PDF can be written as a function of $\bm{\delta}$. Applying a Taylor expansion
of $p^{\delta}(\mathbf{x}_{t})$ with respect to $\bm{\delta}$ in computing the
relative entropy under the tacit assumption that the PDF is differentiable
with respect to the perturbation $\bm{\delta}$ (Majda et al., 2009; Majda and
Wang, 2010; Hairer and Majda, 2010), yields the following leading-order
approximation of the response to the perturbation (Majda and Chen, 2018; Majda
and Gershgorin, 2011),

$\mathcal{P}(p^{\delta}(\mathbf{x}_{t}),p(\mathbf{x}_{t}))=\frac{1}{2}\bm{\delta}\cdot
I(p(\mathbf{x}_{t}))\bm{\delta}+O(|\bm{\delta}|^{3}),$ (5)

where the first term on the right-hand side of (5) is the quadratic form in
$\bm{\delta}$ given by the Fisher information (Williams, 2001; Cover, 1999),

$\bm{\delta}\cdot
I(p(\mathbf{x}_{t}))\bm{\delta}=\int_{\mathbf{x}_{t}}\frac{(\bm{\delta}\cdot\nabla_{\bm{\delta}}p(\mathbf{x}_{t}))^{2}}{p(\mathbf{x}_{t})}\mathrm{d}\mathbf{x}_{t},$
(6)

where the gradients are evaluated at the unperturbed state.

One significant advantage of the quadratic form in (6) is that the most
sensitive perturbation direction, namely the strongest statistical response at
$t$, occurs along the unit direction associated with the largest eigenvalue of
the matrix $I(p(\mathbf{x}_{t}))$. Such an eigenvalue can be easily computed
once the gradients of $p(\mathbf{x}_{t})$ along the directions of the basis
vectors of $\bm{\delta}$ are calculated, which requires only a small number of
evaluations. In contrast, the computationally expensive brute-force search
algorithms, computing the response at all possible directions, have to apply
when the exact formula in (2) or its Gaussian approximation (4) is used to
achieve such a goal. Therefore, the Fisher information provides an efficient
and systematic way in finding the most sensitive perturbation.

####  2.3.3 Fisher information with coarse-grained statistical measurements

In many situations, the observed climatology data are used to compute
unperturbed statistics. Yet, it is worth noting that, due to the limited
amount of data and the possible measurement noise, these data can usually
estimate the first few moments accurately, but the higher-order moments are
very sensitive to the small noise. Therefore, instead of calculating the exact
PDF, the measured leading few moments are typically used to reconstruct the
least biased PDF using the so-called maximum entropy principle (Majda and
Wang, 2006; Bajkova, 1992). It is then used as an approximation of
$p(\mathbf{x}_{t})$ in (6) to compute the most dangerous perturbation
direction. See Majda et al. (2005), Majda and Wang (2006), and Majda and
Gershgorin (2010) for more details.

Denote by
$\mathbf{E}_{L}(\mathbf{x}_{t})=(E_{1}(\mathbf{x}_{t}),\ldots,E_{L}(\mathbf{x}_{t}))$
the $L$ statistical quantities from the observational measurements or model
simulation, for example, the mean, the covariance and up to the $L$-th moment.
The least biased PDF obtained from the maximum entropy principle is denoted by
$p^{{\delta}}_{L}(\mathbf{x}_{t})$. Further denote by
$\overline{\mathbf{E}}_{L}=(\overline{E}_{1},\ldots,\overline{E}_{L})$, where
each component $\overline{{E}}_{l}$ is given by
$\overline{{E}}_{l}=\int_{\mathbf{x}_{t}}{E}_{l}(\mathbf{x}_{t})p^{{\delta}}_{L}(\mathbf{x}_{t})\mathrm{d}\mathbf{x}_{t}$.
The notation ‘overline’ represents the statistical average with respect to the
perturbed PDF, and the resulting $\overline{{E}}_{l}$ is a number that depends
on the perturbed dynamics. Then the quadratic form of $\bm{\delta}\cdot
I(p(\mathbf{x}_{t}))\bm{\delta}$ in (5) can be approximated by (Majda and
Gershgorin, 2010)

$\bm{\delta}\cdot
I(p_{L}(\mathbf{x}_{t}))\bm{\delta}=\bm{\delta}\cdot\left(\nabla_{\bm{\delta}}\overline{\mathbf{E}}_{L}\mathcal{C}^{-1}\nabla_{\bm{\delta}}\overline{\mathbf{E}}_{L}^{\mathtt{\,T}}\right)\bm{\delta},$
(7)

where $\mathcal{C}$ is the $L\times L$ climate correlation matrix

$\mathcal{C}=\overline{(\mathbf{E}_{L}(\mathbf{x}_{t})-\overline{\mathbf{E}}_{L})(\mathbf{E}_{L}(\mathbf{x}_{t})-\overline{\mathbf{E}}_{L})^{\mathtt{T}}},$
(8)

and the Jacobian of the statistical average is evaluated at the unperturbed
state. Note that ${E}_{l}(\mathbf{x}_{t})$ is a function of the state variable
$\mathbf{x}_{t}$ while its statistical average $\overline{{E}}_{l}$ is a
number. For example, if the mean and the variance are adopted as the first two
components of the measurements, then ${E}_{1}(\mathbf{x}_{t})=\mathbf{x}_{t}$
and ${E}_{2}(\mathbf{x}_{t})=(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{2}$.
Correspondingly, $\overline{E}_{1}(\mathbf{x}_{t})=\overline{\mathbf{x}}_{t}$
and
$\overline{E}_{2}(\mathbf{x}_{t})=\overline{(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{2}}$.
The associated first four entries of $\mathcal{C}$ in (8) are given by

$\mathcal{C}_{11}=\overline{(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{2}},\qquad\mathcal{C}_{12}=\mathcal{C}_{21}=\overline{(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{3}},\qquad\mbox{and}\qquad\mathcal{C}_{22}=\overline{(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{4}}-(\overline{(\mathbf{x}_{t}-\overline{\mathbf{x}}_{t})^{2}})^{2}.$

Finally, it is also practically useful to compute the compressed quadratic
form involving fewer measurements, $L^{\prime}\leq L$,

$\bm{\delta}\cdot
I(p_{L}(\mathbf{x}_{t}))\bm{\delta}=\bm{\delta}\cdot\left(\nabla_{\bm{\delta}}\overline{\mathbf{E}}_{L^{\prime}}\mathcal{C}^{-1}\nabla_{\bm{\delta}}\overline{\mathbf{E}}_{L^{\prime}}^{\mathtt{\,T}}\right)\bm{\delta},$
(9)

where
$\mathbf{E}_{L^{\prime}}(\mathbf{x}_{t})=(E_{1}(\mathbf{x}_{t}),\ldots,E_{L^{\prime}}(\mathbf{x}_{t}),0,\ldots,0))$.
That is, the first $L^{\prime}$ entries of
$\mathbf{E}_{L^{\prime}}(\mathbf{x}_{t})$ are the same as those in
$\mathbf{E}_{L}(\mathbf{x}_{t})$ but the remaining entries are zero. The
compressed quadratic form in (9) is relevant in determining the important
practical information regarding whether, for example, changes in the mean
climate statistics alone determine the most sensitive directions of climate
change.

###  2.4 Summary

In the following, the exact formula of computing the statistical response in
(2) and its two approximations, namely the Gaussian approximation (4) and the
approximation via the Fisher information with the quadratic form (5)–(6) will
be applied to compute the statistical response. Three different computational
approaches will be adopted using the Fisher information with the quadratic
form to find the most sensitive perturbation direction. They are the exact PDF
as in (5)–(6) (hereafter “exact quadratic form”), the quadratic form with
measuring only the mean and covariance (7)–(8) (hereafter “quad form w/ mean
and variance”) and the compressed quadratic form involving only the mean (9)
(hereafter “quad form w/ mean only”). In Section 4, the intercomparison
between these five different computational methods will be carried out.

Note that, although the quad form w/ mean and variance utilizes the Gaussian
statistics, it differs from applying the Gaussian approximation (4) to the
exact formula. The quadratic form via the Fisher information (6) is already an
approximation in computing the statistical response since the Taylor expansion
of the relative entropy in (5) is utilized. Building upon this, the first two
moments are adopted to replace the full PDF as a second approximation.
Nevertheless, it is worth highlighting that the methods exploiting the
quadratic form facilitate determining the most sensitive (i.e., dangerous)
perturbation direction by finding the unit eigenvector corresponding to the
largest eigenvalue of the matrix $I(p(\mathbf{x}_{t}))$ or
$I(p_{L}(\mathbf{x}_{t}))$, hereafter named the maximal eigenvector.

The statistical quantities in this work are computed from an ensemble
simulation. It is based on a Monte Carlo simulation with $3000$ ensemble
members. Such an ensemble size is large enough to reproduce the strong non-
Gaussian climatology PDFs for the 6-dimensional conceptual model described in
Section 3. The computational cost of using $3000$ ensemble members remains low
for such a conceptual model. Therefore, it provides an accurate reference
solution of the statistical response using the exact formula. It further
allows us to compare the result with those using the Gaussian approximations
and the quad form w/ mean and variance. Note that a small ensemble size is
usually sufficient to reach reasonably accurate results when applying the
Gaussian approximation or the Gaussian statistics in the quadratic form.
Therefore, these approximate methods can compute the statistical response
using more sophisticated and higher-dimensional operational systems. The
unperturbed initial value is always given by a Gaussian distribution, centered
at the observational values and equipped with a tiny variance of $10^{-4}$
along each direction. Adding such a small uncertainty to the initial value
facilitates the numerical calculation of the relative entropy within a very
short lead time. Other advanced techniques, such as the fluctuation-
dissipation theorem (Majda et al., 2005; Kubo, 1966), can be embedded into the
above methods to compute the statistical response more efficiently. These
methods are helpful for more complicated systems but are not adopted here
since the direct Monte Carlo simulation with $3000$ ensemble members is
sufficient for the 6-dimensional conceptual model. Finally, when calculating
$\nabla_{\bm{\delta}}p(\mathbf{x}_{t})$ in (6), where $p(\mathbf{x}_{t})$ is
the unperturbed distribution, the derivative is approximated by a second-order
accurate centered finite difference. This is achieved by adding a small
numerical perturbation $\bm{\epsilon}$ with
$\|\bm{\epsilon}\|\ll\|\bm{\delta}\|$ to the initial conditions or model
parameters, and then computing the associated PDF at time $t$. Likewise,
compute another PDF at time $t$ corresponding to the perturbation
$-\bm{\epsilon}$. These two PDFs are then used to approximate the derivative
numerically as

$\frac{p^{\epsilon}(\mathbf{x}_{t})-p^{-\epsilon}(\mathbf{x}_{t})}{2\epsilon}.$

The primary issues to be addressed in this work are the following:

  1. 1.

Study the statistical response of different ENSO events, namely the ENSO
complexity, to the perturbations of initial conditions and model parameters.
Find the most sensitive perturbation direction for each type of ENSO event.

  2. 2.

Since the quad form w/ mean only resembles the deterministic trajectory-wise
response, its difference compared with other methods can be utilized to reveal
the crucial role of the uncertainty in affecting the model response.

  3. 3.

Exploit the skill of the computationally efficient methods involving Gaussian
statistics that apply to operational systems.

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

###  3.1 Observational data sets

The monthly ocean temperature and current data are from the GODAS reanalysis
dataset (Behringer and Xue, 2004). The thermocline depth along the equatorial
Pacific is approximated from the potential temperature as the depth of the
$20^{\circ}$C isotherm contour. The analysis period is 36 years, from the
start of 1982 until the end of 2017. Anomalies presented in this study are
calculated by removing the monthly mean climatology of the whole period. In
this work, the Niño 4 ($T_{C}$) and Niño 3 ($T_{E}$) indices are the averages
of the SSTa over the CP ($160^{\circ}$E-$150^{\circ}$W,
$5^{\circ}$S-$5^{\circ}$N) and EP ($150^{\circ}$W-$90^{\circ}$W,
$5^{\circ}$S-$5^{\circ}$N) regions, respectively. The $h_{W}$ index is the
mean thermocline depth anomaly over the western Pacific region
($120^{\circ}$E-$180^{\circ}$, $5^{\circ}$S-$5^{\circ}$N), while the $u$ index
is the mean mixed-layer zonal current in the CP region.

The daily zonal wind data is measured at 850 hPa and is taken from the NCEP-
NCAR reanalysis (Kalnay et al., 1996). It is used to describe the wind bursts
in the intraseasonal scale. Removing the daily mean climatology, the anomalies
are averaged over the WP region to create the wind burst index. Note that the
wind lies in a faster time scale than all other state variables (daily than
monthly). Although a single daily value of the $\tau$ index describing the
wind anomalies has a minor effect on the SST variables, the accumulated wind
over time will modulate the SST variations. Jumping up to the decadal time
scale, Walker circulation strength index data are also included to illustrate
the modulation of the decadal variation on the interannual ENSO characters. It
is defined as the sea level pressure difference over the CP/EP
($160^{\circ}$W-$80^{\circ}$W, $5^{\circ}$S-$5^{\circ}$N) and the Indian
Ocean/WP ($80^{\circ}$E-$160^{\circ}$E, $5^{\circ}$S-$5^{\circ}$N) (Kang et
al., 2020). The monthly zonal SST gradient between the WP and CP region highly
correlates with this Walker circulation strength index (with a simultaneous
Pearson correlation coefficient of around 0.85), suggesting the significance
of the air-sea interactions over the equatorial Pacific. Since the latter is
more directly related to the zonal advective feedback strength over the CP
region, the decadal model state variable ($I$) mainly illustrates this
quantity.

###  3.2 Definition of different types of ENSO events

The definitions of different ENSO events are based on the average SSTa during
boreal winter (DJF). The CP region is defined as
$160^{\circ}$E-$150^{\circ}$W, $5^{\circ}$S-$5^{\circ}$N, with the former
indicating the longitude (with $180^{\circ}$ being the Prime Meridian) and the
latter the latitude (with $0^{\circ}$ being the Equator), and the EP region as
$150^{\circ}$W-$90^{\circ}$W, $5^{\circ}$S-$5^{\circ}$N. Using the definitions
of (Kug et al., 2009), when the EP is warmer than the CP and the EP SSTa,
$T_{E}$, is greater than $0.5^{\circ}$C, it is classified as an EP El Niño.
Based on the classification in (Wang et al., 2019), an extreme EP El Niño
event corresponds to when the maximum of the EP SSTa from April to the
following March is larger than $2.5^{\circ}$C. Accordingly, when the CP is
warmer than the EP and the CP SSTa, $T_{C}$, is larger than $0.5^{\circ}$C, it
is defined as a CP El Niño. Finally, when either the $T_{C}$ or $T_{E}$
anomalies are cooler than $-0.5^{\circ}$C, it is defined as a La Niña event.

###  3.3 The multiscale stochastic dynamical model for the ENSO complexity

The model used to study the statistical response is a recently developed
stochastic conceptual model for the ENSO complexity (Chen et al., 2022):

$\displaystyle\cfrac{\mathrm{d}u}{\mathrm{d}t}$
$\displaystyle=-ru-\delta_{u}\cfrac{T_{C}+T_{E}}{2}+\beta_{u}(I)\tau+\sigma_{u}\dot{W}_{u},$
(10a) $\displaystyle\cfrac{\mathrm{d}h_{W}}{\mathrm{d}t}$
$\displaystyle=-rh_{W}-\delta_{h}\cfrac{T_{C}+T_{E}}{2}+\beta_{h}(I)\tau+\sigma_{h}\dot{W}_{h},$
(10b) $\displaystyle\cfrac{\mathrm{d}T_{C}}{\mathrm{d}t}$
$\displaystyle=(r_{C}-c_{1}(t,T_{C}))T_{C}+\zeta_{C}T_{E}+\gamma_{C}h_{W}+\sigma(I)u+C_{u}+\beta_{C}(I)\tau+\sigma_{C}\dot{W}_{C},$
(10c) $\displaystyle\cfrac{\mathrm{d}T_{E}}{\mathrm{d}t}$
$\displaystyle=(r_{E}-c_{2}(t))T_{E}-\zeta_{E}T_{C}+\gamma_{E}h_{W}+\beta_{E}(I)\tau+\sigma_{E}\dot{W}_{E},$
(10d) $\displaystyle\cfrac{\mathrm{d}\tau}{\mathrm{d}t}$
$\displaystyle=-d_{\tau}\tau+\sigma_{\tau}(T_{C})\dot{W}_{\tau},$ (10e)
$\displaystyle\cfrac{\mathrm{d}I}{\mathrm{d}t}$
$\displaystyle=-\lambda(I-m)+\sigma_{I}(I)\dot{W}_{I}.$ (10f)

The dimensional units and the parameters in the coupled model are summarised
in Table 2 in the Appendix.

The dynamical core of the model is a deterministic three-region interannual
linear model with zonal advective feedback (Fang and Mu, 2018). It extends on
the classical recharge oscillator model (Jin, 1997) and implements the air-sea
interactions over the entire WP, CP, and EP regions. It also incorporates the
ocean content discharge and recharge process controlling the occurrence of El
Niño and La Niña events via the thermal layer and the ocean zonal advection.
In the model, $T_{C}$ and $T_{E}$ are the SSTa in the CP and EP regions,
respectively, while $u$ is the mean zonal current anomaly in the CP region and
$h_{W}$ is the mean thermocline depth anomaly in the WP region. In addition to
these interannual variabilities, two processes, describing the intraseasonal
wind bursts $\tau$ and the decadal variability in the background Walker
circulation $I$, are further incorporated into the model. The intraseasonal
variability $\tau$ accounts for several important atmospheric ENSO triggers,
such as the westerly wind bursts (WWBs), the easterly wind bursts (EWBs), and
the convective envelope of the Madden-Julian Oscillation (MJO) (Chen et al.,
2015; Hu and Fedorov, 2016; Puy et al., 2016; Vecchi et al., 2006). Its
strength is given by a state-dependent (multiplicative) noise that depends on
the SSTa (Jin et al., 2007; Bianucci et al., 2018; Chen and Zhang, 2023),
where a warmer SSTa leads to stronger wind burst activities. As for the
decadal variability in the state variables through $I$, it stems from the
observation that since 1870, through several detailed El Niño-type
classification methods, the EP and CP El Niño events are alternatively
prevalent every 10 to 20 years (Yu and Kim, 2013; Dieppois et al., 2021). This
oscillation between EP-dominant and CP-dominant regimes indicates that the
decadal variability plays an important role in the underlying dynamics, which
is parameterized through a simple linear stochastic differential equation with
multiplicative noise (10f), with no explicit dependence on the state variables
in the faster time scales (Yang et al., 2021). They, together with additional
small Gaussian white noise $\sigma_{u}\dot{W}_{u}$, $\sigma_{h}\dot{W}_{h}$,
$\sigma_{C}\dot{W}_{C}$, $\sigma_{E}\dot{W}_{E}$, characterize the
irregularity and multiscale features of the ENSO complexity (Timmermann et
al., 2018; Fang and Xie, 2020).

It has been shown in the original work (Chen et al., 2022) that the model can
reproduce many observed properties of the ENSO statistics and ENSO diversity.
In terms of statistics, the model can reconstruct the observed power spectrums
in both the CP and EP regions. It perfectly recovers the climatological PDFs
of the SSTa indices with their strong non-Gaussian statistics. It also
captures the observed seasonal phase-locking features. As for the ENSO
complexity with respect to spatiotemporal patterns, the model can reproduce
roughly the same ratio of EP to CP events and the intensity of these events,
including the amplitude and frequency of the extreme ones, as in observations.
The model can also produce delayed super El Niño and mixed CP-EP events.
Furthermore, the model generates multi-year events with more multi-year La
Niña than multi-year El Niño, consistent with observations (Fang and Yu,
2020).

In the following, the non-dimensional form of the model in (10) is utilized to
compute the responses. In the non-dimensional form, the six state variables
have comparable maximum amplitudes. This means the strength of the components
in the 6-dimensional column vector describing the eigenvector of the quadratic
form can be used to intuitively tell the most sensitive direction.
Nevertheless, the results for the evolution of the statistics with or without
the perturbation are shown in the dimensional form that facilitates the
physical explanations.

##  4 Results

###  4.1 Statistical response to the perturbations of initial values

Figure 1 shows the statistical response of $T_{C}$ to the initial value
perturbations for events with distinct starting dates and at different lead
times, where a 30% perturbation is added to the initial value of each of the
six state variables.

The following conclusions can be made from the result using the exact formula
(2) (Panel (a)). Overall, the statistical response at future time instants
heavily relies on the initial state of the system. The initial perturbation
has the most far-reaching impact on the subsequent evolution of the system
when the initial phase is at the peak of strong EP events, where the response
can remain considerable even after two years. This is unsurprising as a strong
initial value can be easily amplified for a chaotic system as time goes on
(Fang and Chen, 2023). The finding implies that the increased strength of El
Niño events under the climate change scenario not only affects the environment
of those years but also has direct subsequent impacts over a long period. In
contrast, if the initial state lies at a La Niña event, the statistical
response is only significant for a very short period. After that, the system
follows the discharge mechanism, which is more predictable (Sharmila et al.,
2023). In such a case, the difference in the initial value, namely the
perturbation, is often damped quickly. Furthermore, if the system is initially
in a neutral state, then the amplitudes of different state variables are all
near zero. Consequently, the perturbation, which is a percentage of the
initial value, is insignificant and the subsequent response becomes
negligible. Qualitatively similar results are found in $T_{E}$, although the
response becomes more significant for $T_{E}$ when the perturbation is imposed
at the initial phases of strong EP events. Finally, the phase-locking
properties are preserved in the statistical response.

In addition to the general conclusion of the statistical response to initial
perturbations, the intercomparison between different methods in Figure 1
provides the following two crucial findings. First, the statistical response
using the Gaussian approximation (Panel (e)) leads to nearly the same results
as that using the exact formula (Panel (a)). At first glance, this may look
controversial as the Gaussian approximation cannot capture the strong non-
Gaussian climatology statistics. However, despite missing the information in
the higher-order moments, the Gaussian approximation may remain accurate in
computing the relative entropy if the information gap between the two non-
Gaussian distributions behaves similarly to that in the low-order Gaussian
statistics. The comparable patterns in Panel (a) and Panel (e) indicate that
the statistical response can be effectively computed using the much cheaper
Gaussian approximation, which facilitates the use of more sophisticated models
in practice. Second, the three approximations with the quadratic form (Panels
(b)–(d)) overall lead to similar results as that using the exact formula. This
is strong evidence indicating that the leading-order expansion of the relative
entropy with the Fisher information in (5)–(6) is appropriate in computing the
statistical response. Yet, the exact quadratic form (Panel (b)), which uses
the full PDF in computing the gradients (6), gives noisier patterns. This is
because taking the numerical gradient of the full PDF can be sensitive to
small errors in the tail of the estimated distribution that affect the
accuracy of the numerical method in Section 2.4. Therefore, an appropriate
practical strategy involves utilizing the low-order moments, e.g., the
Gaussian statistics, as approximations. This becomes especially helpful when
seeking the most sensitive perturbation directions, as shown below. It is
worth noticing that the quad form w/ mean only gives similar response patterns
as the exact formula within a short time range (less than 6 months). It
results in more significant errors at long lead times and for those years with
larger uncertainties (e.g., 1983-1984, 1987-1989, etc). Such a result
indicates that the mean response is the dominant component in the total
response to initial value perturbations for short lead times. This justifies
using the standard trajectory-wise method in studying the response. As the
lead time increases, some differences can be found between these two methods
at medium-range lead times, which implies that the variance and higher-order
moments take over the role of accounting for the statistical response. The
following case studies will demonstrate such a feature. Finally, the
statistical response decays to zero at a longer period as the chaotic system
only has a finite memory length. Note that the findings here are very
different from the scenario with the parameter perturbation (see below), where
the response in the change of climatology is primarily attributed to the
variance and high-order moments.

Figure 1: The statistical response of $T_{C}$ to the initial value
perturbations for events with different starting dates and at different lead
times, where a 30% perturbation is added to the initial value of each of the
six state variables. Different rows show the resulting statistical response
amplitude measured by the relative entropy using different methods. In each
row, the x-axis is the starting date on the first day of each month across the
36 years, and the y-axis is the lead time (months). The horizontal lines above
the x-axis indicate the event type of that year based on the DJF SSTa.

Figure 2 utilizes the principal coordinate direction (PCD) to demonstrate the
most sensitive direction of perturbation at different lead times computed from
the quad form w/ mean and variance. Note that the exact direction is given by
a $6\times 1$ vector consisting of the components for the state variables
$(u,h_{W},T_{C},T_{E},\tau,I)$. Yet, for the convenience of presentation, only
the most significant component (in absolute value) in this 6-dimensional
vector is used. This is named PCD. Later, the full eigenvector will be used to
describe the most sensitive direction in the case studies. As was mentioned at
the end of Section 3.3, the calculation is based on the non-dimensional system
where the amplitudes of all the state variables are roughly the same. Thus,
the largest component in this vector indeed reflects the dominant direction.
Panels (a)–(b) show the PCD for all events across the 37-year observational
period at different lead times when the statistical responses of the $T_{C}$
and $T_{E}$ are evaluated, respectively. Panels (d)–(e) summarize the
schematic structures of the PCD for different events. As shown in Panel (d),
perturbing $T_{C}$ always gives itself the strongest response at a short lead
time, usually within three months. Yet, before the El Niño events, especially
the strong EP events, thermocline depth becomes the predominant component that
triggers the strongest response of $T_{C}$. A stronger thermocline depth
strengthens the recharge mechanism and changes the SST patterns, affecting the
EP and the CP regions. The direction along the thermocline depth is also the
choice for the initial perturbation to maximize the response in the future at
the interannual time scale. Next, if the perturbation is imposed when a La
Niña event transits to an El Niño one, the zonal advection can play an
essential role in the response at around five months lead time. This is
unsurprising as the advection helps accelerate the recharge process and modify
the patterns at such a time scale (An et al., 1999; Tao et al., 2023). The
most complicated scenario is the multi-year La Niña events, for example,
1999-2000, where $T_{C}$ and $h_{W}$ alternate as the most sensible variable
from 5 to 15 months lead time, depending on when the perturbation is imposed.
In contrast, as is shown in Panel (e), when the statistical response of
$T_{E}$ is considered before an EP El Niño event, perturbing $T_{E}$ gives the
strongest response of itself for a short lead time. The PCD becomes $h_{W}$ as
lead time increases. This simple structure is changed at the phase when an EP
El Niño transits to a La Niña, where the advection again becomes important
around a lead time of 5 months. For multi-year La Nina events, $T_{E}$ and
$h_{W}$ alternate at different starting months as the most crucial variable to
perturb that triggers the strongest response of $T_{E}$. Finally, Panel (c)
shows the PCD based on the statistical response of $T_{E}$ but using the quad
form w/ mean only method. The patterns in Panels (b) and (c) are similar for
lead times of less than 6 months. Some differences can be seen around 7 to 9
months. The difference becomes more significant for longer ranges after 14
months, though the amplitude of the response is negligible since the initial
effect almost goes away at such a long range.

Figure 2: The most sensitive direction of perturbation. In Panels (a)–(b), the
directions are computed using the quad form w/ mean and variance. Each point
in the plot represents the variable associated with the largest component in
the 6-dimensional eigenvector, namely the principal coordinate direction
(PCD). It approximates the most sensitive perturbation direction if a
perturbation is imposed on the corresponding starting date (its x-axis value)
that leads to the response at a given lead time (its y-axis value). The two
panels show the cases when the statistics of $T_{C}$ and $T_{E}$ are adopted,
respectively, in computing the relative entropy. Panel (c) shows the PCD based
on the statistical response of $T_{E}$ using the quad form w/ mean only.
Panels (d)–(e) summarize the schematic structures of the most sensitive
direction for different ENSO events corresponding to the findings in Panels
(a)–(b).

Next, Figures 3–4 present the seasonal statistical response by perturbing the
initial conditions along the most sensitive direction of either $T_{E}$
(Figure 3) or $T_{C}$ (Figure 4). Using Figure 3 as an example, the procedure
for generating such a figure is as follows. First, the statistical response is
computed for different dates and lead times. This will give a plot similar to
that in Panel (a) of Figure 1. The difference compared with Figure 1 is that
in Figure 3, the most sensitive direction for $T_{E}$ at the $k$-month lead is
used to determine the initial perturbation direction when computing the
response PDF at such a lead time. In other words, starting from the same date,
different initial perturbations are adopted to calculate the response at
different lead times. Then, the relative entropy between the response PDF and
the unperturbed one is calculated. Finally, the relative entropy values are
averaged over the dates with the same initial month and lead time to reach the
plots. Panels (a) and (d) in these two figures show that similar to the spring
prediction barrier (Lopez and Kirtman, 2014; Duan and Wei, 2013; Zheng and
Zhu, 2010), there is a spring barrier for the response. Such a spring barrier
is consistent when adding the initial perturbation based on the most sensitive
direction of $T_{E}$ or $T_{C}$. The spring barrier is related to the overall
weak initial strength of the signal in the boreal spring. However, compared
with the spring barrier for the standard trajectory-wise prediction, the
spring barrier for the statistical response is less significant. To understand
such a difference, Panels (b)–(c) and (e)–(f) show the response corresponding
to the signal and dispersion components, respectively. It is seen that the
spring barrier is significant in the signal part. Note that the signal part,
defined in (4), is based on the mean time series, which can be regarded as a
surrogate of a trajectory and is thus more consistent with the standard
trajectory-wise-based spring prediction barrier. On the other hand, the
dispersion part, which is based on the variance, shows no apparent spring
barrier. This weakens the overall spring barrier in the statistical response.
The insignificant spring barrier in the variance response is possibly due to
the weak interaction between the mean and variance, so the time inhomogeneous
behavior in the mean response does not affect the variance too much. The
findings here indicate that the model response can behave differently in
different statistical measurements regarding the spring barrier. The spring
barrier affects the path-wise prediction but may not significantly influence
the propagation of the overall uncertainty.

Figure 3: Seasonal statistical response by perturbing the initial conditions
using the most sensitive direction (MSD) based on $T_{E}$. Panels (a)–(c): the
response of $T_{C}$, including the total response and the response in the
signal and the dispersion, respectively. Panels (d)–(f): the response of
$T_{E}$.  Figure 4: Similar as Figure 3 but using the most sensitive direction
(MSD) based on $T_{C}$.

###  4.2 Statistical response to the perturbations of model parameters

Figure 5 shows the PCD and the statistical response to the perturbations of
model parameters. Note that when the parameters are perturbed, they will
remain unchanged afterwards. Thus, the long-term statistics (i.e., the
climatology distribution) will become different from the unperturbed system.
This is distinct from the initial value perturbation that lasts only for a
finite period due to the finite memory of the system.

The vector containing the parameters for perturbation is $10\times 1$
dimension. The following parameters are considered:

  1. 1.

$\delta_{h}$: the average SSTa feedback in the thermocline depth equation.

  2. 2.

$c_{1}(t,T_{C})$: the damping coefficient in the $T_{C}$ equation,
representing the nonlinear parametrization of the subsurface structure as well
as the discharge behavior of the SSTa in the CP region.

  3. 3.

$\zeta_{C}$: the feedback coefficient of $T_{E}$ in the $T_{C}$ equation.

  4. 4.

$\gamma_{C}$: the thermocline feedback coefficient in the $T_{C}$ equation.

  5. 5.

$\sigma(I)$: the decadal variability coupling parameter (zonal ocean current
coupling coefficient).

  6. 6.

$\beta_{C}(I)$: the wind stress coefficient in the $T_{C}$ equation; also
related to the increase or decrease of the MJO or the tropical cyclones in the
CP region.

  7. 7.

$\gamma_{E}$: the thermocline feedback coefficient in the $T_{E}$ equation.

  8. 8.

$c_{2}(t)$: the damping coefficient in the $T_{E}$ equation, representing the
nonlinear parametrization of the subsurface structure as well as the discharge
behavior of the SSTa in the EP region.

  9. 9.

$\zeta_{E}$: the feedback coefficient of $T_{C}$ in the $T_{E}$ equation.

  10. 10.

$\beta_{E}(I)$: the wind stress coefficient in $T_{E}$ dynamics; also related
to the increase or decrease of the MJO or tropical cyclones in the EP region.

It is seen from Panels (a)–(b) that the wind burst coefficients $\beta_{C}$
and $\beta_{E}$ are the dominant parameter that leads to the strongest
response in the statistics of $T_{C}$ and $T_{E}$, respectively. These results
are as expected, since strengthening the wind activities will significantly
increase the SST amplitude and change the resulting statistics. The
coefficient $\gamma_{E}$ in front of the thermocline depth $h_{W}$ and the
coefficient $\zeta_{E}$ in front of the CP SSTa $T_{C}$ in the equation of
$T_{E}$ are also shown to be important intermittently in affecting the
statistical response of $T_{E}$. The parameter $\zeta_{E}$ and the coefficient
$\delta_{h}$ in front of the SST feedback in the $h_{W}$ equation both affect
the response of $T_{C}$. One notable finding by comparing Panel (b) and Panel
(c) is that the most sensitive perturbation direction using the quad form w/
mean only gives a significantly different result than that using the quad form
w/ mean and variance. It reveals that the uncertainty, reflected by the
variance, plays a crucial role in determining the statistical response. The
strengthening of the wind activity may not change the mean response but will
significantly increase the variance of the response. Subsequently, it enhances
the probability of the occurrence of extreme events. Section 4.3 will include
more detailed studies. Finally, as shown in Panels (d)–(e), Similar to the
initial value perturbation, the calculated responses using the exact formula
and the Gaussian approximation are similar, which again justifies utilizing
the latter to improve the computational efficiency in practice.

Figure 5: The most sensitive direction (in the form of the PCD) and the
statistical response to the perturbations of model parameters. Panels (a)–(b)
show the most sensitive directions for the response in $T_{C}$ and $T_{E}$,
respectively, using the quad form w/ mean and variance. Panel (c) shows that
in $T_{E}$ using the quad form w/ mean only. Panels (d)–(e) show statistical
response of $T_{E}$ by perturbing all the variables by 10% using the exact
formula and the Gaussian approximation.

###  4.3 Case studies

####  4.3.1 Response to initial value perturbations for different ENSO events

Figure 6 shows the statistical response of the 1997 extreme EP El Niño event
(Columns (a)–(b)) and the 1987 moderate EP El Niño event (Column (c)). For the
time evolution of the mean and variance responses shown in the first two rows,
the starting date is 4 months before the event peak. The perturbation at the
initial time corresponds to the one that triggers the most sensitive
perturbation at a lead time of 4 months, namely at the event peak, for $T_{C}$
(Column (a)) and $T_{E}$ (Columns (b)–(c)), respectively. The associated
coordinate of the maximal eigenvector is used as the perturbation added to the
initial condition. The third row of this figure shows the most sensitive
perturbation direction at different lead times, where the lead time of 4
months corresponds to the event peak. Similar representations of the results
are adopted in Figures 7, 8, and 10. Note that the results are robust within a
certain range to the choice of the initial perturbations in terms of the lead
time.

For the 1997 extreme event, perturbing $T_{C}$ and $T_{E}$ results in the
strongest responses for themselves at a short lead time, respectively. This is
as expected since the SST variables do not have time to respond to the
perturbation of other variables within such a short time. At long lead times,
the most sensitive direction is predominated by the thermocline feedback,
which is the crucial variable for amplifying the discharge-recharge mechanism
(Jin and An, 1999). The zonal advection $u$ plays a vital role in the medium-
range lead time from 3 to 9 months when considering the statistical response
of $T_{C}$. It does not significantly impact the statistical response of
$T_{E}$. This is consistent with a recent finding in (Tao et al., 2023), which
suggests the zonal current error has the most substantial impact on the
western and central tropical Pacific. Regardless of whether $T_{C}$ or $T_{E}$
is used to assign the initial perturbation, the time evolution of the mean and
the variance responses are similar, as is seen in the first two rows of
Columns (a) and (b). In both cases, the mean response dominates the time
evolution of the statistical response for $T_{C}$. Starting from July 1997,
the mean response is also the main contributor to the statistical response of
$T_{E}$ for a short lead time, but the variance becomes equally essential in
the total explained response from January to May 1998. A stronger response in
the variance indicates a higher probability of triggering extreme events,
which is never seen in the mean response time series that goes toward the
neutral state. The 1987 moderate EP El Niño event has a similar profile of the
statistical response. The only difference is that the mean response quickly
disappears, but the variance response lasts much longer until the follow-up La
Niña event. This means the perturbation of the precursor of the El Niño event
can have a far-reaching impact on the subsequent years and affect the
associated discharge phase of the cold event.

Figure 6: Statistical responses of the 1997 extreme EP El Niño event (Columns
(a)–(b)) and the 1987 moderate EP El Niño event (Column (c)). The first two
rows show the time evolution of the mean and variance. The starting date is 4
months before the event peak, where the event peak is marked by a black arrow.
The perturbation at the initial time corresponds to the one that triggers the
most sensitive perturbation at a lead time of 4 months, namely at the event
peak, for $T_{C}$ (Column (a)) and $T_{E}$ (Columns (b)–(c)), respectively.
The blue and red curves show the time evolution of the statistics associated
with the unperturbed and the perturbed initial conditions, respectively. The
dashed black curve in the panels of the mean time series shows the single true
trajectory. The third row shows the most sensitive perturbation direction at
different lead times, where the starting date is the same as the first two
rows.

Figure 7 shows the statistical response of the 2014-2015 delayed El Niño event
(Allan et al., 2020; Ludescher et al., 2014; Santoso et al., 2017). The
initial perturbation is imposed on June 2014, 4 months before the first peak
of the event (first and third rows), and on July 2015, 4 months before the
second and the strong peak of the El Niño event (second and fourth rows). The
initial perturbation for the time series in the first two rows corresponds to
the strongest response at a lead time of 4 months using $T_{E}$ while that in
the last two rows uses the response of $T_{C}$. The following conclusions can
be made by comparing the responses in the first and the third rows. Starting
from June 2014, along the most sensitive perturbation direction of $T_{E}$
(first row), the response in the mean of $T_{C}$ and $T_{E}$ dominates the
total response while the response in the variance is negligible. In contrast,
the overall statistical response is relatively weak along the most sensitive
perturbation direction of $T_{C}$ (third row). Next, it can be seen from the
second and the fourth rows that the situation becomes very different when the
starting date becomes July 2015. Both the mean and the variance will
contribute to the statistical response regardless of using $T_{C}$ or $T_{E}$
as the variable for determining the most sensitive perturbation direction.
Notably, when the direction causing the strongest response of $T_{C}$ at the
lead time of 4 months is utilized as the initial perturbation direction (the
fourth row), the time evolution of the statistical response is significant for
both the SST variables, especially in the variance response. This is similar
to the 1997 extreme event shown in Figure 6. The results imply that their
statistical response remains the same despite the intrinsic differences in the
formation mechanisms of these two strong El Niño events. It is also worth
remarking that wind bursts are crucial for each realization of the events
(Chen et al., 2015; Hu and Fedorov, 2016; Levine and McPhaden, 2016; Chen et
al., 2017; Thual et al., 2019, 2016), not only in observations but also in the
model utilized here. Nevertheless, the wind burst activity does not appear to
trigger strong statistical responses of the SST variables. This is because the
intraseasonal wind bursts occur in a much faster time scale than the SST
variables. Therefore, their statistical contribution is averaged out. This
also illustrates a fundamental difference in studying the El Niño events using
trajectory-wise and statistical methods.

Figure 7: Statistical responses of the 2014-2015 delayed El Niño event. The
first four columns show the time evolution of the mean of $T_{C}$, variance of
$T_{C}$, mean of $T_{E}$, and variance of $T_{E}$, respectively. The last
column shows the PCD (as a surrogate for the most sensitive perturbation
direction) at different lead times, where the starting date is the same as the
corresponding first four rows. The first and the third rows show the results
when the initial perturbation is imposed on June 2014, 4 months before the
first peak of the event. The second and the fourth rows show the results when
the initial perturbation is imposed on July 2015, 4 months before the second
(and the strong) peak of the El Niño event. The initial perturbation for the
time series in the first two rows corresponds to the strongest response at a
lead time of 4 months using $T_{E}$ while that in the last two rows uses the
response of $T_{C}$.

Figure 8 includes the cases of the 1992 mixed CP-EP event (Columns (a)–(b))
and the 1995 CP El Niño event (Column (c)). For the mixed CP-EP event, the
first two rows in Column (a) show the time evolution of the mean and variance
for the unperturbed and perturbed initial conditions. The initial perturbation
corresponds to the direction of the strongest response of $T_{C}$ at the lead
time of 4 months. With the first 5 months, the mean response of $T_{C}$
outweighs the variance response while the response in $T_{E}$ is negligible.
However, at the range of the 6- to 14-month lead times, both the mean and
variance responses become significant for $T_{C}$ and $T_{E}$. Note that the
unperturbed system predicts a 50% possibility of a La Niña event in October
1992 (since the mean value is $-0.5^{o}$C), though the true signal does not
exceed this threshold. With the perturbation, the statistical forecast
indicates that a La Niña event will almost surely occur. Thus, such an initial
perturbation is critical since it drives the dynamics from a possibly neutral
state to a La Niña event at a lead time of around one year. Column (b) shows
that if the initial perturbation is imposed based on the strongest response of
$T_{E}$ at the lead time of 4 months, then the perturbation has almost no
influence on the statistical evolution of $T_{C}$ and the dominant component
of the response in $T_{E}$ is in the mean. The comparison between these two
columns implies the role of different initial perturbations in modulating the
dynamics. In addition, the mean and variance responses are not always
synchronized, especially for long lead times, which further indicates the
importance of utilizing the entire PDF or at least considering the variance in
calculating the response. Finally, by imposing the initial perturbation along
the direction of the strongest response of $T_{C}$ at the lead time of 4
months, the time evolution of the statistical response of the 1995 CP El Niño
in Column (c) is similar to the CP-EP event in Column (a). That is, both the
mean and the variance responses are significant for lead times up to 14
months.

Figure 8: Statistical responses of the 1992 mixed CP-EP event (Columns
(a)–(b)) and the 1995 CP El Niño event (Column (c)). The caption description
is similar to Figure 6.

Figure 9 demonstrates the time evolution of the entire PDF and that of the
skewness and kurtosis of $T_{C}$ of the 1992 mixed CP-EP event. The initial
perturbation corresponds to the direction of the strongest response of $T_{C}$
at the lead time of 4 months. This figure supplements Column (a) of Figure 8.
A strong response in the PDF of $T_{C}$ happens at the end of 1992, which is
consistent with that in the leading two moments in Figure 8. Notably, the PDF
associated with the perturbed initial condition can show even more non-
Gaussian features than the one with the unperturbed initial condition.
Nevertheless, since the skewness and kurtosis peak at the same time as the
variance, the result justifies that the Gaussian approximation is appropriate
for finding the most sensitive perturbation direction.

Figure 9: Time evolution of the statistics of $T_{C}$ of the 1992 mixed CP-EP
event when the perturbation is imposed on December 1991. The initial
perturbation corresponds to the direction of the strongest response of $T_{C}$
at the lead time of 4 months. This figure supplements Column (a) of Figure 8.
The first two rows show the time evolution of the PDFs associated with the
unperturbed (blue) and perturbed (red) initial conditions. The third row shows
the time evolution of the third and fourth moments, namely the skewness and
kurtosis, representing the non-Gaussian features.

Finally, Figure 10 shows the case studies of the 2010-2011 multi-year La Niña
event (Iwakiri and Watanabe, 2021; Luo et al., 2017) and the 1988 single-year
La Niña event. For the 2010-2011 multi-year La Niña event, initial
perturbations are imposed on two different dates, October 2010 and November
2011. The two dates are 4 months before the peak of the first and second year
La Niña. For the single-year La Niña event, the perturbation is imposed 4
months before the peak time. As in the previous figures, the initial
perturbation corresponds to the direction of the strongest response of $T_{C}$
or $T_{E}$ at the lead time of 4 months. For these La Niña events, the mean
response is the dominant component of the statistical response. This indicates
that the initial perturbation is dissipated over the discharge phase, which
explains the longer predictability of the La Niña events.

Figure 10: Statistical responses of the 2010-2011 multi-year La Niña event and
the 1988 single-year La Niña event. The caption description is similar to
Figure 7.

####  4.3.2 Response to model parameter perturbations

Figure 11 shows three case studies about the statistical response to model
parameter perturbations. In all cases, perturbation is imposed 4 months before
the event’s peak. The direction of the perturbation is determined by
triggering the largest perturbation of either $T_{E}$ (Column (a)) or $T_{C}$
(Columns (b)–(c)) at the event peak using the quad form w/ mean and variance
method.

Column (a) shows the case of the 1987 moderate EP El Niño. The most sensitive
perturbation direction in terms of the statistics of $T_{E}$ at the lead time
of 4 months completely follows the direction of perturbing the wind stress
coefficient $\beta_{E}$ in the $T_{E}$ equation. Since the wind stress
coefficient only changes the amplitude of the wind but not the preferences of
the WWB or the EWB, the mean response is zero. However, because of the
strengthening of the wind forcing, the response in the variance of $T_{E}$
significantly increases. Consequently, the probability of the occurrence of
extreme events becomes large. Note that despite the coupling between $T_{C}$
and $T_{E}$, the resulting variance response of $T_{C}$ is only significant
for about a year after the initial perturbation.

Column (b) shows the case of the 1995 CP El Niño, where the most sensitive
perturbation direction is determined by maximizing the statistical response of
$T_{C}$ at the lead time of 4 months. As a symmetry to Column (a), the most
sensitive direction is entirely given by the wind stress coefficient
$\beta_{C}$ in the $T_{C}$ equation. Again, the perturbed statistics are
mainly due to the change of the variance rather than the mean, which triggers
more variabilities in the subsequent SST.

Column (c) shows the case of the 1998 La Niña. Since the wind burst is very
weak during the La Niña phase. The wind stress coefficient is no longer the
major contributor to a strong statistical response. The perturbation of the
parameter $\zeta_{C}$, representing the feedback from $T_{E}$ in the $T_{C}$
equation, is the one that triggers the strongest response in the dynamics of
$T_{C}$ in the medium range within 4 to 7 months. In contrast, the feedback
coefficient $\delta_{h}$ accounting for the feedback from $T_{C}$ and $T_{E}$
to the dynamics of $h_{W}$ is the one that triggers the most significant
response for the 8 to 13 months lead time. Finally, the wind stress
coefficient $\beta_{C}$ becomes the dominant factor for a longer range after
the system leaves the La Niña phase.

It is worth highlighting that when computing the most sensitive perturbations
with parameter perturbations, the response in the variance is much more
significant than that in the mean. A strong response in the variance and
higher-order statistical moments implies the increased probability of extreme
events and the enhancement of the uncertainty from the model output. This
indicates the fundamental difference between statistical and traditional
trajectory-wise methods in studying the model response and seeking the most
sensitive direction.

Figure 11: Statistical response to parameter perturbations. The first two rows
show the time evolutions of the mean and variance using the unperturbed (blue)
and perturbed (red) initial condition. The perturbation is imposed 4 months
ahead of the event peak. The direction of the perturbation is determined by
triggering the largest perturbation of either $T_{E}$ (Column (a)) or $T_{C}$
(Columns (b)–(c)) at the event peak using the quad form w/ mean and variance
method. The third row shows the PCD at different lead months starting 4 months
ahead of the event peak. Columns (a)–(c) show the cases of the 1987 moderate
EP El Niño, the 1995 CP El Niño, and the 1988 La Niña events, respectively.

Finally, Table 1 shows the occurrence frequency of different ENSO events in
observations (first row), original system (second row), and the perturbed
systems by perturbing each listed parameter by 30% (third to seventh rows).
The observational period is from 1950 to 2020, totaling 71 years. When
perturbing $\beta_{E}$ and $\beta_{C}$, the EP and CP El Niño events will
become more frequent. Note that the number of extreme EP events is
significantly increased when the wind stress or the wind amplitude is
amplified, corresponding to an enhanced $\beta_{E}$. It also triggers more
multi-year La Niña events. Increased occurrences of the strong El Niño and
multi-year La Niña are related to the global warming scenario (Geng et al.,
2023; Cai et al., 2015). Next, when perturbing $\zeta_{C}$, the feedback
coefficient of $T_{E}$ in the $T_{C}$ equation, the system tends to increase
the number of CP events and decrease the number of EP events, pushing more
activities toward the CP region due to the strengthening of the feedback. A
decrease of $\delta_{h}$, the average SSTa feedback in the thermocline depth
equation, will increase the chance of multi-year El Niño events. Finally,
recall that $c_{1}$ is the damping coefficient in the TC equation,
representing the nonlinear parametrization of the subsurface structure as well
as the discharge behavior of the SSTa in the CP region. If such a parameter is
decreased, then the number of CP events will dramatically increase due to the
more active response in the CP region (Zhao et al., 2021). It will also
trigger more multi-year El Niño events but fewer multi-year La Niña events.

| El Niño (EN) | EP EN | CP EN | Extreme EP | Multi-year EN | La Niña (LN) | Multi-year LN
---|---|---|---|---|---|---|---
Observations | 24 | 14 | 10 | 4 | 5 | 24 | 8
Original system |  22.6$\pm$2.8 |  13.4$\pm$2.4 |  09.2$\pm$2.1 |  4.5$\pm$1.7 |  4.0$\pm$1.5 |  29.9$\pm$3.0 |  07.5$\pm$1.8
Perturbing $\beta_{E}(I)$ |  24.4$\pm$2.8 |  17.0$\pm$2.9 |  07.4$\pm$2.7 |  9.3$\pm$2.4 |  5.4$\pm$1.5 |  33.8$\pm$3.0 |  10.0$\pm$2.1
Perturbing $\beta_{C}(I)$ |  22.2$\pm$3.3 |  12.2$\pm$2.4 |  10.0$\pm$2.6 |  4.3$\pm$2.0 |  3.8$\pm$1.9 |  31.0$\pm$3.0 |  08.2$\pm$2.3
Perturbing $\zeta_{C}$ |  22.3$\pm$2.8 |  11.6$\pm$1.8 |  10.7$\pm$1.9 |  4.1$\pm$1.6 |  3.9$\pm$1.8 |  30.8$\pm$3.4 |  07.5$\pm$2.5
Perturbing $\delta_{h}$ |  22.5$\pm$3.0 |  13.7$\pm$2.9 |  08.8$\pm$2.8 |  4.4$\pm$2.1 |  4.6$\pm$1.9 |  29.8$\pm$4.2 |  07.8$\pm$2.4
Perturbing $c_{1}(t,T_{C})$ |  26.8$\pm$2.9 |  12.6$\pm$2.0 |  14.2$\pm$2.4 |  4.1$\pm$1.8 |  5.3$\pm$1.7 |  31.1$\pm$2.5 |  06.4$\pm$2.0
Table 1: The occurrence frequency of different ENSO events in observations
(first row), original system (second row), and the perturbed systems by
perturbing each listed parameter by 30% (third to seventh rows). The
observational period is from 1950 to 2020, totaling 71 years. A long
simulation is carried out for each case. It is then divided into 30 non-
overlapped segments, each being 71 years as being consistent with the length
of the observations. The mean occurrence frequency numbers plus and minus the
corresponding standard errors of these segments are shown.

##  5 Conclusion

In this paper, a mathematical framework for computing the statistical response
of a complex system is developed, where information theory is utilized to
measure the strength of the response. The method is then applied to study the
response of different ENSO events to the perturbations of initial conditions
and model parameters. It is also utilized to find the most sensitive
perturbation direction for each ENSO event. The main conclusions are
summarized as follows.

  * •

Depending on the initial phase and the time horizon, different state variables
contribute to the most sensitive perturbation direction. While initial
perturbations in SST and thermocline depth usually lead to the most
significant response of SST at short- and long-range, respectively, initial
adjustment of the zonal advection can be crucial to trigger strong statistical
responses at medium-range around 5 to 7 months, especially at the transient
phases between El Niño and La Niña.

  * •

Despite the mean response dominating the total response with a short range for
the initial value perturbation, the variance and higher-order moments
contribute to the response at medium-range lead times.

  * •

The spring barrier in the statistical response is overall weaker than that in
the standard trajectory-wise prediction. Notably, the spring barrier is only
significant in the signal part of the response PDF (corresponding to the mean)
but is not apparent in the dispersion part (corresponding to the variance or
uncertainty).

  * •

The response in the variance triggered by external random forcing
perturbations, such as the wind bursts, often dominates the mean response at
long range in the parameter perturbation scenario, making the resulting most
sensitive direction very different from the trajectory-wise methods.

  * •

Despite the strong non-Gaussian climatology distributions, using Gaussian
approximations in the information theory is efficient and accurate to compute
the statistical response, allowing the method to be applied to more
sophisticated systems, such as the intermediate coupled models in (Chen and
Fang, 2023; Geng and Jin, 2022) or operational systems.

###### Acknowledgements.

The research of N.C. is funded by ARO W911NF-23-1-0118. M.A. is supported as a
graduate research assistant under this grant. The monthly ocean temperature
and current data were downloaded from GODAS
(https://www.esrl.noaa.gov/psd/data/gridded/data.godas.html). The daily zonal
wind data at 850 hPa were downloaded from the NCEP–NCAR reanalysis
(https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html). The code used
in the analysis and generation of the figures in this work can be provided
upon contact with the corresponding author.

### Parameter values and characteristic scales

Here we include the model parameters in the standard run (i.e., unperturbed
system), and the characteristic scales appearing in the three-region
multiscale stochastic model (10a)–(10f).

$\varrho$ (scaling factor) | $0.65$
---|---
$b_{0}$ | $2.5$ | $\mu$ | $0.5$
$\alpha_{1}$ | $0.0625\varrho$ | $\alpha_{2}$ | $0.125\varrho$
$[u]$ |  $1.5$ms-1 | $[h]$ |  $150$m
$[T]$ |  $7.5^{\circ}$C | $[\tau]$ |  $5$ms-1
$[t]$ |  $2$ months | $d_{\tau}$ | $2$
$r$ | $0.25\varrho$ | $r_{C}$ | $0.75b_{0}\mu\varrho/2$
$r_{E}$ | $3r_{C}=2.25b_{0}\mu\varrho/2$ | $\lambda$ | $2/60$
$\delta_{u}$ | $\alpha_{1}b_{0}\mu$ | $\delta_{h}$ | $2\delta_{u}=\alpha_{2}b_{0}\mu$
$\zeta_{C}$ | $0.75b_{0}\mu\varrho/2$ | $\zeta_{E}$ | $0.75b_{0}\mu\varrho/2$
$\gamma_{C}$ | $0.75\varrho$ | $\gamma_{E}$ | $0.75\varrho$
$C_{u}$ | $0.03\varrho$ | $m$ | $2$
$\sigma(I)$ | $I\varrho/5$ | $\beta_{E}(I)$ | $0.15(1-0.2I)\sqrt{\varrho}$
$\beta_{u}(I)$ | $-0.2\beta_{E}(I)$ | $\beta_{h}(I)$ | $-0.4\beta_{E}(I)$
$\beta_{C}(I)$ | $0.8\beta_{E}(I)$ | $\sigma_{u}$ | $0.04\sqrt{\varrho}$
$\sigma_{h}$ | $0.02\sqrt{\varrho}$ | $\sigma_{C}$ | $0.04\sqrt{\varrho}$
$\sigma_{E}$ | $0$ | $\sigma_{I}(I)$ | $\sqrt{\lambda(4-I)I}$
$\sigma_{\tau}(t,T_{C})$ | $0.9[\tanh(7.5T_{C})+1]\times\left[1+0.3\cos\left(\frac{2\pi}{6}t+\frac{2\pi}{6}\right)\right]$
$c_{1}(t,T_{C})$ | $\varrho\left[25\left(T_{C}+\frac{0.75}{7.5}\right)^{2}+0.9\right]\times\left[1+0.3\sin\left(\frac{2\pi}{6}t-\frac{\pi}{6}\right)\right]$
$c_{2}(t)$ | $1.4\varrho\left[1+0.3\sin\left(\frac{2\pi}{6}t+\frac{2\pi}{6}\right)+0.25\sin\left(\frac{2\pi}{3}t+\frac{2\pi}{6}\right)\right]$
Table 2: Parameters of the stochastic conceptual model for the ENSO complexity
(10).

## References

  * Allan et al. (2020) Allan, R. J., J. Gergis, and R. D. D’Arrigo, 2020: Placing the AD 2014–2016 ‘protracted’ El Niño episode into a long-term context. The Holocene, 30 (1), 90–105.
  * An et al. (1999) An, S.-I., F.-F. Jin, and I.-S. Kang, 1999: The role of zonal advection feedback in phase transition and growth of ENSO in the cane-zebiak model. Journal of the Meteorological Society of Japan. Ser. II, 77 (6), 1151–1160.
  * Ashok et al. (2007) Ashok, K., S. K. Behera, S. A. Rao, H. Weng, and T. Yamagata, 2007: El Niño Modoki and its possible teleconnection. Journal of Geophysical Research: Oceans, 112 (C11).
  * Bajkova (1992) Bajkova, A. T., 1992: The generalization of maximum entropy method for reconstruction of complex functions. Astronomical and Astrophysical Transactions, 1 (4), 313–320.
  * Behringer and Xue (2004) Behringer, D., and Y. Xue, 2004: Evaluation of the global ocean data assimilation system at NCEP: The Pacific Ocean. Proc. Eighth Symp. on Integrated Observing and Assimilation Systems for Atmosphere, Oceans, and Land Surface.
  * Berner et al. (2017) Berner, J., and Coauthors, 2017: Stochastic parameterization: Toward a new view of weather and climate models. Bulletin of the American Meteorological Society, 98 (3), 565–588.
  * Bianucci et al. (2018) Bianucci, M., A. Capotondi, S. Merlino, and R. Mannella, 2018: Estimate of the average timing for strong El Niño events using the recharge oscillator model with a multiplicative perturbation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28 (10).
  * Billingsley (2017) Billingsley, P., 2017: Probability and measure. John Wiley & Sons.
  * Boucharel et al. (2021) Boucharel, J., R. Almar, E. Kestenare, and F.-F. Jin, 2021: On the influence of ENSO complexity on Pan-Pacific coastal wave extremes. Proceedings of the National Academy of Sciences, 118 (47).
  * Branicki et al. (2013) Branicki, M., N. Chen, and A. J. Majda, 2013: Non-Gaussian test models for prediction and state estimation with model errors. Chinese Annals of Mathematics, Series B, 34 (1), 29–64.
  * Branicki and Majda (2014) Branicki, M., and A. Majda, 2014: Quantifying Bayesian filter performance for turbulent dynamical systems through information theory. Communications in Mathematical Sciences, 12 (5), 901–978.
  * Branicki and Majda (2012) Branicki, M., and A. J. Majda, 2012: Quantifying uncertainty for predictions with model error in non-Gaussian systems with intermittency. Nonlinearity, 25 (9), 2543.
  * Branstator and Teng (2010) Branstator, G., and H. Teng, 2010: Two limits of initial-value decadal predictability in a CGCM. Journal of Climate, 23 (23), 6292–6311.
  * Cai et al. (2020) Cai, W., A. Santoso, G. Wang, L. Wu, M. Collins, M. Lengaigne, S. Power, and A. Timmermann, 2020: ENSO response to greenhouse forcing. El Niño southern oscillation in a changing climate, 289–307.
  * Cai et al. (2015) Cai, W., and Coauthors, 2015: Increased frequency of extreme La Niña events under greenhouse warming. Nature Climate Change, 5 (2), 132–137.
  * Callahan et al. (2021) Callahan, C. W., C. Chen, M. Rugenstein, J. Bloch-Johnson, S. Yang, and E. J. Moyer, 2021: Robust decrease in El Niño/southern oscillation amplitude under long-term warming. Nature Climate Change, 11 (9), 752–757.
  * Capotondi et al. (2018) Capotondi, A., P. D. Sardeshmukh, and L. Ricciardulli, 2018: The nature of the stochastic wind forcing of ENSO. Journal of Climate, 31 (19), 8081–8099.
  * Capotondi et al. (2015) Capotondi, A., and Coauthors, 2015: Understanding ENSO diversity. Bulletin of the American Meteorological Society, 96 (6), 921–938.
  * Chen et al. (2015) Chen, D., and Coauthors, 2015: Strong influence of westerly wind bursts on El Niño diversity. Nature Geoscience, 8 (5), 339–345.
  * Chen et al. (2017) Chen, L., T. Li, B. Wang, and L. Wang, 2017: Formation mechanism for 2015/16 super El Niño. Scientific reports, 7 (1), 1–10.
  * Chen (2023) Chen, N., 2023: Stochastic Methods for Modeling and Predicting Complex Dynamical Systems. Springer International Publishing, 10.1007/978-3-031-22249-8, URL https://doi.org/10.1007/978-3-031-22249-8.
  * Chen and Fang (2023) Chen, N., and X. Fang, 2023: A simple multiscale intermediate coupled stochastic model for El Niño diversity and complexity. Journal of Advances in Modeling Earth Systems, 15 (4), e2022MS003 469.
  * Chen et al. (2022) Chen, N., X. Fang, and J.-Y. Yu, 2022: A multiscale model for El Niño complexity. npj Climate and Atmospheric Science, 5 (1), 1–13.
  * Chen and Zhang (2023) Chen, N., and Y. Zhang, 2023: Rigorous derivation of stochastic conceptual models for the El Niño-southern oscillation from a spatially-extended dynamical system. Physica D: Nonlinear Phenomena, 453, 133 842.
  * Cover (1999) Cover, T. M., 1999: Elements of information theory. John Wiley & Sons.
  * DelSole (2004) DelSole, T., 2004: Predictability and information theory. Part I: Measures of predictability. Journal of the Atmospheric Sciences, 61 (20), 2425–2440.
  * DelSole and Tippett (2007) DelSole, T., and M. K. Tippett, 2007: Predictability: Recent insights from information theory. Reviews of Geophysics, 45 (4).
  * Dieppois et al. (2021) Dieppois, B., A. Capotondi, B. Pohl, K. P. Chun, P.-A. Monerie, and J. Eden, 2021: ENSO diversity shows robust decadal variations that must be captured for accurate future projections. Communications Earth & Environment, 2 (1), 1–13.
  * Duan and Mu (2018) Duan, W., and M. Mu, 2018: Predictability of el niño-southern oscillation events. Oxford University Press, URL https://oxfordre.com/climatescience/view/10.1093/acrefore/9780190228620.001.0001/acrefore-9780190228620-e-80, 10.1093/acrefore/9780190228620.013.80.
  * Duan and Wei (2013) Duan, W., and C. Wei, 2013: The ‘spring predictability barrier’ for ENSO predictions and its possible mechanism: results from a fully coupled model. International Journal of Climatology, 33 (5), 1280–1292.
  * Fang and Yu (2020) Fang, S.-W., and J.-Y. Yu, 2020: Contrasting transition complexity between El Niño and La Niña: Observations and cmip5/6 models. Geophysical Research Letters, 47 (16), e2020GL088 926.
  * Fang and Chen (2023) Fang, X., and N. Chen, 2023: Quantifying the predictability of ENSO complexity using a statistically accurate multiscale stochastic model and information theory. Journal of Climate, 36 (8), 2681–2702.
  * Fang and Xie (2020) Fang, X., and R. Xie, 2020: A brief review of ENSO theories and prediction. Science China Earth Sciences, 63 (4), 476–491.
  * Fang and Mu (2018) Fang, X.-H., and M. Mu, 2018: A three-region conceptual model for central Pacific El Niño including zonal advective feedback. Journal of Climate, 31 (13), 4965–4979.
  * Franzke et al. (2015) Franzke, C. L., T. J. O’Kane, J. Berner, P. D. Williams, and V. Lucarini, 2015: Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change, 6 (1), 63–78.
  * Gardiner et al. (1985) Gardiner, C. W., and Coauthors, 1985: Handbook of stochastic methods, Vol. 3. springer Berlin.
  * Geng and Jin (2022) Geng, L., and F.-F. Jin, 2022: ENSO diversity simulated in a revised Cane-Zebiak model. Frontiers in Earth Science, 10, 899 323.
  * Geng et al. (2023) Geng, T., F. Jia, W. Cai, L. Wu, B. Gan, Z. Jing, S. Li, and M. J. McPhaden, 2023: Increased occurrences of consecutive La Niña events under global warming. Nature, 619 (7971), 774–781.
  * Hairer and Majda (2010) Hairer, M., and A. J. Majda, 2010: A simple framework to justify linear response theory. Nonlinearity, 23 (4), 909–922, 10.1088/0951-7715/23/4/008.
  * Hayashi and Watanabe (2017) Hayashi, M., and M. Watanabe, 2017: ENSO complexity induced by state dependence of westerly wind events. Journal of Climate, 30 (9), 3401–3420.
  * Hu and Fedorov (2016) Hu, S., and A. V. Fedorov, 2016: Exceptionally strong easterly wind burst stalling El Niño of 2014. Proceedings of the National Academy of Sciences, 113 (8), 2005–2010.
  * Hu and Fedorov (2017) Hu, S., and A. V. Fedorov, 2017: The extreme El Niño of 2015–2016 and the end of global warming hiatus. Geophysical Research Letters, 44 (8), 3816–3824.
  * Iwakiri and Watanabe (2021) Iwakiri, T., and M. Watanabe, 2021: Mechanisms linking multi-year La Niña with preceding strong El Niño. Scientific reports, 11 (1), 17 465.
  * Jin (1997) Jin, F.-F., 1997: An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model. Journal of the Atmospheric Sciences, 54 (7), 811–829.
  * Jin and An (1999) Jin, F.-F., and S.-I. An, 1999: Thermocline and zonal advective feedbacks within the equatorial ocean recharge oscillator model for ENSO. Geophysical research letters, 26 (19), 2989–2992.
  * Jin et al. (2007) Jin, F.-F., L. Lin, A. Timmermann, and J. Zhao, 2007: Ensemble-mean dynamics of the ENSO recharge oscillator under state-dependent stochastic forcing. Geophysical Research Letters, 34 (3).
  * Kalnay et al. (1996) Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society, 77 (3), 437–472.
  * Kang et al. (2020) Kang, S. M., S.-P. Xie, Y. Shin, H. Kim, Y.-T. Hwang, M. F. Stuecker, B. Xiang, and M. Hawcroft, 2020: Walker circulation response to extratropical radiative forcing. Science Advances, 6 (47), eabd3021.
  * Kao and Yu (2009) Kao, H.-Y., and J.-Y. Yu, 2009: Contrasting eastern-Pacific and central-Pacific types of ENSO. Journal of Climate, 22 (3), 615–632.
  * Kleeman (2002) Kleeman, R., 2002: Measuring dynamical prediction utility using relative entropy. Journal of the Atmospheric Sciences, 59 (13), 2057–2072.
  * Kleeman (2011) Kleeman, R., 2011: Information theory and dynamical system predictability. Entropy, 13 (3), 612–649.
  * Klein et al. (1999) Klein, S. A., B. J. Soden, and N.-C. Lau, 1999: Remote sea surface temperature variations during ENSO: Evidence for a tropical atmospheric bridge. Journal of Climate, 12 (4), 917–932.
  * Kubo (1966) Kubo, R., 1966: The fluctuation-dissipation theorem. Reports on progress in physics, 29 (1), 255.
  * Kug et al. (2009) Kug, J.-S., F.-F. Jin, and S.-I. An, 2009: Two types of El Niño events: cold tongue El Niño and warm pool El Niño. Journal of Climate, 22 (6), 1499–1515.
  * Kullback (1959) Kullback, S., 1959: Statistics and information theory. J Wiley Sons, New York.
  * Kullback (1987) Kullback, S., 1987: Letter to the editor: The kullback-leibler distance. AMERICAN STATISTICIAN.
  * Kullback and Leibler (1951) Kullback, S., and R. A. Leibler, 1951: On information and sufficiency. The Annals of Mathematical Statistics, 22 (1), 79–86.
  * Kuo and Sloan (2005) Kuo, F. Y., and I. H. Sloan, 2005: Lifting the curse of dimensionality. Notices of the AMS, 52 (11), 1320–1328.
  * Levine et al. (2016) Levine, A., F. F. Jin, and M. J. McPhaden, 2016: Extreme noise–extreme El Niño: How state-dependent noise forcing creates El Niño–La Niña asymmetry. Journal of Climate, 29 (15), 5483–5499.
  * Levine and McPhaden (2016) Levine, A. F., and M. J. McPhaden, 2016: How the July 2014 easterly wind burst gave the 2015–2016 El Niño a head start. Geophysical research letters, 43 (12), 6503–6510.
  * Lopez and Kirtman (2014) Lopez, H., and B. P. Kirtman, 2014: WWBs, enso predictability, the spring barrier and extreme events. Journal of Geophysical Research: Atmospheres, 119 (17), 10–114.
  * Lorenz (1965) Lorenz, E. N., 1965: A study of the predictability of a 28-variable atmospheric model. Tellus, 17 (3), 321–333.
  * Ludescher et al. (2014) Ludescher, J., A. Gozolchiani, M. I. Bogachev, A. Bunde, S. Havlin, and H. J. Schellnhuber, 2014: Very early warning of next El Niño. Proceedings of the National Academy of Sciences, 111 (6), 2064–2066.
  * Luo et al. (2017) Luo, J.-J., G. Liu, H. Hendon, O. Alves, and T. Yamagata, 2017: Inter-basin sources for two-year predictability of the multi-year La Niña event in 2010–2012. Scientific reports, 7 (1), 2276.
  * Majda et al. (2005) Majda, A., R. V. Abramov, and M. J. Grote, 2005: Information theory and stochastics for multiscale nonlinear systems, Vol. 25. American Mathematical Soc.
  * Majda and Wang (2006) Majda, A., and X. Wang, 2006: Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows. Cambridge University Press, 10.1017/CBO9780511616778.
  * Majda and Wang (2010) Majda, A., and X. Wang, 2010: Linear response theory for statistical ensembles in complex systems with time-periodic forcing. Communications in Mathematical Sciences, 8 (1), 145 – 172.
  * Majda et al. (2009) Majda, A. J., R. Abramov, and B. Gershgorin, 2009: High skill in low-frequency climate response through fluctuation dissipation theorems despite structural instability. Proceedings of the National Academy of Sciences, 107 (2), 581–586, 10.1073/pnas.0912997107.
  * Majda and Branicki (2012) Majda, A. J., and M. Branicki, 2012: Lessons in uncertainty quantification for turbulent dynamical systems. Discrete and Continuous Dynamical Systems, 32 (9), 3133.
  * Majda and Chen (2018) Majda, A. J., and N. Chen, 2018: Model error, information barriers, state estimation and prediction in complex multiscale systems. Entropy, 20 (9), 644.
  * Majda and Gershgorin (2010) Majda, A. J., and B. Gershgorin, 2010: Quantifying uncertainty in climate change science through empirical information theory. Proceedings of the National Academy of Sciences, 107 (34), 14 958–14 963.
  * Majda and Gershgorin (2011) Majda, A. J., and B. Gershgorin, 2011: Link between statistical equilibrium fidelity and forecasting skill for complex systems with model error. Proceedings of the National Academy of Sciences, 108 (31), 12 599–12 604.
  * Majda and Qi (2018) Majda, A. J., and D. Qi, 2018: Strategies for reduced-order models for predicting the statistical responses and uncertainty quantification in complex turbulent dynamical systems. SIAM Review, 60 (3), 491–549.
  * Majda and Qi (2019) Majda, A. J., and D. Qi, 2019: Linear and nonlinear statistical response theories with prototype applications to sensitivity analysis and statistical control of complex turbulent dynamical systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29 (10), 103 131, 10.1063/1.5118690, https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.5118690/14624899/103131“˙1“˙online.pdf.
  * Mayer et al. (2016) Mayer, M., J. T. Fasullo, K. E. Trenberth, and L. Haimberger, 2016: ENSO-driven energy budget perturbations in observations and CMIP models. Climate Dynamics, 47, 4009–4029.
  * McPhaden et al. (2006) McPhaden, M. J., S. E. Zebiak, and M. H. Glantz, 2006: ENSO as an integrating concept in earth science. Science, 314 (5806), 1740–1745.
  * Mu (2000) Mu, M., 2000: Nonlinear singular vectors and nonlinear singular values. Science in China Series D: Earth Sciences, 43, 375–385.
  * Mu and Duan (2003) Mu, M., and W. Duan, 2003: A new approach to studying ENSO predictability: Conditional nonlinear optimal perturbation. Chinese Science Bulletin, 48, 1045–1047.
  * Mu and Wang (2001) Mu, M., and J. Wang, 2001: Nonlinear fastest growing perturbation and the first kind of predictability. Science in China Series D: Earth Sciences, 44, 1128–1139.
  * Palmer et al. (2009) Palmer, T. N., R. Buizza, F. Doblas-Reyes, T. Jung, M. Leutbecher, G. J. Shutts, M. Steinheimer, and A. Weisheimer, 2009: Stochastic parametrization and model uncertainty. ECMWF Technical Memorandum, 598.
  * Philander (1983) Philander, S. G. H., 1983: El Niño Southern Oscillation phenomena. Nature, 302 (5906), 295–301.
  * Puy et al. (2016) Puy, M., J. Vialard, M. Lengaigne, and E. Guilyardi, 2016: Modulation of equatorial Pacific westerly/easterly wind events by the Madden–Julian oscillation and convectively-coupled Rossby waves. Climate Dynamics, 46 (7-8), 2155–2178.
  * Robert et al. (2010) Robert, C. P., G. Casella, and G. Casella, 2010: Introducing monte carlo methods with r, Vol. 18. Springer.
  * Ropelewski and Halpert (1987) Ropelewski, C. F., and M. S. Halpert, 1987: Global and regional scale precipitation patterns associated with the El Niño/Southern Oscillation. Monthly Weather Review, 115 (8), 1606–1626.
  * Samelson and Tziperman (2001) Samelson, R. M., and E. Tziperman, 2001: Instability of the chaotic ENSO: The growth-phase predictability barrier. Journal of the atmospheric sciences, 58 (23), 3613–3625.
  * Santoso et al. (2017) Santoso, A., M. J. Mcphaden, and W. Cai, 2017: The defining characteristics of ENSO extremes and the strong 2015/2016 el niño. Reviews of Geophysics, 55 (4), 1079–1129.
  * Sharmila et al. (2023) Sharmila, S., H. Hendon, O. Alves, A. Weisheimer, and M. Balmaseda, 2023: Contrasting El Niño–La Niña predictability and prediction skill in 2-year reforecasts of the twentieth century. Journal of Climate, 36 (5), 1269–1285.
  * Sun and Yu (2009) Sun, F., and J.-Y. Yu, 2009: A 10–15-yr modulation cycle of ENSO intensity. Journal of Climate, 22 (7), 1718–1735.
  * Tao et al. (2023) Tao, L., M. Mu, L. Wang, X. Fang, W. Duan, and R.-H. Zhang, 2023: Impacts of initial zonal current errors on the predictions of two types of El Niño events. Journal of Geophysical Research: Oceans, 128 (6), e2023JC019 833.
  * Thual et al. (2019) Thual, S., A. J. Majda, and N. Chen, 2019: 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.
  * Thual et al. (2016) Thual, S., A. J. Majda, N. Chen, and S. N. Stechmann, 2016: Simple stochastic model for El Niño with westerly wind bursts. Proceedings of the National Academy of Sciences, 113 (37), 10 245–10 250.
  * Timmermann et al. (2018) Timmermann, A., and Coauthors, 2018: El Niño-Southern Oscillation complexity. Nature, 559 (7715), 535–545.
  * Toniazzo et al. (2008) Toniazzo, T., M. Collins, and J. Brown, 2008: The variation of ENSO characteristics associated with atmospheric parameter perturbations in a coupled model. Climate dynamics, 30, 643–656.
  * Vecchi et al. (2006) Vecchi, G. A., A. Wittenberg, and A. Rosati, 2006: Reassessing the role of stochastic forcing in the 1997–1998 El Niño. Geophysical Research Letters, 33 (1).
  * Wang et al. (2019) Wang, B., X. Luo, Y.-M. Yang, W. Sun, M. A. Cane, W. Cai, S.-W. Yeh, and J. Liu, 2019: 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), 22 512–22 517.
  * Williams (2001) Williams, D., 2001: Weighing the odds: a course in probability and statistics. Cambridge University Press.
  * Yang et al. (2021) Yang, Q., A. J. Majda, and N. Chen, 2021: ENSO diversity in a tropical stochastic skeleton model for the MJO, El Niño, and dynamic Walker circulation. Journal of Climate, 1–56.
  * Yu and Fang (2018) Yu, J.-Y., and S.-W. Fang, 2018: The distinct contributions of the seasonal footprinting and charged-discharged mechanisms to enso complexity. Geophysical Research Letters, 45 (13), 6611–6618.
  * Yu and Kao (2007) Yu, J.-Y., and H.-Y. Kao, 2007: Decadal changes of ENSO persistence barrier in SST and ocean heat content indices: 1958–2001. Journal of Geophysical Research: Atmospheres, 112 (D13).
  * Yu and Kim (2013) Yu, J.-Y., and S. T. Kim, 2013: Identifying the types of major El Niño events since 1870. International Journal of Climatology, 33 (8), 2105–2112.
  * Zhao et al. (2021) Zhao, S., F.-F. Jin, X. Long, and M. A. Cane, 2021: On the breakdown of ENSO’s relationship with thermocline depth in the central-equatorial Pacific. Geophysical Research Letters, 48 (9), e2020GL092 335.
  * Zheng and Zhu (2010) Zheng, F., and J. Zhu, 2010: Spring predictability barrier of ENSO events from the perspective of an ensemble prediction system. Global and Planetary Change, 72 (3), 108–117.