Rigorous Derivation of Stochastic Conceptual Models for the El Niño-Southern Oscillation from a Spatially-Extended Dynamical System

The simple spatially-extended stochastic model was originally developed in [13]. The starting system is a coupled atmosphere-ocean-SST model, which is given by

Since the most predominant behavior of the ENSO is around the equator and along the zonal direction, it is natural to apply a meridional expansion, which serves as the separation of variables, and then implement a meridional truncation to the first basis in both atmosphere $\phi_{0}(Y)$ and ocean $\psi_{0}(y)$ to reduce the complexity of the system. The basis functions in the meridional direction are given by the parabolic cylinder functions [42]. See Panels (a)–(b) in Figure 1 for the meridional basis functions, where the first one has a Gaussian profile. The meridional truncation triggers atmosphere Kelvin, Rossby waves $K_{A},R_{A}$, and ocean Kelvin, Rossby waves $K_{O},R_{O}$. Therefore, after removing the meridional dependence from (1)–(3), the coupled system (with dependence only on time $t$ and zonal coordinate $x$) yields.

In addition to the interannual state variables, the intraseasonal variabilities are also part of the indispensable ENSO dynamics. Although these intraseasonal variabilities behave like random noise in the interannual time scale, they play vital roles in generating the ENSO events and increasing the ENSO complexity. In particular, atmosphere wind bursts in the intraseasonal time scale, such as the westerly wind bursts (WWBs) [43, 44, 45], easterly wind bursts (EWBs) [46, 47] and the MJO [48, 49], are essential in triggering and terminating the El Niño events. Due to the fast dynamics of wind activities, it is natural to adopt a simple stochastic process to effectively characterize their time evolutions. With the stochastic wind bursts, the zonal wind stress $\tau$ now has two components, $\tau=\gamma(u+u_{p})$, where $u$ is directly from the atmosphere model (4) while $u_{p}$ is the contribution from the stochastic wind bursts. Here,

where for simplicity $s_{p}(x)$ is a fixed spatial basis function being localized in the western Pacific (see Panel (c) of Figure 1), since most of the active wind bursts are observed there. The stochastic amplitude $a_{p}$ is given by the following stochastic differential equation

where $\dot{W}(t)$ is a standard white noise [50] and $d_{p}$ is a damping term characterizing the temporal correlation of the wind bursts at the intraseasonal time scale.

Note that the noise coefficient ${\sigma_{p}(T_{W})}$ is state-dependent (namely, the multiplicative noise) as a function of the averaged SST in the western Pacific. This is because a warmer SST induces more convective activities, which in turn triggers more wind bursts [44, 43]. Here, the state dependence of ${\sigma_{p}(T_{W})}$ is given by a two-state Markov jump process with one active state and one quiescent state:

The transition rate from the quiescent to the active state $\mu_{01}$ and that from the active to the quiescent state $\mu_{10}$ are shown in Panel (d) of Figure 1. The explicit equations of them are in the following:

Here, transition rates $\mu_{01}$ and $\mu_{10}$ are correlated and anti-correlated with $T_{W}$, respectively, which reflect the nature of the multiplicative noise. Note that the wind burst activity from the model does not favor specifically westerly or easterly. The WWBs and EWBs generated from the model are entirely random. Only the amplitude of the wind bursts is determined by the transitions related to the SST. The two-state Markov jump process for parameterizing the wind activities is utilized here for simplicity, categorizing the winds by only the active and quiescent phases. A continuous dependence of the wind strength on the SST can easily be incorporated, which leads to a similar SST behavior in the ENSO dynamics [51, 52].

As a final remark, the coupled model (4)–(9) focuses on the eastern Pacific El Niño events, which is the primary interest of this study. The model does not include the mechanism of generating the central Pacific events [53, 54, 55]. Incorporating the central Pacific events into the model will be briefly discussed in Section 6.

The dynamical core of the system (4)–(6) is deterministic, linear and stable. In other words, the solution is a linear damped regular oscillator. The random wind bursts serve as external forcing to trigger the El Niño events and introduce irregularity in the solution. The subsequent La Niña is a consequence of the decaying phase of the ENSO cycle.

