## Abstract

El Niño–Southern Oscillation (ENSO) has global effects on the hydrological cycle, agriculture, ecosystems, health, and society. We present a novel nonhomogeneous hidden Markov model (NHMM) for studying the underlying dynamics of sea surface temperature anomalies (SSTA) over the region 15°N–15°S, 150°E–80°W from January 1856 to December 2019, using the monthly SSTA data from the Kaplan extended SST v2 product. This nonparametric machine learning scheme dynamically simulates and predicts the spatiotemporal evolution of ENSO patterns, including their asymmetry, long-term trends, persistence, and seasonal evolution. The model identifies five hidden states whose spatial SSTA patterns are similar to the so-called ENSO flavors in the literature. From the fitted NHMM, the model shows that there are systematic trends in the frequency and persistence of the regimes over the last 160 years that may be related to changes in the mean state of basin temperature and/or global warming. We evaluated the ability of NHMM to make out-of-sample probabilistic predictions of the spatial structure of temperature anomalies for the period 1995–2016 using a training period from January 1856 to December 1994. The results show that NHMMs can simulate the behavior of the Niño-3.4 and Niño-1.2 regions quite well. The NHMM results over this period are comparable or superior to the commonly available ENSO prediction models, with the additional advantage of directly providing insights as to the space patterns, seasonal, and longer-term trends of the SSTA in the equatorial Pacific region.

## 1. Introduction

El Niño–Southern Oscillation (ENSO) is the most important interannual climate mode, with significant effects on the planetary hydrological cycle, agriculture, ecosystems, health, and society. The name “El Niño” was proposed at the end of the nineteenth century by Peruvian geographers to describe the warming of the surface waters of the Pacific Ocean around the coasts of Peru and Equator during Christmas (Carranza 1891). Later Sir Gilbert Walker described the existence of the Southern Oscillation following a severe drought in India between 1918 and 1919 (Walker 1923, 1925; Walker and Bliss 1932). Bjerknes (1969) recognized that El Niño (EN) and the Southern Oscillation (SO) represent two different aspects of the same system. Moreover, he clarified the critical role of the positive feedback between the ocean and the atmosphere for the warming of the tropical Pacific Ocean (Lindzen and Nigam 1987), Bjerknes called this feedback ENSO and described the existence of a cold phase, which would later be called “La Niña” (LN) (Philander 1990).

The occurrence of the 1982/83 El Niño event made ENSO a recent focus of research on tropical climatology. Important advances on its understanding, characterization, and modeling are presented by Rasmusson and Carpenter (1982), Cane and Zebiak (1985), Cane et al. (1986), Zebiak and Cane (1987), Enfield (1989), Philander (1990), McCreary and Anderson (1991), Battisti and Sarachik (1995), Trenberth (1997), and McPhaden (2012).

The spatial pattern of sea surface temperature (SST) is key to the evolution and impacts of each ENSO event. Trenberth and Stepaniak (2001) introduced the trans-Niño index (TNI), defined as the difference between normalized SSTA between the regions Niño-4 and Niño-1 + 2, to capture the different “flavors” of ENSO. Using the TNI, Larkin and Harrison (2005) postulated the existence of two types of El Niño, the first followed the “conventional” definition of simultaneous heating in the eastern and central regions of the Pacific Ocean, and the second, related only to warming in the vicinity of the date line. This latter flavor began to be of particular interest to the scientific community because of the El Niño episode of 2004. In the literature, there are several names to these patterns, including El Niño Modoki (Ashok et al. 2007), central El Niño (CP), and east Pacific (EP) events (Kao and Yu 2009). Presently, the term used to denote the variety of spatial patterns associated with warm and cold events is ENSO “diversity” (Johnson 2013; Capotondi et al. 2015). Although, Takahashi et al. (2011) propose that the different “variability modes” of ENSO may represent a continuum.

A fundamental aspect of the observed ENSO is the positive asymmetry between its warm and cold events. Generally, strong EN events tend to reach a higher intensity than strong LN events. Major EN events are centered in the Eastern Pacific, while the center of strong LN events is often in the central Pacific (Rodgers et al. 2004; Sun and Yu 2009). The temporal dynamics of ENSO may be changing over time. A particular observation is that episodes associated with central heating have appeared more frequently and persistently during the last few decades (Fedorov and Philander 2000; Ashok et al. 2007; Kao and Yu 2009; Kug et al. 2009; Lee and McPhaden 2010; Kohyama and Hartmann 2017; Kohyama et al. 2018). This is a sign of nonlinearity (Kohyama and Hartmann 2017).

### a. ENSO trends and projections

There is weak consensus on whether observed ENSO trends can be attributed to human influence or natural variability (Capotondi et al. 2015). Most climate models project a weakening of Walker’s cell as a result of global warming (Vecchi and Wittenberg 2010), with a consequent decrease in the zonal slope of the equatorial thermocline. McPhaden (2012) show that such changes may explain the increase in CP events, and Yeh et al. (2009) suggest that anthropogenic global warming could cause these changes. Therefore it is expected that the frequency of EP events may increase over the next decades (Cai et al. 2014). However, based on multimillennial modeling, Yeh et al. (2012), suggest that such a thermocline inclination could also be due to natural variability, consistent with the results of Newman et al. (2011) and McPhaden et al. (2011).

Tsonis et al. (2005) analyze the relationship between global temperature and ENSO. They conclude that it is a two-way relationship. It is a well-known fact that a warm event corresponds to an increase in global temperature. However, they find that ENSO is itself affected by the global temperature, that positive (negative) global temperature tendency tends to trigger an El Niño (La Niña). Thus, in a warming climate, El Niño events (based on the current definitions) may be more frequent than La Niña events.

