An improved noise model for representing westerly wind bursts in the recharge oscillator model of ENSO

Westerly wind bursts (WWBs) have long been known to have a major impact on the development of El Niño events. In particular, they amplify these events, with stronger events associated with a higher number and stronger WWBs. We consider here a noise-driven recharge oscillator model of ENSO. Commonly, WWBs are represented by a state-dependent Gaussian noise that naturally reproduces the amplification of warm events. However, we show that many properties of WWBs and their effects on sea surface temperature (SST) are better captured by a conditional additive and multiplicative (CAM) noise, which presents a promising alternative to represent WWBs. In addition to recovering the sporadic nature of WWBs, CAM noise leads to an asymmetry between El Niño and La Niña events without the need for deterministic nonlinearities. Furthermore, CAM noise generates SST dynamics with a higher frequency of WWBs accompanying the largest events. This suggests that extreme warm events are better modelled by CAM noise. To cover the full spectrum of warm events, we propose a conditional noise model in which the wind stress is modelled by additive Gaussian noise for sufficiently small SSTs and by additive CAM noise once the SST exceeds a certain threshold. We show that this conditional noise model captures observed bulk statistical properties of ENSO equally well as the commonly used multiplicative Gaussian red noise model, but additionally better reproduces dynamical signatures such as the increased number of WWBs preceding large El Niño events.

1 Introduction

Westerly wind bursts (WWBs) in the equatorial Pacific last a week or two, have a longitudinal scale of a thousand km, and have preceded and amplified every major El Niño event in the observed record (Harrison and Vecchi, 1997; Lengaigne et al., 2004; Yu and Fedorov, 2022). WWBs trigger ocean Kelvin waves that accelerate the East Pacific warming (Kessler et al., 1995). While these are weather events, they are not completely random, and tend to occur when the equatorial SST begins to warm (Tziperman and Yu, 2007), therefore amplifying developing El Ni no events (Eisenman et al., 2005; Gebbie et al., 2007). Further, WWBs occur more frequently during the active phase of the MJO (Chiodi et al., 2014; Liang and Fedorov, 2021a). Because of this implied state-dependency, WWB events are commonly represented as a multiplicative Gaussian noise term in simplified ENSO models studying their role in El Niño events (e.g., Perez et al., 2005; Vialard et al., 2025). In this formulation, the red noise amplitude is multiplied by a Heaviside function of the temperature and by the East Pacific temperature itself, restricting the events to warmer than normal temperatures, and making sure the stochastic noise amplitude increases with the SST anomaly, as motivated by observational analysis (Tziperman and Yu, 2007). This allows for an excellent fit to the observed ENSO SST characteristics, including the spectrum and PDF (Vialard et al., 2025; Han et al., 2026).

Figure 1 WWBs from a climate model and NINO3 from observations. (a) An inferred measure of the time-integrated wind stress due to westerly wind bursts in a global climate model (CESM2, Danabasoglu et al., 2020). Shown is $\mathrm{\Delta }{t}_{\mathrm{WWB}}{U}_{\mathrm{WWB}}^{\mathrm{2}}$, where the WWB duration is given by Δt_WWB (in days) and the surface wind speed strength is given by U_WWB (in m s⁻¹). (b) Observed NINO3 index from 1870 until 2026 (Rayner et al., 2003). A moving average over 4 months was applied to observations. The horizontal line demarcates an index of 1.5, which we use to separate large El Niño events from normal ones. The vertical lines denote the years 1997 and 2015, for which the large El Niño events were accompanied by strong WWBs.

