A Stochastic Conceptual Model for the Coupled ENSO and MJO

The interannual ocean model (2.1) is a generalized extension of the classical recharge oscillator model (Jin, 1997), applied over the WP, CP and EP to capture both types of ENSO events and their diversity (Chen et al., 2022b). This model incorporates the ocean’s heat content discharge-recharge processes and zonal advection. In the model, $T_{C}$ and $T_{E}$ represent the SST anomalies in the CP and EP, respectively, $U$ denotes the zonal ocean current in the CP, and $h_{W}$ corresponds to the thermocline depth in the WP. The model reads:

	$\displaystyle\frac{{\,\rm d}U}{{\,\rm d}t}$	$\displaystyle=-rU-\frac{\alpha_{1}b_{0}\mu}{2}(T_{C}+T_{E})+{\beta_{u}(I)\tau},$		(2.1a)
	$\displaystyle\frac{{\,\rm d}h_{W}}{{\,\rm d}t}$	$\displaystyle=-rh_{W}-\frac{\alpha_{2}b_{0}\mu}{2}(T_{C}+T_{E})+{\beta_{h}(I)\tau},$		(2.1b)
	$\displaystyle\frac{{\,\rm d}T_{C}}{{\,\rm d}t}$	$\displaystyle=\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}+\rho IU+C_{U}+{\beta_{C}(I)\tau},$		(2.1c)
	$\displaystyle\frac{{\,\rm d}T_{E}}{{\,\rm d}t}$	$\displaystyle=\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.$		(2.1d)

Two additional variables appear in the interannual ocean model, which will be discussed in the following subsections. The first is the decadal variable $I$, appearing in (2.1c). This variable accounts for the long-term modulation of Walker circulation’s influence on ENSO, allowing the model to favor the development of either EP or CP events. It also represents the zonal SST gradient between the WP and CP, affecting interannual variability by modulating the efficiency of zonal advection through the $IU$ term. The second key variable is $\tau$, which represents atmospheric wind in the WP region and incorporates both intraseasonal and interannual components. In addition, the influence of seasonality has been incorporated into the two damping coefficients, $c_{1}$ and $c_{2}$, to more accurately capture the seasonal phase-locking behavior of ENSO. This behavior is characterized by the tendency of ENSO events to peak during boreal winter (Tziperman et al., 1997, Stein et al., 2014, Fang and Zheng, 2021).

It is worth explaining the model mechanisms for triggering CP and EP El Niño events. In the absence of the stochastic wind forcing, the model is a deterministic and nearly linear system if the decadal variable $I$ is held constant. Note that there is a non-linearity in the damping term of (2.1c), but it is weak and does not play a significant role when examining the function of the strong non-linear term $IU$. The model exhibits different oscillatory patterns depending on the fixed value of $I$ (Fang and Mu, 2018). When $I$ is small, corresponding to weak advection, the model generates an EP El Niño cycle. In contrast, when $I$ is large, the strengthened Walker circulation enhances the role of advection, leading to a dominant CP El Niño oscillatory pattern. The time-varying $I$, combined with the stochastic forcing, allows for the occurrence of different El Niño events over time. Depending on the value of $I$, one type of event becomes more likely during a given period, thereby reproducing the observed alternation between EP- and CP-dominant phases (Yu and Kim, 2013).

The decadal variability in the conceptual model is primarily used to modulate the background state over the equatorial Pacific (Capotondi et al., 2023, Power et al., 2021), favoring the generation of more EP or CP events over specific time periods. Developing a comprehensive model describing the decadal variation is not the main focus. To this end, a simple linear stochastic process is employed as the governing equation for the decadal variable $I$,

$$\frac{{\,\rm d}I}{{\,\rm d}t}=-\lambda(I-\bar{I})+\sigma_{I}(I)\dot{W}_{I},$$

(2.2)

where the damping coefficient $\lambda$ is chosen such that $I$ varies in the decadal timescale and $\bar{I}$ is the mean state. It is important to note that the trade winds in the lower level of the Walker circulation on decadal timescales are predominantly easterly. This implies that the sign of $I$ remains constant over time, resulting in a non-Gaussian distribution of $I$. This feature can be easily integrated into the linear stochastic process of $I$ by using a state-dependent noise coefficient (Averina and Artemiev, 1988, Yang et al., 2021). Given the limited observational data and the principle of deriving the least biased maximum entropy solution for a distribution (Chen, 2023), a uniform distribution is adopted for $I$ in this work. Here, a larger value of $I$ corresponds to stronger easterly trade winds.

The intraseasonal model generates intraseasonal wind, a key trigger for El Niño events. Unlike prior studies that used a single process or noise to represent wind in ENSO models (Chen et al., 2022b, Geng et al., 2020, Jin et al., 2007), this model employs simple atmospheric equations to capture more detailed dynamics, including the MJO. This allows for studying MJO-ENSO coupling.