The amplitude response of ENSO to global warming is a fundamental issue that remains open despite numerous analyses (Chen et al. 2017). For instance, Cai et al. (2015) project a Walker circulation slowdown, which weakens equatorial Pacific Ocean currents, boosting the occurrence of eastward propagating warm surface anomalies that characterize observed extreme El Niño events, whereas Kohyama et al. (2018) argue for a weakening of strong events. This lack of consensus reflects the limits of the theory and modeling of tropical Pacific dynamics and ENSO.

Held and Soden (2006) and Vecchi and Soden (2007) argue that the hydrologic and energetic constraints robustly require the global-mean circulation to be weakened by global warming. As a consequence, a weakened Walker circulation is expected to yield an El Niño–like response via reduced upwelling. However, Kohyama and Hartmann (2017) argue that a strengthened Walker circulation does not violate the requirement of the weakening global-mean circulation and hence, their argument does not rule out a La Niña–like mean-state change with increased upwelling. On the other hand, a realistic nonlinear subsurface temperature response to waves could imply an ocean that is more stable where the warmer water in the upper ocean is less likely to be vertically mixed with the colder water in the deeper ocean. The suppressed vertical heat exchange further stabilizes the system. Kohyama and Hartmann (2017) refer to this mechanism as NEWS, the Nonlinear ENSO Warming Suppression, which implies a weakening of the ENSO amplitude under global warming. Besides, this asymmetric weakening response can rectify the mean-state SST to become La Niña-like (Kohyama et al. 2018).

Kim et al. (2017) show that ENSO drives the interannual variation in the global carbon cycle. They state that the sensitivity of the terrestrial carbon flux to ENSO increases under greenhouse warming by 44% ± 15%, indicating a future amplification of carbon–climate interactions. These arguments suggest that hydroclimate changes caused by anthropogenic forcing may enhance the ENSO-related carbon cycle. Kim et al. (2014) predict an increasing trend in the present century ENSO amplitude.

### b. ENSO and general circulation models

Current state-of-the-art coupled general circulation models (CGCM) still have problems simulating important features of ENSO variability including biases in amplitude, spatial structure, frequency, seasonal phase locking, asymmetry between El Niño and La Niña or the strength of feedbacks (Wang and Picaut 2004; Guilyardi et al. 2009; Bellenger et al. 2014; Atwood et al. 2017). Also, future projections about ENSO remain highly uncertain, as ENSO representation strongly varies among current CGCMs (Kim et al. 2014; Bayr et al. 2018). The problems relate to the stability—or lack of it—of the atmosphere–ocean interactions in the Bjerknes feedback (Tang et al. 2018).

Various analyses claim that some GCM models do a relatively good job simulating and/or predicting ENSO. For instance Tang et al. (2018) say that various coupled ocean–atmosphere models with hierarchical complexity have improved predictive skill to a level at which it is possible to make successful real-time ENSO prediction at lead times of seasons and longer. Also, Karamperidou et al. (2014) find the presence of rich ENSO variability in the long unforced simulation of GFDL’s CM2.1. But Atwood et al. (2017) find that like most coupled general circulation models, the CM2.1 suffers from large biases in its ENSO simulation, including ENSO variance that is nearly twice that seen in the last 50 years of observations. Also, Bayr et al. (2018) find that many state-of-the-art climate models exhibit an equatorial cold SST bias in the Niño-4 region. As a consequence, they simulate a too westward located rising branch of the Walker Circulation and only a weak convective response.

The CGCMs could be more “chaotic” than the real world, due to an overactive thermocline feedback, and deficient damping from evaporation and cloud shading, likely related to an equatorial cold tongue bias. Alternately, the noise level, associated with tropical weather, could be higher in the GCM than in nature (Karamperidou et al. 2014). This discrepancy may arise from external radiative forcings present in the historical reconstruction. Comparison to observations showed that the GCMs have a strong, positive feedback signal over the central Pacific, which is not present in the observations (Dessler 2013) and a seemingly too positive tropical mean cloud radiative feedback (Kolly and Huang 2018). Kohyama et al. (2018) demonstrated that the ENSO amplitude responses are associated with the nonlinear subsurface temperature response to oceanic waves and many GCMs do not reproduce this kind of nonlinearity.

### c. ENSO models

Deterministic and stochastic models have been used specifically to address the dynamics of ENSO. Low-order deterministic models with the capacity to reproduce some essential characteristics of the ENSO include: the delayed oscillator (Battisti and Sarachik 1995; Suarez and Schopf 1988), the western Pacific oscillator (Weisberg and Wang 1997), the recharge oscillator (Jin 1997), and the advective-reflective oscillator (Picaut et al. 1996), all of them unified in the model of Wang (2001). Also, Zebiak and Cane (1987) developed a coupled atmosphere–ocean model that reproduced certain key features including the recurrence of warm events at irregular intervals with no anomalous external forcing.

Delay differential equations (DDE) models relevant for ENSO phenomenology (Ghil and Zaliapin 2015) are also presented. These models reproduce the intraseasonal activity reminiscent of westerly wind bursts, spontaneous interdecadal oscillations, the observed quasi-biennial oscillation in tropical Pacific SSTs, a phase locking of the solutions’ local extrema to the seasonal cycle, and the irregularity of occurrence of strong El Niño events (Mukhin et al. 2015).

