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

Yinling Zhang, Nan Chen, Jérôme Vialard, Xianghui Fang

► To cite this version:

Yinling Zhang, Nan Chen, Jérôme Vialard, Xianghui Fang. A Physics-Informed Auto-Learning Framework for Developing Stochastic Conceptual Models for ENSO Diversity. Journal of Climate, 2024, 37, pp.6323-6347. 10.1175/JCLI-D-24-0092.1 . insu-04821260

HAL Id: insu-04821260

https://insu.hal.science/insu-04821260v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Physics-Informed Auto-Learning Framework for Developing Stochastic Conceptual Models for ENSO Diversity

YINLING ZHANG<sup>a</sup>,<sup>a</sup> NAN CHEN<sup>a</sup>,<sup>a</sup> JÉRÔME VIALARD,<sup>b</sup> AND XIANGHUI FANG<sup>c</sup>

(Manuscript received 8 February 2024, in final form 2 August 2024, accepted 5 September 2024)

**ABSTRACT:** Understanding El Niño–Southern Oscillation (ENSO) dynamics has tremendously improved over the past few decades. The ENSO diversity in spatial pattern, peak intensity, and temporal evolution is, however, still poorly represented in conceptual ENSO models. In this paper, a physics-informed auto-learning framework is applied to derive ENSO stochastic conceptual models with varying degrees of freedom. The framework is computationally efficient and easy to apply. Once the state vector of the target model is set, causal inference is exploited to build the right-hand side of the equations based on a mathematical function library. Fundamentally different from standard nonlinear regression, the auto-learning framework provides a parsimonious model by retaining only terms that improve the dynamical consistency with observations. It can also identify crucial latent variables and provide physical explanations. This methodology successfully reconstructs the equations of a realistic six-dimensional reference ENSO model based on the recharge oscillator theory from its data. A hierarchy of lower-dimensional models is derived, and their representation of ENSO (including its diversity) is systematically assessed. The minimum model that represents ENSO diversity is four-dimensional, with three interannual variables describing the western Pacific thermocline depth, the eastern and central Pacific sea surface temperatures (SSTs), and one intraseasonal variable for westerly wind events. Without the intraseasonal variable, the resulting three-dimensional model underestimates extreme events and is too regular. A limited number of weak nonlinearities in the model are essential in reproducing the observed extreme El Niño events and the observed nonlinear relationship between eastern and western Pacific SSTs.

**SIGNIFICANCE STATEMENT:** This study develops a physics-informed auto-learning approach to improve the modeling and understanding of El Niño–Southern Oscillation (ENSO), a major climate phenomenon influencing global weather and climate. The auto-learning framework explores the causality between key processes to systematically produce stochastic conceptual models with different climate factors that simulate the diversity of observed ENSO events. The key finding is that the minimal model for characterizing the ENSO diversity (with the least number of climate factors) is a four-variable model capturing thermocline depth, sea surface temperatures, and wind bursts that can reproduce intensity and spatial pattern variation. This advancement provides an interpretable tool to identify the minimum sufficient processes governing ENSO behavior for improved predictability.

**KEYWORDS:** El Niño; ENSO; Wind stress; Ocean models; Stochastic models; Machine learning

## 1. Introduction

El Niño–Southern Oscillation (ENSO) dominates the planet's interannual variability with significant impacts on global weather and climate through atmospheric teleconnections (Philander 1983; Ropelewski and Halpert 1987; Klein et al. 1999; McPhaden et al. 2006; Dai and Wigley 2000). In the traditional viewpoint, ENSO consists of alternating periods of anomalously warm El Niño conditions and cold La Niña conditions that peak in the equatorial eastern Pacific (EP). However, in recent decades, many El Niño events have been observed to peak in the central Pacific (CP) (Ashok et al. 2007; Kao and Yu 2009; Kim et al. 2012). The discovery of this new type of ENSO events leads to the El Niño diversity concept (Capotondi et al. 2015), which underlines that El Niño events vary continuously between the CP and EP ideal types (Larkin and Harrison 2005; Yu and Kao 2007; Ashok et al. 2007; Kao and Yu 2009; Kug

et al. 2009; Jin 2022). The shift of the warming center can cause significant differences in the air–sea coupling over the equatorial Pacific, changing how ENSO affects the global climate (Chen and Cane 2008; Jin et al. 2008; Barnston et al. 2012; Hu et al. 2012; Zheng et al. 2014; Fang et al. 2015; Sohn et al. 2016; Santoso et al. 2019; Stuecker et al. 2013). In addition to these two major categories, individual ENSO events further exhibit diverse characteristics in spatial pattern, peak intensity, and temporal evolution. This is known as the ENSO complexity (Timmermann et al. 2018).

Conceptual models of ENSO are useful due to their low dimensionality, which makes them physically explainable and mathematically tractable, and then further provide the testbeds for examining different physical hypotheses and interdependence between state variables (Jin 1997a,b). Several low-order conceptual models have been developed, including the recharge–discharge oscillator (Jin 1997a), the delayed oscillator (Suarez and Schopf 1988; Battisti and Hirst 1989), the western Pacific oscillator (Weisberg and Wang 1997), and the advective–reflective oscillator (Picaut et al. 1997). Later, a unified ENSO

<footnote><sup>a</sup> Department of Mathematics, University of Wisconsin–Madison, Madison, Wisconsin</footnote>

<footnote><sup>b</sup> LOCEAN-IPSL, CNRS/IRD/MNHN/Sorbonne Université, Paris, France</footnote>

<footnote><sup>c</sup> Department of Atmospheric and Oceanic Sciences and Institute of Atmospheric Sciences, Fudan University, Shanghai, China</footnote>

<footnote>Corresponding author: Nan Chen, chennan@math.wisc.edu</footnote>
oscillator motivated by the dynamics and thermodynamics of Zebiak and Cane's coupled ocean-atmosphere model has been built (Wang 2001). Some other conceptual models have also been constructed to emphasize the nonlinear dynamics of ENSO and explore the ENSO features in the decadal time scale (Timmermann et al. 2003; Roberts et al. 2016; Timmermann and Jin 2002). These models were mainly proposed based on physical arguments and highlighted one or two specific dynamical features of ENSO as building blocks. They primarily focused on EP El Niños and led to many successes in understanding and predicting their key characteristics.

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

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

To address these gaps, a physics-informed auto-learning framework is developed in this paper. This framework can automatically derive stochastic conceptual models with different complexities that focus on capturing ENSO diversity.

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

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

## 2. The physics-informed auto-learning framework

### a. Overview

The physics-informed auto-learning framework aims to extract the known physical knowledge from existing models with the help of additional data to develop new models automatically. It determines the model structure based on causal inference and provides a physics-informed parsimonious model. An overview of the framework is as follows:
1) Either an observational time series or an appropriate existing model is needed to provide a long time series for the auto-learning framework. In the latter case, such a model is predetermined by physical knowledge and is calibrated by observations. The model simulation is expected to meet as many desirable validation criteria as possible. In contrast to short observational data, the long model simulation provides robust learning outcomes.

2) The quantities of interests, namely, the state variables, are predetermined to derive a new model using the framework.

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

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

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

Note that the observational data for many natural phenomena are not long enough to provide robust results. Therefore, the framework allows the use of an appropriate predetermined model as a surrogate to generate a long time series as input. In such a situation, the observations can still be valuable as the outcome from the auto-learning framework can be validated by the dynamical features and statistics computed from the observations.

It is worthwhile highlighting several significant features of this framework. First, the framework differs fundamentally from nonlinear regression in determining model structure. Typical nonlinear regression retains many nonphysical terms with small coefficients to minimize error. Sparse regression can eliminate some terms (Santosa and Symes 1986; Tibshirani 1996) but lacks physical constraints. In contrast, the causal inference exploits the causal relationship, which reflects the underlying physics, to automatically determine the model structure that provides a physically explainable parsimonious model (AlMomani et al. 2020; Fish et al. 2021; Kim et al. 2017; AlMomani and Bollt 2020; Elinger 2021; Elinger and Rogers 2021; Chen and Zhang 2023). The model identified from the auto-learning framework can naturally include additional terms that are not present in the reference model. This allows the framework to start from a general class of models, which is not limited to the form of the reference model. The causal inference can determine if these terms contribute to the underlying dynamics.

Second, the framework utilizes relationships and insights from the established models as a starting point rather than fitting models directly to short observational data. Candidate functional forms for the model structure are constrained by incorporating physical dependencies learned from the existing models. Unlike the pure data-driven method, the framework

imposes known dynamical couplings from the known knowledge to inform model derivation. Physical insights from the established models guide the exploration of model structures and dependencies.

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

The details of steps 2–5 are provided below, where a schematic illustration is shown in Fig. 1. The established model used in step 1, also known as the reference model, will be described at the end of this section.

## b. Determining model structure using causal inference

### 1) DETERMINING THE STATE VARIABLES

The state variables of the model to be derived from the auto-learning framework are predetermined. These variables are denoted by an n-dimensional column vector $U = (u_1, ..., u_n)^T$, where each $u_i$ is one state variable. To characterize ENSO diversity, the average SST in the eastern Pacific $T_E$ and that in the central Pacific $T_C$ are the two variables appearing in all the models to be derived. Depending on the interests, additional variables are included in different models, such as the zonal current in the central Pacific to represent advection processes and the extra heat content variables to account for more ocean adjustment time scales.

### 2) DEVELOPING A FUNCTION LIBRARY

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

$$f = \{f_1, ..., f_{m-1}, f_m, f_{m+1}, ..., f_M\}. \quad (2.1)$$

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

### 3) DETERMINING MODEL STRUCTURE USING CAUSATION INFERENCE

Next, causal inference is utilized to evaluate the relevance of the candidate functions to the underlying dynamics and thus determine the model structure. To this end, the so-called causation entropy $C_{f_m \to \hat{u}_i || [f \setminus f_m]}$ is computed to detect if the
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic illustration of the auto-learning framework.</figcaption></figure>

candidate function $f_m$ contributes to the right-hand side of the equation for each $u_i$, namely, $du_i/dt := \dot{u}_i$. The causation entropy is given by (Elinger and Rogers 2021; AlMomani et al. 2020; Fish et al. 2021)

$$C_{f_m \to \dot{u}_i} |_{[f_m]} = H(\dot{u}_i | [f_m]) - H(\dot{u}_i | \mathbf{f}), \quad (2.2)$$

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

$$\mathbf{f}f_m = \{f_1, \dots, f_{m-1}, f_{m+1}, \dots, f_M\}. \quad (2.3)$$

The term $H(\cdot|\cdot)$ is the conditional entropy from information theory, which is related to Shannon's entropy $H(\cdot)$ and the joint entropy $H(\cdot, \cdot)$. Two multidimensional random variables $\mathbf{X}$ and $\mathbf{Y}$ (with the corresponding states being $\mathbf{x}$ and $\mathbf{y}$) are defined as (Cover 1999)

$$\begin{align*}
H(\mathbf{X}) &= -\int_{\mathbf{x}} p(\mathbf{x}) \ln[p(\mathbf{x})] d\mathbf{x}, \\
H(\mathbf{Y}|\mathbf{X}) &= -\int_{\mathbf{x}} \int_{\mathbf{y}} p(\mathbf{x}, \mathbf{y}) \ln[p(\mathbf{y}|\mathbf{x})] d\mathbf{y} d\mathbf{x}, \\
H(\mathbf{X}, \mathbf{Y}) &= -\int_{\mathbf{x}} \int_{\mathbf{y}} p(\mathbf{x}, \mathbf{y}) \ln[p(\mathbf{x}, \mathbf{y})] d\mathbf{y} d\mathbf{x}, \tag{2.4}
\end{align*}$$

where the term $\ln[p(\mathbf{x})]$ [and other log functions in (2.4)] represents the (lack of) information, originating from the coding theory (Shannon 1948). Intuitively, when binary digits are

utilized to convey a message, the number of digits needed to fully characterize $m$ different messages (events) is $\log_2(m)$. Mimicking this logic and generalizing the situation to any continuous function, if an event occurs with probability $p(\mathbf{x})$, then $\log[1/p(\mathbf{x})]$ bits of information are needed to fully characterize it. On the right-hand side of (2.2), the difference between the two conditional entropies indicates the information in $\dot{u}_i$ contributed by the specific function $f_m$ given the contributions from all the other functions. Thus, it tells if $f_m$ provides additional information to $\dot{u}_i$ conditioned on the other potential terms in the dynamics. It is worthwhile to highlight that the causation entropy in (2.2) is fundamentally different from directly computing the correlation between $\dot{u}_i$ and $f_m$, as the causation entropy also considers the influence of the other library functions. If both $\dot{u}_i$ and $f_m$ are caused by a common factor $f'_m$ in the function library, then $\dot{u}_i$ and $f_m$ can be highly correlated. Yet, in such a case, the causation entropy $C_{f_m \to \dot{u}_i} |_{[f_m]}$ will be zero as $f_m$ is not the causation of $\dot{u}_i$.