The governing equations are based on the stochastic skeleton model for the MJO (Thual et al., 2014). However, only the modes up to the third zonal wavenumber are retained, simplifying the model while still capturing large-scale MJO features (Majda and Stechmann, 2009, 2011, Ling et al., 2017). The governing equations, written in terms of Fourier modes, are as follows,

	$\displaystyle\frac{{\,\rm d}\hat{K}_{k}}{{\,\rm d}t}$	$\displaystyle=-d_{k}\hat{K}_{k}-ik\hat{K}_{k}-\overline{H}\hat{A}_{k}/2,$		(2.3a)
	$\displaystyle\frac{{\,\rm d}\hat{R}_{k}}{{\,\rm d}t}$	$\displaystyle=-d_{k}\hat{R}_{k}+ik\hat{R}_{k}/3-\overline{H}\hat{A}_{k}/3,$		(2.3b)
	$\displaystyle\frac{{\,\rm d}\hat{Q}_{k}}{{\,\rm d}t}$	$\displaystyle=-d_{k}\hat{Q}_{k}-ik\overline{Q}(\hat{K}_{k}-\hat{R}_{k}/3)+\overline{H}\hat{A}_{k}(\overline{Q}/6-1)+\sigma_{\hat{Q}_{k}}(E_{q})\dot{W}_{\hat{Q}_{k}},$		(2.3c)
	$\displaystyle\frac{{\,\rm d}\hat{A}_{k}}{{\,\rm d}t}$	$\displaystyle=-\lambda_{A}\hat{A}_{k}+\sum_{j}\Gamma\hat{Q}_{j}(\hat{A}_{k-j}+\bar{A}_{k-j})+\sigma_{\hat{A}_{k}}(Q,A)\dot{W}_{\hat{A}_{k}},$		(2.3d)

where $k=0,\pm 1,\pm 2$ and $\pm 3$. The four model variables (indicated with hats) on the left-hand side of (2.3) are the Fourier coefficients of the Kelvin wave ($K$), Rossby wave ($R$), moisture ($Q$), and convective activity ($A$), which are all anomalies. Note that $\bar{A}_{k}$ represents the Fourier coefficients of $\bar{A}$, which is the background convective activity computed by

$$H\bar{A}=E_{q}+S^{q}-\bar{Q}S^{\theta}/(1-\bar{Q}),$$

(2.4)

where $S^{\theta}$ and $S^{q}$ are the external source of cooling and moistening, respectively. Both $S^{\theta}$ and $S^{q}$ are pre-determined, peaking in the warm pool region. In (2.4), $E_{q}$ is the latent heat that is affected by the SST. In this context, the latent heat $E_{q}$ is approximated by a function that is proportional to the SST in the CP, serving as a bulk representation,

$$E_{q}=\alpha_{q}T_{C}.$$

(2.5)

As the SST increases, it is expected that convection and moisture will also increase, which modulates the atmosphere. The other two factors $\overline{H}$ and $\overline{Q}$ are also constants (Thual et al., 2014). A linear combination of $\hat{K}_{k}$, $\hat{R}_{k}$, $\hat{Q}_{k}$, and $\hat{A}_{k}$ can also be utilized to reconstruct the wind and the MJO, where a temporal filter within the band 30 to 90 days is further carried out for reconstructing the MJO. Please refer to Appendix 6 for further details.

The skeleton model for the MJO (Majda and Stechmann, 2009) captures several key observational features: (i) the eastward propagation speed of 5 m/s, (ii) a dispersion rate of ${\,\rm d}\omega/{\,\rm d}k\approx 0$, and (iii) a horizontal quadrupole structure. The stochastic version (Thual et al., 2014) further simulates (iv) the intermittent generation of MJO events and (v) the organization of MJO events into wave trains that exhibit growth and dissipation. These properties are retained in the low-order representation provided in (2.3).

Unlike the box model used for the ENSO component in (2.1), the Fourier representation in (2.3) enables a more accurate depiction of the eastward propagation of the MJO’s large-scale features. Since the modes $k=-1,-2$, and $-3$ are complex conjugates of those $k=1,2$, and $3$, only 16 equations (four for each wavenumber, together with the mode $k=0$) are needed for characterizing the intraseasonal atmospheric component.

The original stochastic skeleton model employed a stochastic birth-death process (Gardiner, 2009, Thual et al., 2014) to capture the irregularity of MJO events. It has been demonstrated that, in the limit of infinitesimal jumps, a continuous stochastic differential equation driven by white noise $\dot{W}_{A}$ with a state-dependent coefficient can be derived (Chen and Majda, 2016). This formulation facilitates the spectral decomposition of the model and is the one utilized here. Although the system (2.3) is expressed in Fourier space, the state-dependent noises are easier to interpret when the system is presented in physical space. In physical space, the state-dependent noises in the equations for convective activity and moisture take the following form:

	$\displaystyle\sigma_{A}$	$\displaystyle=\sqrt{\nu|Q|(A+\bar{A})},$		(2.6a)
	$\displaystyle\sigma_{Q}$	$\displaystyle=\max\left(\tilde{\sigma}_{Q}\left(1-e^{-c_{q}E_{q}}\right),{\sigma_{Q}^{\min}}\right),$		(2.6b)

where $\nu$, $\tilde{\sigma}_{Q}$, $c_{q}$, and $\sigma_{Q}^{\min}$ are all positive constants. The noise coefficient $\sigma_{A}$ in (2.6a), derived from the limit of infinitesimal jumps in the original stochastic skeleton model (Chen and Majda, 2016), indicates that the strength of the convective activity is dependent on the moisture level. This relationship is consistent with the deterministic part of the model, $\Gamma Q(A+\bar{A})$, where moisture influences the tendency of the SST temporal variation. On the other hand, the noise coefficient $\sigma_{Q}$ links the ocean and the atmosphere via the latent heat variable $E_{q}$. Higher SST causes an increase in the latent heat, which strengthens the statistical coupling between the atmosphere and ocean. A minimum noise in the moisture process, denoted as $\sigma_{Q}^{\min}$, is prescribed to ensure that the noise coefficient remains non-negative. Since only the leading three Fourier modes are retained, the stochastic noises not only capture the explainable physics but also compensate for contributions from smaller scales.