There are also various stochastic models of ENSO, including the Penland and Sardeshmukh (1995) linear model that identifies an optimal pattern of sea surface anomalies for the development of warm events. Also, Gehne et al. (2014) use a two-dimensional linear oscillator forced stochastically to estimate ENSO period and decaying time. Before, Trenberth and Hoar (1996) have used a linear system with noise. Navarra et al. (2013) developed a nonlinear system with noise. They suggest that a discrete state sequence can describe ENSO. Hassan et al. (2016) used Markov chains with four states to study the relationships between ENSO and sunspots, and more recently Conti et al. (2017) analyzed the transition probabilities for four prespecified states of ENSO. Some previous works on ENSO dynamics with homogeneous and nonhomogeneous Markov models include Rajagopalan et al. (1997) and Conti et al. (2017). Machine learning approaches have also been developed as noted in Lima et al. (2009, 2015), Tangang et al. (1998), Hsieh (2001), and Ruiz et al. (2005).

### d. This paper’s structure

ENSO is a spatiotemporal dynamical system, where the underlying dynamics represent air–sea interactions or feedbacks that lead to chaotic or quasiperiodic spatiotemporal behavior. There is an open debate about various critical issues: nonlinear versus linear, deterministic versus stochastic, trends in the underlying dynamics, and the role of phase locking to the seasonal cycle. Given this background, we considered the use of the nonhomogeneous hidden Markov model (NHMM) with an automatic selection of hidden states and other parameters to be appropriate for exploring the tropical Pacific SST data and its underlying dynamics. The HMM has been used extensively in machine learning applications and also for modeling regional precipitation dynamics, with latent states that depend on atmospheric circulation regimes. A similar idea is applied here with the notion that regime dynamics and their transitions could be identified from the spatiotemporal SST data. If these regimes are similar to those identified by others, then their transition probabilities, and trends would provide useful insights as to the potential predictability of the different ENSO flavors. The paper proceeds with a brief overview of the SST data that we use, followed by a description and application of the HMM to identify the associated latent states or regimes. We then consider the seasonal phase locking of ENSO through a nonhomogeneous hidden Markov model, using a predictor that describes the seasonal cycle to inform the transition probabilities of the hidden states. A discussion of these results and the resulting temporal expression or trends in the frequency of occurrence of the hidden states follows. The skill of the fitted NHMM for multiseason forecasts of the tropical Pacific SSTA field and ENSO is then assessed for the 1995–2016 period using the model fit to the prior data. We end with a summary of the key findings and some thoughts as to applications and future work.

## 2. Tropical Pacific sea surface temperature data

We use the SSTA dataset derived from the latest version of the Kaplan Extended V2 SST anomaly (SSTA) data product, which combines the work of Kaplan et al. (1998) with the dataset of Smith and Reynolds optimally interpolated (OI) SST (Parker et al. 1994). The method to generate SSTA is described in Kaplan et al. (1997). This dataset is aimed at analyses of large-scale climate variability over the long term and no trend was removed. Monthly data from January 1856 to December 2019 with a spatial resolution of 5° were used. These translated into 156 time series of SSTA distributed along the tropical Pacific Ocean between 15°N–15°S and 150°E–80°W as shown in Fig. 1.

Additional data, like zonal wind and thermocline depth, may help improve a forecast model, but they have a much shorter record. So they do not contribute to the identification of long-term trends, which is an objective of our paper. For the same reason, no trend was removed from the SSTA record.

## 3. The nonhomogeneous hidden Markov model

A hidden Markov model considers a finite set of latent states with Markov transitions over time. Each latent state is identified with a partition of the state space of the random field modeled. The process is similar to that of cluster analysis, but the clusters have a Markovian dynamics in time, and the spatial field at each time step may correspond to a certain latent state with some probability. Using the transition probabilities for the states, each time period can also be assigned to a specific state in a way that maximizes the likelihood of the assignment over the time period of the data. An external predictor can also be used to inform the change of the transition probabilities over time. In this case the model is considered to be a nonhomogeneous hidden Markov model.

Let $Rt=(Rt1,\u2009\u2026,RtM)$ be an *M*-dimensional vector of SSTA at *M* places at time *t*. Let **R**_{1:T} = (**R**_{1}, …, **R**_{T}) denote a sequence of *T* vectors corresponding to the successive monthly observations since the beginning of the record. The sequence **R**_{1:T} is assumed to be characterized by a Markov chain with hidden (latent) states **S**_{t} = (*S*_{1}, …, *S*_{T}) taking on values from 1 to an optimum hidden state *k*. A HMM defines a joint distribution of **R**_{1:T} and **S**_{1:T} using two conditional stochastic independence assumptions.

The first assumption is that the state sequence **S**_{1:T} is a stationary Markovian process. Therefore the joint distribution function can be computed for any *i*_{t} = 1, …, *k*, *t* = 1, …, *T* in terms of the initial probability function and the transition probability matrix:

where *π*_{i} = Pr{*S*_{1} = *i*}, *i* = 1, …, *k* is the initial probability function and *γ*_{i,j} = Pr{**S**_{t} = *j*|**S**_{t−1} = *i*} is the *k* × *k* state-transition probability matrix.

The second assumption is that given the state **S**_{t} at time *t*, individual observations of **R**_{t} are conditionally independent of all other variables in the model:

Therefore the joint probability of the data and the hidden states can then be written as

For the NHMM, the **R**_{t} is considered to be dependent on a set of exogenous variables **X**_{t}. The NHMM structure is visualized as edges in a directed graph of the HMM as show in Fig. 2, with **X**_{1:T} = (**X**_{1}, …, **X**_{T}) is a sequence of exogenous input vectors.

Here, **R**_{t} represents the spatially distributed SSTA dataset with *M* = 156 locations in the tropical Pacific. The Gamma distribution is considered as a model for an individual SSTA time series in a specific location **R**_{t}. The SSTA field is then modeled as a spatial Gamma process with hidden states **S**_{t} that may be related to ENSO flavors. For the NHMM, we also consider a sequence of *D*-dimensional input column vectors **X**_{1:T} = (**X**_{1}, …, **X**_{T}) such that the transition probabilities of the hidden states also depend on the value of the corresponding input vector **X**_{t}. Specifically, the hidden state on month *t*, depends both on the predictor vector **X**_{t} for month *t* and on the value of the hidden state **S**_{t−1} on month *t* − 1. Because **X**_{t} can vary in time, this implies that the transition probabilities between states (**Γ** matrix) can vary in time in response of multivariate predictor **X**_{1:T}. The NHMM can be written as