The causation entropy is computed from each of the candidate functions in $\mathbf{f}$ to each $\dot{u}_i$. Thus, there are in total $NM$ causation entropies, which can be written as an $N \times M$ matrix, called the causation entropy matrix. Note that the dimension of $\mathbf{X}$ in (2.4) is $M$ when it is applied to compute the difference between two terms on the right-hand side of the causation entropy in (2.2). This implies that the direct calculation of the entropies in (2.4) involves a high-dimensional numerical integration, which is a well-known computationally challenging issue (Bellman 1961). To circumvent the direct numerical integration, the entropy calculation approximates all the joint and marginal distributions as Gaussians. In such a way, the causation entropy can be approximated by
$$ 
\begin{aligned}
C_{X \to X|Y} &= H(X|Y) - H(X|Y, X) \\
&= H(X, Y) - H(Y) - H(X, Y, X) + H(Y, X) \\
&\approx \frac{1}{2}\ln[\det(\mathbf{R}_{XY})] - \frac{1}{2}\ln[\det(\mathbf{R}_Y)] \\
&\quad - \frac{1}{2}\ln[\det(\mathbf{R}_{XY})] + \frac{1}{2}\ln[\det(\mathbf{R}_{YX})],
\end{aligned}
\quad (2.5) 
$$

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

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

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

### 4) PARAMETER ESTIMATION

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

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

where $\mathbf{\Phi}(\mathbf{U})$ is an $n$-dimensional column vector with each entry representing the summation of the remaining nonlinear

candidate functions from the causal inference for the corresponding component of $\mathbf{U}$. The matrix $\boldsymbol{\sigma} \in \mathbb{R}^{n \times d}$ is the noise amplitude, and $\mathbf{\dot{W}}(t) \in \mathbb{R}^{d \times 1}$ is a column vector with each entry being a white noise. Note that in the most general case, $n$ does not necessarily equal $d$ although in many systems, including those used in this work, $d = n$ is adopted. The matrix $\mathbf{\Theta} \in \mathbb{R}^s$ denotes a column vector containing all the parameters in the model with $s \gg n$ in a typical situation. The first term on the right-hand side of (2.6) can be written as

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

where $\mathbf{M} \in \mathbb{R}^{n \times s}$ is a matrix, while $\mathbf{Q} \in \mathbb{R}^n$ is a column vector. The entries of $\mathbf{M}$ and $\mathbf{Q}$ are functions of $\mathbf{U}$. Only $\mathbf{M}$ is multiplied by the parameter $\mathbf{\Theta}$. Note that $\mathbf{M}$ and $\mathbf{Q}$ can be functions of time, and therefore, the framework here covers the situation with seasonality and other functions with explicit time dependency. The parameter $\mathbf{\Theta}$ can be easily determined using a maximum likelihood estimator. Meanwhile, the noise coefficients $\boldsymbol{\sigma}$ are computed based on the residual via the least squares as part of the maximum likelihood estimator. Analytic solutions are available for estimating both $\mathbf{\Theta}$ and $\boldsymbol{\sigma}$. See Chen (2020) for the technical details. Notably, the entire parameter estimation can be solved via closed analytic formulas, making the procedure efficient and accurate.

### c. Auto-learning with latent variables

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

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

In the presence of latent variables, the auto-learning framework needs to be slightly modified by introducing a joint learning procedure to derive the equations for both groups of variables. Note that the time series is only available for the predetermined state variables (also called observed variables in the following context). Therefore, an iterative method is used to alternate between inferring the time series of the latent variables and deriving the entire model structure. The model structure is determined using the causal inference discussed in section 2b when the inferred time series of the latent variables is updated. The recovered time series of the latent variables, based on the current structure of the latent processes, is obtained by a Bayesian sampling technique (Chen 2020). The iterative procedure terminates when the model structure and parameters converge, which is ensured by the
Available Information

<figure><img src="image_3.png" /><figcaption>FIG. 2. Overview of the auto-learning framework with latent variables.</figcaption></figure>

<figure><img src="image_4.png" /></figure>

information revealed from observed variables. The initial value of the unobserved variable is assigned arbitrarily, for example, by taking the mean value of the process. It only takes a short period (typically a few months here) for the time series to be adjusted to eliminate the error in the initial guess. See Fig. 2 for the illustration of incorporating this iterative process into the auto-learning framework. The technical details of this iterative method can be found in Chen and Zhang (2023).

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

### 3. Reference model and observational data

#### a. Reference model

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

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

$$ \frac{du}{dt} = -ru - \frac{a_1 b_0 \mu}{2} (T_C + T_E) + \beta_u(I)\tau + \sigma_u \dot{W}_u, \quad (3.1a) $$

$$ \frac{dh_W}{dt} = -rh_W - \frac{a_2 b_0 \mu}{2} (T_C + T_E) + \beta_h(I)\tau + \sigma_h \dot{W}_h, \quad (3.1b) $$

$$ \frac{dT_C}{dt} = \left[ \frac{\gamma b_0 \mu}{2} - c_1(T_C, t) \right] T_C + \frac{\gamma b_0 \mu}{2} T_E + \gamma h_W + c_I I u + C_u + \beta_C(I) \tau + \sigma_C \dot{W}_C, \quad (3.1c) $$

$$ \frac{dT_E}{dt} = \gamma h_W + \left[ \frac{3\gamma b_0 \mu}{2} - c_2(t) \right] T_E - \frac{\gamma b_0 \mu}{2} T_C + \beta_E(I)\tau + \sigma_E \dot{W}_E, \quad (3.1d) $$

$$ \frac{d\tau}{dt} = -d_{\tau}\tau + \sigma_{\tau}(T_C, t)\dot{W}_{\tau}, \quad (3.1e) $$

$$ \frac{dI}{dt} = -\lambda(I - m) + \sigma_I(I)\dot{W}_I, \quad (3.1f) $$

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

In this model, stochasticity plays a crucial role in coupling variables at different time scales and parameterizing the unresolved features in the model. First, intraseasonal wind bursts $\tau$ are modeled by a stochastic differential equation with state-dependent $\sigma_{\tau}(T_C)$. The wind burst $\tau$ is then coupled to the
<table><caption>TABLE 1. Summary of the nondimensional units and reference model parameters. The square brackets around the h, u, τ, T, and t mean the unit of each variable. Note that the time unit is 2 months. Therefore, the annual cycle is represented by cos[(2π/6)t] and sin[(2π/6)t].</caption><tr><td>[h<sub>w</sub>]</td><td>150 m</td><td>[T<sub>E</sub>], [T<sub>C</sub>]</td><td>7.5°C</td></tr><tr><td>[u]</td><td>1.5 m s<sup>-1</sup></td><td>[t]</td><td>2 months</td></tr><tr><td>[τ]</td><td>5 m s<sup>-1</sup></td><td>d<sub>τ</sub></td><td>2 (1 month<sup>-1</sup>)</td></tr><tr><td>γ</td><td>0.45</td><td>r</td><td>0.15</td></tr><tr><td>α<sub>1</sub></td><td>0.0375</td><td>α<sub>2</sub></td><td>0.075</td></tr><tr><td>b<sub>0</sub></td><td>2.5</td><td>μ</td><td>0.5</td></tr><tr><td>σ</td><td>0.12I</td><td>λ</td><td>2/60 (5 yr<sup>-1</sup>)</td></tr><tr><td>p(I)</td><td>0.25 in I ∈ (0, 4)</td><td>c<sub>U</sub></td><td>0.018</td></tr><tr><td>Φ(x)</td><td>∫<sub>0</sub><sup>x</sup>(y - m)p(y)dy</td><td>σ<sub>I</sub>(I)</td><td>√[2/p(I)][-λΦ(I)]</td></tr><tr><td>β<sub>E</sub>(I)</td><td>0.1239 × (2 - 0.2I)</td><td>β<sub>u</sub>(I)</td><td>-0.2β<sub>E</sub>(I)</td></tr><tr><td>β<sub>h</sub>(I)</td><td>-0.4β<sub>E</sub>(I)</td><td>β<sub>C</sub>(I)</td><td>0.8β<sub>E</sub>(I)</td></tr><tr><td>σ<sub>u</sub></td><td>0.0310</td><td>σ<sub>h</sub></td><td>0.0155</td></tr><tr><td>σ<sub>C</sub></td><td>0.0310</td><td>σ<sub>E</sub></td><td>0.0232</td></tr><tr><td>σ<sub>τ</sub>(T<sub>C</sub>, t)</td><td>0.9[tanh(4.5T<sub>C</sub>) + 1] [1 + 0.25 cos(2π/6 t)]</td><td></td><td></td></tr><tr><td>c<sub>1</sub>(T<sub>C</sub>, t)</td><td>[15.6(T<sub>C</sub> + 0.75/7.5)<sup>2</sup> + 0.57][1 + 0.4 sin(2π/6 t)]</td><td></td><td></td></tr><tr><td>c<sub>2</sub>(t)</td><td>0.9[1 + 0.4 sin(2π/6 t + 2π/6) + 0.2 sin(2π/3 t + 2π/6)]</td><td></td><td></td></tr></table>

processes of the interannual part serving as an external forcing. In addition, four Gaussian random noises $\sigma_u \dot{W}_u$, $\sigma_h \dot{W}_h$, $\sigma_C \dot{W}_C$, and $\sigma_E \dot{W}_E$ are further added to the processes, which effectively parameterize the untracked factors, such as the subtropical atmospheric forcing at the Pacific Ocean. More broadly, the noise increases variability for statistics to match observations (Palmer et al. 2009). Second, since the details of the background Walker circulation consist of uncertainties and randomness (Chen et al. 2015), a simple but effective stochastic process is used to describe the temporal evolution of the decadal variability I (Yang et al. 2021). The multiplicative noise in the process of I is aimed at guaranteeing its positivity because the long-term average of the background Walker circulation is nonnegative. It is worth highlighting that the theoretical interpretation of ENSO is mainly categorized into two frameworks. One is a damped oscillatory mode interacting with multiplicative noise from the intraseasonal variabilities, and the other is a fully nonlinear oscillator. The model here belongs mainly to the former category. The model contains a minimum number of nonlinearities. The extreme El Niño events and non-Gaussian statistics of $T_E$ are mainly triggered by the multiplicative noise in the wind burst noise. The model can also be understood as having two fundamental oscillators: one for the CP events and the other for the EP events (Fang and Mu 2018). These two oscillators are linked by the decadal variability I that determines which one plays a dominant role alternatively. Finally, the effects of seasonality are added to both the wind activity $\sigma_\tau$ and the collective damping $c_1$ and $c_2$ to depict the seasonal phase-locking characteristics. The former accounts for the active phase of the MJO in boreal winter, while the latter is due to the seasonal migration of the intertropical convergence zone (ITCZ), which modulates the

strength of the upwelling and horizontal advection processes to influence the evolution of the SST. Seasonality manifests as the tendency of ENSO events to peak during boreal winter (Tziperman et al. 1997; Stein et al. 2014; Fang and Zheng 2021).

The model parameters are listed in Table 1.

## b. Observational datasets

The following datasets are utilized in this study to provide some estimates of the observed ENSO properties, including their diversity. For the assessment of the model's dynamical properties, monthly ocean temperatures, zonal current, thermocline depth, and daily zonal winds at 850 hPa are used. For this, the oceanic variables are derived from the GODAS dataset (Behringer and Xue 2004) (https://www.cpc.ncep.noaa.gov/products/GODAS), with thermocline depth calculated from potential temperature as the depth of the 20°C isotherm, while the atmospheric variable is obtained from the NCEP-NCAR reanalysis (Kalnay et al. 1996) (http://www.esrl.noaa.gov/psd/). The analysis period is from 1982 to 2019. This period is used because it is during the satellite era, and its observations are more reliable than those before. Anomalies are calculated by removing the monthly mean climatology of the whole period. For the assessment of the model performance on simulating different types of ENSO events, a longer period for SST (i.e., 1950–2020) is also used, which is obtained from the Extended Reconstructed SST, version 5 (Huang et al. 2017) (https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html). The Niño-4 ($T_C$ in the model) and Niño-3 ($T_E$ in the model) indices are the averages of SST anomalies over the regions 160°E–150°W, 5°S–5°N (CP) and 150°–90°W, 5°S–5°N (EP), respectively. The $h_W$ index in the model is the mean thermocline depth anomaly over the
<figure><img src="image_5.png" /><figcaption>FIG. 3. Summary of validation metrics for assessing the model performance.</figcaption></figure>