In this work, we examine an alternative to this stochastic forcing formulation, which provides an improved representation of the effects of WWBs on ENSO. We show that while Gaussian noise reproduces the global statistical features of ENSO, an alternative formulation is required to better represent the character of WWBs. To motivate our proposed noise model, consider Fig. 1a, which shows an inferred measure of the time-integrated wind stress associated with WWBs from a numerical simulation of a global configuration of the Community Atmospheric Model (Conley et al., 2012). Lian et al. (2018) show that even an earlier version of the Community Atmospheric Model simulates WWBs quite realistically. Specifically, we show $\mathrm{\Delta }{t}_{\mathrm{WWB}}{U}_{\mathrm{WWB}}^{\mathrm{2}}$, where Δt_WWB is the WWB duration and U_WWB is the surface wind speed strength. This measure attempts to represent the effect of these events on the ocean, which depends on both their amplitude and their duration. The observed time series in Fig. 1a with its sporadic high-amplitude peaks resembles that of so-called correlated additive and multiplicative (CAM) noise, which has been used to study non-Gaussian oceanic SST anomalies and atmospheric dynamics (Sura and Sardeshmukh, 2008; Sardeshmukh and Sura, 2009; Penland and Sardeshmukh, 2012; Sardeshmukh and Penland, 2015). Numerical simulations suggest that resolving the sporadic peaky nature of WWBs is crucial in modelling ENSO (Yu and Fedorov, 2020). This suggests replacing the commonly employed Gaussian stochastic forcing used to represent WWBs in idealized ENSO models with non-Gaussian CAM noise. As we will show, such CAM noise also naturally reproduces the observed asymmetry between El Niño and La Niña events, without the inclusion of any deterministic nonlinearities designed to promote such asymmetry. Using CAM noise to explain the observed asymmetry in amplitude and persistence of El Niño and La Niña events has also been suggested by Martinez-Villalobos et al. (2019), who showed that the asymmetry can be captured by a multivariate linear model driven by CAM noise. Moreover, we will show that CAM noise can generate characteristic dynamical signatures of major El Niño events. In particular, major El Niño events are the result of a sustained driving by several large amplitude WWBs (Chiodi et al., 2014; Liang and Fedorov, 2021a). We will show that a recharge oscillator (RO) model driven by CAM noise well reproduces this dynamical feature, whereas standard multiplicative Gaussian noise fails to capture it.

We wish to make the case that CAM noise is a more appropriate noise model to model WWBs that are known to occur with a developing El Niño warming, further amplifying large El Niño events (Yu et al., 2003; Tziperman and Yu, 2007; Liang and Fedorov, 2021b). Stochastic wind forcing over cooler oceanic environments, on the other hand, may be more appropriately captured by Gaussian noise. We therefore propose a conditional noise model, that we will coin CON, in which the noise is Gaussian unless the SST exceeds a certain threshold when the noise switches to CAM noise. This conditional noise model will be shown to well reproduce the observed signatures of large El Niño events, while preserving the overall observed statistical properties of ENSO.

The paper is organized as follows. In Sect. 2.1 we introduce the recharge oscillator model (Jin, 1997; Burgers et al., 2005; Vialard et al., 2025), including a deterministic nonlinearity representing an asymmetric response of SST to changes in the thermocline. Section 2.2 introduces the stochastic forcing. We discuss both additive and multiplicative noise and introduce a coloured Gaussian noise model as well as CAM noise. We demonstrate that coloured Gaussian noise gives rise to Brownian motion when integrated, whereas the CAM noise gives rise to Lévy noise with abrupt jumps, which are more similar to WWBs and have the potential to trigger large El Niño events. Section 3 numerically explores the effect of three stochastic forcing scenarios for representing WWBs in an RO model: (1) a standard Gaussian noise used as multiplicative noise based on the East Pacific SST, representing what is the current consensus in the ENSO literature (Perez et al., 2005; Vialard et al., 2025); (2) a CAM noise, and finally, (3) a conditional noise that alternates between Gaussian noise for cold SST anomalies and CAM noise for positive anomalies. We argue that the conditional CAM noise may better capture the dynamics of large El Niño events. We conclude in Sect. 4 with a discussion of our results.

2 Methods

2.1 The recharge oscillator model

We consider the recharge oscillator (RO) model (Jin, 1997; Burgers et al., 2005; Jin et al., 2007; Vialard et al., 2025) for the Western Pacific thermocline depth anomaly h_w and the Eastern Pacific sea-surface temperature (SST) anomaly T_e,

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\dot{T}}_{\mathrm{e}}=-r\left(T_{\mathrm{e}}-\mathit{\gamma }{h}_{\mathrm{e}}\left(t\right)\right)\text{(2)}& {\displaystyle }& {\displaystyle }{\dot{h}}_{\mathrm{w}}=-\mathit{\epsilon }\left(h_{\mathrm{w}}+a\mathit{\tau }\left(t\right)\right),\end{array}$$

where the Eastern Pacific thermocline depth is expressed as

$$\begin{array}{}\text{(3)}& h_{\mathrm{e}}=h_{\mathrm{w}}+d\mathit{\tau },\end{array}$$

with non-dimensional wind stress

$$\begin{array}{}\text{(4)}& \mathit{\tau }=b{T}_{\mathrm{e}}+\mathit{\nu }\mathit{\xi }\left(t\right).\end{array}$$