The total wind $\tau$ appearing in (2.1) has two components:

$$\tau={u}_{W}+{\bar{u}_{M}},$$

(2.7)

where $u_{W}$ represents the intraseasonal wind $u$ reconstructed by $K$ and $R$ averaged over the western Pacific, the region where wind activity is most pronounced. The second component of the total wind is the interannual component, $\bar{u}_{M}$, which requires a separate governing equation.

The interannual wind averaged over the western Pacific is driven by a simple process:

$$\frac{{\,\rm d}\bar{u}_{M}}{{\,\rm d}t}=-\lambda_{M}\bar{u}_{M}+{\alpha_{3}E_{q}}+\sigma_{M}\dot{W}_{M}.$$

(2.8)

In this equation, the latent heat $E_{q}$ from (2.5) acts as a forcing for the wind velocity. Together with the wind forcing in the interannual ocean model (2.1), this creates a feedback mechanism for air-sea interaction. The additional damping and stochastic terms introduce variability to the interannual wind that is not solely determined by latent heat. It is important to note that the fluctuations in interannual wind resulting from this noise are significantly weaker than those from the intraseasonal atmospheric model, allowing the latter to remain the dominant contributor to intraseasonal variability. See Figure 6.6 in Appendix 6. Nevertheless, the noise in (2.8) enhances the irregularity of the model simulation. It also establishes a more realistic dependence between latent heat and interannual wind, which prevents a purely deterministic relationship.

In the conceptual model developed here, ENSO is characterized by two SST indices, $T_{E}$ and $T_{C}$, as described in (2.1). Similarly, the MJO is represented by the leading three Fourier modes in (2.3). Although both phenomena are represented in a low-dimensional form using time series, a spatiotemporal reconstruction based on these time series enables us to visualize the spatiotemporal patterns and spatial propagation of the observed features. The details and validation of the spatiotemporal reconstructions can be found in Appendix 6.

The monthly SST data used in this study is sourced from the GODAS dataset (Behringer and Xue, 2004). The SST anomalies are calculated by subtracting the monthly mean climatology over the entire analysis period. The Niño 4 and Niño 3 indices represent the average SST anomalies over the zonal regions $160^{\circ}$E–$150^{\circ}$W and $150^{\circ}$W–$90^{\circ}$W, respectively, with a meridional average between $5^{\circ}$S and $5^{\circ}$N. The reanalysis period spans from 1982 to 2019, corresponding to the satellite era when observations are more reliable. This period is used to study ENSO-MJO spatiotemporal patterns. For statistical comparisons, a longer SST dataset (1950–2019) is used, obtained from the Extended Reconstructed Sea Surface Temperature version 5 (Huang et al., 2017).

The classification of different El Niño and La Niña events, which helps examine ENSO variability, is based on the criteria outlined in (Kug et al., 2009, Wang et al., 2019). These classifications rely on the average SST anomalies during boreal winter (December–January–February; DJF). An event is classified as an EP El Niño if the EP is warmer than the CP, and the EP SST anomaly exceeds $0.5^{o}$C. An extreme El Niño occurs when the maximum EP SST anomaly between April and the following March exceeds $2.5^{o}$C. A CP El Niño is identified when the CP is warmer than the EP, with a CP SST anomaly above $0.5^{o}$C. Lastly, a La Niña event is defined when the SST anomaly in either the CP or EP falls below $-0.5^{o}$C.

To represent convective activity ($a$), we use the National Oceanic and Atmospheric Administration’s (NOAA) interpolated outgoing longwave radiation (OLR) dataset (Liebmann and Smith, 1996). While several proxies for convective activity exist, OLR is chosen here as an initial, straightforward option for analysis. For other variables such as zonal wind, geopotential height, and specific humidity, we rely on the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis data (Kalnay et al., 1996). Both datasets offer a horizontal spatial resolution of 2.5${}^{o}\times$ 2.5^o and are available at daily temporal intervals. The period used for this study spans from 1982 to 2019 (Stechmann and Majda, 2015).