WP region (120°E–180°, 5°S–5°N). The *u* index in the model is the mean mixed layer zonal current in the CP region.

### c. Definitions of different types of the ENSO events

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

### 4. Validation metrics for assessing the model performance

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

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

To examine the performance of the models in reproducing fundamental ENSO properties, a time series of 20 000 years from each model is generated. Omitting the first 1000 years as possibly the burning period, the remaining time series of 19 000 years is divided into 271 nonoverlapping segments, each having a length of 70 years as the observational time series. These different segments allow us to estimate the 95% confidence interval when computing these statistics. Note that the burning period is when the model transits from a specific initial condition to statistical equilibrium, which is about a few decorrelation times of the system. A period of 1000 years is chosen as a simple cutoff, especially since the model can easily run for a few thousand years.

#### a. Metric for assessing the fundamental dynamical properties

##### NONLINEAR DEPENDENCE BETWEEN STATE VARIABLES

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

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

$$p[u(t)|v^{\text{obs}}(t)] \sim p[u(t)] \times p[v^{\text{obs}}(t)|u(t)], \quad (4.1)$$

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

## b. Metrics for fundamental properties of ENSO

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

### 1) SEASONAL SYNCHRONIZATION

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

### 2) PDFS

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

### 3) ACFs

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

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

## c. Metrics for ENSO diversity

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

### 1) BIVARIATE DIAGRAM OF LOCATION AND STRENGTH OF ENSO EVENTS

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

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

$$\text{SST}(x, t) = r_C(x)T_C(t) + r_E(x)T_E(t), \quad (4.2)$$

where $r_C$ and $r_E$ are the two regression coefficients between the observational time series at the longitude $x$ and the $T_C$ and $T_E$ time series, respectively.

### 2) OCCURRENCE FREQUENCY OF DIFFERENT TYPES OF EVENTS

The second criterion is the occurrence frequency of different ENSO events. This is a refined metric that counts the number of events in multiple subcategories, including all the El Niño events, the EP El Niños, the CP El Niños, the extreme El Niños, the multiyear El Niños, all the La Niña events, and the multiyear La Niñas (Chen et al. 2022) according to the definitions mentioned in section 4c.

## 5. Hierarchy of stochastic conceptual models for ENSO diversity

### a. Overview of the models

The auto-learning framework is applied to derive five models with different state variables and degrees of freedom, summarized in Table 2. The detailed expressions of these models
<table><caption>TABLE 2. Summary of the five models derived from the auto-learning framework.</caption><thead><tr><th>Model name</th><th>Acronym</th><th>Dimension</th><th>State variables</th></tr></thead><tbody><tr><td>Interannual intraseasonal decadal model with advection</td><td>I<sub>a</sub>I<sub>s</sub>DMA</td><td>6</td><td>u, h<sub>W</sub>, T<sub>C</sub>, T<sub>E</sub>, τ, I</td></tr><tr><td>Interannual intraseasonal model with advection</td><td>I<sub>a</sub>I<sub>s</sub>MA</td><td>5</td><td>u, h<sub>W</sub>, T<sub>C</sub>, T<sub>E</sub>, τ</td></tr><tr><td>Interannual intraseasonal decadal model</td><td>I<sub>a</sub>I<sub>s</sub>DM</td><td>5</td><td>h<sub>W</sub>, T<sub>C</sub>, T<sub>E</sub>, τ, I</td></tr><tr><td>Interannual intraseasonal model</td><td>I<sub>a</sub>I<sub>s</sub>M</td><td>4</td><td>h<sub>W</sub>, T<sub>C</sub>, T<sub>E</sub>, τ</td></tr><tr><td>Interannual model</td><td>I<sub>a</sub>M</td><td>3</td><td>h<sub>W</sub>, T<sub>C</sub>, T<sub>E</sub></td></tr></tbody></table>

are included in the appendix. Depending on the prescribed state variables in the model, the corresponding time series generated from the reference model is utilized for model identification in the auto-learning framework. The goal here is twofold. One is to discover a parsimonious model given the dimension and state variables. The other goal is to reveal the minimum model that captures the ENSO diversity based on the criteria proposed here.

### 1) I<sub>a</sub>I<sub>s</sub>DMA

The interannual intraseasonal decadal model with advection (I<sub>a</sub>I<sub>s</sub>DMA) is a six-dimensional model that involves three time scales, consistent with the multiscale system in Chen et al. (2022).

### 2) I<sub>a</sub>I<sub>s</sub>MA

The interannual intraseasonal model with advection (I<sub>a</sub>I<sub>s</sub>MA) is a five-dimensional model that excludes the decadal variable *I* and only contains interactions between the interannual and intraseasonal variables. This model will be used to examine if the decadal modulation significantly impacts the dynamics and statistics of ENSO diversity.

### 3) I<sub>a</sub>I<sub>s</sub>DM

The interannual intraseasonal decadal model (I<sub>a</sub>I<sub>s</sub>DM) is a five-dimensional model that involves three time scales but excludes the variable describing the ocean zonal current. This current variable was introduced to resolve advection, which plays a strong role in the central Pacific (Vialard et al. 2001). However, Jin and An (1999) pointed out that advection is implicitly included in the standard recharge oscillator framework, suggesting that the *u* variable may not be necessary.

### 4) I<sub>a</sub>I<sub>s</sub>M

The interannual intraseasonal model (I<sub>a</sub>I<sub>s</sub>M) is a four-dimensional model that excludes the decadal variability and the ocean zonal current. It aims to use a minimum number of state variables to characterize the interactions between interannual and intraseasonal variabilities.

### 5) I<sub>a</sub>M

The interannual model (I<sub>a</sub>M) is a three-dimensional model that further excludes the intraseasonal wind variable. Therefore, it reduces to a single time-scale model describing the interannual variabilities. Note that this model still includes stochastic forcing. However, the intraseasonal variability, which is characterized by a state-dependent red noise, is

omitted in this model. This model will be compared with the I<sub>a</sub>I<sub>s</sub>M to reveal the crucial role of intraseasonal variability in contributing to ENSO diversity (Chen et al. 2015; Eisenman et al. 2005; Jin et al. 2007). Despite the lack of skill to fully characterize ENSO diversity, this I<sub>a</sub>M will serve as the basis for the auto-learning framework to include latent variables for model improvement.

Note that the time scales mentioned in the above models are those explicitly appearing when designing the model. Due to the seasonal modulation of the collective damping terms and fast interactions between state variables and noise terms (Stuecker et al. 2013; Stein et al. 2014), the model simulation will exhibit features across different time scales and result in combination tones.

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

$$ 
\begin{aligned}
& u, h_W, T_C, T_E, \tau, I, u^2, h_W^2, T_C^2, T_E^2, I^2, uh_W, uT_C, \\
& uT_E, u\tau, uI, h_W T_C, h_W T_E, h_W \tau, h_W I, T_C T_E, \\
& T_C \tau, T_C I, T_E \tau, T_E I, \tau I, u^3, h_W^3, T_C^3, T_E^3, \tau^3, I^3, 1.
\end{aligned} 
\quad (5.1) 
$$

Despite the general form of the candidate functions, the choice of this function library contains some physical justification. Linear functions are natural choices that include damping mechanisms and linear oscillations. Quadratic functions are the mechanisms for nonlinear advections and oscillations. Cubic functions usually represent the nonlinear damping terms. In addition, since many existing ENSO conceptual models only use polynomials as the building block, such knowledge is utilized in developing the above function library. Note that for models with only a subset of the six possible variables, the corresponding candidate functions involving the missing state variables are removed from the library.

It is worth mentioning that the derivation of the six-dimensional I<sub>a</sub>I<sub>s</sub>DMA is also considered a proof-of-concept twin experiment to validate the auto-learning framework. As seen in the following (e.g., Fig. 9), the identified model has almost an identical structure to the reference model. The simulations from these two models also resemble each other. This evidence supports the auto-learning framework.

As a remark, since the auto-learning framework utilizes a relatively long time series, which contains a length of 2000 years, to discover the models, the threshold for selecting the terms on the right-hand side of the equations can be easily seen from the values in the causation entropy matrix. The model structure is robust with respect to the choice of the
<figure><img src="image_6.png" /><figcaption>FIG. 4. Comparison between different models using criterion 1a—assessing the fundamental dynamical properties via Bayesian inference. The correlations between the recovered mean time series and observations are used to quantify the performance of these models in the right box chart. (a) Correlations of the results for recovery of *u* and *h*<sub>W</sub> by observing *T*<sub>C</sub>, *T*<sub>E</sub>, and *I*. (b) Correlations of the results for recovery of *T*<sub>C</sub> and *T*<sub>E</sub> by observing other variables. On the left side, there is an example showing the recovered time series of the reference model in two experiments. The red lines represent the observation, the solid blue lines represent the ensemble mean simulation, and the blue shading areas are the confidence intervals of the ensemble simulations. The correlation *r* between the truth (red line) and the recovered ensemble mean time series (solid blue line) for each state variable is also listed.</figcaption></figure>

threshold. Therefore, the threshold, the only tuning parame-
ter in the method, will not be a crucial factor for the model
development.

b. *Intercomparison between different models*

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

1) CRITERION CATEGORY 1: ASSESSMENT OF THE
FUNDAMENTAL DYNAMICAL PROPERTIES

Figure 4 shows the results using the Bayesian inference to assess the fundamental dynamical properties. The correlation between the recovered mean time series and observations is used to quantify the performance of different models. Figures 4a and 4b illustrate the results of two experiments. In the first experiment, the SST variables *T*<sub>C</sub>, *T*<sub>E</sub>, and decadal variability *I* are observed to recover the time series of the remaining variables in each model. Note that *I* is the slowest variable and is always treated as observations. The results show that all models, except the three-dimensional interannual model, provide comparably high correlations in recovering *h*<sub>W</sub>, which are around 0.7. In contrast, the three-dimensional interannual model only has a correlation of around 0.5. One interesting finding in Fig. 4a is that recovery of *u* seems quite challenging because of the weak coupling relationship between the SST and current *u* in the models. This also means that *u* is less balanced with the

SST variables compared with thermocline fluctuation and wind
stress on the interannual time scale. Such a finding suggests
that the process of *u* may not be entirely necessary for building
the minimum conceptual model to capture the observed statis-
tical features of the ENSO diversity.

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

Two dashed boxes on the left side in Fig. 4 show the direct
results from Bayesian inference, including the ensemble mean
time series (blue), the confidence interval, and the observa-
tions (red) of the reference model as an example. The refer-
ence model has a good performance in characterizing the
nonlinear relationship between the state variables. Specifically,
<figure><img src="image_7.png" /><figcaption>FIG. 5. Comparison between different models using criteria 2a–2c including the seasonal variation of the variance, PDFs, and ACFs of $T_C$ and $T_E$.</figcaption></figure>

the recovered $T_E$ time series shows consistency with the observation in depicting the events around 1982–83 and 1997–98 and the 2015–16 super El Niños. The other models, except the three-dimensional interannual model, have comparable performance.

## 2) CRITERION CATEGORY 2: QUANTIFICATION OF FUNDAMENTAL ENSO PROPERTIES

The first row in Fig. 5 shows the seasonal dependence of the $T_C$ and $T_E$ amplitude (criterion 2a). Only the results of three models are shown to avoid redundancy—the first four models (from the six-dimensional interannual intraseasonal decadal model with advection to the four-dimensional interannual intraseasonal model) perform comparably similarly. All the models have the skill to reproduce seasonal synchronization.

Next, the comparison of the non-Gaussian PDFs (criterion 2b) is included in the second row of Fig. 5. All the models have similar profiles regarding these PDFs, especially the sign of the skewness for $T_C$ and $T_E$. The only noticeable difference among different models is that the PDF of $T_C$ associated with the three-dimensional interannual model is slightly biased compared with other models. Yet, it is worth remarking that the higher-order moments that characterize non-Gaussian features may not be easily seen from the standard PDF plots. Therefore, Fig. 6 shows the boxplots for the variance, skewness, and kurtosis from the model based on the 271 simulations. Despite the signs of the skewness being correctly recovered by all these models, the three-dimensional interannual model has significant biases in reproducing the skewness and kurtosis values. For $T_C$, both the skewness and kurtosis are stronger than the observations. In particular, the kurtosis exceeds three, which means the probability tail is fatter than a Gaussian distribution that will cause more occurrences of extreme CP La Niña events in such a model. In contrast, both the skewness and kurtosis are weaker than observations for