We set d=1 m${}^{\mathrm{3}}/N$ throughout. The parameters r and ε control the characteristic decay time scales for the SST and for the thermocline depth, respectively, and a, b are positive constants. The noise ξ(t) with amplitude ν>0 represents the effect of westerly wind bursts (Vialard et al., 2025). For westerly wind anomalies, ξ>0, the noise acts to reduce the West Pacific thermocline depth h_w and amplifies the East Pacific SST T_e, and vice versa for ξ<0. The RO model can be concisely written as

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\dot{T}}_{\mathrm{e}}=R{T}_{\mathrm{e}}+F_{\mathrm{1}}{h}_{\mathrm{w}}+{\mathit{\sigma }}_{\mathrm{T}}\mathit{\xi }\left(t\right)\text{(6)}& {\displaystyle }& {\displaystyle }{\dot{h}}_{\mathrm{w}}=-F_{\mathrm{2}}{T}_{\mathrm{e}}-\mathit{\epsilon }{h}_{\mathrm{w}}-{\mathit{\sigma }}_{\mathrm{h}}\mathit{\xi }\left(t\right),\end{array}$$

where $R=r(\mathit{\gamma }bd-\mathrm{1})$ represents the Bjerknes feedback and F₁=rγ and F₂=abε are constants, and σ_T=F₁dν and ${\mathit{\sigma }}_{\mathrm{h}}=F_{\mathrm{1}}{F}_{\mathrm{2}}d\mathit{\nu }/(R+r)$ denote the amplitudes of the noise (Vialard et al., 2025). Without stochastic driving (i.e., ν=0) the dynamics of Eqs. (1)–(2) is that of a linear oscillator with frequency ${\mathit{\omega }}^{\mathrm{2}}=F_{\mathrm{1}}{F}_{\mathrm{2}}-(\mathit{\epsilon }+R{)}^{\mathrm{2}}/\mathrm{4}$ and a growth/damping rate $\mathit{\lambda }=(\mathit{\epsilon }-R)/\mathrm{2}$.

We choose parameters such that one time unit in Eqs. (5)–(6) corresponds to 1 month. We remark that an oversimplified fit of the RO model to match observations or climate model simulation of ENSO leads to a biased set of parameters that shows ENSO to be damped (Weeks and Tziperman, 2025) even if it is self-sustained (positive growth rate), although that should not affect our focus here.

