A Non-Gaussian Stochastic Model for the El Niño Southern Oscillation

Air-sea interactions in the tropical Pacific drive the largest interannual variability process in climate called the El Niño-Southern Oscillation (ENSO), with an influence that reaches the higher latitudes via atmospheric and oceanic teleconnections McPhaden2006 ; GhilRMP . Due to its spatiotemporal impacts, estimating the seasonal to interannual state of ENSO is essential for predicting a wide range of regional and global climate phenomena Battisti1995 ; Havlin2015 . However, the complexity of ENSO, involving as it does stochastic forcing from atmospheric transients (e.g., Kleeman1997, ) and nonlinear air-sea interactions (e.g., An2004, ), poses great challenges in both modeling and prediction.

There are two traditional methods of modeling ENSO (see e.g., Tang2018, , for a review). One is to use the state-of-art coupled general circulation models (GCMs), which solve the governing conservation (energy, mass, momentum) equations for the rotating planet with many sub-grid-scale parameterizations of unresolved small-scale processes. The other is to build low-order statistical models. The former treat the main physical aspects of ENSO and thus provide reasonable representations of the spatiotemporal patterns and air-sea interactions from intraseasonal to interannual time scales Jan2005 . The latter, commonly designed to describe and predict particular ENSO indices, characterize the large-scale features of the system and are much cheaper computationally than GCMs Petrova2020 . Nonetheless, despite great progress in developing these methods, there remain limitations to advancing analysis and prediction. On the one hand, GCMs incur systematic errors originating from inaccurate state estimations and incomplete sub-grid-scale parameterizations associated with the ocean thermocline, tropical instability waves, and the Madden-Julian Oscillation (MJO) Chen2008 ; Barnston2015 . On the other hand, the statistical models based on multivariate regressions cannot provide detailed physical information regarding the phase and intensity of ENSO Tippet2019 . One direction for progress is to utilize a hybrid strategy that simultaneously maximizes the advantages of both physically-oriented and statistical models.

The basic theoretical understanding of the tropical air-sea interactions underlying ENSO Cane1985 ; Jin1993 has lead to a hierarchy of simple models of the main processes that control the sea surface temperature (SST) in the tropical ocean (e.g., Wang2017, ). A prominent approach is the recharge-discharge model Jin1 , and its extensions and generalizations, derived from the forced shallow water equations using the two-strip approximation. Despite it being a two-dimensional linear model, it (and its generalizations) can capture inter-annual variability, seasonal variability with time-periodic coefficients Stein2010 ; Levine2015 , and a range of stochastic forcing representing weather and intraseasonal processes, such as westerly wind bursts (WWBs) Chen2017 .

Here we describe a two-stage stochastic model approach that captures the large-scale dynamical and statistical characteristics of ENSO. Stochastic recharge-discharge oscillator models have used observations to recover statistics on inter-annual time-scales Burgers1999 ; Burgers2005 , and have incorporated the seasonal cycle of the Bjerknes feedback to treat ENSO phase-locking and the spring predictability barrier Stein2010 ; Stein2011 ; Levine2015 . However, while our approach builds upon the recharge-discharge model, we incorporate a slowly-varying low-order deterministic component that systematically treats time-varying coefficients and multiplicative noise to accurately describe the central intraseasonal to interannual features of ENSO. In particular, our framework reproduces ENSO’s observed seasonal phase locking, uncertainty and highly non-Gaussian statistics. However, we clearly depart from multivariate regression statistical models that require observations of all state variables. Rather, our non-Gaussian model incorporates a subset of observations from which we can infer unobserved quantities and quantify their uncertainties using an efficient and exact data assimilation scheme, and can thereby accurately recover the difficult to observe thermocline depth and wind bursts. Finally, by systematically determining time-dependent model coefficients, we precisely simulate the seasonal to inter-annual ENSO statistics.

We begin with the following coupled model,

in which $x(t)$ is the averaged SST anomaly in the equatorial central-eastern Pacific (the so-called Niño 3.4 region), and $h(t)$ is the thermocline depth, which is a surrogate of the heat content, averaged over the western tropical Pacific. The time-dependent function, $a(t)$, represents the seasonal evolution of the Bjerknes feedback and the noise amplitude, $N(t)$, captures the role of the relatively short time-scale forcing of the SST anomaly $x(t)$. Both $\xi_{x}(t)$ and $\xi_{h}(t)$ are independent Gaussian white noise processes. Finally, the other time-dependent function, $\omega(t)$, controls the coupling of tropical atmosphere and the upper ocean, which is one of the key parameters that determines the quasi-oscillatory behavior of the dominant ENSO modes, and in traditional low-order models is treated as a constant Stein2010 . Here, however, treating it as time-dependent facilitates a more accurate description of the intraseasonal and interannual statistics.