Panel (e) in Figure 1 shows one simulation of the coupled model (4)–(9) for different fields. The model can reproduce the quasi-regular oscillation behavior of the SST anomaly. When the SST in the western Pacific increases, the wind activities become more significant. If the strong wind activity is dominated by the WWBs, then El Niño events are likely to be triggered. During this process, the center of the atmosphere wind convergence moves towards the east, and the strong eastward ocean zonal current is observed in the western Pacific. The thermocline depth becomes deeper in the eastern Pacific and induces the El Niño events. The model can generate not only the moderate El Niño events, but also the extreme events, known as the super El Niños, that have strong amplitudes (e.g., the event at $t=490$). In addition to reproducing the super El Niño that mimics the observed 1997-1998 event (e.g., at $t=490$), the model succeeds in generating the so-called delayed super El Niño [56], similar to the observed 2014-2016 event [57, 47]. One such example is the event around $t=474$-$476$, which is due to the peculiar wind structures. A series of WWBs induce the onset of an El Niño event. However, a strong EWB occurs immediately to present the development of such an event to a super El Niño. After a year, another period of strong WWBs happen to trigger the extreme event.

Averaging the SST over the eastern Pacific, the resulting time series is denoted by $T_{E}$. The model can reproduce very similar non-Gaussian probability density function (PDF) of $T_{E}$ as the observations [13]. This non-Gaussian feature is reproduced with the contribution from the multiplicative noise, without which the model is linear and the associated PDF is Gaussian. Likewise, the spectrum of the model has the most significant power within the band of 3 to 7 years, which is also consistent with observations.

Without the stochastic wind bursts, the coupled atmosphere-ocean-SST model (4)–(6) is a linear, deterministic and stable model. Since $u$ (and therefore the wind stress $\tau$) is a function of $K_{A}$ and $R_{A}$, the coupled system can be rewritten as

where a small damping $d_{A}=10^{-8}$ is added to guarantee the solution is unique [13], and the new constants introduced in (12) are the following ones:

The domain of $K_{O},R_{O}$ and $T$ is the equatorial Pacific ocean $[0,L_{O}]$ while the domain of $K_{A}$ and $R_{A}$ is the entire equatorial band $[0,L_{A}]$.

Define a vector $\mathbf{u}$ that contains all the prognostic variables at discrete equal-partitioned grid points,

where $N_{O}$ is the total number of grid points over the equatorial Pacific region $[0,L_{O}]$ and $\cdot^{\mathtt{T}}$ is the vector transpose. Using an upwind scheme, the discrete form of the linear system (12) is given by:

The detailed form of $\mathbf{M}$ is included in the Appendix. With the stochastic wind burst, the above equation can be formally written as:

where $\mathbf{F}$ is associated with the stochastic forcing. Denote by $\mathbf{s}_{\mathbf{u}}=(\chi_{O}c_{1}\gamma s_{p}/2,\\ -\chi_{O}c_{1}\gamma s_{p}/3,\mathbf{0})^{\mathtt{T}}$ and then $\mathbf{F}=\mathbf{s}_{\mathbf{u}}a_{p}$, where $s_{p}$ is written in the discrete form at grid points $1,\ldots,N_{O}$ and $\mathbf{0}$ is a $1\times N_{O}$ vector with all zero entries. This is the wind burst forcing acting on the three components $\mathbf{u}=(\mathbf{K}_{O},\mathbf{R}_{O},\mathbf{T})$ (ocean Kelvin wave, ocean Rossby wave, and SST) of the state variables. See Panel (a) of Figure 2 for the profile of $\mathbf{s}_{\mathbf{u}}$. As expected, a positive response is on the ocean Kelvin waves and a negative one is on the ocean Rossby waves, both in the western Pacific. The forcing is not directly imposed on SST; therefore, the third segment of $\mathbf{s}_{\mathbf{u}}$, corresponding to SST components, is zero.