To account for the observed asymmetry that El Niño events are typically of a larger magnitude than La Niña events, the RO model can be extended to include nonlinear terms in the temperature equation (see Vialard et al., 2025, for a recent review). Here, we use in one of our three numerical experiments below a nonlinear response of the SST to the thermocline depth, motivated by the parameterization of subsurface temperature in the CZ model (Zebiak and Cane, 1987). This involves employing a nonconstant response γ of the SST to changes in the thermocline with

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{\gamma }=\left\{\begin{array}{cl}{\mathit{\gamma }}_{+}& \mathrm{if}{h}_{\mathrm{e}}\ge \mathrm{0}\\ {\mathit{\gamma }}_{-}& \mathrm{if}{h}_{\mathrm{e}}<\mathrm{0},\end{array}\right..\end{array}$$

With ${\mathit{\gamma }}_{+}>{\mathit{\gamma }}_{-}$, Eq. (1) implies a stronger effect for a positive thermocline depth anomaly on the SST than for a negative thermocline depth anomaly. This enhances El Niño amplitudes relative to those of La Niña events, consistent with observations. Alternatively, one could add quadratic nonlinearities such as βT², with β>0, in the temperature equation (Eq. 1) to represent physical processes that favor the growth of El Niño relative to La Niña (see, e.g., Vialard et al., 2025; Liu et al., 2024, for a systematic investigation of the impact of different nonlinear terms). Besides such deterministic nonlinearities, state-dependent stochastic drivers (Jin et al., 2007; Levine et al., 2016; Levine and Jin, 2017) and non-Gaussian noise (Bianucci, 2016; Bianucci et al., 2018) were proposed as further dynamic ingredients to account for the observed asymmetry. In the following, we introduce several prototypical noise models. In Sect. 3 we will discuss their effect on the dynamics of ENSO, including its asymmetry, within the RO model, and then combine them to generate a conditional noise model that exhibits more realistic signatures consistent with observations.

2.2 Noise models

We consider both additive and multiplicative noise with

$$\begin{array}{}\text{(8)}& \mathit{\nu }={\mathit{\nu }}_{\mathrm{0}}+{\mathit{\nu }}_{\mathrm{1}}\max \left(\mathrm{0},T_{\mathrm{e}}\right).\end{array}$$

This form of the multiplicative noise, with ν₁≠0, takes into account that atmospheric noise, such as WWBs, occurs more frequently and with higher amplitude over a warmer equatorial ocean. Such multiplicative noise naturally introduces an asymmetry with larger El Niño events compared to La Niña events (Jin et al., 2007; Levine et al., 2016; Han et al., 2026).

Additive noise with ν₁=0 cannot generate the desired ENSO asymmetry for a Gaussian driving noise such as coloured Ornstein–Uhlenbeck (OU) processes. However, as we will show, certain non-Gaussian noise models, such as correlated additive and multiplicative (CAM) noise, can generate the asymmetry even when only applied additively. CAM noise is defined as the stochastic process governed by (e.g., Thompson et al., 2017),

$$\begin{array}{}\text{(9)}& \mathrm{d}\mathit{\xi }=c_{\mathrm{1}}\mathit{\xi }\mathrm{d}t+\left(c_{\mathrm{2}}\mathit{\xi }+c_{\mathrm{3}}\right)\circ \mathrm{d}{W}_{\mathrm{1}}+c_{\mathrm{4}}\mathrm{d}{W}_{\mathrm{2}}.\end{array}$$

Here, W_1,2 are independent Brownian motions (Wiener noise), with independent normally distributed random increments dW_1,2, and the ∘-symbol denotes Stratonovich noise. This noise model naturally appears when modelling the effect of fast dynamic processes onto slower ones and has found numerous applications in atmospheric and climate dynamics (Sura and Sardeshmukh, 2008; Sardeshmukh and Sura, 2009; Majda et al., 2009; Penland and Sardeshmukh, 2012; Sardeshmukh and Penland, 2015; Gottwald et al., 2017; Gottwald, 2021) and in particular, in the context of ENSO (Martinez-Villalobos et al., 2019).

The noise model (Eq. 9) contains both a purely Gaussian noise and a non-Gaussian noise as particular limits, depending on its parameters. For $c_{\mathrm{2}}=c_{\mathrm{3}}=\mathrm{0}$ and c₁<0, we have

$$\begin{array}{}\text{(10)}& \mathrm{d}\mathit{\xi }=-\left|c_{\mathrm{1}}\right|\mathit{\xi }\mathrm{d}t+c_{\mathrm{4}}\mathrm{d}{W}_{\mathrm{2}},\end{array}$$

and the CAM noise reduces to an Ornstein–Uhlenbeck noise with an asymptotic Gaussian distribution and zero mean (Pavliotis and Stuart, 2008). The central limit theorem ensures that if a Gaussian process is integrated in time, the resulting random variables are also distributed according to a Gaussian distribution.

However, for c₁<0 and c₂≠0 CAM noise is non-Gaussian and lies in the domain of attraction of α-stable processes (Kuske and Keller, 2001). This means that if such a process is integrated in time, it generates random variables that are drawn from an α-stable distribution. Such α-stable processes $L_{\mathit{\alpha },\mathit{\beta },\mathit{\eta }}$ are characterized by discrete jumps similarly to Lévy processes, and are parametrized by three parameters α, β and η. The stability parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{2}]$ determines the occurrence and size of the jumps. For α=2 we obtain a continuous Gaussian process without any discrete jumps. For α<2, however, the variance of such a process is not defined, as discrete jumps of arbitrary size have non-vanishing probability, and for α<1 even the mean ceases to exist. The skewness parameter $\mathit{\beta }\in [-\mathrm{1},\mathrm{1}]$ controls the direction of the jumps with β=1 allowing for only positive jumps, $\mathit{\beta }=-\mathrm{1}$ allowing for only negative jumps, β=0 allowing, on average, for as many positive as negative jumps, and values of β in between quantifying the probability of having positive or negative jumps. The scale parameter η reduces to the variance for the Gaussian case with α=2. For more details on α-stable processes, we refer the reader to Applebaum (2009) and Chechkin et al. (2008).

For c₄≠0 the mean of ξ is well-defined, and one has explicit expressions for the parameters of the resulting Lévy process α, β and η as functions of the parameters of the CAM process (Kuske and Keller, 2001; Thompson et al., 2017). For integrated CAM noise, the stability parameter α of the resulting α-stable process $L_{\mathit{\alpha },\mathit{\eta },\mathit{\beta }}$ is given by