Figure 4.1 presents the model simulations and associated statistics compared with observations. Panel (a) demonstrates that the simulated SST time series qualitatively match the observations. Notably, several large positive values are observed in $T_{E}$ (e.g., at $t=141$ and $145$), corresponding to strong EP El Niño events. Similarly, there are periods (e.g., at $t=129$ and $133$) where both $T_{E}$ and $T_{C}$ exceed 0.5^oC, with $T_{C}>T_{E}$, indicating CP events. Panels (b)–(d) compare key statistical measures. First, the model simulations exhibit power spectra very similar to the observations, reflecting the average oscillation period and its irregularity. Additionally, the model accurately reproduces the non-Gaussian probability density functions (PDFs), capturing the positive skewness and one-sided fat tail of the Niño 3 SST distribution. This allows the model to simulate extreme EP El Niño events and the El Niño-La Niña asymmetry. Likewise, the model captures the negative skewness of the Niño 4 SST distribution, preventing extreme events from occurring in the CP region. Furthermore, the model successfully reproduces the seasonal variation of SST variance, indicating a higher likelihood of event occurrence during boreal winter. Panel (e) shows the bivariate distribution of DJF SST peaks. Similar to observations, the strength of the SST peaks increases as events shift eastward in the Pacific, where CP events generally have weaker amplitudes than EP events in the model. Finally, Panel (f) displays the frequency of different ENSO events, which is a useful indicator of ENSO complexity. Overall, the model succeeds in reproducing the number of different ENSO events, consistent with nature. Appendix 6 contains the statistics of other variables, the role of the decadal variability $I$, and the reconstructed spatiotemporal fields. Particularly, the wind statistics also closely match the observations, capturing key features like positive skewness and a one-sided fat tail, which correspond to the westerly wind bursts, an important precursor for El Niño events.

The atmospheric zonal velocity ($u$) and convective activity ($A$) are important variables for characterizing and reconstructing the MJO. The power spectra of these two variables are natural indicators to statistically describe the dominant spatiotemporal signal. The model simulation shows a high spectral density within the intraseasonal band (30 to 90 days) for modes $k=1,2$, and $3$, indicating the dominant eastward-propagating MJO signal. Likewise, the westward-propagating moisture Rossby waves are captured by modes $k=-1,-2$, and $-3$. The power spectra of the model within this intraseasonal band closely match observations. The results indicate that the model succeeds in capturing these crucial intraseasonal variabilities. See Figure 6.3 in the Appendix 6 for more details.

Figure 4.2 shows the Hovmoller diagrams of the spatially reconstructed SST and MJO from the conceptual model. These diagrams illustrate the diversity and complexity of ENSO, as well as the coupled relationship between ENSO and the MJO. A similar figure based on observational data can be found in Appendix 6.

First, the model simulates different types of ENSO events with frequencies similar to observations (see also Panel (f) in Figure 4.1). A few examples of each type, with corresponding years of observed events in parentheses, are as follows: (i) moderate EP El Niño (1983): $t=$149; (ii) extreme EP El Niño (1998): $t=$141; (iii) delayed super El Niño (2014-2015): $t=$118-119; (iv) multi-year EP El Niño (1987-1988): $t=$123-126; (v) CP El Niño (2010): $t=$137; (vi) multi-year CP El Niño (1993-1994): $t=$144-145; (vii) La Niña (2011): $t=$127; (viii) multi-year La Niña (1999-2000): $t=$116-117. Additional examples can be found in Figure 6.5.

Second, the model successfully reproduces strong interactions between the MJO and ENSO. Typically, intense MJO events are observed either preceding or during the El Niño phase, for both EP and CP El Niño events. The MJO signals extend further towards the EP, particularly during strong EP El Niño events. However, some MJO events are observed during La Niña phases (e.g., at $t=$144) or do not trigger ENSO events at all (e.g., at $t=$131). These findings are consistent with observations, as detailed in Appendix 6.

Figure 4.3 shows the lagged correlation between MJOI, an index representing the MJO averaged over the WP, or its strength $|$MJOI$|$, and an SST index, either $T_{C}$ or $T_{E}$. This correlation is calculated specifically for years when El Niño events occur, aiming to highlight the statistical dependence between ENSO and the MJO. The analysis explores the lagged correlation across all El Niño events, including both EP and CP types (Panels (a)–(d)), as well as conditioning separately on EP events (Panels (e)–(f)) and CP events (Panels (g)–(h)).

Overall, the lagged correlations from the model simulations align well with observations. First, neither the simulations nor the observations provide clear evidence that MJOI statistically leads or lags SST (Thual et al., 2018), likely due to the faster oscillation and occurrence of the MJO on much shorter timescales. Consequently, the randomness of different MJO events tends to cancel out any consistent correlation with ENSO. However, the amplitude of the MJOI ($|$MJOI$|$) does show a correlation with SST, with the MJO becoming more intense prior to or during El Niño phases. Notably, the lagged correlation with $T_{E}$ during all El Niño phases exhibits strong asymmetry, with a longer correlation observed before an El Niño event (Panel (a)). Such an extended correlation range reflects the process of MJO buildup and the triggering of El Niño events. After the El Niño peak, the MJO remains active for a few months until the SST anomaly dissipates. The model accurately captures this behavior, showing a similar lagged correlation pattern (Panel (b)). For $T_{C}$, the lagged correlation during all El Niño events resembles that with $T_{E}$ in observations, though the model produces a slightly more symmetric correlation band. When conditioned only on CP events, the lagged correlation with $T_{C}$ is more symmetric in both observations and model results (Panels (g)–(h)), with significant correlations spanning a shorter window, from -5 to 5 months, compared to conditioning on all El Niño events. Finally, Panels (e)–(f) depict the lagged correlation with $T_{E}$ conditioned on EP events. Although the model and observations have similar patterns, a stronger correlation is observed using the model outcomes. However, caution is needed when interpreting such a result since the number of EP events in the observational period is quite limited.