It can be numerically validated that the $3N_{O}$ eigenvalues all have distinct values and are non-zero. In addition, the $i$-th left eigenvector and the $j$-th right eigenvector of the matrix $\mathbf{M}$ are orthogonal with each other for $i\neq j$. Therefore, in light of the matrix similarity, the following result is arrived at:

or equivalently

where $\mathbf{\Lambda}$ is a diagonal matrix with diagonal entries being eigenvalues and $\mathbf{X}$ contains the $3N_{O}$ right eigenvectors. Since the right eigenvectors $\mathbf{X}$ are complete, the matrix inverse exists.

Panel (a) of Figure 3 shows the eigenvalues associated with the $3N_{O}\times 3N_{O}$ matrix $\mathbf{M}$, where $N_{O}=56$ is utilized here and in the subsequent numerical simulations. These eigenvalues are complex numbers. The real and imaginary parts are shown in the y- and x-axis, representing the growth rate and the frequency, respectively. Since the starting system is stable, the growth rates of all the eigenvalues are negative, representing the decay of the solution. These eigenvalues are ranked according to the growth rate. Therefore, the leading two eigenmodes (marked by red dots) are the dominant modes of the system since they have the slowest decaying rates. These two modes appear in pairs; both the eigenvalues and the eigenvectors are complex conjugates. Their eigenvalues have a frequency of $4.6$ years and a decay rate of $1.5$ years, which are consistent with the observed features of ENSO. The reconstructed spatiotemporal patterns associated with these two eigenmodes (by setting the growth rate to be zero for the plotting purpose), as shown in Panel (b) of (3), mimic the large-scale ENSO oscillation solution. Notably, the leading two eigenmodes are very robust with respect to the spatial resolution $N_{O}$ as long as $N_{O}$ is not too small. Therefore, the analysis here provides a theoretical justification for developing a low-order model based on these leading two eigenmodes.

Before developing a reduced-order conceptual model based on the eigenvalue analysis, as shown above, it is essential to study the similarity and the gap between the full solution from the spatially-extended stochastic system (4)–(9) and that from the low-dimensional representation consisting of only the leading two eigenmodes, namely the projected solution.

Figure 4 includes a comparison of the SST between these two solutions. The projected solution (Panel (b)) is qualitatively similar to the full solution (Panel (a)) with an almost identical pattern and a nearly negligible residual error (Panel (c)). This is not surprising since the decaying rates of the remaining eigenmodes are more significant than those of the leading two components. Thus, the time evolution of the signal is overall explained by the two leading modes. In addition to the spatiotemporal pattern, the two key statistics, namely the power spectrum and the PDF, of the averaged SST over the eastern Pacific $T_{E}$ are almost indistinguishable between the full and the projected solutions. In particular, the power spectrum peaks between 3 and 7 years, which is consistent with the observed frequency of the ENSO cycles. The projected solution also recovers the non-Gaussian PDF of $T_{E}$, which has almost the same skewness and kurtosis as the full solution. Non-Gaussianity is a crucial feature of the ENSO statistics, which indicates a non-symmetry between the El Niño and the La Niña. The fat tail on the positive side corresponds to the super El Niños, which remains in the projected solution.

Figure 5 includes a more detailed comparison between the full and the projected solutions by showing different physical and wave variables. A delayed super El Niño event (from $t=474$ to $t=476$) occurs during this period. Overall, the projected solution captures the spatiotemporal patterns of the full one. Especially, the delayed super El Niño is reproduced in the projected solution. Yet, the small-scale wave propagation is missed in the projected solution. This is particularly unambiguous in the panels of the ocean Kelvin wave $K_{O}$ and the ocean Rossby wave $R_{O}$. As a result, the values across the western Pacific or the eastern Pacific remain constant for different fields. This is perhaps the most significant difference between the full and the projected solutions. Small-scale wave propagations are expected to appear at least with three degrees of freedom (e.g., either three variables or three regions across Pacific) [9]. Nevertheless, the projected solution clearly captures the large-scale eastward moving patterns, representing the ENSO development.