$$\begin{array}{}\text{(11)}& \mathit{\alpha }=-\mathrm{2}{c}_{\mathrm{1}}/c_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

the skewness parameter is given by,

$$\begin{array}{}\text{(12)}& \mathit{\beta }=\tanh \left({\displaystyle \frac{\mathit{\pi }{c}_{\mathrm{3}}(\mathit{\alpha }-\mathrm{1})}{\mathrm{2}{c}_{\mathrm{4}}}}\right),\end{array}$$

and the scale parameter η is given by

$$\mathit{\eta }={\left({\displaystyle \frac{\mathrm{2}\cosh \left(\frac{\mathit{\pi }{c}_{\mathrm{3}}(\mathit{\alpha }-\mathrm{1})}{\mathrm{2}{c}_{\mathrm{4}}}\right)}{c_{\mathrm{2}}^{\mathit{\alpha }+\mathrm{1}}\mathit{\alpha }N}}\mathrm{\Gamma }(\mathrm{1}-\mathit{\alpha })\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\mathit{\alpha }\right)\right)}^{\frac{\mathrm{1}}{\mathit{\alpha }}},$$

with

$$N=\mathrm{2}\mathit{\pi }(\mathrm{2}{c}_{\mathrm{4}}{)}^{-\mathit{\alpha }}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathit{\alpha }\right)}{c_{\mathrm{2}}\mathrm{\Gamma }\left(z\right)\mathrm{\Gamma }\left(\stackrel{\mathrm{‾}}{z}\right)}},$$

where the bar denotes the complex conjugate and

$$z={\displaystyle \frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}+i{\displaystyle \frac{c_{\mathrm{3}}(\mathit{\alpha }-\mathrm{1})}{c_{\mathrm{4}}}}.$$

Figure 2 shows examples of CAM processes ξ(t) and their integrals $\mathrm{\Xi }={\int }^t\mathit{\xi }\left(s\right)\mathrm{d}s$. Figure 2a depicts a Gaussian OU process with $c_{\mathrm{1}}=-\mathrm{2}/\mathrm{3}$, c₄=0.7, $c_{\mathrm{2}}=c_{\mathrm{3}}=\mathrm{0}$, which when integrated yields Brownian motion as shown in Fig. 2c, in accordance with the central limit theorem. Figure 2b shows non-Gaussian CAM noise with intermittent unbounded peaks with c₁=1.22, c₂=1.14, c₃=0.65 and c₄=0.8, which, when integrated, leads to α-stable noise with jumps, as shown in Fig. 2d. In the following, we will call noise obtained from Eq. (9) OU noise if $c_{\mathrm{2}}=c_{\mathrm{3}}=\mathrm{0}$ and c₁<0. We will use the term CAM noise only for those processes (Eq. 9) that are not OU processes. For parameters of the CAM process that lead to α-stable noise with only positive jumps, i.e., β=1, we see that jumps are caused by sporadic peaks of varying sizes of the CAM noise. These sporadic large-amplitude peaks will constitute our prototypical noise representation of WWBs (cf. Fig. 1a). In our application, we need to limit the magnitude of the CAM noise for WWBs not to have an arbitrarily large amplitude. We will therefore replace below the output of the CAM noise (Eq. 9) by max(ξ,θ), with θ=3 unless stated otherwise.

Figure 2 Time series and integrated time series of OU and CAM noise processes. (a) OU process with a characteristic decorrelation decay time of 4.9 d ($c_{\mathrm{1}}=-\mathrm{1.4}$, c₄=0.7, $c_{\mathrm{2}}=c_{\mathrm{3}}=\mathrm{0}$). (b) CAM noise process with intermittent peaks (c₁=1.22, c₂=1.14, c₃=0.65 and c₄=0.8 with corresponding α=1.88 and β=0.81). (c) and (d) show the integrals of (a) and (b), correspondingly.

3 Results

To examine the effect on ENSO of the noise model as a representation of WWBs, we consider three different noise models:

1. A multiplicative OU noise, a commonly used representation of WWBs (Jin et al., 2007; Levine et al., 2016; Vialard et al., 2025), denoted OU.

2. An additive noise model with non-Gaussian CAM noise, denoted CAM.

3. A conditional noise model, denoted CON, where we employ additive OU noise for negative SST anomaly, T_e<0, and then allow for more intense WWBS modeled by employing additive CAM noise for T_e>0.