An important aspect to explore further is the role of the state-dependent noises, $\sigma_{\hat{Q}_{\mathbf{k}}}(E_{q})\dot{W}_{\hat{Q}_{k}}$ in (2.3c) and $\sigma_{\hat{A}_{\mathbf{k}}}(Q,A)\dot{W}_{\hat{A}_{k}}$ in (2.3d). These noises are designed to capture statistical feedback from latent heat and generate intermittent MJO events, respectively. Without random noise in the moisture equation ($\sigma_{\hat{Q}_{\mathbf{k}}}=0$), the MJO is less affected by SST, resulting in fewer strong wind events and a weakened correlation between ENSO and MJO. Consequently, the model simulation will underestimate extreme El Niño occurrences and reduce non-Gaussian features. Conversely, when $\sigma_{\hat{A}_{\mathbf{k}}}=0$, MJO events become more regular, and the relationship between MJO and ENSO becomes more deterministic, resulting in a much larger correlation compared to observations. See Figures 6.7–6.9 in Appendix 6.

The MJO skeleton model is a nonlinear oscillator model for the MJO skeleton as a neutrally stable wave (Majda and Stechmann, 2009, 2011). The core mechanism behind this oscillation is the interaction between (i) planetary-scale lower-tropospheric moisture anomalies, denoted by $q$, and (ii) subplanetary-scale convection and wave activity, represented by $a$. The planetary envelope, $a\geq 0$, characterizes the collective influence of unresolved subplanetary convection and wave activity. A critical aspect of the $q$-$a$ interaction lies in the assumption that moisture anomalies $q$ directly affect the rate of change of $a$, expressed by the relation $a_{t}=\Gamma qa$, where $\Gamma>0$ is a constant governing the strength of this interaction.

In the skeleton model, the interaction between moisture anomalies ($q$) and convection/wave activity ($a$) is further coupled with the linearized primitive equations projected onto the first vertical baroclinic mode. The non-dimensionalized system of equations is given as follows:

	$\displaystyle u_{t}-yv-\theta_{x}$	$\displaystyle=0,$		(6.1a)
	$\displaystyle yu-\theta_{y}$	$\displaystyle=0,$		(6.1b)
	$\displaystyle\theta_{t}-u_{x}-v_{y}$	$\displaystyle=\bar{H}a-s^{\theta},$		(6.1c)
	$\displaystyle q_{t}+\bar{Q}(u_{x}+v_{y})$	$\displaystyle=-\bar{H}a+s^{q},$		(6.1d)
	$\displaystyle a_{t}$	$\displaystyle=\Gamma qa,$		(6.1e)

with periodic boundary conditions along the equatorial belt and planetary-scale equatorial long-wave scaling. In the dry dynamics ((6.1a)–(6.1c)), $u$, $v$, and $\theta$ represent zonal velocity, meridional velocity, and potential temperature, respectively, while (6.1d) governs the evolution of low-level moisture ($q$). All variables are anomalies from a radiative-convective equilibrium, except for $a$. The skeleton model incorporates a minimal set of parameters: $\bar{Q}$ represents the background vertical moisture gradient, and $\Gamma$ is a proportionality constant for the $q$-$a$ interaction. While $\bar{H}$ is not dynamically influential, it defines a cooling/drying rate $\bar{H}a$ in dimensional terms. External cooling and moistening effects, $s^{\theta}$ and $s^{q}$, must be prescribed to complete the system.

Next, the system (6.1) is projected onto the first Hermite function in the meridional direction, where the fields are represented as $a(x,y,t)=\tilde{A}(x,t)\phi_{0},q=Q\phi_{0},s^{q}=S^{q}\phi_{0},s^{\theta}=S^{\theta}\phi_{0}$, where $\phi_{0}(y)=\sqrt{2}(4\pi)^{-1/4}\exp(-y^{2}/2)$. This choice of meridional heating structure excites only Kelvin waves and the first symmetric equatorial Rossby waves. The resulting meridionally truncated system is:

	$\displaystyle K_{t}+K_{x}$	$\displaystyle=(S^{\theta}-\bar{H}\tilde{A})/2,$		(6.2a)
	$\displaystyle R_{t}-R_{x}/3$	$\displaystyle=(S^{\theta}-\bar{H}\tilde{A})/3,$		(6.2b)
	$\displaystyle Q_{t}+\bar{Q}(K_{x}-R_{x}/3)$	$\displaystyle=(\bar{H}\tilde{A}-S^{q})(\bar{Q}/6-1),$		(6.2c)
	$\displaystyle\tilde{A}_{t}$	$\displaystyle=\Gamma Q\tilde{A},$		(6.2d)

provided that $S^{\theta}=S^{q}$. Note that $\tilde{A}$ represents the total convective activity, namely the summation of the anomaly $A$ and the background state $\bar{A}$. The dry dynamics components can be reconstructed as follows:

$$\begin{split}u&=(K-R)\phi_{0}+R\phi_{2}/\sqrt{2},\\ v&=(4\partial_{x}R-\bar{H}\tilde{A})\phi_{1}/3\sqrt{2},\\ \theta&=-(K+R)\phi_{0}-R\phi_{2}/\sqrt{2}.\end{split}$$

(6.3)

Here, the higher-order Hermite functions $\phi_{1}(y)=2y(4\pi)^{-1/4}\exp(-y^{2}/2)$ and $\phi_{2}(y)=(2y^{2}-1)(4\pi)^{-1/4}\exp(-y^{2}/2)$, though irrelevant to the dynamics, are necessary for retrieving the MJO’s quadruple structure (Majda and Stechmann, 2009).