We emphasize the role of the seasonality of background fields in the tropical Pacific by making all the time dependent coefficients in (1) periodic on the annual cycle. We determine the functions $a(t)$ and $N(t)$ using the monthly statistics of the Niño 3.4 index Moon2017 . The coefficient $a(t)$ is related to the relationship between the monthly variance $\langle x^{2}(t)\rangle$ and the covariance between two adjacent months $\langle x(t)x(t+\Delta t)\rangle$. We use the derived $a(t)$ to estimate $N(t)$ from the statistics of the derived data; $y=x(t+\Delta t)-x(t)-a(t)x(t)\Delta t$. It has been shown previously that $a(t)$ and $N(t)$ can be accurately estimated independent of the other coefficients (see supplementary material in Moon2017 ). The coupling function, $\omega(t)$, and the constants $\lambda$ and $\sigma$, are determined by an expectation-maximization algorithm Dempster1977 ; Sundberg1974 ; Sundberg1976 , which is a useful and efficient statistical estimation method involving incomplete data. We note that only the Niño 3.4 index data is used as the partial observation for model calibration, whereas the thermocline depth, $h(t)$, is assumed to be unknown since it is not directly available from satellite observations. Therefore, as detailed in Appendix B, C, we infer $h(t)$ together with unknown coefficients $\omega(t)$, $\lambda$, and $\sigma$ chen2020learning .

We evaluate the quality of the model by comparing the model statistics with those of the observed Niño 3.4 index in Fig. 5. A qualitative comparison between a random realization of the model and the observed data show comparable variability in Fig. 5(a) and quantitative comparisons of the probability density functions (PDFs) and the autocorrelation functions (ACFs) of the SST anomaly shown in Figs. 5(b) and (c) demonstrate that the model is able to accurately simulate the interannual quasi-oscillatory behavior of ENSO. Moreover, as seen in Fig. 5(d), the model reproduces monthly standard deviations of the Niño 3.4 index data with high precision. Therefore, two important seasonal features of ENSO, “phase locking” to the seasonal cycle Tziperman1998 and the “spring predictability barrier” Levine2015 , are accurately captured in this two dimensional model. Phase locking is related to the activation of the positive Bjerknes feedback from summer to early winter, which leads to the concentration of abnormal ENSO events during November and December, as is clear from the maximum of the standard deviation near the end of the year. This phase locking feature is clearly seen from the model trajectory (Fig. 5a), where the peaks of nearly all of the major ENSO events occur in the boreal winter, consistent with the observations. The spring predictability barrier is related to the increase of the noise magnitude, $N(t)$, during the spring, which increases model uncertainties during the boreal summer thereby causing the loss of predictability Moon2017 . Indeed, this is reflected by the structural similarity between the model and the observations in the persistence forecast diagrams shown in Figs. 5(e) and (f).

Despite its simplicity, the model system of Eqs. (1) recovers the monthly standard deviation and ACFs of the Niño 3.4 index with higher accuracy than coupled GCMs Bellenger2014 , linear and nonlinear statistical models Kondrashov2005 , and simpler recharge-discharge based approaches Stein2010 . Our systematic statistical estimation of coefficients enable this stochastic model to regenerate the statistics of the Niño 3.4 index. In particular, as shown in Fig. 5(c), the precise representation of the ACFs beyond 2 years suggests that the overall dynamical features of the thermocline depth, $h(t)$, are well captured in this model. Variations in $h(t)$, and thus of the equatorial warm water volume (WWV) in the western tropical Pacific, influence the SST through vertical advection of temperature anomalies by mean upwelling, a process known as the “thermocline feedback”. Thus, the thermocline depth acts as a precursor for the inter-annual prediction of SST anomalies Meinen2000 .

Although the linear model described above successfully captures the seasonal variability and the interannual large-scale dynamical features of ENSO, it only allows Gaussian PDFs, whereas the observed statistics are highly non-Gaussian. Importantly, a model must include atmospheric wind burst forcing, which plays a critical role in generating the extreme ENSO events responsible for non-Gaussianity Thual2016 ; Chen2017 . Therefore, we include wind bursts in the stochastic model of Eqs. (1) as follows,