ENSO events normally begin their evolution in the spring [March–May (MAM)] or summer [June–August (JJA)] and tend to reach the maximum phase during the autumn [September–November (SON)] or winter [December–February (DJF)]. Overall, ENSO has a strong tendency to be coupled with the annual cycle (Rasmusson and Carpenter 1982; Stein et al. 2010, 2014; Levine and McPhaden 2015). To take that coupling into consideration we choose the annual cycle represented with the following periodic deterministic signals as a conditioning variable for the NHMM:

Under the NHMM model, the conditional log-likelihood LL(Θ)of the observed SSTA data, given the inputs, is

We use a multinomial regression framework to parameterize the hidden state transitions in Eq. (7) by using

Given a fixed number of hidden states *K*, we find the parameters Θ of the NHMM by searching for parameters that best fit the observed data by maximum likelihood. The set of parameters Θ that maximize LL(Θ) can be obtained using the Baum–Welch algorithm (Rabiner 1989), a variation of an iterative expectation–maximization (EM) algorithm (Dempster et al. 1977) for obtaining maximum likelihood parameter estimates for models with hidden variables and/or missing data. For the NHMM case, we start the EM algorithm 10 times from different random starting positions in parameter space, then run the EM algorithm to convergence for each starting point, and choose as our solution the parameter vector **Θ** that has the maximum likelihood over all 10 runs (this helps avoid poor local maxima). For full details on the specific EM procedure used in this paper for NHMM parameter estimation see Robertson et al. (2003) and Kirshner (2005).

## 4. Underlying dynamics of the tropical Pacific SSTA

### a. Hidden states of tropical Pacific sea surface temperature anomalies from NHMM

In this section, we identify the hidden states of ENSO dynamics using the SSTA data. For this purpose, we used the Bayesian nonhomogeneous Markov and mixture models (multiple time series) R package (Holsclaw et al. 2017). The results of the algorithm include the Viterbi sequence with an assignment of a hidden state to each time instance, the monthly transition probability matrix, as a function of the exogenous variable **X**_{t} and all model parameters. The asymmetry in the SSTA series can produce a nonzero residual effect in its temporal mean, for this reason, it is advisable to use non-Gaussian probability distribution functions. In this work, we use the gamma distribution, this requires a positive time series, for that we add a constant value to all the data that is below the minimum of all recorded data for each gridbox time series.

NHMMs were run five times from different initial conditions, applying 40 iterations for the parameters estimation. We use cross validation to evaluate the quality of the fitted NHMMs as a function of *K*, the number of hidden states. The resulting normalized out-of-sample values of the cross-validated log-likelihood for each model are given in Table 1 for *k* = 1, 2, 3, 4, 5, 6, together with a normalized Bayes information criterion (BIC) and Akaike information criterion (AIC) for each model.

The lowest value of the BIC criterion occurs for *k* = 2, while for the AIC criterion it is for *k* = 3. However, the AIC values are essentially the same for *k* = 3, 4, 5, and 6. For five states, the classification leads to spatial composites of SST that are quite similar to the patterns identified in previous work. Since the AIC and the BIC penalize bias and variance of estimation in slightly different ways, choosing 5 states will likely lead to lower bias and higher variance, as reflected by these criteria. We chose 5 states on this basis, since they approximately match the 5 states identified in the literature. The selection of *k* = 6 did not seem to identify a distinct additional state. The lack of clear distinction may actually support the notion of a continuum in the dynamics, rather than these distinct states. Figure 3 shows the composite of the SSTA fields constructed from the most likely states derived from the HMM with *K* = 5.

The states identified can be interpreted from the coldest to the warmest: state one corresponds to episodes of strong La Niña, the second corresponds to moderate (weak) episodes of La Niña, the third corresponds to a neutral state, the fourth corresponds to CP El Niño events, and the fifth one represents episodes of EP El Niño. The spatial distribution of the five states shows that the strongest events of EN present greater anomalies in the eastern Pacific, while strong La Niña events tend to be located within the central Pacific. Thus, the five-state NHMM can reproduce the spatial asymmetry of the SSTA described in the literature. Figure 3 shows the composites fields of SSTA for each state. The use of the annual cycle as an exogenous variable within the NHHM allowed an identification of the ENSO diversity without needing to emphasize specific regions of the Tropical Pacific Ocean.

### b. Temporal evolution of ENSO states

Figure 4 shows the Viterbi matrix of the ENSO system using the NHHM; one can see that the dynamic system suggests a smooth transition from the coldest state to the warmest, for example, the process only arrives to state 5 if it previously was in a state 4. We emphasize that the process may persist in intermediate states such as 2 and 4 without evolving to a stronger state. This suggests that the central El Niño (state 4) and the weak La Niña (state 2) are perhaps worthy of being identified as “regimes” of the ENSO dynamical system. We emphasize that given the AIC/BIC values for the number of states and the continuum in the underlying spatial representation, we expect that there is an underlying continuous dynamics of spatiotemporal evolution that is being approximated here by five discrete states with a monthly transition probability matrix. We do not attach a deeper meaning to these empirically derived states.