$T_E$. This corresponds to the underestimation of the extreme EP El Niño events. These findings indicate that the state-dependent red noise is crucial in reproducing these strong non-Gaussian features. The simple additive white noise is insufficient to capture the observed features in nature.

Finally, the third row in Fig. 5 compares different models in recovering the ACFs (criterion 2c). Similar to the variances and PDFs, only the three-dimensional interannual model demonstrates different features compared with other models. The ACFs of the three-dimensional interannual model have stronger oscillation patterns for both $T_C$ and $T_E$. Such a feature is related to the strong peak in the power spectrums. Due to the simplest structure among all the models, the three-dimensional model is dominated by such an oscillation pattern. In contrast, other models have more intrinsic time scales and mechanisms that help trigger more irregular behavior.

## 3) CRITERION CATEGORY 3: EVALUATION OF ENSO DIVERSITY

Figure 7 shows the bivariate diagram of the location and strength of El Niño events (criterion 3a). Again, the first four derived models have similar performance. Thus, only the four-dimensional interannual intraseasonal model and three-dimensional interannual model are shown. The events occur at nearly all longitudes across the equatorial Pacific. In general, the CP events always exhibit weak-to-moderate amplitudes, while amounts of the strong events can be found in the EP region, which is consistent with observations (Capotondi et al. 2015). Among different models, the three-dimensional interannual model is the only one that underestimates the number of extreme events in EP. This is related to the inaccuracy in the skewness and kurtosis shown in Fig. 6.

Figure 8 presents the occurrence frequency of different types of events (criterion 3b). All the models, except the three-dimensional interannual model, can capture the observed occurrence frequency
<figure><img src="image_8.png" /><figcaption>FIG. 6. Additional supporting information of criterion 2b. The second, third, and fourth moments, i.e., the variance, skewness, and kurtosis, respectively, of $T_C$ and $T_E$ for each model. These statistics can be computed using the PDFs in Fig. 5. In each boxplot, the central mark on each box denotes the median, while the bottom and top edges of the box signify the 5th and 95th percentiles, respectively. The horizontal dashed line in each panel shows the true value using observations.</figcaption></figure>

of all the seven types of ENSO events. This guarantees the skill of these models in simulating ENSO diversity and its complex features. Note that the three-dimensional interannual model nevertheless reproduces the occurrence frequency of most ENSO events. The only two exceptions are the underestimation of extreme events, which has already been discussed in Figs. 6 and 7, and the overestimation of the number of multiyear El Niño events due to the longer memory shown in the third row of Fig. 5.

## 6. Mechanisms and the minimum model for ENSO diversity

Figure 9 illustrates the main structure of each model determined by the auto-learning framework. Note that the same

number of columns is utilized for all the models to keep a unified format in the figure. The columns with dark blue shadings indicate the functions not built in the library for the corresponding model. Additional additive noise is added to the equations of $u$, $h_W$, $T_E$, and $T_C$, while multiplicative noise is imposed on $\tau$ and $I$. The details of the models can be found in the appendix.

### a. Minimum model

The four-dimensional $I_s I_a M$ is the minimum model that succeeds in capturing the target features of ENSO diversity characterized by the validation criteria in Figs. 4–8. A further reduction of the model dimension to three by eliminating the

<figure><img src="image_9.png" /><figcaption>FIG. 7. Comparison between different models using criterion 3a: bivariate diagram of location and strength of El Niño events. Each panel shows the distribution of equatorial Pacific SSTA maxima for El Niño events of the standard simulations. For each of the El Niño events, the winter-mean SST anomalies are averaged over the equatorial zone, and then, the Pacific zonal maximum is located. For the presentation purpose, the blue dots show the El Niño events based on the time series of only the first 4000 years of simulation in each panel.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Comparison between different models using criterion 3b: occurrence frequency of different types of events. In the boxplot, the central mark on each box denotes the median, while the bottom and top edges of the box signify the 5th and 95th percentiles, respectively. The horizontal dashed line in each panel shows the true value using observations.</figcaption></figure>

intraseasonal variability will significantly decrease the model's capacity to reproduce many observed features.

The model, according to Fig. 9e or (EA4), also has a relatively simple structure. It is almost a linear model containing the interactions between the three interannual variabilities. The only nonlinear terms are the cubic damping in the process of $T_C$, which have been shown to be crucial in

reproducing the non-Gaussian features in $T_C$ (Chen et al. 2022). Similarly, the multiplicative noise in the wind burst process contributes to the extreme events and non-Gaussian features in $T_E$ (Thual et al. 2016; Levine and Jin 2017).

Figure 10 shows the time series and the spatiotemporal pattern of the four-dimensional $I_a I_a M$. The time series in Fig. 10a illustrates a clear relationship between the wind $\tau$ and the EP

<figure><img src="image_11.png" /><figcaption>FIG. 9. A summary of the model structures for (a) the reference model and (b)–(f) the five derived models from the auto-learning framework. The yellow grids indicate the terms that appear in the model, while the blue ones mean the functions do not exist on the right-hand side of the derived model. The columns with light purple shadings indicate the functions not built in the library for the corresponding model. The pink boxes mark the terms that are present in the derived models but not in the reference model.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Simulation of the four-dimensional $I_3I_8M$. (a) Time series of the model (blue), which is compared with observations (red) in a qualitative way. (b) Hovmöller diagram of the spatiotemporal pattern of the model simulation. Several different types of ENSO events are marked by dashed boxes, where the orange ones stand for El Niño events and the blue ones represent the La Niña events. (c) Kernel density estimate of joint PDF (blue shading) of $T_E$ and $h_W$ from the model. The contour plot overlaid on top of the shading depicts the joint PDF of $T_E$ and $h_W$ derived from observations.</figcaption></figure>

events in the model, for example, around $t = 230, 236$, and 254. The coupled relationship between the three interannual variabilities in the model is also quite similar to those in the observations. The spatiotemporal patterns of different types of ENSO events are presented in Fig. 10b, including EP El Niño ($t = 202$), CP El Niño ($t = 138$), La Niña ($t = 203$), extreme El Niño ($t = 147$), multiyear El Niño ($t = 206–207$), and multiyear La Niña ($t = 215–216$). In addition, Fig. 10c shows the joint PDF of $T_E$ and $h_W$ from the model compared with the one from the observations. The two PDFs resemble each other, showing the dependence between $T_E$ and $h_W$.

b. Improving the model performance by incorporating latent variables

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

To improve the performance of the three-dimensional interannual model, one latent variable is added to the system. The auto-learning framework then involves a joint learning

procedure to derive the equations for the three predetermined state variables ($h_W, T_C$, and $T_E$) and the latent variable from the synthetic data generated from the reference model. Notably, the time series is only available for the three predetermined state variables, as the meaning of the latent variable remains unclear at this stage. Therefore, the learning procedure here fundamentally differs from deriving the four-dimensional $I_3I_8M$, in which all the four state variables have observations. An iterative method is applied to alternate between inferring the time series of the latent variables and deriving the model structure (Chen and Zhang 2023). The output of this method is a new model involving existing state variables and the latent variable after iteration, as shown in Fig. 2. For simplicity, only additive noise is assumed in the latent process since there is no prior knowledge about the potential dependence of the noise coefficient on the state variable.

To figure out if the latent variable mimics the behavior of a certain unobserved state variable, a comparison between the latent variable and this unobserved state variable in a pathwise way is necessary. Here, synthetic data generated from the reference model and real observations are taken as truth in the following two experiments, respectively. Based on the newly derived four-dimensional model from the auto-learning framework and true trajectories of existing state variables, the time series of the latent variable could be sampled via the conditional sampling method, which is computed from the Bayesian formula in (4.1). Figure 11 shows the observed time series of the three interannual variables ($h_W, T_C$, and $T_E$) and one sampled time series of the latent variable using the ensemble Kalman smoother data assimilation approach (see also section 4).
<figure><img src="image_13.png" /><figcaption>FIG. 11. Comparison of observed time series for interannual variables ($h_W$, $T_C$, and $T_E$) and a sampled latent variable trajectory derived from the auto-learning framework. (a) Trajectories (pink) of four variables ($h_W$, $T_C$, $T_E$, and $\tau$) from the six-dimensional reference model and a sampled trajectory (blue) of the latent variable $\tau$ based on the new four-dimensional model from the learning framework when the time series of $h_W$, $T_C$, and $T_E$ from the reference model is given. (b) As in (a), but for the time series from the reference model that is replaced by the observational data (red).</figcaption></figure>

Figure 11a shows a sampled trajectory based on the time series of $h_W$, $T_C$, and $T_E$ generated from the reference model, while Fig. 11b utilizes the observational dataset as the input. A significant finding is that the sampled trajectory of the latent variable has a relatively strong correlation with the atmospheric wind at the interannual time scale in both cases. It also contains obvious intraseasonal variabilities. However, recovery of individual intraseasonal wind events is challenging due to its strong uncertainty and randomness. Therefore, it is not expected that the time series of the latent variable will match the observations on such a fast time scale. Overall, the auto-learning framework succeeds in identifying the dominant component missing in the three-dimensional interannual model. While concluding that the latent variable being the atmospheric wind based on the strong correlation might be too subjective, the finding here implies some potential connections between them. Hereafter, the latent variable is denoted by $\tau$. As a remark, the difference between the identified latent variable and the wind time series can be used to compute the statistics, which allows us to determine the multiplicative noise coefficient in the latent process for further improvement of the model. Note that the amplitude of the sampled trajectories differs from that of the actual observed wind burst $\tau$ (see the last row of Fig. 11). This is not surprising, as $\tau$ multiplied by the wind stress coefficient is the actual forcing to the SST processes; for example, the forcing is given by $L_{34}\tau$ in (A7c). Neither $L_{34}$ nor $\tau$ is uniquely determined, but the actual wind forcing $L_{34}\tau$ can be accurately learned. If the wind stress coefficient is known, then the latent variable is naturally the wind burst.

Figure 12 displays the structure of the resulting new four-dimensional model (Fig. 12b) and compares it with the six-dimensional reference model (Fig. 12a). The latent variable is denoted by $\tau$ in this figure since it has been identified as a surrogate for the atmospheric wind. The resulting new four-dimensional model provides a very similar structure to the reference model in terms of the dynamics ($h_W$, $T_C$, and $T_E$), which are dominated by linear terms but contain a cubic damping for the process of $T_C$. It is worth highlighting that $T_C\tau$ appears on the right-hand side of the processes for the three interannual variables. Since only additive noise is used in the latent process, this term $T_C\tau$ plays a similar role as the multiplicative noise in the original (3.1), which represents the red noise and is consistent with the original stochastic re-charge-discharge oscillator (Jin et al. 2007).

Figure 13 includes the dynamical and statistical features of this new four-dimensional model. These features are mostly missed by the three-dimensional interannual model (see Figs. 4–8) but are generally accurately recovered by the new four-dimensional model. Figure 13a shows the results of two experiments utilizing Bayesian inference. In the first experiment, the observed time series of $T_C$ and $T_E$ are used to estimate the state of $h_W$ and the latent variable (the surrogate for the wind burst $\tau$ in the observational data), where the amplitude of the recovered wind time series has been adjusted to match the observations. In the second experiment, the observed time series of $h_W$ and $\tau$ are used to recover the state of $T_C$ and $T_E$. Among the two experiments, the recovered time series of $h_W$ is more accurate than that of the three-dimensional model, indicating that the new four-dimensional model captures a better oscillation mechanism between SST and
<figure><img src="image_14.png" /><figcaption>FIG. 12. Comparison of the model structures of (a) the six-dimensional reference model and (b) the resulting new four-dimensional model containing one latent variable $\tau$. (c) Causation entropy of the new model. The light purple shadings show the terms not included in the predetermined function library.</figcaption></figure>

thermocline with the contribution of the latent variable. Next, Fig. 13b presents the bivariate plot of the location and strength of El Niño events. The new four-dimensional model captures the feature that the average strength of CP events is lower than that of EP events. Besides, the model can generate a sufficient number of extreme events. This is also supported by the boxplot of the occurrence frequency of different ENSO events in Fig. 13c. The occurrence frequencies of different ENSO events from observational data (red diamonds) are in the corresponding box. Similarly, the model performs well in reproducing long-term statistics, especially PDFs, which can be seen in Fig. 13d. Recall that a noticeable bias was found in the PDF of $T_c$ of the three-dimensional model in Fig. 5c. With the addition of the latent variable, the new four-dimensional model can accurately capture the non-Gaussian features. Further validation is shown in Fig. 13e, which contains boxplots of variance, skewness, and kurtosis. It illustrates similar profiles among these aspects with observations (black dashed line). In conclusion, the latent variable that serves as the wind process provides essential feedback to the interannual variables. This allows the new model to reproduce the observed features of ENSO diversity.