Figure 3 Time series of multiplicative Gaussian noise and of CON noise over 50 years. (a) Multiplicative Gaussian noise model with ${\mathit{\xi }}_T={\mathit{\sigma }}_T(\mathrm{1}+B\max (T_{\mathrm{e}},\mathrm{0}\left)\right){\mathit{\zeta }}_t+T_{\mathrm{e}}^{\mathrm{2}}+c{T}_{\mathrm{e}}^{\mathrm{3}}$ with parameters from Vialard et al. (2025) and Gaussian noise. (b) CON noise model with parameters as in Table A1.

Note that in CON, the noise is state-dependent, and can therefore be considered multiplicative noise, even if its amplitude is not proportional to the temperature as in OU.

Figure 1b shows the NINO3 index from 1870 until 2026 from the NOAA data set (Rayner et al., 2003), where we subtracted the seasonal cycle and employed a moving average over a 4-month window. The focus of this work is on large El Niño events, which we define to be those events with an NINO3 index larger than 1.5 °C. We expect such large events to be preceded by strong WWBs, consistent with observed large El Niño events such as those of 1997 and 2015 (marked by vertical gray lines in Fig. 1b), which were accompanied by multiple strong WWBs (Puy et al., 2019).

Figure 4 Comparison of typical SST time series for different noise models. Time series of 60 years of the NINO3 index (a) are shown together with the SST T_e obtained from the RO model (Eqs. 1–2) when driven by (b) OU, (c) CAM and (d) CON.

To illustrate to what degree the conditional noise model CON generates forcing signals that are more consistent with the observed WWB time series as shown in Fig. 1a, we show in Fig. 3 a random realization of the CON model and of the multiplicative Gaussian noise model with ${\mathit{\xi }}_{\mathrm{T}}={\mathit{\sigma }}_{\mathrm{T}}(\mathrm{1}+B\max (T_{\mathrm{e}},\mathrm{0}\left)\right){\mathit{\zeta }}_{\mathrm{t}}+T_{\mathrm{e}}^{\mathrm{2}}+c{T}_{\mathrm{e}}^{\mathrm{3}}$ with Gaussian noise ζ_t that is used in Vialard et al. (2025), representing a commonly used choice in RO models. The additive CON noise ξ nicely produces sporadic bursts that are clearly separated from a noisy background. In contrast, the large amplitude signals supported by the multiplicative Gaussian noise model are less pronounced.

Our aim is to reproduce the statistical behaviour of both WWBs and the resulting observed signatures of ENSO, including those associated with extreme El Niño events. Each of our three noise models is calibrated to reproduce the empirical histogram of the observed NINO3 time series, its power spectrum, and its variance and skewness, all based on the 153-year-long period from 1870 until 2026 (i.e., 1872 months). We summarize the model parameters in Table A1. For the conditional noise CON, we examine the three month average of the SST to determine the persistence of a positive or negative SST anomaly, based on which the noise is switched between OU and CAM. For computational ease, in the RO model (Eqs. 1–2) the average over the past three months is calculated as the average over the last three monthly snapshots, rather than as an average over all time steps occurring during the past three months. To match the empirical histogram for the conditional noise model CON, we find we need to use the nonlinear SST response (Eq. 7), which in effect shifts the histograms to higher values of the SST, enhancing its asymmetry.

Typical time series of the East Equatorial Pacific SST, T_e, are shown for the three noise models in Fig. 4. We show some statistical characterizations of these time series below. Yet one might argue just based on these snapshots that multiplicative OU (Fig. 4b) and the CON noise (Fig. 4d) seem to produce similar and reasonable ENSO simulations, as also supported by the statistical analyses below. The large sporadic peaks of the CAM noise as seen in Fig. 2b give rise to large El Niño events and hence to a higher degree of asymmetry between El Niño and La Niña and a more skewed distribution compared to the OU process. We recall that we do not use the stochastic signal obtained from Eq. (9) directly but cap the signal to be bounded and not to exceed a threshold of θ=3.

Figure 5 shows a contour plot of the empirical 2D histogram of the variance and skewness for each 1872-month-long segment of a 10⁸ months-long simulation of the SST T_e for the RO model (Eqs. 1–2), together with the mean of the variance and skewness on the observed NINO3 index. Note that the nonzero skewness for the additive CAM noise is entirely generated by the CAM noise process, which favours positive amplitudes with β=0.8075. Hence, CAM noise is capable of generating ENSO's asymmetry of having stronger El Niño events than La Niña events without any multiplicative noise or deterministic nonlinearities.