The original stochastic skeleton model (Thual et al., 2014) introduces a simple stochastic parameterization to represent synoptic-scale processes. In this model, the amplitude equation (6.1e) or (6.2d) is replaced by a stochastic birth-death process, which allows for intermittent fluctuations in the synoptic activity envelope (Gardiner, 2009). Let $\tilde{A}$ be a random variable taking discrete values $\tilde{A}=\Delta\tilde{A}\eta$, where $\eta$ is a positive integer. The transition probabilities between states over a time step $\Delta t$ are given by:

$$\begin{split}&P\{\eta(t+\Delta t)=\eta(t)+1\}=\lambda\Delta t+o(\Delta t),\\ &P\{\eta(t+\Delta t)=\eta(t)-1\}=\mu\Delta t+o(\Delta t),\\ &P\{\eta(t+\Delta t)=\eta(t)\}=1-(\lambda+\mu)\Delta t+o(\Delta t),\\ &P\{\eta(t+\Delta t)\neq\eta(t)+1,\eta(t),\eta(t)+1\}=o(\Delta t),\end{split}$$

(6.4)

where $\lambda$ and $\mu$ are the birth and death rates defined as:

$$\lambda=\begin{cases}\Gamma|q|\eta+\delta_{\eta 0}\qquad&\mbox{if~{}}q\geq 0\\ \delta_{\eta 0}\qquad&\mbox{if~{}}q<0\end{cases}\qquad\mbox{and}\qquad\mu=\begin{cases}0\qquad&\mbox{if~{}}q\geq 0\\ \Gamma|q|\eta\qquad&\mbox{if~{}}q<0\end{cases}$$

(6.5)

with $\delta_{\eta 0}$ being the Kronecker delta operator. These transition rates ensure that $\partial_{t}E(\tilde{A})=\Gamma E(Q\tilde{A})$ for $\Delta\tilde{A}$ small, recovering the average $Q$-$\tilde{A}$ interaction described (6.1).

In (Chen and Majda, 2016), a continuous stochastic differential equation (SDE) for convective activity is derived for small $\Delta a$:

$$\frac{{\,\rm d}\tilde{A}}{{\,\rm d}t}=\Gamma Q\tilde{A}+\sqrt{\Delta\tilde{A}\Gamma|Q|\tilde{A}}\dot{W}_{A}.$$

(6.6)

where $\dot{W}_{A}$ represents a standard Wiener process, capturing the stochastic fluctuations in convective activity.

The leading three Fourier modes from the system (6.2a)–(6.2c), along with the stochastic equation (6.6), are employed to construct a conceptual model for the atmospheric intraseasonal component. It aims to capture the essential large-scale dynamics and provide a simplified yet effective representation of intraseasonal atmospheric processes.

The structure of the MJO can be characterized by the eigenvector associated with the linearized form of (6.2). Let $\mathbf{U}=(K,R,Q,A)^{\mathtt{T}}$ represent the collection of state variables in (6.2). By assuming a plane-wave ansatz for $\mathbf{U}$ and substituting it into (6.2), we derive an eigenvalue problem for each wavenumber $k$. This results in a four-dimensional system that features four eigenmodes: dry Kelvin, moisture Kelvin, dry Rossby, and MJO. Each of these waves can be expressed as a linear combination of $K,R,Q$ and $A$.

Let the eigenvalue corresponding to the MJO for wavenumber $k$ be denoted as follows:

$$\omega_{k}:=\omega_{\mbox{\tiny MJO}}(k).$$

The evolution of the MJO for wavenumber $k$ is then represented by the time series of the Fourier coefficients $\hat{K}_{k},\hat{R}_{k},\hat{Q}_{k}$, and $\hat{A}_{k}$ projected onto the eigenvector associated with $\omega_{k}$. The numerical values of the eigenvalues and eigenvectors pertaining to the MJO modes are presented in Table 1.

It is worth noting that the frequencies of the first three modes remain nearly constant. The final step in constructing the MJO signal involves summing the modes for $k=1,2$, and $k=3$, followed by the application of a temporal filter to retain only the signals within the intraseasonal band, specifically between 30 and 90 days.

The processing of observational data for the intraseasonal atmospheric model follows the procedure outlined in (Stechmann and Majda, 2015).

$$\theta:=\theta_{BC,1}=-Z_{BC,1}=-\frac{Z(850\mathrm{hPa})-Z(200\mathrm{hPa})}{2\sqrt{2}}.$$

(6.8)

The reference scales $\theta$ and $Z$ are approximately $\bar{\alpha}\approx 15.6\mathrm{~{}K}$ and $c^{2}/g\approx 265\mathrm{~{}m}$, respectively. To achieve non-dimensionalization, the geopotential height data $Z$ is divided by $c^{2}/g$ for non-dimensionalization.

Note that in the datasets for $q$, data at the 725 hPa level is unavailable; therefore, the data from the 700 hPa level is used instead. Additionally, $Q$ has been non-dimensionalized using the natural reference scale $L_{v}/c_{p}\tilde{\alpha}$, where $\tilde{\alpha}$ represents the reference potential temperature scale.