Observations from Fig. 4 provide a first approximation to the temporal evolution of the dynamic system. From 1856 to 1920, the cold conditions were quite frequent with a higher recurrence of state 1, whose average duration was over 12 months, whereas state 5 was less frequent and with an average duration of shorter than 6 months. In contrast, from 1968 to 2019 warm events have been the most recurrent, with state 4 being the most abundant, and with state 1 persistence less than six months while the state 5 has average duration over 12 months. These observations are coincident with a general warming trend, but do reflect distinct spatial patterns as recurrent and persistent, rather than a generalized basin wide warming. Consequently, it is tempting to ascribe these changes to a response to the generalized warming. However, based on just our data analysis, and the underlying stationarity assumption, the notion that these changes could still be a part of the natural system variability cannot be ruled out.

### c. Persistence and transition between states

Table 2 presents the average transition probability matrix for the 165 years considered in the present study. The highest values are on the diagonal; this implies a high persistence of the states. The highest value of the diagonal corresponds to remaining in state four whereas the lowest is to stay in state 5, this means that CP like events tend to persist over time while EP like events have historically been of shorter duration.

Notice that entries next to the diagonal are second in magnitude. This observation implies that a transition to a near state is more probable than to farther away states. This reinforces the idea of continuity in the transitions. For instance, going from state 1 to state 2 has an 8.4% probability (i.e., leaving strong La Niña it is more likely to go to a weak La Niña). From state 2 the probability of a transition to state 3 is 7.4% and to state 1 is 4.4%. This lack of symmetry for the cold events indicates a tendency to a shorter duration.

The probability of leaving state 3 to state 2 is 5.4% and to state 4 is 7.0%. However, form state 4 it is more likely to return to state 3 (7.9%) than to state 5 (3.4%). These comparisons could explain the frequency of purely CP events in the ENSO system because warming in the tropical Pacific Ocean does not necessarily evolve to EP events. Finally, leaving from stage 5 a return to state 4 (13.3%) is more likely, therefore the tendency of EP episodes to become CP events once their cycle is over. Recall that a Gamma distribution is used to parameterize the SST for each grid box for each state. The parameters of this distribution then provide insight as to how strong a CP or EP event may be. In the current application, we do not use an exogenous predictor for the parameters of this distribution. Consequently, we have the same probability distribution apply to a particular state. As an NHMM simulation progresses, the predictive probability distribution of the next time period will then be a mixture of Gamma distributions, where the probability of being in a particular state at the next time acts as a weight for the associated Gamma distribution of SSTA. Thus, depending on the pathway taken one could generate simulations that have different strengths for different CP, EP or La Niña events. However, this will be potentially damped in amplitude relative to a model that could inform the amplitude of the SSTA. A possibility here is to include dependence on the prior SSTA value for that grid in addition to the latent state. We have not yet explored this strategy.

Thus, the transition matrix in Table 2 shows that it is unlikely to move from a very cold state to a very warm state directly without any transition. CP events are highly persistent and have a low probability of evolving toward EP events. An EP event has a high probability of going back to a CP event. The self-persistence, as well as low probability of transition to EP events, favor the development of CP type events. As a result, intensification in frequency of CP events is one of the most exciting characteristics of the dynamic system obtained from the NHHM for the ENSO.

Given that the exogenous variable used to force the transition matrix in the NHMM is a sinusoidal function associated to each of the months of the year, it is possible to calculate the annual cycle of the transition matrix. This computation is done averaging the monthly values of each entry. Figure 5 shows the annual cycle of the transition matrix scaled by the annual average value. Again, one notes that the probability of remaining in each of the states is higher and remains almost constant during the year while there are clear cycles for the rest of the transition probabilities.

The following are observations from an analysis of the annual cycle of the transition probabilities to near-neighbor states:

Higher probabilities of going from state 1 to state 2 occur during March, April, and May, with significantly lower values during September, October, and November. Therefore, a transition from a strong La Niña event to a weak La Niña is more likely during the boreal spring.

A transition from state 2 to state 3 is more likely during spring, but the transition from state 2 to state 1 is more likely during the fall; so weak La Niña events tend become neutral events during spring, while weak La Niña events tend to evolve to a strong La Niña events during fall.

Transitions from state 3 to state 4, or from state 3 to state 2 are more likely during the fall, so one can conclude that going from a neutral event to either of the two adjacent states (the weak La Niña or CP) is more likely during September, October, and November.

Transition probabilities from state 4 to state 5 are greater during the fall, while from state 4 to state 3 are higher during spring, so the greater probability that an event CP El Niño evolves into EP is during September, October, November, and December. The weakening of the central El Niño to a neutral condition occurs more frequently during spring.

Finally moving from state 5 to state 4 is more likely during the spring, so usually, an event EP El Niño may become a CP event during March, April, and May.

These results suggest that ENSO dynamics has two predictability barriers, periods when the probability of transition to different states increases. The first one is the spring barrier, where the different events weaken, and the second one is the autumn barrier that usually initiates the different events from the neutral condition or strengthens them until they evolve to strong warm or cold episodes. Thus, the NHMM analysis is effective in adding some insight as to the seasonal dynamics.

### d. Comparison with other methods

To better understand the implications of our findings on ENSO diversity, we make a comparison with several indices that were introduced to identify different event types, with emphasis on the warm (El Niño) ENSO phase. Recall that those indices use SST. Since to reach state 5 the process needs to pass first by state 4, we assume that events type EP-EN are those in which a state 4 goes to state 5. Besides, events CP-EN correspond to those where the process remains in state four before going to a neutral state. We compare our NHMM classifications of El Niño types with four other methods: the central location method (i.e., the NINO method hereafter) used in Kug et al. (2009) and Yeh et al. (2009); the EMI method used in Ashok et al. (2007); the CT and WP index method (CT/WP method hereafter) proposed by Ren and Jin (2011); and the pattern correlation method (PTN) used in Yu and Kim (2013) as shown in Table 3.