Figure 5 Histograms of the variance and skewness of the SST T_e for different noise models. The red dot demarcates the observed variance V[T_e] and skewness Skew[T_e] from the NINO3 index. The black dot is the average of a 10⁸ month-long simulation of the recharge oscillator model (Eqs. 1–2) when driven by (a) OU, (b) CAM, (c) CON.

Figure 6 Comparison of histograms and power spectra of observations and the RO model driven by the different noise models. Top panels: histograms of the observed NINO3 index from 1870 until 2026 with a monthly resolution and with the seasonal cycle removed (blue bars) and of the solution of the RO model (Eqs. 1–2) shown by continuous curves, for the three noise models, (a) OU, (b) CAM, (c) CON. We show the 5th and 95th percentile confidence levels for the NINO3 index data as bars (beige) estimated from 10 000 bootstrapped samples. Bottom panels: power spectra S(T_e) of the observed NINO3 index from 1870 until 2026 with a monthly resolution and with the seasonal cycle removed (red) and of the solution of the RO model (Eqs. 1–2) for the three noise models (blue), (d) OU, (e) CAM, (f) CON. A Welch window of 16 years was employed.

Figure 7 Occurrence of WWBs before peak events for the different noise models. Empirical histogram of the number of large noise amplitude events with $\mathit{\xi }>\mathrm{10}/{\mathit{\nu }}_{\mathrm{0}}$ occurring in the 12 months preceding a warm event in the RO model (Eqs. 1–2). Shown are the empirical histograms for the largest and the smallest fifths of El Niño events (with a duration of at least 4 months). Parameters as in Fig. 5. (a) OU, (b) CAM, (c) CON.

Figure 6 shows the corresponding histograms for the SST T_e and the power spectrum together with the corresponding curves for the observed NINO3 index. It is seen that all three noise models are capable of reproducing the global observed ENSO statistics reasonably well. To quantify to what degree each noise model can reproduce these statistics, we have compared several diagnostics. We performed a two-sample Kolmogorov–Smirnov test to check if any of the RO models generates SST data that could have been drawn from the distribution implied by the observational NINO3 index data. The null hypothesis was rejected for all three noise models, with CAM performing worst, as to be expected, as CAM noise is designed to generate only extreme warming events, rather than a wide range of events. The p value of the OU model and the CON model are comparable with p=0.0114 and p=0.0103, respectively. We further estimated the Wasserstein distance between the synthetic T_e model output and the NINO3 index, which measures the difference between their respective distributions. The OU model and the CON model again exhibit similar Wasserstein distances from observations, of 0.0552 and 0.0569, respectively. The data generated by the CAM noise exhibits again a larger Wasserstein distance from observations, of 0.1234, which is to be expected as these data are heavy tailed with several extreme events. To quantify the differences in the generated power spectra, we used a simple mean-square error. Here all three models show roughly the same error with 0.1344, 0.1091 and 0.115 for OU, CAM and CON, respectively. We remark that we did not optimize the equation and noise parameters to produce optimal fits with the observations but rather used a visual assessment.

In summary, all three noise models are capable of producing the global bulk statistics of the observations, measured by the histograms, power spectra, and the variance and skewness. The CAM noise model is designed to generate extreme events and hence performs worse on the histograms. However, our aim for the noise model here is to not only fit the global statistical features of ENSO, but also to reproduce some of the dynamical features and signatures associated with WWBs.

To better distinguish the capability of the respective noise models to reproduce the effect of large WWBs on the dynamics, we now seek more fine-grained dynamical signatures associated with large-amplitude El Niño events. Observed strong El Niño events are accompanied by a sequence of WWBs (Chiodi et al., 2014; Liang and Fedorov, 2021a). This is part of the positive feedback of warmer SSTs promoting the probability of the occurrence of WWBs, and WWBs intensifying the SST warming (e.g., Tziperman and Yu, 2007; Eisenman et al., 2005; Gebbie et al., 2007). Figure 7 shows the number of months in which WWBs occurred during the 12-month period preceding an El Niño event. Strong WWBs are defined as large noise events with $\mathit{\xi }>\mathrm{10}/{\mathit{\nu }}_{\mathrm{0}}$; for CON we choose the smaller of the two values for ν₀. We show results for the 20 % largest and for the 20 % smallest El Niño events (out of a long simulation of a total of 10⁶ months). Consistent with the observation of strong El Niño events co-occurring with several WWBs, the temperature T_e of the strongest El Niño events is accompanied by a much larger number of WWBs for CAM and CON noise compared to small El Niño events (Fig. 7b and c). This effect is still observable to some degree for the multiplicative OU process (Fig. 7a), albeit to a much smaller degree.