where we adopt a simple linear process with multiplicative noise to generate the stochastic wind bursts $\tau(t)$ that drive the SST and thermocline dynamics. The parameters in the wind burst equation are chosen to accurately reproduce the observed wind stress (Fig. 6 (d)-(f)), while the coefficients of Eqs. (1) were slightly modified to take into account the additional process in the model. Despite the non-Gaussianity of the model, the thermocline depth and the wind bursts remain as Gaussian processes conditioned on the observed SST. Moreover, by using the partial observations that involve only the Niño SST data, we have closed analytic formulae for solving the conditional distribution of the unobserved thermocline depth and the wind bursts (see chen2018conditional, , and Appendix B),

where $\mu_{f}(t)$ and $R_{f}(t)$ are the so-called filter mean and filter covariance in data assimilation.

The filter estimate (3) of the non-Gaussian model (2) can be used to carry out an online time series reconstruction of the unobserved variables and analyze the associated uncertainty, which can also serve as the forecast initialization. Note that the observational data for the equatorial heat content (the integrated WWV above the 20^oC isotherm between 5N-5S, 120E to 80W Meinen2000 , which is equivalent to the thermocline depth $h(t)$ in (2)) is only available since 1983 ¹¹1The WWV values are derived from the temperature analyses of the Bureau National Operations Centre (BNOC) at the Australian Bureau of Meteorology Smith1995 , whereas SST data is available from the pre-industrial era. Therefore, the available SST data can be incorporated into the model (2) to reconstruct the thermocline data over a much longer period.

We first validate the accuracy of the recovered $h(t)$ and $\tau(t)$ in the period after 1983. Fig. 6 (a) shows the filter mean time series of the thermocline depth, $h(t)$, (red) from the non-Gaussian model (2), which compares favorably to the observational value (blue) with a pattern correlation that is nearly $0.9$. Moreover, the small error between the observed and the recovered time series lies approximately within one standard deviation of the filter estimated uncertainty (pink shaded region). In particular, as expected from the physical role of the WWBs, the extreme events are accurately recovered. Similar conclusions are reached for the wind burst time series (Fig. 6e). In addition to a skillful state estimation, the reconstructed wind burst time series has the further influence of filtering out the independent white noise in observations. As shown in Figs. 6(b, f), the model PDFs capture those of the observations for both the thermocline depth and the wind bursts, particularly the significant non-Gaussian skewness and kurtosis. Moreover, similar accuracy is found in ACFs and the power spectra as shown Figs. 6(c, d, g, h). The success in recovering unobserved processes is due to the confluence of the model containing the essential physics and properly constrained stochastic forcing and coefficients.

As shown Fig. 6, the accurate nonlinear filter estimates of Eq. (3) provide an online estimation method of the thermocline depth $h(t)$ and the wind bursts $\tau(t)$ relevant for model initialization. On the other hand, reconstructing historical data can be achieved by solving a slightly different conditional distribution

where $t\in[0,T]$. This is the so-called smoother estimate. It differs from the filter in the sense that it exploits all the available information up to the current time instant. One desirable feature of the nonlinear smoother estimate in Eq. (4) associated with the non-Gaussian system in Eqs. (2) is that it can be solved via closed analytic formulae (see Appendix B). We use the smoother estimate to recover the historical realizations of the two hidden variables conditioned on the entire observational period of the Niño 3.4 index from 1870 to present. The results are shown in Fig. 3. Since the thermocline depth and the wind bursts are the most important precursors for predicting the state of ENSO, these recovered variables can be used to evaluate and improve prediction models and provide a guideline for improving parameterizations in coupled GCMs.

Understanding the complexity of climate dynamics presents a dilemma of completeness, between using solely statistical regression models lack the essential physics, to global climate models, which have more variables and parameterizations than observables. Our approach presents a compromise through hybridization. We combine a slowly-varying low-order dynamical system with an observationally faithful statistical representation of short-time processes. In consequence we find a realistic description of ENSO as measured by the success in (i) reproducing the observed seasonal phase locking and its uncertainty; (ii) capturing the spring predictability barrier; (iii) simulating the observed highly non-Gaussian statistics, and (iv) accurately recovering the intermittent time series of the hidden processes via a nonlinear but analytically solvable data assimilation scheme. Furthermore, we suggest that extending the approach described here to examine climate variability on a wide range of time scales GhilRMP ; PinkPRL , wherein proxy data are available, will provide important insights into mechanisms and feedbacks in the climate system. Finally, given the generality of the stochastic dynamical systems in physics, it is hoped that our approach can be adopted not only for other problems in climate, but far more broadly in the physical sciences.