Figure 6 compares the full observed SST using the reanalysis data (Panel (a)) with the projected data using the leading two eigenmodes computed from the model. Panel (b) shows the results by directly projecting the SST to these two modes. The projected solution resembles the full observations in the eastern Pacific with a pattern correlation above $0.85$ at all longitude grid points. In contrast, the projected solution in the central-western Pacific differs from the full solution. This is not surprising since the starting model and the low-order presentation here focus only on the eastern Pacific El Niño events and do not have the mechanism to create the central Pacific events. The projected solution in the western Pacific captures the overall patterns of the full observations, though the details are somewhat missed. In addition to using the SST data, the atmosphere wind data can also be used to reconstruct the low-order representation of the SST. See Panel (c). Here, the atmosphere wind data is first projected to the leading two eigenvectors corresponding to the atmosphere components to compute the projection coefficients. Then these coefficients are multiplied by the eigenvectors corresponding to the SST components for reconstruction. The reconstructed low-order representation of SST using the atmosphere reanalysis data has a similar pattern as that based on the direct projection using the SST in the eastern Pacific. Finally, Panel (d) shows the reconstructed SST data using a direct data-driven regression method. This is achieved by first computing the correlation coefficient $r(x)$ between the SST at each longitude with the Niño 3 SST index using the observational data and then carrying out the spatiotemporal reconstruction using the formula $\mbox{SST}(x,t)=r(x)\mbox{Ni\~{n}o 3}(t)$. Note that this procedure uses the entire observational data and can be regarded as one of the optimal linear reconstruction methods. It serves as a benchmark solution. The results show that, similar to the low-dimensional projected solutions in Panels (b)–(c), the reconstruction in the eastern Pacific resembles the truth while that in the central-western Pacific contains a significant error. The comparison between the results in Panels (b)–(c) and Panel (d) implies that the reconstruction using the model-induced basis functions is skillful in reproducing the eastern Pacific El Niño events. It also indicates that both the wind and the SST data can be used for the reconstruction, allowing large freedom in choosing the available data sets.

Plugging the relationship of the matrix decomposition $\mathbf{M}=\mathbf{X}\mathbf{\Lambda}\mathbf{X}^{-1}$ in (17) into the full system (15) yields

where $\mathbf{\Lambda}$ is a diagonal matrix. Then multiplying (18) by $\mathbf{X}^{-1}$ gives $3N_{O}$ independent equations

where

Once $\mathbf{v}$ is solved, $\mathbf{u}$ can be obtained by

Denote by $\mathbf{z}^{\mathtt{T}}$ and $\bar{\mathbf{z}}^{\mathtt{T}}$ the first two rows of $\mathbf{X}^{-1}$. See Panel (b) of Figure 2 for the profile of $\mathbf{z}$. Likewise, denote by $\mathbf{x}$ and $\bar{\mathbf{x}}$ the first two columns of $\mathbf{X}$. The corresponding leading two eigenvalues are denoted by $\lambda_{o}$ and $\bar{\lambda}_{o}$ (the two red dots in Panel (a) of Figure 3), which are written as $-d_{o}\pm i\omega_{o}$ with $-d_{o}<0$ and $\omega_{o}$ being the growth rate and the frequency, respectively. Here $\cdot^{\mathtt{T}}$ is the transpose, and $\bar{\cdot}$ is the complex conjugate. Clearly, $\mathbf{z}^{\mathtt{T}}\mathbf{x}=\bar{\mathbf{z}}^{\mathtt{T}}\bar{\mathbf{x}}=1$.

Figure 7 shows the eigenvector associated with $\lambda_{o}$ for different model components.

It is worthwhile to highlight that, as was shown in (13), the solution $\widetilde{\mathbf{u}}$ consists of the ocean Kelvin wave $K_{O}$, the ocean Rossby wave $R_{O}$, and the SST $T$. Recall in (7) that the ocean current $U$ and the thermocline depth $H$ are the linear combinations of $K_{O}$ and $R_{O}$. In addition, the atmosphere wind $u$ and the potential temperature $\theta$ can be solved by exploiting the linear relationship with the ocean and SST variables according to (12) (see Appendix for details). Therefore, the two-dimensional low-order conceptual model (22) provides the solutions for all the atmosphere, ocean, and SST fields with detailed spatiotemporal solutions. This contrasts with many existing conceptual models, which focus only on a small set of specified state variables and often merely involve the averaged quantities over a large domain.