Recall that, contrary to the OU and the CAM noise model, the CON noise model involves an asymmetric response of the SST with γ₁≠γ₂. An obvious question is, if the observed behaviour depicted in Fig. 7c is due to this asymmetric response rather than to the proposed conditional CAM noise model. To clarify this, we show in Fig. 8 the empirical histogram of the number of large noise amplitude events occurring during warm events for a conditional noise with a symmetric response (and adjusted noise amplitudes ν₀ to match the observed power spectrum as well as the implied variance and skewness of T_e). It is clearly seen that the conditional noise model with a symmetric thermocline/SST response also exhibits the characteristic behaviour observed for real WWBs. The asymmetric response is still required, however, to allow for a better approximation of the empirical histogram of the temperature (cf. Fig. 6c); a symmetric response leads to a histogram shifted to smaller temperatures (not shown).

Figure 8 Showing that the asymmetric thermocline response is not responsible for the increased number of WWBs prior to extreme El Niño events for CON (Fig. 7c). Empirical histogram of the number of large noise amplitude events with $\mathit{\xi }>\mathrm{10}/{\mathit{\nu }}_{\mathrm{0}}$ occurring in the 12 months preceding a warm event for the RO model (Eqs. 1–2) driven by a conditional noise model. All parameters are as for CON but with a symmetric response ${\mathit{\gamma }}_{\mathrm{1}}={\mathit{\gamma }}_{\mathrm{2}}=\mathrm{0.08}$ and adjusted noise amplitudes ν₀=17.5 for the OU noise component and ν₀=5.4 for the CAM component.

4 Conclusions

We explored three different noise models for representing the effects of WWBs on ENSO in an RO model. One is multiplicative OU noise, a choice often used in the literature for this purpose (Vialard et al., 2025). Another is based on correlated additive and multiplicative (CAM) noise (Sura and Sardeshmukh, 2008; Sardeshmukh and Sura, 2009; Penland and Sardeshmukh, 2012; Sardeshmukh and Penland, 2015) noise, characterized by the occurrence of sporadic large peaks, resembling realistic WWB events (cf. Fig. 1a). And a final one, CON, that conditionally switches from OU for negative NINO3 to CAM for positive NINO3, where the conditional switching is motivated by the observation that WWBs respond to the SST (Yu et al., 2003; Tziperman and Yu, 2007).

All three models were able to explain the global ENSO statistics such as the empirical histogram, the power spectrum, and the variance and the skewness of the observed NINO3 time series. The various noise models, however, give rise to different noise characteristics that are not all appropriate for representing WWBs. The point we are making here is that, from an understanding point of view, it is not sufficient for the noise model to produce ENSO's overall statistical characteristics, but the noise itself should resemble the observed WWB characteristics. Furthermore, characteristics such as having more WWBs preceding large El Niño events should be accounted for as well. Our results suggest OU noise may be an appropriate forcing leading to normal El Niño events of small and moderate amplitudes. Whereas we find that CAM noise is better suited to generate high-amplitude events, accounting for a larger number of WWBs preceding them. Our proposed conditional noise model CON combines OU and CAM noise models to construct a noise model that accounts for both small and large amplitude El Niño events, allowing for a better dynamical representation of WWBs. In particular, CON noise generates a noise forcing that more faithfully represents the sporadic nature of WWBs when compared to additive and multiplicative Gaussian noise models. Moreover, CON noise reproduces the observed dynamical signature that extreme El Niño events are preceded by several WWB events in the 12 months preceding the peak.

Further exploration of this idea using more realistic climate models seems an appropriate future direction.

Appendix A: Model parameters

Table A1 Parameters for the RO model (Eqs. 1–2) and the noise model (Eq. 9). Here the units are as follows: $\left[r\right]=\mathrm{1}/\mathrm{month}$, $\left[\mathit{\epsilon }\right]=\mathrm{1}/\mathrm{month}$, $\left[a\right]=m^{\mathrm{3}}/N$, $\left[b\right]=\mathrm{1}/K$, $\left[{\mathit{\gamma }}_{\pm }\right]=K/m$, [ν₀]=1, $\left[{\mathit{\nu }}_{\mathrm{1}}\right]=\mathrm{1}/K$.