# 7. Further discussion and conclusions

## a. Summary

In this paper, a physics-informed auto-learning framework is built to derive a hierarchy of stochastic conceptual models with different state variables and degrees of freedom for ENSO diversity. These models are intercompared based on validation criteria. Except for the three-dimensional interannual model, all other models have comparably skillful results

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

## b. Multiscale nature in modeling ENSO diversity

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

On the other hand, the role of decadal variability in contributing to the ENSO diversity can be understood by comparing the four-dimensional $I_aI_sM$ with the five-dimensional $I_aI_sDM$ and the six-dimensional $I_aI_sDMA$. It turns out that
<figure><img src="image_15.png" /><figcaption>FIG. 13. Validation tests for the new four-dimensional model derived from the auto-learning framework. (a) Assessing the fundamental dynamical properties. Note that the wind time series has been plotted as a monthly average for illustration purposes. (b) Bivariate diagram of the locations and strengths of different ENSO events. (c) Occurrence frequency of different ENSO events. In the boxplot, the central mark on each box denotes the median, and the bottom and top edges present the 5th and 95th percentiles, respectively. (d) Long-time statistics including PDFs and ACFs of $T_C$ and $T_E$. (e) Additional supporting information of PDFs in (d). It presents the statistical features of SST, containing the variance, skewness, and kurtosis of $T_C$ and $T_E$.</figcaption></figure>

the results based on the validation metrics do not have a sig-
nificant change when the decadal variable is excluded from
the models. However, caution is needed to explain such a re-
sult. All the results are based on 70-yr segments of data as
same as the length of observations, which could be too short
for the random modulation of the decadal variability to play a
significant role in changing the model statistics. In other
words, its role in modifying the statistics may not be distin-
guishable from certain random noise. It is worth highlighting
that decadal variability is essential (Dieppois et al. 2021) in
situations with regime switching, such as the climate change
scenario (Yu and Kao 2007; Yu and Kim 2013). This was the
guideline for designing the multiscale model in Chen et al.
(2022), where the decadal variability determines the CP- or
EP-favored regimes. If the mean state of the decadal variabil-
ity is changed, then it will affect the occurrence frequency of
the CP and EP events (Chen et al. 2022), which cannot be eas-
ily reproduced by simple random noises. The model perfor-
mance for such a sensitivity experiment is not included in the
model validation criteria adopted here. These criteria are pri-
marily designed to examine the skill of the models in repro-
ducing the current climate.

### c. Applicability to CMIP and paleo databases

The auto-learning framework utilizes information theory to
identify model structures requiring long-time datasets. While
real-world observations are limited, model simulations and pa-
leo reconstructions provide ample data for training conceptual

models. Coupled Model Intercomparison Project (CMIP)
model outputs contain multicentury projections enabling robust
structure extraction (Eyring et al. 2016). Paleo-proxy records
also have long climate time series that could inform minimal
model development for different climate epochs (Brown et al.
2008). Therefore, applying the auto-learning approach with
CMIP and paleo databases is an intriguing future direction. The
framework may extract underlying dynamics from complex
model projections. Paleo-informed models could also give in-
sights into ENSO evolution under past climate regimes (Ford
et al. 2012). The feasibility of extracting structures from such al-
ternative datasets merits in-depth exploration.

### d. Relation to recharge oscillator models

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

### e. Potential future directions

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

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

*Data availability statement.* The datasets analyzed during this study are derived from publicly available sources. Oceanic variables, including monthly ocean temperatures, zonal currents, and thermocline depth, were obtained from the Global Ocean Data Assimilation System (GODAS) dataset, available at https://www.cpc.ncep.noaa.gov/products/GODAS, with thermocline depth calculated based on the depth of the 20°C isotherm following Behringer and Xue (2004). Atmospheric variables, specifically daily zonal winds at 850 hPa, were sourced from the NCEP-NCAR reanalysis project, accessible at http://www.esrl.noaa.gov/psd/, as described by Kalnay et al. (1996). Sea surface temperature (SST) data for the period 1950–2020, used in simulating different types of

ENSO events, were derived from the Extended Reconstructed SST, version 5 (ERSST.v5), available at https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html, as documented by Huang et al. (2017). For further inquiries or access to the processed data and modeling scripts used in this study, including the calculation of anomalies and indices, please contact the corresponding author.

## APPENDIX

### Details of the Stochastic Conceptual Models and Parameters Resulting from the Auto-Learning Framework

#### a. Interannual intraseasonal decadal model with advection

The interannual intraseasonal decadal model with advection ($I_aI_sDMA$) is as follows:

$$ \frac{du}{dt} = L_{11}u + L_{12}\tau + L_{13} + \sigma_u \dot{W}_u, \quad (A1a) $$

$$ \frac{dh_W}{dt} = L_{21}h_W + L_{22}T_E + L_{23}\tau + L_{24}\tau I + L_{25} + \sigma_h \dot{W}_h, \quad (A1b) $$

$$ \frac{dT_C}{dt} = L_{31}h_W + L_{32}(t)T_C + L_{33}T_E + L_{34}\tau + L_{35}(t)T_C^2 + L_{36}uI + L_{37}\tau I + L_{38}(t)T_C^3 + L_{39} + \sigma_C \dot{W}_C, \quad (A1c) $$

$$ \frac{dT_E}{dt} = L_{41}h_W + L_{42}T_C + L_{43}(t)T_E + L_{44}\tau + L_{45}\tau I + L_{46} + \sigma_E \dot{W}_E, \quad (A1d) $$

$$ \frac{d\tau}{dt} = L_{51}\tau + L_{52} + \sigma_{\tau}(T_C)\dot{W}_{\tau}, \quad (A1e) $$

$$ \frac{dI}{dt} = L_{61}I + L_{62} + \sigma_I(I)\dot{W}_I, \quad (A1f) $$

where the parameters are listed in Table A1. Note that the functions of multiplicative noises are not included in the

<table><caption>TABLE A1. Parameters of the $I_aI_sDMA$.</caption><tr><td>L<sub>11</sub></td><td>-0.1400</td><td>L<sub>31</sub></td><td>0.4540</td><td>L<sub>41</sub></td><td>0.4493</td></tr><tr><td>L<sub>12</sub></td><td>-0.0428</td><td>L<sub>32</sub></td><td>-0.4505 - 0.2989s<sub>1</sub></td><td>L<sub>42</sub></td><td>-0.2830</td></tr><tr><td>L<sub>13</sub></td><td>0.0000</td><td>L<sub>33</sub></td><td>0.2850</td><td>L<sub>43</sub></td><td>-0.0558 - 0.3618s<sub>2</sub> - 0.1788s<sub>3</sub></td></tr><tr><td>L<sub>21</sub></td><td>-0.1663</td><td>L<sub>34</sub></td><td>0.1983</td><td>L<sub>44</sub></td><td>0.2470</td></tr><tr><td>L<sub>22</sub></td><td>-0.0694</td><td>L<sub>35</sub></td><td>-3.1142 - 1.2090s<sub>1</sub></td><td>L<sub>45</sub></td><td>-0.0245</td></tr><tr><td>L<sub>23</sub></td><td>-0.1007</td><td>L<sub>36</sub></td><td>0.1218</td><td>L<sub>46</sub></td><td>0.0001</td></tr><tr><td>L<sub>24</sub></td><td>0.0100</td><td>L<sub>37</sub></td><td>-0.0196</td><td>L<sub>51</sub></td><td>-1.9942</td></tr><tr><td>L<sub>25</sub></td><td>0.0001</td><td>L<sub>38</sub></td><td>-15.6559 - 6.2024s<sub>1</sub></td><td>L<sub>52</sub></td><td>0.0045</td></tr><tr><td>σ<sub>u</sub></td><td>0.0310</td><td>L<sub>39</sub></td><td>0.0177</td><td>L<sub>61</sub></td><td>-0.0323</td></tr><tr><td>σ<sub>h</sub></td><td>0.0155</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td><td>L<sub>62</sub></td><td>0.0639</td></tr><tr><td>σ<sub>C</sub></td><td>0.0310</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2πt/6)</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2πt/6)</td></tr><tr><td>σ<sub>E</sub></td><td>0.0232</td><td>σ<sub>τ</sub></td><td>0.8999[tanh(4.5T<sub>C</sub>) + 1][1 + 0.25 cos(2πt/6)]</td><td></td><td></td></tr><tr><td>Φ(x)</td><td>∫<sub>b</sub><sup>x</sup>(y - m)p(y)dy</td><td>σ<sub>I</sub>(I)</td><td>√(2/p(I))[-λΦ(I)]</td><td></td><td></td></tr></table>
<table><caption>TABLE A2. Parameters of the I<sub>a</sub>I<sub>s</sub>MA.</caption><tbody><tr><td>L<sub>11</sub></td><td>-0.1400</td><td>L<sub>31</sub></td><td>0.2183</td><td>L<sub>41</sub></td><td>0.4443</td></tr><tr><td>L<sub>12</sub></td><td>-0.0428</td><td>L<sub>32</sub></td><td>0.4222</td><td>L<sub>42</sub></td><td>-0.2741</td></tr><tr><td>L<sub>13</sub></td><td>0.0000</td><td>L<sub>33</sub></td><td>-0.4061 - 0.2989s<sub>1</sub></td><td>L<sub>43</sub></td><td>-0.0553 - 0.3597s<sub>2</sub> - 0.1791s<sub>3</sub></td></tr><tr><td>L<sub>21</sub></td><td>-0.1481</td><td>L<sub>34</sub></td><td>0.2455</td><td>L<sub>44</sub></td><td>0.1976</td></tr><tr><td>L<sub>22</sub></td><td>-0.0501</td><td>L<sub>35</sub></td><td>0.1568</td><td>L<sub>45</sub></td><td>0.0001</td></tr><tr><td>L<sub>23</sub></td><td>-0.0456</td><td>L<sub>36</sub></td><td>-2.8721 - 1.1140s<sub>1</sub></td><td>L<sub>51</sub></td><td>-1.9942</td></tr><tr><td>L<sub>24</sub></td><td>-0.0793</td><td>L<sub>37</sub></td><td>-14.5992 - 5.6033s<sub>1</sub></td><td>L<sub>52</sub></td><td>0.0045</td></tr><tr><td>L<sub>25</sub></td><td>0.0000</td><td>L<sub>38</sub></td><td>0.0166</td><td></td><td></td></tr><tr><td>σ<sub>u</sub></td><td>0.0310</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td><td></td><td></td></tr><tr><td>σ<sub>h</sub></td><td>0.0155</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td><td></td><td></td></tr><tr><td>σ<sub>C</sub></td><td>0.0310</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td><td></td><td></td></tr><tr><td>σ<sub>E</sub></td><td>0.0232</td><td>σ<sub>τ</sub></td><td>0.8999[tanh(4.5T<sub>C</sub>) + 1][1 + 0.25 cos(2πt/6)]</td><td></td><td></td></tr></tbody></table>

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

### b. Interannual intraseasonal model with advection

The interannual intraseasonal model with advection (I<sub>a</sub>I<sub>s</sub>MA) is as follows:

$$ \frac{du}{dt} = L_{11}u + L_{12}\tau + L_{13} + \sigma_u \dot{W}_u, \quad (A2a) $$

$$ \frac{dh_W}{dt} = L_{21}h_W + L_{22}T_C + L_{23}T_E + L_{24}\tau + L_{25}\tau + L_{25} + \sigma_h \dot{W}_h, \quad (A2b) $$

$$ \frac{dT_C}{dt} = L_{31}u + L_{32}h_W + L_{33}(t)T_C + L_{34}T_E + L_{35}\tau + L_{36}(t)T_C^2 + L_{37}(t)T_C^3 + L_{38} + \sigma_C \dot{W}_C, \quad (A2c) $$

$$ \frac{dT_E}{dt} = L_{41}h_W + L_{42}T_C + L_{43}(t)T_E + L_{44}\tau + L_{45} + \sigma_E \dot{W}_E, \quad (A2d) $$

$$ \frac{d\tau}{dt} = L_{51}\tau + L_{52} + \sigma_{\tau}(T_C)\dot{W}_{\tau}, \quad (A2e) $$

where the parameters are listed in Table A2.

### c. Interannual intraseasonal decadal model