The general framework in (22) also allows the development of the reduced-order models with arbitrarily specified quantities (in addition to the direct model output $u$, $\theta$, $H$, $U$ and $T$) as state variables. This can be achieved by introducing two new real-valued state variables $\widetilde{a}$ and $\widetilde{b}$, which are the linear combinations of $v$ and $\bar{v}$ by angle rotations and amplifications,

Here $\alpha$ and $\beta$ are rotating angles to be in phase with the two new physical variables of interest while the real variables $A_{a}$ and $A_{b}$ are the amplitudes. Here, it is required that $\alpha,\beta\in[-\pi,\pi]$ and $\beta\neq\alpha$ to rule out the degenerated case. The two variables $\widetilde{a}$ and $\widetilde{b}$ can be defined as any variables that can be represented by a linear combination of the model state variables.

The above framework is adapted to the cases where the main interest lies in certain averaged quantities, as in many of the classical conceptual models for the ENSO. For example, $\widetilde{a}$ can be the SST in the eastern Pacific and $\widetilde{b}$ the thermocline in the western Pacific. They can also be the atmosphere wind and ocean current averaged over the entire Pacific or a specific region. In such situations, the amplitudes $A_{a}$ and $A_{b}$ can be determined by two scalars $x_{\tilde{a}}$ and $x_{\tilde{b}}$, which are the averaged values of $\mathbf{x}_{\tilde{a}}$ and $\mathbf{x}_{\tilde{b}}$ over a certain domain. This allows the development of a two-dimensional system for $\widetilde{a}$ and $\widetilde{b}$ that is linked with the traditional conceptual models.

See the Appendix for the proofs.

The discharge-recharge oscillator argues that discharge and recharge of equatorial heat content cause the coupled system to oscillate [7]. Define $\widetilde{a}:=T_{E}$ and $\tilde{b}:=H_{W}$ in (29), representing the averaged SST in the eastern Pacific and the averaged thermocline depth in the western Pacific. The corresponding scalars $x_{T_{E}}$ and $x_{H_{W}}$ are given by

which are computed by averaging the eigenvectors associated with $T$ and $H$ over the eastern Pacific and western Pacific, respectively. See Figure 8.

The other coefficients in (29) can be computed using the eigenvalues of the leading two modes $\lambda_{o}$ as well as $x_{T_{E}}$ and $x_{H_{W}}$ in (31). They are:

and

The remaining thing is to determine the multiplicative noise coefficient $\sigma_{p}(T_{W})$ of the stochastic amplitude $a_{p}$ in (9), which is a function of the averaged SST in the western Pacific $T_{W}$. Yet, $T_{W}$ is not the state variable in the two-dimensional discharge-recharge oscillator model. Nevertheless, $T_{W}$ can be written as a linear function of $T_{E}$ and $H_{W}$. Recall that

which can be used to solve $\mbox{Re}(v)$ and $\mbox{Im}(v)$,

Once the basis $v$ is determined, the value of $T_{W}$ can be computed,

which are functions of $H_{W}$ and $T_{E}$, and $x_{T_{W}}$ is given by

The direct representation of the $T_{W}$ by $T_{E}$ also justifies and unifies the use of SST in different regions (averaged values over the eastern Pacific, over the western Pacific or basin averaged value) as multiplicative noise coefficients in various conceptual models [58, 59, 60, 61].

After obtaining the value of $T_{W}$, the multiplicative noise coefficient $\sigma_{p}(T_{W})$ is determined. Note that the $\sigma_{p}(T_{W})$ is slightly modified to $\sigma_{p0}=0.2$ and $\sigma_{p1}=1.8$ to compensate for the approximation error from the low-order representation.