Using the meridional basis functions, the definitions of the Kelvin wave $K$ and the first symmetric equatorial Rossby wave $R$ are given by:

		$\displaystyle K=\frac{1}{2}\left(u_{0}-\theta_{0}\right)\text{ and }$		(6.11)
		$\displaystyle R=-\frac{1}{4}\left(u_{0}+\theta_{0}\right)+\frac{\sqrt{2}}{4}\left(u_{2}-\theta_{2}\right),$

where $u_{m}$ and $\theta_{m}$ represent the meridional projections.

The values and definitions of the constant parameters in the conceptual model are listed in Table 2. In addition to these, the model uses functions of time and state variables for some parameters. The wind coefficients,

	$\displaystyle\beta_{E}(I)$	$\displaystyle=\frac{\sqrt{0.6}\left(8-\frac{4I}{5}\right)}{3},$		(6.12a)
	$\displaystyle\beta_{C}(I)$	$\displaystyle=0.8\beta_{E}(I),$		(6.12b)
	$\displaystyle\beta_{u}(I)$	$\displaystyle=-0.4\beta_{E}(I),$		(6.12c)
	$\displaystyle\beta_{h}(I)$	$\displaystyle=-0.2\beta_{E}(I),$		(6.12d)

which appear in (2.1), are some of the variable parameters in the model. Such functions incorporate the decreased sensitivity to changes in wind during periods of strong Walker circulation.

Other parameters that are non-constant include the damping coefficients $c_{1}(T_{C},t)$ and $c_{2}(t)$ in (2.1). Both are sinusoidal functions of time in order to reproduce the observed seasonal phase locking of SST in the CP and EP. In addition, $c_{1}$ is a quadratic function of $T_{C}$ which allows the model to produce the negatively skewed PDF for CP SST anomalies. Specifically, they are expressed as

	$\displaystyle c_{1}(T_{C},t)$	$\displaystyle=0.78\left[25\left(T_{C}+0.1\right)^{2}+0.6\right]\left[1+0.6\sin\left(\frac{2\pi}{6}t\right)\right],$		(6.13a)
	$\displaystyle c_{2}(t)$	$\displaystyle=0.84\left[1+0.4\sin\left(\frac{2\pi}{6}t+\frac{2\pi}{6}\right)+0.4\sin\left(\frac{4\pi}{6}t+\frac{2\pi}{6}\right)\right].$		(6.13b)

The state dependent noise coefficient $\sigma_{I}(I)$ in (2.2) can be derived using the Fokker-Planck equation evaluated at the stationary solution of the PDF for I, $p(I)$ (Averina and Artemiev, 1988). The resulting noise coefficient is

$$\sigma_{I}(I)=Re\left(\sqrt{-8\lambda\int_{0}^{I}\lambda(y-\overline{I})dy}\right).$$

(6.14)

The model uses non-dimensional variables. The dimensional units of each variable are listed in Table 3.

Recall that the stochastic conceptual model includes two SST variables, $T_{C}$ and $T_{E}$, as defined in (2.1). These variables can be used to approximately reconstruct the entire SST field across the equatorial Pacific. The method employed here is a bivariate regression,

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

(6.15)

where $a_{E}(x)$ and $a_{C}(x)$ are the spatial basis functions that depends on the location $x$.

For each longitude $x^{*}$, the two scalar regression coefficients, $a_{E}(x^{*})$ and $a_{C}(x^{*})$, are computed using (6.15) based on observational SST anomaly data from 1980 to 2020. Repeating this process across all longitude grids yields the spatially dependent functions $a_{E}(x)$ and $a_{C}(x)$. These functions are displayed in Panel (a) of Figure 6.1. The regression coefficient function $a_{C}(x)$ ($a_{E}(x)$) peaks in the CP (EP) region, as the SST at those longitudes is closely correlated with $T_{C}$ ($T_{E}$).

Given that the spatiotemporal reconstruction based on this bivariate regression is applied to study the coupled ENSO-MJO phenomena, it is crucial to validate the accuracy of the reconstructed field. Panels (b) and (c) of Figure 6.1 compare the actual SST field with the reconstructed one, while Panel (d) presents their difference. Overall, the reconstructed field almost perfectly captures the exact SST patterns, especially in the EP and CP regions, with only minor biases in the WP and along the eastern boundary.

The reconstruction of the intraseasonal atmospheric variables, which can further describe the spatiotemporal patterns of the MJO, is given by the Fourier summation. For instance, the reconstruction of the moisture variable $Q$ in (2.3) in physical space is expressed as:

$$Q(x,t)=\sum_{k=\pm 1,\pm 2,\pm 3}\hat{Q}(t)e^{ikx}.$$

(6.16)

Figure 6.2 presents Hovmoller diagrams of ENSO and MJO patterns during the observational period from 1982 to 2018. Similar to Figure 4.2, the SST spans the equatorial Pacific, while the MJO also extends into the Indian Ocean. The red vertical line in the MJO panels marks the boundary of the Western Pacific (WP) at 120°E. In the SST panels, the averaged atmospheric wind over the WP is overlaid on the SST, with red and blue indicating westerly and easterly wind bursts, respectively. The black curve represents the interannual wind. This figure is used to qualitatively validate the coupled relationship between ENSO and MJO as illustrated in the stochastic conceptual model shown in Figure 4.2.