The interannual intraseasonal decadal model (I<sub>a</sub>I<sub>s</sub>DM) is as follows:

$$ \frac{dh_W}{dt} = L_{11}h_W + L_{12}T_E + L_{13}\tau + L_{14}\tau I + L_{15} + \sigma_h \dot{W}_h, \quad (A3a) $$

$$ \frac{dT_C}{dt} = L_{21}h_W + L_{22}(t)T_C + L_{23}T_E + L_{24}\tau + L_{25}(t)T_C^2 + L_{26}T_C I + L_{27}T_E I + L_{28}\tau I + L_{29}(t)T_C^3 + L_{2,10} + \sigma_C \dot{W}_C, \quad (A3b) $$

$$ \frac{dT_E}{dt} = L_{31}h_W + L_{32}T_C + L_{33}(t)T_E + L_{34}\tau + L_{35}\tau I + L_{36} + \sigma_E \dot{W}_E, \quad (A3c) $$

$$ \frac{d\tau}{dt} = L_{41}\tau + L_{42} + \sigma_{\tau}(T_C)\dot{W}_{\tau}, \quad (A3d) $$

$$ \frac{dI}{dt} = L_{51}I + L_{52} + \sigma_I(I)\dot{W}_I, \quad (A3e) $$

where the parameters are listed in Table A3.

### d. Interannual intraseasonal model

The interannual intraseasonal model (I<sub>a</sub>I<sub>s</sub>M) is given as follows:

$$ \frac{dh_W}{dt} = L_{11}h_W + L_{12}T_C + L_{13}T_E + L_{14}\tau + L_{15} + \sigma_h \dot{W}_h, \quad (A4a) $$

<table><caption>TABLE A3. Parameters of the I<sub>a</sub>I<sub>s</sub>DM.</caption><tbody><tr><td>L<sub>11</sub></td><td>-0.1663</td><td>L<sub>27</sub></td><td>-0.0563</td><td>L<sub>41</sub></td><td>-1.9942</td></tr><tr><td>L<sub>12</sub></td><td>-0.0694</td><td>L<sub>28</sub></td><td>-0.0203</td><td>L<sub>42</sub></td><td>0.0045</td></tr><tr><td>L<sub>13</sub></td><td>-0.1007</td><td>L<sub>29</sub></td><td>-15.3784 - 6.2253s<sub>1</sub></td><td>L<sub>51</sub></td><td>-0.0323</td></tr><tr><td>L<sub>14</sub></td><td>0.0100</td><td>L<sub>2,10</sub></td><td>0.0166</td><td>L<sub>52</sub></td><td>0.0639</td></tr><tr><td>L<sub>15</sub></td><td>0.0001</td><td>L<sub>31</sub></td><td>0.4493</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td></tr><tr><td>L<sub>21</sub></td><td>0.5030</td><td>L<sub>32</sub></td><td>-0.2830</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td></tr><tr><td>L<sub>22</sub></td><td>-0.4785 - 0.2993s<sub>1</sub></td><td>L<sub>33</sub></td><td>-0.0558 - 0.3518s<sub>2</sub> - 0.1788s<sub>3</sub></td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td></tr><tr><td>L<sub>23</sub></td><td>0.3142</td><td>L<sub>34</sub></td><td>0.2470</td><td>σ<sub>h</sub></td><td>0.0155</td></tr><tr><td>L<sub>24</sub></td><td>0.1981</td><td>L<sub>35</sub></td><td>-0.0245</td><td>σ<sub>C</sub></td><td>0.0310</td></tr><tr><td>L<sub>25</sub></td><td>-2.9602 - 1.1545s<sub>1</sub></td><td>L<sub>36</sub></td><td>0.0001</td><td>σ<sub>E</sub></td><td>0.0232</td></tr><tr><td>L<sub>26</sub></td><td>0.0555</td><td>σ<sub>τ</sub></td><td>0.8999[tanh(4.5T<sub>C</sub>) + 1][1 + 0.25 cos(2πt/6)]</td><td></td><td></td></tr></tbody></table>
<table><caption>TABLE A4. Parameters of the I<sub>a</sub>I<sub>s</sub>M.</caption><tr><td>L<sub>11</sub></td><td>-0.1481</td><td>L<sub>25</sub></td><td>-2.7519 - 1.0865s<sub>1</sub></td><td>L<sub>41</sub></td><td>-1.9942</td></tr><tr><td>L<sub>12</sub></td><td>-0.0501</td><td>L<sub>26</sub></td><td>-14.2956 - 5.7679s<sub>1</sub></td><td>L<sub>42</sub></td><td>0.0045</td></tr><tr><td>L<sub>13</sub></td><td>-0.0456</td><td>L<sub>27</sub></td><td>0.0158</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td></tr><tr><td>L<sub>14</sub></td><td>-0.0793</td><td>L<sub>31</sub></td><td>0.4446</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td></tr><tr><td>L<sub>15</sub></td><td>0.0000</td><td>L<sub>32</sub></td><td>-0.2741</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td></tr><tr><td>L<sub>21</sub></td><td>0.4861</td><td>L<sub>33</sub></td><td>-0.0553 - 0.3597s<sub>2</sub> - 0.1791s<sub>3</sub></td><td>σ<sub>h</sub></td><td>0.0155</td></tr><tr><td>L<sub>22</sub></td><td>-0.3443 - 0.3014s<sub>1</sub></td><td>L<sub>34</sub></td><td>0.1976</td><td>σ<sub>C</sub></td><td>0.0310</td></tr><tr><td>L<sub>23</sub></td><td>0.1813</td><td>L<sub>35</sub></td><td>0.0000</td><td>σ<sub>E</sub></td><td>0.0232</td></tr><tr><td>L<sub>24</sub></td><td>0.1558</td><td>σ<sub>τ</sub></td><td>0.8999[tanh(4.5T<sub>C</sub>) + 1][1 + 0.25 cos(2πt/6)]</td><td></td><td></td></tr></table>

$$ \frac{dT_C}{dt} = L_{21}h_W + L_{22}(t)T_C + L_{23}T_E + L_{24}\tau \\ + L_{25}(t)T_C^2 + L_{26}(t)T_C^3 + L_{27} + \sigma_C \dot{W}_C, \quad (A4b) $$

$$ \frac{dT_E}{dt} = L_{31}h_W + L_{32}T_C + L_{33}(t)T_E \\ + L_{34}\tau + L_{35} + \sigma_E \dot{W}_E, \quad (A4c) $$

$$ \frac{d\tau}{dt} = L_{41}\tau + L_{42} + \sigma_{\tau}(T_C)\dot{W}_{\tau}, \quad (A4d) $$

where the parameters are listed in Table A4.

### e. Interannual model

The interannual model (I<sub>a</sub>M) is given as follows:

$$ \frac{dh_W}{dt} = L_{11}h_W + L_{12}T_C + L_{13}T_E + L_{14}T_C^2 \\ + L_{15} + \sigma_h \dot{W}_h, \quad (A5a) $$

$$ \frac{dT_C}{dt} = L_{21}h_W + L_{22}(t)T_C + L_{23}T_E + L_{24}h_W T_C \\ + L_{25}T_C T_E + L_{26} + \sigma_C \dot{W}_C, \quad (A5b) $$

$$ \frac{dT_E}{dt} = L_{31}h_W + L_{32}(t)T_E + L_{33}T_C^2 + L_{34}(t)T_E^2 \\ + L_{35}(t)T_E^3 + L_{36} + \sigma_E \dot{W}_E, \quad (A5c) $$

where the parameters are listed in Table A5.

### f. Linear model

The linear six-dimensional model is as follows:

$$ \frac{du}{dt} = L_{11}u + L_{12}\tau + L_{13} + \sigma_u \dot{W}_u, \quad (A6a) $$

$$ \frac{dh_W}{dt} = L_{21}h_W + L_{22}T_E + L_{23}\tau \\ + L_{24}\tau I + L_{25} + \sigma_h \dot{W}_h, \quad (A6b) $$

$$ \frac{dT_C}{dt} = L_{31}h_W + L_{32}(t)T_C + L_{33}T_E + L_{34}\tau \\ + L_{35}uI + L_{36}\tau I + L_{37} + \sigma_C \dot{W}_C, \quad (A6c) $$

$$ \frac{dT_E}{dt} = L_{41}h_W + L_{42}T_C + L_{43}(t)T_E \\ + L_{44}\tau + L_{45}\tau I + L_{46} + \sigma_E \dot{W}_E, \quad (A6d) $$

$$ \frac{d\tau}{dt} = L_{51}\tau + L_{52} + \sigma_{\tau}(T_C)\dot{W}_{\tau}, \quad (A6e) $$

$$ \frac{dI}{dt} = L_{61}I + L_{62} + \sigma_I(I)\dot{W}_I, \quad (A6f) $$

where the parameters are listed in Table A6.

### g. New four-dimensional model with the latent variable

The new four-dimensional model with the latent variable is given as follows:

$$ \frac{dh_W}{dt} = L_{11}h_W + L_{12}T_C + L_{13}\hat{\tau} \\ + L_{14}T_C \hat{\tau} + L_{15} + \sigma_h \dot{W}_h, \quad (A7a) $$

$$ \frac{dT_C}{dt} = L_{21}h_W + L_{22}(t)T_C + L_{23}\hat{\tau} + L_{24}(t)T_C^2 \\ + L_{25}T_C \hat{\tau} + L_{26}(t)T_C^3 + L_{27} + \sigma_C \dot{W}_C, \quad (A7b) $$

$$ \frac{dT_E}{dt} = L_{31}h_W + L_{32}T_C + L_{33}(t)T_E \\ + L_{34}\hat{\tau} + L_{35}T_C \hat{\tau} + L_{36} + \sigma_E \dot{W}_E, \quad (A7c) $$

$$ \frac{d\hat{\tau}}{dt} = L_{41}\hat{\tau} + L_{42} + \sigma_{\tau} \dot{W}_{\tau}, \quad (A7d) $$

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

<table><caption>TABLE A5. Parameters of the I<sub>a</sub>M.</caption><tr><td>L<sub>11</sub></td><td>-0.0678</td><td>L<sub>24</sub></td><td>0.6729</td><td>L<sub>36</sub></td><td>-0.0078</td></tr><tr><td>L<sub>12</sub></td><td>-0.1927</td><td>L<sub>25</sub></td><td>-1.1864</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td></tr><tr><td>L<sub>13</sub></td><td>-0.0593</td><td>L<sub>26</sub></td><td>0.0070</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td></tr><tr><td>L<sub>14</sub></td><td>-0.6729</td><td>L<sub>31</sub></td><td>0.2719</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td></tr><tr><td>L<sub>15</sub></td><td>0.0049</td><td>L<sub>32</sub></td><td>0.0523 - 0.2549s<sub>2</sub> - 0.1792s<sub>3</sub></td><td>σ<sub>h</sub></td><td>0.0156</td></tr><tr><td>L<sub>21</sub></td><td>0.3146</td><td>L<sub>33</sub></td><td>1.1864</td><td>σ<sub>C</sub></td><td>0.0311</td></tr><tr><td>L<sub>22</sub></td><td>-0.2553 - 0.1799s<sub>1</sub></td><td>L<sub>34</sub></td><td>-0.0134 + 0.0284s<sub>2</sub> + 0.0381s<sub>3</sub></td><td>σ<sub>E</sub></td><td>0.0236</td></tr><tr><td>L<sub>23</sub></td><td>0.2340</td><td>L<sub>35</sub></td><td>-0.2611 - 0.3158s<sub>2</sub> + 0.1207s<sub>3</sub></td><td></td><td></td></tr></table>
<table><caption>TABLE A6. Parameters of the linear six-dimensional model.</caption><tr><td>L<sub>11</sub></td><td>-0.1400</td><td>L<sub>31</sub></td><td>0.4461</td><td>L<sub>41</sub></td><td>0.4493</td></tr><tr><td>L<sub>12</sub></td><td>-0.0428</td><td>L<sub>32</sub></td><td>-0.5564 - 0.2799s<sub>1</sub></td><td>L<sub>42</sub></td><td>-0.2830</td></tr><tr><td>L<sub>13</sub></td><td>0.0000</td><td>L<sub>33</sub></td><td>0.1990</td><td>L<sub>43</sub></td><td>-0.0558 - 0.3618s<sub>2</sub> - 0.1788s<sub>3</sub></td></tr><tr><td>L<sub>21</sub></td><td>-0.1663</td><td>L<sub>34</sub></td><td>0.1766</td><td>L<sub>44</sub></td><td>0.2470</td></tr><tr><td>L<sub>22</sub></td><td>-0.0694</td><td>L<sub>35</sub></td><td>0.1000</td><td>L<sub>45</sub></td><td>-0.0245</td></tr><tr><td>L<sub>23</sub></td><td>-0.1007</td><td>L<sub>36</sub></td><td>-0.0162</td><td>L<sub>46</sub></td><td>0.0001</td></tr><tr><td>L<sub>24</sub></td><td>0.0100</td><td>L<sub>37</sub></td><td>-0.0007</td><td>L<sub>51</sub></td><td>-1.9942</td></tr><tr><td>L<sub>25</sub></td><td>0.0001</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td><td>L<sub>52</sub></td><td>0.0045</td></tr><tr><td>σ<sub>u</sub></td><td>0.0310</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td><td>L<sub>61</sub></td><td>-0.0323</td></tr><tr><td>σ<sub>h</sub></td><td>0.0155</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td><td>L<sub>62</sub></td><td>0.0639</td></tr><tr><td>σ<sub>C</sub></td><td>0.0310</td><td>σ<sub>τ</sub></td><td>0.8999[tanh(4.5T<sub>C</sub>) + 1][1 + 0.25 cos(2πt/6)]</td><td></td><td></td></tr><tr><td>σ<sub>E</sub></td><td>0.0232</td><td></td><td></td><td></td><td></td></tr></table>