According to (32), the coefficient $\widetilde{c}_{11}$ and $\widetilde{c}_{12}$ associated with $T_{E}$ equations are both positive. The coefficient $\widetilde{c}_{11}>0$ represents the Bjerknes positive feedback while $\widetilde{c}_{12}>0$ is the recharge mechanism. On the other hand, a negative coefficient $\widetilde{c}_{21}$ appears in the equation of $H_{W}$, which indicates the discharge mechanism. Another finding is that, although the wind bursts directly impact the SST in heuristic thinking, the stochastic forcing appears in both the equations of $T_{E}$ and $H_{W}$ for the completeness of the stochastic discharge-recharge diagram. The opposite signs of $\alpha_{T_{E}}$ and $\alpha_{H_{W}}$ imply that the increase of SST in the eastern Pacific by the WWBs will simultaneously cause a decrease of the thermocline in the western Pacific.

Figure 9 shows the time series and phase plot of both the deterministic and stochastic versions of the discharge-recharge models developed above. The deterministic model (by setting the damping to be zero and ignoring the stochastic terms) leads to a regular oscillator while the time series of $H_{W}$ and $T_{E}$ generated from the stochastic model mimic those in observations. It is worth highlighting that the phase difference between $H_{W}$ and $T_{E}$ is $122^{o}$ (not $\pi/2$ or $90^{o}$), which is also confirmed by the lag between the two time series. Such a phase difference also resembles that from the observations. In addition, the observational time series and the stochastic model simulation of $T_{E}$ are compared. Due to the random noise, the trajectories from the stochastic model and the observations do not expect to have a one-to-one point-wise match between each other. Nevertheless, the stochastic model simulation succeeds in capturing the EP events. Finally, the PDF and the autocorrelation function (ACF) of $T_{E}$ from the stochastic model have very similar behavior to the observations. The PDF is skewed with a one-sided fat tail, which indicates that the model can reproduce a sufficient number of the super El Niño events. Similarly, matching the ACF with observations implies that the time series of $T_{E}$ from the model has, on average, the same temporal correlation (namely the time scale of the memory effect) as the observations. These findings suggest that the low-order model can capture the observed dynamical and statistical features of the time series $T_{E}$, which is an important indicator for the eastern Pacific El Niño.

The two-dimensional stochastic oscillator (29) can be further reduced to a one-dimensional oscillator with a delayed effect, which leads to the delayed oscillator of the ENSO.

See the Appendix for the proofs. Numerically, these constants and coefficients are given by

The detailed expression of the stochastic forcing with the memory effect derived here is an essential supplement to the original delayed oscillator model [29, 16], which was deterministic. Several conclusions can be drawn from this delayed oscillator. First, In front of the delayed term $T_{E}\left(t-\varphi/\omega_{o}\right)$, there is a factor with an exponential decay $e^{-d_{o}\varphi/\omega_{o}}$. This is a crucial feature, which indicates that the contribution of the delayed signal will be affected by the damping effect and becomes weaker as time goes on. Second, there are two ways that stochastic noise contributes to the system. The term $\alpha_{T_{E}}a_{p}(t)$ provides an instantaneous response to the system, while the stochastic forcing inside the integral represents the contribution from the noise accumulation. The latter also has a delayed contribution. Notably, although the damping time of the noise $d_{p}$ is intraseasonal, the noise interacting with the coupled system can affect the interannual time scale. Third, similar to the discharge-recharge oscillator, the delay here is $\varphi=2.1287$ ($122^{o}$) instead of $\pi/2$ ($90^{o}$). Finally, the sign of the delayed term depends on the sign of $\widetilde{c}_{12}$, which can be either positive or negative, depending on the state variables used for the model. The parameter $\widetilde{c}_{12}$ is positive in the above model since $\varphi=2.1287>\pi/2$. Yet, when the angle $\varphi$ is less than $\pi/2$, then $\widetilde{c}_{12}<0$ representing a negative delay, which is consistent with the discussions in [16]. In particular, if the average of the thermocline over the whole Pacific is adopted as was in [16], then a negative delay (about $\pi/6$) appears.

This section will provide details for parameter values and dimensional values.