Figure 6.3 presents the power spectra of atmospheric zonal velocity ($u$) and convective activity ($A$), two key variables for reconstructing the MJO. The x-axis spans wavenumbers from $k=-3$ to $k=3$, which are the wavenumbers used in the model. The y-axis represents frequency in cycles per day (cpd). The focus is on the intraseasonal band, highlighted between the two horizontal dashed lines. Black dots mark the dispersion curves from linear analysis of the MJO skeleton model. The high density within this band for modes $k=1,2$, and $3$ indicates the dominant eastward-propagating MJO signal. Conversely, the westward-propagating moisture Rossby waves are captured by modes $k=-1,-2$, and $-3$. The power spectra of the model within the intraseasonal band resemble observations. On the other hand, the model underestimates the power density at lower frequencies (where cpd approaches zero). This is because the observations contain information across all temporal scales, whereas this model is primarily designed to capture intraseasonal variabilities. Further, the model density patterns at higher frequencies for modes $k=1,2$, and $3$ are pretty consistent with observations. The modes with negative wavenumbers in the model have stronger spectral density, possibly due to the stochastic noise. Since the model only focuses on coupling MJO and ENSO, matching the higher frequency signals does not impact the model’s function. Note that the primary mechanism through which SST drives the intraseasonal model is the state-dependent noise coefficient on moisture, representing statistical feedback that does not directly affect the mean state of intraseasonal variability.

Figure 6.4 presents the standard deviation of MJO as a function of longitude conditioned on CP and EP El Niño events. The longitudes span from the prime meridian to $70^{\circ}W$. The warm pool region, which includes the Indian Ocean and the WP, is highlighted by dashed lines. The high MJO activity in this region indicates that the warm SSTs fuel the MJO. Additionally, the standard deviation of MJO during EP El Niño events is higher than that of CP El Niño events in the eastern Pacific while the opposite is true in the warm pool region for both the model simulation and observations. Note that the amplitude of the MJO signal in the model is slightly weaker than that of observations since the model is coarse grained and thus averages the high frequency larger amplitudes.

Figure 6.5 presents additional model results. While the main text focuses on the two SST variables, this figure includes time series and statistics for other variables in the stochastic conceptual model. The thermocline in the WP, $h_{W}$, exhibits statistics comparable to observations. The wind statistics closely match the observations, capturing key features like positive skewness and a one-sided fat tail, which correspond to the westerly wind bursts, an important precursor for El Niño events. Accompanying the time series, Hovmoller diagrams of SST and MJO during the same period are also shown. Notably, the decadal variability $I$ is more pronounced after $t=510$, resulting in a stronger zonal ocean advection effect (see (2.1c)). Consequently, more CP El Niño events are observed during this period compared to the prior half of the simulation. In contrast, when $I$ approaches zero before $t=510$, several strong EP El Niño events occur, along with more active MJO events.

Figure 6.6 displays $T_{E}$ alongside the two wind components in the WP: the interannual wind ($u_{W}$) and the intraseasonal wind ($\bar{u}_{M}$). The interannual wind is closely synchronized with $T_{E}$, showing only minor fluctuations due to randomness. However, the amplitude of this randomness is much smaller than that of the intraseasonal wind, which is linked to the MJO. As a result, variability in the interannual wind does not significantly overshadow the signals from the intraseasonal model.

Figure 6.7 compares three model configurations: the full model (Panel (a)), the model without state-dependent noise in the convective activity ($\sigma_{\hat{A}_{k}}=0$ in (2.3d); Panel (b)), and the model without state-dependent noise in the moisture process ($\sigma_{\hat{Q}_{k}}=0$ in (2.3c); Panel (c)). When state-dependent noise is removed from the convective activity process (2.3d) (Panel (b)), the intermittent behavior of the MJO vanishes, and the occurrences of MJO and ENSO become more synchronized, diverging from natural behavior. This also introduces a significant bias in the lagged regression between MJO and ENSO, as seen in Figure (6.8). In contrast, when state-dependent noise is removed from the moisture process (2.3c) (Panel (c)), the intermittency is restored. However, in the absence of SST feedback, which is an external driving force for the MJO, strong MJO events are severely underestimated. Consequently, extreme El Niño events are harder to trigger, resulting in biases in reproducing the non-Gaussian, fat-tailed statistics of ENSO. Additionally, the lack of SST feedback drastically weakens the correlation between SST and the MJO as seen in Figure 6.9.

Figure 6.8 presents the lagged correlation between MJO and ENSO when the state-dependent noise in convective activity is removed ($\sigma_{\hat{A}_{\mathbf{k}}}=0$), based on the MJOI and SST time series. Compared with the full model results, as shown in Figure 4.3, the lagged correlation without this noise is significantly higher than observed. As a result, the relationship between MJO and ENSO becomes much more deterministic, deviating from the natural system. The findings here highlight the crucial role of state-dependent noise in convective activity, which is essential for generating the many intermittent MJO events seen in nature.

Similarly, figure 6.9 depicts the lagged correlation between MJOI and SST time series for the model with no noise in the low-level moisture equation. In comparison with observations, and the full model results (Figure 4.3), without this noise the link between MJO and SST is severely weakened. Removing this state dependent noise breaks the feedback loop between the atmospheric and oceanic components. Ultimately wind bursts are significantly reduced thus inhibiting the formation of extreme El Niño events and in turn reducing MJO activity.