<table><caption>TABLE A7. Parameters of the new four-dimensional model with one latent variable.</caption><tr><td>L<sub>11</sub></td><td>-0.1003</td><td>L<sub>25</sub></td><td>0.1834</td><td>L<sub>41</sub></td><td>-1.5815</td></tr><tr><td>L<sub>12</sub></td><td>-0.0819</td><td>L<sub>26</sub></td><td>-10.1500 - 1.8734s<sub>1</sub></td><td>L<sub>42</sub></td><td>-0.1297</td></tr><tr><td>L<sub>13</sub></td><td>-0.0289</td><td>L<sub>27</sub></td><td>0.0122</td><td>s<sub>1</sub></td><td>sin(2πt/6)</td></tr><tr><td>L<sub>14</sub></td><td>-0.1069</td><td>L<sub>31</sub></td><td>0.3222</td><td>s<sub>2</sub></td><td>sin(2πt/6 + 2π/6)</td></tr><tr><td>L<sub>15</sub></td><td>0.0012</td><td>L<sub>32</sub></td><td>-0.1863</td><td>s<sub>3</sub></td><td>sin(2πt/3 + 2π/6)</td></tr><tr><td>L<sub>21</sub></td><td>0.3470</td><td>L<sub>33</sub></td><td>-0.2110 - 0.1402s<sub>2</sub> - 0.0768s<sub>3</sub></td><td>σ<sub>h</sub></td><td>0.0155</td></tr><tr><td>L<sub>22</sub></td><td>-0.2554 - 0.2584s<sub>1</sub></td><td>L<sub>34</sub></td><td>0.0782</td><td>σ<sub>C</sub></td><td>0.0310</td></tr><tr><td>L<sub>23</sub></td><td>0.0589</td><td>L<sub>35</sub></td><td>0.2580</td><td>σ<sub>E</sub></td><td>0.0232</td></tr><tr><td>L<sub>24</sub></td><td>-2.1976 - 0.5077s<sub>1</sub></td><td>L<sub>36</sub></td><td>-0.0017</td><td>σ<sub>τ</sub></td><td>2.2034</td></tr></table>

REFERENCES

AlMomani, A. A., and E. Bollt, 2020: ERFit: Entropic regression fit MATLAB package, for data-driven system identification of underlying dynamic equations. arXiv, 2010.02411v1, https://doi.org/10.48550/arXiv.2010.02411.

—, J. Sun, and E. Bollt, 2020: How entropic regression beats the outliers problem in nonlinear system identification. *Chaos*, **30**, 013107, https://doi.org/10.1063/1.5133386.

An, S.-I., and F.-F. Jin, 2004: Nonlinearity and asymmetry of ENSO. *J. Climate*, **17**, 2399–2412, https://doi.org/10.1175/1520-0442(2004)017<2399:NAAOE>2.0.CO;2.

Ashok, K., S. K. Behera, S. A. Rao, H. Weng, and T. Yamagata, 2007: El Niño Modoki and its possible teleconnection. *J. Geophys. Res.*, **112**, C11007, https://doi.org/10.1029/2006JC003798.

Barnston, A. G., M. K. Tippett, M. L. L'Heureux, S. Li, and D. G. DeWitt, 2012: Skill of real-time seasonal ENSO model predictions during 2002–11: Is our capability increasing? *Bull. Amer. Meteor. Soc.*, **93**, 631–651, https://doi.org/10.1175/BAMS-D-11-00111.1.

Battisti, D. S., and A. C. Hirst, 1989: Interannual variability in a tropical atmosphere–ocean model: Influence of the basic state, ocean geometry and nonlinearity. *J. Atmos. Sci.*, **46**, 1687–1712, https://doi.org/10.1175/1520-0469(1989)046<1687:IVIATA>2.0.CO;2.

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*, Seattle, WA, Amer. Meteor. Soc., 2.3, https://ams.confex.com/ams/84Annual/techprogram/paper_70720.htm.

Bellman, R. E., 1961: *Dynamic Programming Treatment of the Traveling Salesman Problem*. RAND Corporation, 7 pp.

Branicki, M., and A. J. Majda, 2012: Quantifying uncertainty for predictions with model error in non-Gaussian systems with

intermittency. *Nonlinearity*, **25**, 2543–2578, https://doi.org/10.1088/0951-7715/25/9/2543.

Brown, J., A. W. Tudhope, M. Collins, and H. V. McGregor, 2008: Mid-holocene ENSO: Issues in quantitative model proxy data comparisons. *Paleoceanography Paleoclimatol.*, **23**, PA3202, https://doi.org/10.1029/2007PA001512.

Brunton, S. L., J. L. Proctor, and J. N. Kutz, 2016: Discovering governing equations from data by sparse identification of nonlinear dynamical systems. *Proc. Natl. Acad. Sci. USA*, **113**, 3932–3937, https://doi.org/10.1073/pnas.1517384113.

Capotondi, A., and Coauthors, 2015: Understanding ENSO diversity. *Bull. Amer. Meteor. Soc.*, **96**, 921–938, https://doi.org/10.1175/BAMS-D-13-00117.1.

Chen, D., and M. A. Cane, 2008: El Niño prediction and predictability. *J. Comput. Phys.*, **227**, 3625–3640, https://doi.org/10.1016/j.jcp.2007.05.014.

—, and Coauthors, 2015: Strong influence of westerly wind bursts on El Niño diversity. *Nat. Geosci.*, **8**, 339–345, https://doi.org/10.1038/ngeo2399.

Chen, L., T. Li, B. Wang, and L. Wang, 2017: Formation mechanism for 2015/16 super El Niño. *Sci. Rep.*, **7**, 2975, https://doi.org/10.1038/s41598-017-02926-3.

Chen, N., 2020: Learning nonlinear turbulent dynamics from partial observations via analytically solvable conditional statistics. *J. Comput. Phys.*, **418**, 109635, https://doi.org/10.1016/j.jcp.2020.109635.

—, 2023: *Stochastic Methods for Modeling and Predicting Complex Dynamical Systems: Uncertainty Quantification, State Estimation, and Reduced-Order Models*. Springer Nature, 199 pp.

—, and A. J. Majda, 2017: Simple stochastic dynamical models capturing the statistical diversity of El Niño Southern Oscillation. *Proc. Natl. Acad. Sci. USA*, **114**, 1468–1473, https://doi.org/10.1073/pnas.1620766114.
—, and X. Fang, 2023: A simple multiscale intermediate coupled stochastic model for El Niño diversity and complexity. *J. Adv. Model. Earth Syst.*, **15**, e2022MS003469, https://doi.org/10.1029/2022MS003469.

—, and Y. Zhang, 2023: A causality-based learning approach for discovering the underlying dynamics of complex systems from partial observations with stochastic parameterization. *Physica D*, **449**, 133743, https://doi.org/10.1016/j.physd.2023.133743.

—, X. Fang, and J.-Y. Yu, 2022: A multiscale model for El Niño complexity. *npj Climate Atmos. Sci.*, **5**, 16, https://doi.org/10.1038/s41612-022-00241-x.

Cover, T. M., 1999: *Elements of Information Theory*. John Wiley & Sons, 542 pp.

Dai, A., and T. M. L. Wigley, 2000: Global patterns of ENSO-induced precipitation. *Geophys. Res. Lett.*, **27**, 1283–1286, https://doi.org/10.1029/1999GL011140.

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. *Commun. Earth Environ.*, **2**, 212, https://doi.org/10.1038/s43247-021-00285-6.

Eisenman, I., L. Yu, and E. Tziperman, 2005: Westerly wind bursts: ENSO's tail rather than the dog? *J. Climate*, **18**, 5224–5238, https://doi.org/10.1175/JCLI3588.1.

Elinger, J., 2021: Information theoretic causality measures for parameter estimation and system identification. Ph.D. dissertation, Georgia Institute of Technology, 172 pp.

—, and J. Rogers, 2021: Causation entropy method for covariate selection in dynamic models. *2021 American Control Conf. (ACC)*, New Orleans, LA, Institute of Electrical and Electronics Engineers, 2842–2847, https://doi.org/10.23919/ACC50511.2021.9483371.

Emerick, A. A., and A. C. Reynolds, 2013: Ensemble smoother with multiple data assimilation. *Comput. Geosci.*, **55**, 3–15, https://doi.org/10.1016/j.cageo.2012.03.011.

Evensen, G., 2003: The ensemble Kalman filter: Theoretical formulation and practical implementation. *Ocean Dyn.*, **53**, 343–367, https://doi.org/10.1007/s10236-003-0036-9.

Eyring, V., S. Bony, G. A. Meehl, C. A. Senior, B. Stevens, R. J. Stouffer, and K. E. Taylor, 2016: Overview of the Coupled Model Intercomparison Project phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.*, **9**, 1937–1958, https://doi.org/10.5194/gmd-9-1937-2016.

Fang, X., and N. Chen, 2023: Quantifying the predictability of ENSO complexity using a statistically accurate multiscale stochastic model and information theory. *J. Climate*, **36**, 2681–2702, https://doi.org/10.1175/JCLI-D-22-0151.1.

Fang, X.-H., and M. Mu, 2018: A three-region conceptual model for central Pacific El Niño including zonal advective feedback. *J. Climate*, **31**, 4965–4979, https://doi.org/10.1175/JCLI-D-17-0633.1.

—, and F. Zheng, 2021: Effect of the air–sea coupled system change on the ENSO evolution from boreal spring. *Climate Dyn.*, **57**, 109–120, https://doi.org/10.1007/s00382-021-05697-w.

—, —, and J. Zhu, 2015: The cloud-radiative effect when simulating strength asymmetry in two types of El Niño events using CMIP5 models. *J. Geophys. Res. Oceans*, **120**, 4357–4369, https://doi.org/10.1002/2014JC010683.

Fish, J., A. DeWitt, A. A. R. AlMomani, P. J. Laurienti, and E. Bollt, 2021: Entropic regression with neurologically motivated applications. *Chaos*, **31**, 113105, https://doi.org/10.1063/5.0039333.

Ford, H. L., A. C. Ravelo, and S. Hovan, 2012: A deep Eastern Equatorial Pacific thermocline during the early Pliocene warm period. *Earth Planet. Sci. Lett.*, **355–356**, 152–161, https://doi.org/10.1016/j.epsl.2012.08.027.

Geng, L., and F.-F. Jin, 2022: ENSO diversity simulated in a revised Cane-Zebiak model. *Front. Earth Sci.*, **10**, 899323, https://doi.org/10.3389/feart.2022.899323.

Geng, T., W. Cai, and L. Wu, 2020: Two types of ENSO varying in tandem facilitated by nonlinear atmospheric convection. *Geophys. Res. Lett.*, **47**, e2020GL088784, https://doi.org/10.1029/2020GL088784.

Gubner, J. A., 2006: *Probability and Random Processes for Electrical and Computer Engineers*. Cambridge University Press, 639 pp.

Hameed, S. N., D. Jin, and V. Thilakan, 2018: A model for super El Niños. *Nat. Commun.*, **9**, 2528, https://doi.org/10.1038/s41467-018-04803-7.

Herrnstein, R. J., and C. Murray, 2010: *The Bell Curve: Intelligence and Class Structure in American Life*. Simon and Schuster, 912 pp.