Table 4 presents a way of comparing the different methods used to classify the warm ENSO events between CP and EP events (we take MIX events as EP events for this comparison). The accuracy between a pair of models is the ratio of the number of same predictions to the total number of predictions. As the last row shows, the NHMM model is more in concordance with the EMI method. The closest pair of methods is the NINO-CP/WP, with 90% accuracy. Minimum accuracy is 54% for NHMM with both NINO and CP/WP classification methods. Among all the methods the average accuracy is 70%. It is clear that there is not a correct method. There are other indexes used to compare models, but they make more sense when there is a correct standard to compare. For instance, if one assumes that each of the other methods (PTN, NINO, EMI, CP/WP) is the true one, the NHMM is precise 75%, 100%, 82%, and 88% of the time predicting CP events, respectively. One can conclude that the classification methods are not in concordance all the time.

## 5. ENSO evolution and trends

To analyze the temporal evolution of the probability of occurrence of each of the five states associated with the NHMM, logistic regression was performed to estimate the trend of the probability of occurrence of each state during the period 1856–2019. We use the R package for local regression, locfit (https://cran.r-project.org/web/packages/locfit/locfit.pdf) with a logistic link function using default value of the nearest-neighbor fraction (alpha) of 0.7 to choose the neighborhood, with the Viterbi sequence of states identified from the NHMM. Figure 6 presents the temporal evolution of the probability of occurrence of the different states associated with NHMM. This figure suggests that even though the NHMM was fit assuming underlying stationary dynamics, the variability of state occurrence that leads to the maximum likelihood solution does change over time. Again, this figure may also account to some degree for the nonrobustness in the selection of the number of states. Since in the past, some states occur infrequently; their inclusion could potentially increase the variance of estimation but reduce the bias when they do occur, leading to an inability of the AIC/BIC metrics to choose a clear winner under the stationary assumption. One could indeed include global temperature or CO_{2} or time as a covariate in the NHMM to address these transitions. However, as we discussed earlier, ENSO and the global temperature have feedback in their expression, and so this may not lead to a definite diagnosis of how the dynamics may have changed. In effect, by assuming no trends in the dynamics, fitting a model for trends and finding that there are trends in the state occurrence we have a more explicit indication of change. Thus, at this stage, our analysis and model fitting is diagnostic rather than the best predictive model of the near-term evolution of the system.

Figures 6a and 6b together show that the cold events, associated with states 1 and 2, respectively, strong and weak La Niña events, have a decreasing trend in their probability since 1856. Figure 6a shows the time-dependent probability of occurrence for state 1. It varies nonmonotonically with a maximum of 40% to a minimum of 2% at present. The same decreasing trend in Fig. 6b is for the probability of state 2, which goes from 38% in the middle of the nineteenth century to 8% at present. On the other hand, the probability of events associated with neutral conditions, state 3, has a significant increasing trend from 1850 to 1940, when it reached a maximum of 50%, since then there is a decrease to 30%, as shown in Fig. 6c. The decreasing trend in the probability of occurrence of cold and neutral events contrasts with the rapid increase of the warm event probabilities of central El Niño events, shown in Fig. 6d. This probability has increased steeply since 1940, going from 20% to 60% now, explaining the increase in the frequency of this type of events. There is also an increase in the likelihood of occurrence of events associated with EP El Niño, this probability has doubled since 1940 from 4% to 8% as shown in Fig. 6e.

These results indicate significant shifts in the tropical Pacific Ocean conditions, which modulate the probability of occurrence of the different ENSO states, decreasing the occurrence of cold events and considerably increasing the probability of occurrence of warm events. Today central El Niño events are more likely, but it is important to highlight that no extrapolation is valid from those local regressions. Figure 7 shows the evolution of global temperature anomalies from 1880 to 2019 compared to the probability of occurrence of events associated with central El Niño events. As shown, since 1960 the global warming trend follows the increasing trend in the probability of occurrence of central El Niño events. This coincidence suggests that global climate change may have some relation to the probability of occurrence of warm ENSO’s central events.

## 6. ENSO probabilistic forecast

### a. SST anomaly forecast

In this section, we analyze the ability of NHMM to predict the spatial structure of temperature anomalies in the tropical Pacific Ocean. For this purpose, we fit the NHMM model using a training period from January 1856 to December 1994, then perform a probabilistic prediction month by month from January 1995 to December 2016, for multiple lead times. For each lead time, we forecast an ensemble of SSTAs for all the 156 grid boxes.

From the simulated ensemble, at each grid point and for each month, we compute the empirical distribution function *F*(*x*) of the forecasts and then evaluate its value at the actual observed value *x*′. So if the forecast is good, *F*(*x*′) should be approximately 0.5 (i.e., the median of the forecast ensemble). The function *F*(*x*′) lies between 0 and 1. A 0 (1) means that the observation is lower (higher) than all the values in the forecast ensemble. As an aggregate measure of performance, it is useful to measure departures of *F*(*x*′) from the median value of 0.5. This is done by recording *F*(*x*)* = min{*F*(*x*′), [1 − *F*(*x*′)]} at each grid box and time step and then reporting its average for each lead time forecast, for each grid. For a 9-month lead time forecast, the average of *F*(*x*)* is reported just for the ninth month.

Figure 8 presents the spatial distribution of the average of the *F*(*x*)* for +1, +3, +6, and +9 lead times. Continuous lines represent values of *F*(*x*)* greater than 0.3. The results indicate that the NHHM adequately represent the spatial distribution of the SSTA for the entire Pacific Ocean for +1 and +3 lead times, while for +6 and +9 some of the spatial forecasting ability is lost in the eastern Pacific, but the predictions in the central Pacific and the portion closest to the coasts of South America remain acceptable. Therefore, it is a good model for predicting ENSO conditions, as well as the certainty in the estimation of temperature anomalies in the region closest to southern America, even for +6 and +9 leads. These capabilities mean that the model can estimate the temperature contrast between the different regions in the Pacific Ocean and consequently the ENSO flavors.

We also compute a pseudo-log-likelihood across all grid boxes as the product of ln[*F*(*x*)*] for each time. The results presented in Fig. 9 show the forecasting ability for the La Niña event 2011/12 and El Niño 2015/16.

### b. Niño-3.4 SSTA forecast

We aggregated the ensemble forecasts across the grids that comprise the Niño-3.4 region, to examine the predictive ability of the model relative to selected ENSO forecast models that are publicly available. Then we proceed to estimate the pdf of those ensemble forecasts to see whether the 5%–95% interval of the forecast covered the observed historical value at each lead time approximately 90% of the time. Figure 10 shows that the historical values are mostly (more than 90% of the time) within the range 5%–95% for the different forecast leads. This coverage of observed values by the probabilistic forecast indicates the ability of the model to reproduce the magnitude of the temperature anomalies in the Niño-3.4 region. We also estimated the correlation between the ensemble mean and the historical data. These correlations are often higher than 0.65 for the different lead times, and the corresponding RMSEs between the mean ensemble forecast and the observations are typically low. These high correlations and low RMSEs confirm the previous analysis of the *F*(*x*)*.

Finally a comparison between the NHMM approach and 18 different CPC-IRI dynamical (D) and statistical (S) model forecasts, including the persistence forecast, of the 3-month running mean of SSTA in the Niño-3.4 region is presented in Table 5. The NHHM ENSO approach is competitive in a short time (lags 1 and 3) and much better for the longer lead times (lags 6 and 9). Recall that these are out-of-sample forecasts for the NHMM over the 2010–16 period.

To provide probabilistic forecasts of ENSO events consistent with those reported by other modelers, we developed categorical probabilities, using the following thresholds: La Niña when the SST < −0.5°C and El Niño when the SSTA > 0.5°C. As an example, Fig. 11 presents the monthly probabilistic forecast from August 2015 to July 2016 provided by the NHMM associated with the strong event of El Niño 2015/16. The probability of occurrence of each type is provided by the bars while the lines represent the historical probability for each event. In general, for EN event, the NHMM ENSO model indicates a higher chance of EN event occurring relative to the historical average from August 2015 to July 2016. Given these high probability anomalies, an extreme EN event is likely. Moreover, one of the NHMM model properties is the ability to predict the ENSO flavor using the forecast state probability for each month and each lead time. In this case, the NHMM model predicts an EP-EN case (the warmest state 5), from July 2015 to April 2016, and afterward transitions to a CP event.

We also computed verification scores in order to show skills of NHMM forecasts issued for El Niño and La Niña events as defined previously using the probabilistic forecasts from January 2010 to December 2016. Table 6 presents the results of rate of return, likelihood score, likelihood skill score, hit proportion, relative operating characteristic (ROC) area below −0.5, and ROC area above +0.5. Results for most indicators are good until +6 lead time. Figure 12 presents the ROC for the ENSO probabilistic forecast using NHMM. The area under the curve for warm and cold episodes is higher than 0.6 for all leads.

### c. ENSO flavor forecast

As discussed earlier, an especially interesting and novel aspect of NHMM is that we can provide a forecast probability of the whole SST state or for ENSO flavors. Since every hidden state has a specific spatial structure, one can estimate for a given observed SST field the maximum likelihood classification that matches it for each month. It is interesting to see if the forecast state probability is high for the actual classified state (i.e., more specifically for the 2015 and 2016 event, when the central heating areas were not the same as in the past El Niño events). We get a high probability that the state that we forecast matches the heating patterns observed—of course adjoining states (similar ones) may also have a high probability for the NHMM with five states, we can estimate the annual cycle of the marginal probability for each state from Viterbi sequence *p*(*k*, *m*) where *k* denotes the state and *m* the month as is presented in Table 7.

With the specific Viterbi sequence for 2010–16 (Fig. 4) we estimate the marginal likelihood for the historical data during a specific forecast with horizon *l* as $\u220ft=1lp\u2061(t|k,\u2009m)$. Using the NHMM forecast we can also compute the marginal probability for each state at time *t*, *q*(*t*, *k*), and the marginal likelihood for each forecast as $\u220ft=1lq\u2061(t|k)$. Finally the log-likelihood ratio for each forecast at time *t* with horizon *l* is

When LL(*t*, *l*) is higher than zero, the SSTA NHMM forecasts are skillful at reproducing the ENSO flavors. Sometimes when $\u220ft=1lq\u2061(t|k)$ is zero and LL(*t*, *l*) undefined, it means that the marginal probability is zero for at least one of flavors in the *l* forecast sequence, and when LL(*t*, *l*) < 0 the historical marginal probability is better than NHMM to explain the ENSO flavors. Figure 13 shows the log-likelihood ratio for +1, +3, +6, and +9 leads. LL(*t*, *l*) is above zero 78% of the time for lead +1, 74% for lead+3, 72% for lead 6%, and 68% for lead 9; most of the time the NHMM is reproducing ENSO flavors. Again, this is an out-of-sample test. Since we associated each state with a particular ENSO flavor, we could go beyond forecasting Niño-3.4 SSTA to include probabilities of the different flavors. The result is an enhanced probability of the specific ENSO states for each lead time, which also leverages the NHMM method much better than a Niño-3.4 SSTA forecast.

## 7. Conclusions

We find that the NHMM can provide useful diagnostics as to the underlying fluctuations in the dynamics of SSTA associated with the ENSO phenomenon. A novel description of the regimes associated with the spatiotemporal evolution of the SSTA in the tropical Pacific and their persistence characteristics is provided using NHMM. The spatial expression of the regimes is similar to the patterns others have identified through clustering or other classification mechanisms, but with the additional requirement that the seasonal evolution as well as the longer-term evolution admits and models persistence, regime transition dynamics, and spatiotemporal asymmetry.

It is relevant to highlight some additional conclusions:

Although it is well established that ENSO dynamics are the result of ocean–atmosphere interactions, our model shows that SSTs provide significant redundancy and elements to learn from the slow variable in the system with an intelligent embedding.

While we modeled the SSTA field as a stationary stochastic process, with modulation by the seasonal cycle, we found that there are systematic trends in the frequency and persistence of the regimes over the last 160 years that may be related to changes in the mean state of basin temperature and global warming. A caveat is that the data quality is not uniform over time, and we have not addressed that in our analysis. However, knowing how the older data were extrapolated using recent spatiotemporal correlation patterns, one would expect that the amplitudes of events in the past may be muted, but the spatial patterns would be more similar to recent spatial patterns. The finding that the persistence and duration of the events with specific, identified spatial patterns, corresponding to the different flavors of ENSO phenomenon do change over time in a direction that is consistent with what some expect for a global warming response is consequently interesting. One could indeed include global temperature or CO2 as a covariate in the NHMM and argue that the significance of this predictor could provide causal evidence of change. We chose not to present such results here, since we recognize that ENSO dynamics also influence global temperature and this may confound the model if the interaction of these dynamics was not properly specified. In the future, we expect to explore a systematic model development that would account for such mutual feedbacks between the tropical Pacific and the global climate system.

For the short lead times from 1 to 9 months, we found that the NHMM gives comparable to superior results relative to existing models for out-of-sample prediction over the 2010–16 time period. This period has been challenging for most models. By assuming stationarity over the last 150 years over which the NHMM was fit, we have potentially compromised some predictive capacity that would account for changing conditions. On the other hand, we would increase the variance of parameter estimates by adding more parameters. A central purpose in our using an established machine learning tool (e.g., the NHMM) was that given its structure and the detailed characterization of ENSO dynamics in the literature, it was a promising tool for useful diagnostics as to the underlying fluctuations in the dynamics of SSTA. We find that the HMM can identify the ENSO flavors and draw some insights as to their temporal evolution, including their dependence on the seasonal cycle, and long-term trends. Broadly speaking the insights are consistent with those presented in the literature. However, it is important to show that a single model can encode the spatiotemporal structure of the regime dynamics, provide insights into the spring and fall “barriers” or transitions in predictability, provide insights into the temporal trends in the dynamics (regime transition probabilities), and also provide useful near-term forecasts for ENSO, the regimes, and the full spatial field’s evolution.

## REFERENCES

*J. Geophys. Res.*

_{2. 1}coupled GCM

*Climate Dyn.*

*Rev. Geophys.*

*Climate Dyn.*

_{3}to CMIP

_{5}

*Climate Dyn.*

*Mon. Wea. Rev.*

*Nat. Climate Change*

*Nat. Climate Change*

*Science*

*Nature*

*Bull. Amer. Meteor. Soc.*

*Bol. Soc. Geogr. Lima.*

*Climate Dyn.*

*J. Climate*

*J. Roy. Stat. Soc.*

*J. Climate*

*Rev. Geophys.*

*Science*

*Climate Dyn.*

*Bull. Amer. Meteor. Soc.*

*J. Atmos. Sol.-Terr. Phys.*

*J. Climate*

*Ann. Appl. Stat.*

*J. Climate*

*J. Atmos. Sci.*

*J. Climate*

*J. Climate*

*J. Geophys. Res.*

*J. Geophys. Res*

*Climate Dyn.*

*Nat. Commun.*

*Nat. Climate Change*

*J. Climate*

*Geophys. Res. Lett.*

*J. Geophys. Res. Atmos.*

*J. Climate*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*J. Geophys. Res. Atmos.*

*Geophys. Res. Lett.*

*J. Climate*

*Machine Learning and Data Mining Approaches to Climate Science*, V. Lakshmanan, Ed., Springer, 13–21

*J. Atmos. Sci.*

*J. Geophys. Res.*

*Geophys. Res. Lett.*

**39**,

*Geophys. Res. Lett.*

*J. Climate*

*J. Climate*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*J. Climate*

*El Niño, La Niña, and the Southern Oscillation*

*Science*

*Proc. IEEE*

*J. Climate*

*Mon. Wea. Rev.*

*Geophys. Res. Lett.*

*J. Climate*

*J. Climate*

*J. Climate*

*J. Climate*

*J. Atmos. Sci.*

*J. Climate*

*Geophys. Res. Lett.*

*Natl. Sci. Rev.*

*J. Climate*

*Bull. Amer. Meteor. Soc.*

*Geophys. Res. Lett.*

*J. Climate*

*Geophys. Res. Lett.*

*J. Climate*

*Wiley Interdiscip. Rev.: Climate Change*

*Memoirs of the India Meteorological Department*, Vol. 24, Meteorological Office, 75–131

*Mon. Wea. Rev.*

*Mem. Roy. Meteor. Soc.*

*J. Climate*

*Earth’s Climate: The Ocean–Atmosphere Interaction*, C. Wang, S.-P. Xie, and J. A. Carton, Eds., Amer. Geophys. Union, 21–48

*Geophys. Res. Lett.*

*Nature*

_{3}to CMIP

_{5}and its implication of ENSO

*J. Climate*

*Int. J. Climatol.*

*Mon. Wea. Rev.*

Extreme Events: Observations, Modeling and Economics, M. Chavez, M. Ghil, and J. Urrutia-Fucugauchi, Eds., Amer. Geophys. Union