Hu, Z.-Z., A. Kumar, B. Jha, W. Wang, B. Huang, and B. Huang, 2012: An analysis of warm pool and cold tongue El Niños: Air–sea coupling processes, global influences, and recent trends. *Climate Dyn.*, **38**, 2017–2035, https://doi.org/10.1007/s00382-011-1224-9.

Huang, B., and Coauthors, 2017: Extended Reconstructed Sea Surface Temperature, version 5 (ERSSTv5): Upgrades, valuations, and intercomparisons. *J. Climate*, **30**, 8179–8205, https://doi.org/10.1175/JCLI-D-16-0836.1.

Jin, E. K., and Coauthors, 2008: Current status of ENSO prediction skill in coupled ocean–atmosphere models. *Climate Dyn.*, **31**, 647–664, https://doi.org/10.1007/s00382-008-0397-3.

Jin, F.-F., 1997a: An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model. *J. Atmos. Sci.*, **54**, 811–829, https://doi.org/10.1175/1520-0469(1997)054<0811:AEORPF>2.0.CO;2.

—, 1997b: An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-down coupled model. *J. Atmos. Sci.*, **54**, 830–847, https://doi.org/10.1175/1520-0469(1997)054<0830:AEORPF>2.0.CO;2.

—, 2022: Toward understanding El Niño Southern-Oscillation’s spatiotemporal pattern diversity. *Front. Earth Sci.*, **10**, 899139, https://doi.org/10.3389/feart.2022.899139.

—, and S.-I. An, 1999: Thermocline and zonal advective feedbacks within the equatorial ocean recharge oscillator model for ENSO. *Geophys. Res. Lett.*, **26**, 2989–2992, https://doi.org/10.1029/1999GL002297.

—, L. Lin, A. Timmermann, and J. Zhao, 2007: Ensemble-mean dynamics of the ENSO recharge oscillator under state-dependent stochastic forcing. *Geophys. Res. Lett.*, **34**, L03807, https://doi.org/10.1029/2006GL027372.

—, H.-C. Chen, S. Zhao, M. Hayashi, C. Karamperidou, M. F. Stuecker, R. Xie, and L. Geng, 2020: Simple ENSO models. *El Niño Southern Oscillation in a Changing Climate, Geophys. Monogr.*, Vol. 253, Amer. Geophys. Union, 119–151, https://doi.org/10.1002/9781119548164.ch6.

Kalanay, E., and Coauthors, 1996: The NCEP/NCAR 40-Year Reanalysis Project. *Bull. Amer. Meteor. Soc.*, **77**, 437–472, https://doi.org/10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2.

Kao, H.-Y., and J.-Y. Yu, 2009: Contrasting eastern-Pacific and central-Pacific types of ENSO. *J. Climate*, **22**, 615–632, https://doi.org/10.1175/2008JCLI2309.1.
Kay, S. M., 1988: *Modern Spectral Estimation: Theory and Application*. Prentice Hall, 543 pp.

Kim, J.-S., K.-Y. Kim, and S.-W. Yeh, 2012: Statistical evidence for the natural variation of the central Pacific El Niño. *J. Geophys. Res.*, **117**, C06014, https://doi.org/10.1029/2012JC008003.

Kim, P., J. Rogers, J. Sun, and E. Bollt, 2017: Causation entropy identifies sparsity structure for parameter estimation of dynamic systems. *J. Comput. Nonlinear Dyn.*, **12**, 011008, https://doi.org/10.1115/1.4034126.

Kleeman, R., 2011: Information theory and dynamical system predictability. *Entropy*, **13**, 612–649, https://doi.org/10.3390/ent13030612.

Klein, S. A., B. J. Soden, and N.-C. Lau, 1999: Remote sea surface temperature variations during ENSO: Evidence for a tropical atmospheric bridge. *J. Climate*, **12**, 917–932, https://doi.org/10.1175/1520-0442(1999)012<0917:RSSTVD>2.0.CO;2.

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. *J. Climate*, **22**, 1499–1515, https://doi.org/10.1175/2008JCLI2624.1.

Larkin, N. K., and D. E. Harrison, 2005: Global seasonal temperature and precipitation anomalies during El Niño autumn and winter. *Geophys. Res. Lett.*, **32**, L16705, https://doi.org/10.1029/2005GL022860.

Levine, A. F. Z., and F. F. Jin, 2017: A simple approach to quantifying the noise–ENSO interaction. Part I: Deducing the state-dependency of the windstress forcing using monthly mean data. *Climate Dyn.*, **48** (1–2), 1–18, https://doi.org/10.1007/s00382-015-2748-1.

Majda, A. J., and N. Chen, 2018: Model error, information barriers, state estimation and prediction in complex multiscale systems. *Entropy*, **20**, 644, https://doi.org/10.3390/e20090644.

McPhaden, M. J., S. E. Zebiak, and M. H. Glantz, 2006: ENSO as an integrating concept in Earth science. *Science*, **314**, 1740–1745, https://doi.org/10.1126/science.1132588.

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 Tech. Memo. 598, 44 pp., https://www.ecmfw.int/sites/default/files/elibrary/2009/11577-stochastic-parametrization-and-model-uncertainty.pdf.

Percival, D. B., and A. T. Walden, 1993: *Spectral Analysis for Physical Applications*. Cambridge University Press, 583 pp.

Philander, S. G. H., 1983: El Niño Southern Oscillation phenomena. *Nature*, **302**, 295–301, https://doi.org/10.1038/302295a0.

Picaut, J., F. Masia, and Y. Du Penhoat, 1997: An advective-reflective conceptual model for the oscillatory nature of the ENSO. *Science*, **277**, 663–666, https://doi.org/10.1126/science.277.5326.663.

Priestley, M. B., 1981: *Spectral Analysis and Time Series*. Academic Press, 706 pp.

Ren, H.-L., and F.-F. Jin, 2013: Recharge oscillator mechanisms in two types of ENSO. *J. Climate*, **26**, 6506–6523, https://doi.org/10.1175/JCLI-D-12-00601.1.

Roberts, A., J. Guckenheimer, E. Widiasih, A. Timmermann, and C. K. R. T. Jones, 2016: Mixed-mode oscillations of El Niño–Southern Oscillation. *J. Atmos. Sci.*, **73**, 1755–1766, https://doi.org/10.1175/JAS-D-15-0191.1.

Ropelewski, C. F., and M. S. Halpert, 1987: Global and regional scale precipitation patterns associated with the El Niño/Southern Oscillation. *Mon. Wea. Rev.*, **115**, 1606–1626, https://doi.org/10.1175/1520-0493(1987)115<1606:GARSPP>2.0.CO;2.

Santosa, F., and W. W. Symes, 1986: Linear inversion of band-limited reflection seismograms. *SIAM J. Sci. Stat. Comput.*, **7**, 1307–1330, https://doi.org/10.1137/0907087.

Santoso, A., and Coauthors, 2019: Dynamics and predictability of El Niño–Southern Oscillation: An Australian perspective on progress and challenges. *Bull. Amer. Meteor. Soc.*, **100**, 403–420, https://doi.org/10.1175/BAMS-D-18-0057.1.

Shannon, C. E., 1948: A mathematical theory of communication. *Bell Syst. Tech. J.*, **27**, 379–423, https://doi.org/10.1002/j.1538-7305.1948.tb01338.x.

Sohn, S.-J., C.-Y. Tam, and H.-I. Jeong, 2016: How do the strength and type of ENSO affect SST predictability in coupled models. *Sci. Rep.*, **6**, 33790, https://doi.org/10.1038/srep33790.

Stein, K., A. Timmermann, N. Schneider, F.-F. Jin, and M. F. Stuecker, 2014: ENSO seasonal synchronization theory. *J. Climate*, **27**, 5285–5310, https://doi.org/10.1175/JCLI-D-13-00525.1.

Stuecker, M. F., A. Timmermann, F.-F. Jin, S. McGregor, and H.-L. Ren, 2013: A combination mode of the annual cycle and the El Niño/Southern Oscillation. *Nat. Geosci.*, **6**, 540–544, https://doi.org/10.1038/ngeo1826.

Suarez, M. J., and P. S. Schopf, 1988: A delayed action oscillator for ENSO. *J. Atmos. Sci.*, **45**, 3283–3287, https://doi.org/10.1175/1520-0469(1988)045<3283:ADAOFE>2.0.CO;2.

Thual, S., and B. Dewitte, 2023: ENSO complexity controlled by zonal shifts in the walker circulation. *Nat. Geosci.*, **16**, 328–332, https://doi.org/10.1038/s41561-023-01154-x.

—, A. J. Majda, N. Chen, and S. N. Stechmann, 2016: Simple stochastic model for El Niño with westerly wind bursts. *Proc. Natl. Acad. Sci. USA*, **113**, 10 245–10 250, https://doi.org/10.1073/pnas.1612002113.

—, —, and —, 2019: Statistical occurrence and mechanisms of the 2014–2016 delayed super El Niño captured by a simple dynamical model. *Climate Dyn.*, **52**, 2351–2366, https://doi.org/10.1007/s00382-018-4265-5.

Tibshirani, R., 1996: Regression shrinkage and selection via the lasso. *J. Roy. Stat. Soc.*, **58B**, 267–288, https://doi.org/10.1111/j.2517-6161.1996.tb02080.x.

Timmermann, A., 2003: Decadal ENSO amplitude modulations: A nonlinear paradigm. *Global Planet. Change*, **37**, 135–156, https://doi.org/10.1016/S0921-8181(02)00194-7.

—, and F.-F. Jin, 2002: A nonlinear mechanism for decadal El Niño amplitude changes. *Geophys. Res. Lett.*, **29**, 1003, https://doi.org/10.1029/2001GL013369.

—, —, and J. Abshagen, 2003: A nonlinear theory for El Niño bursting. *J. Atmos. Sci.*, **60**, 152–165, https://doi.org/10.1175/1520-0469(2003)060<0152:ANTFEN>2.0.CO;2.

—, and Coauthors, 2018: El Niño–Southern Oscillation complexity. *Nature*, **559**, 535–545, https://doi.org/10.1038/s41586-018-0252-6.

Tippett, M. K., R. Kleeman, and Y. Tang, 2004: Measuring the potential utility of seasonal climate predictions. *Geophys. Res. Lett.*, **31**, L22201, https://doi.org/10.1029/2004GL021575.

Tziperman, E., S. E. Zebiak, and M. A. Cane, 1997: Mechanisms of seasonal–ENSO interaction. *J. Atmos. Sci.*, **54**, 61–71, https://doi.org/10.1175/1520-0469(1997)054<0061:MOSEI>2.0.CO;2.

Vialard, J., C. Menkes, J.-P. Boulanger, P. Delecluse, E. Guilyardi, M. J. McPhaden, and G. Madec, 2001: A model study of oceanic mechanisms affecting equatorial Pacific sea surface temperature during the 1997–98 El Niño. *J. Phys. Oceanogr.*, **31**,
1649–1675, https://doi.org/10.1175/1520-0485(2001)031<1649: AMSOOM>2.0.CO;2.

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. *Proc. Natl. Acad. Sci. USA*, **116**, 22512–22517, https://doi.org/ 10.1073/pnas.1911130116.

Wang, C., 2001: A unified oscillator model for the El Niño–Southern Oscillation. *J. Climate*, **14**, 98–115, https://doi.org/10.1175/1520-0442(2001)014<0098:AUOMFT>2.0.CO;2.

Weisberg, R. H., and C. Wang, 1997: A western Pacific oscillator paradigm for the El Niño-Southern Oscillation. *Geophys. Res. Lett.*, **24**, 779–782, https://doi.org/10.1029/97GL00689.

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. *J. Climate*, **34**, 3481–3501, https:// doi.org/10.1175/JCLI-D-20-0447.1.

Ying, X., 2019: An overview of overfitting and its solutions. *J. Phys.*, **1168**, 022022, https://doi.org/10.1088/1742-6596/1168/ 2/022022.

Yu, J.-Y., and H.-Y. Kao, 2007: Decadal changes of ENSO persis-
tence barrier in SST and ocean heat content indices: 1958–
2001. *J. Geophys. Res.*, **112**, D13106, https://doi.org/10.1029/ 2006JD007654.

—, and S. T. Kim, 2013: Identifying the types of major El Niño events since 1870. *Int. J. Climatol.*, **33**, 2105–2112, https://doi.org/10.1002/joc.3575.

Zheng, F., X.-H. Fang, J.-Y. Yu, and J. Zhu, 2014: Asymmetry of the Bjerknes positive feedback between the two types of El Niño. *Geophys. Res. Lett.*, **41**, 7651–7657, https://doi.org/10.1002/2014GL062125.