## Abstract

One of the key characteristics of El Niño–Southern Oscillation (ENSO) is its synchronization to the annual cycle, which manifests in the tendency of ENSO events to peak during boreal winter. Current theory offers two possible mechanisms to account the for ENSO synchronization: frequency locking of ENSO to periodic forcing by the annual cycle, or the effect of the seasonally varying background state of the equatorial Pacific on ENSO’s coupled stability. Using a parametric recharge oscillator (PRO) model of ENSO, the authors test which of these scenarios provides a better explanation of the observed ENSO synchronization.

Analytical solutions of the PRO model show that the annual modulation of the growth rate parameter results directly in ENSO’s seasonal variance, amplitude modulation, and 2:1 phase synchronization to the annual cycle. The solutions are shown to be applicable to the long-term behavior of the damped model excited by stochastic noise, which produces synchronization characteristics that agree with the observations and can account for the variety of ENSO synchronization behavior in state-of-the-art coupled general circulation models. The model also predicts spectral peaks at “combination tones” between ENSO and the annual cycle that exist in the observations and many coupled models. In contrast, the nonlinear frequency entrainment scenario predicts the existence of a spectral peak at the biennial frequency corresponding to the observed 2:1 phase synchronization. Such a peak does not exist in the observed ENSO spectrum. Hence, it can be concluded that the seasonal modulation of the coupled stability is responsible for the synchronization of ENSO events to the annual cycle.

## 1. Introduction

El Niño–Southern Oscillation (ENSO) is the largest global climate signal on interannual time scales (Neelin et al. 1998); strong ENSO events cause changes in the tropical Pacific climate that are large enough to influence the global atmospheric circulation (Trenberth et al. 1998), leading to significant environmental and socioeconomic impacts that occur in areas throughout the world (McPhaden et al. 2006). ENSO events occur irregularly, with 2–7-yr spans between them, but they each follow a similar pattern of developing during boreal summer and peaking during boreal winter (Rasmusson and Carpenter 1982; Larkin and Harrison 2002). Such seasonal synchronization is a defining characteristic of ENSO, and understanding the cause is of central importance to ENSO predictions (Balmaseda et al. 1995; Torrence and Webster 1998). The exact mechanism responsible for the synchronization of ENSO to the annual cycle has not yet been determined, although current ENSO theory offers two possible candidates: frequency locking of ENSO to periodic forcing by the annual cycle (Jin et al. 1994; Tziperman et al. 1994), or the modulation of ENSO’s coupled stability due to the seasonal variation of the background state of the equatorial Pacific (Philander et al. 1984; Hirst 1986). The goal of this study is to determine which of these two synchronization mechanisms best explains the observed seasonal characteristics of ENSO.

Evidence supporting frequency locking of ENSO to the annual cycle as a synchronization mechanism comes from investigating the behavior of ENSO models under the variation of relevant model parameters, in particular the simple delay oscillator model (Suarez and Schopf 1988) and the intermediate complexity Zebiak–Cane model (hereafter ZC model; Zebiak and Cane 1987). Generally, a parameter related to the amplitude–growth rate of the model’s ENSO is varied along with a parameter related to the strength of the seasonal forcing (e.g., Tziperman et al. 1995) or the intrinsic frequency of the ENSO mode (e.g., Jin et al. 1996). Within the resultant model parameter space, frequency-locked solutions are a common feature. This is due to the fact that the various ENSO models used in such studies, although differing in details, each follow the quasiperiodic route to chaos and thus admit the same suite of possible model solutions: quasiperiodic solutions that include both the annual cycle and ENSO frequencies, frequency-locked solutions where the ENSO frequency is a rational multiple of the annual cycle, and chaotic solutions that result from the overlapping of multiple frequency-locked solutions within the model parameter space (Tziperman et al. 1995; Jin et al. 1996). This behavior has been demonstrated within the periodically forced delay oscillator (Tziperman et al. 1994; Liu 2002), a two-equation dynamic system model of ENSO (Wang and Fang 1996), the ZC model (Tziperman et al. 1995; Pan et al. 2005), and variations of the ZC model that include coupling the atmosphere to total SST (Chang et al. 1994) and reducing the ocean component to zonal equatorial strip with fixed meridional structure (Jin et al. 1994, 1996).

As the various ENSO models each follow the quasiperiodic route to chaos, one can investigate the model results simultaneously by examining the relevance of each type of model solution to the observed ENSO synchronization. Quasiperiodic solutions do not reproduce the observed ENSO seasonal synchronization and frequency-locked solutions do not reproduce the observed ENSO irregularity, so realistic solutions must be either chaotic or frequency-locked solutions that are perturbed by high-frequency atmospheric forcing. Chaotic solutions are relatively rare compared to quasiperiodic and frequency-locked solutions (Jin et al. 1996), so the most likely realistic solution is a frequency-locked solution perturbed by stochastic noise. Moreover, both chaotic and stochastically forced solutions retain the *subharmonic* peaks in the ENSO spectra that are characteristic of the frequency-locked solutions (Jin et al. 1996). Thus, the results from the various studies can be fairly said to show that ENSO synchronization within the idealized models results from frequency locking of ENSO to the annual cycle, a mechanism that is characterized by subharmonic peaks in the ENSO spectrum.

Alternatively, several physical processes have been proposed whereby the annual cycle could affect the growth rate of ENSO anomalies, beginning with the idea that the seasonal movement of the intertropical convergence zone (ITCZ) should have a strong effect on the coupled instability of the equatorial Pacific ocean–atmosphere system because of its influence on atmospheric heating (Philander 1983). Analytical results based on a suite of four linear coupled models showed that the unstable modes allowed by the models were highly dependent on the parameterization of SST anomalies (SSTAs) and large-scale latent heating, and that the observed climatological background state did not permit the growth of ENSO-like instabilities (Hirst 1986). However, ENSO events could be initialized during more favorable conditions, including high SST, a shallow thermocline, a large zonal SST gradient, and strong surface winds (Hirst 1986). Isolating the effect of individual variables within the ZC model indicated that the seasonality in the wind divergence (Tziperman et al. 1997) and SST (Yan and Wu 2007) fields is most critical for the synchronization of ENSO events. Sensitivity analysis of a hybrid coupled model, used to capture the structure of the mixed layer and thermocline, found that the seasonal outcropping of the thermocline increased the coupled instability of the model by linking thermocline anomalies to the surface (Galanti et al. 2002). Lastly, at the end of the calendar year the location of ENSO-associated western Pacific wind anomalies shifts from along the equator to the Southern Hemisphere, which forces oceanic equatorial Kelvin waves that act to reduce or reverse the eastern equatorial Pacific SST anomalies (Harrison and Vecchi 1999). The shift in wind anomalies has been linked to the southward displacement of highest SST in boreal winter (Lengaigne et al. 2006), which is associated with increased convection and minimal surface momentum damping of wind anomalies (McGregor et al. 2012). The wind shifts have also been associated with a recently identified climate mode with energy at *combination tone* frequencies that emerges through an atmospheric nonlinear interaction between ENSO and the annual cycle (Stuecker et al. 2013).

The effect of the seasonal cycle on ENSO variance has been confirmed statistically by studies that examine the optimal perturbation growth around a seasonally varying background state within both ENSO model output and observations. For example, singular vector decomposition of a linearized version of the ZC model (Thompson and Battisti 2000), as well as the ZC forward tangent model along a trajectory in reduced EOF space (Xue et al. 1997), results in singular values that have a strong seasonal dependence, with growth of the singular vectors peaking in boreal winter. Similarly, cyclic Markov models derived from the ZC model (Pasmanter and Timmermann 2003), an anomaly coupled GCM (Kallummal and Kirtman 2008), and observations (Johnson et al. 2000) reveal a strong seasonality in the internal dynamics of the equatorial Pacific coupled ocean–atmosphere system. It has been suggested that this internal seasonality is sufficient to produce the observed the ENSO seasonal variance, without the need for nonlinear dynamics or seasonality in the noise forcing (Thompson and Battisti 2000; Kallummal and Kirtman 2008; Stein et al. 2010). Moreover, the Markov models can be explicitly related to Floquet analysis (Pasmanter and Timmermann 2003), which can be used to show that the dynamics of the most unstable mode of the ZC model with a seasonally varying background are the same as in the annual average case (Jin et al. 1996; Thompson and Battisti, 2000), which forms the basis of our dynamical understanding of ENSO (Philander et al. 1984; Hirst 1986; Neelin and Jin 1993a,b).

In this study, a parametric recharge oscillator (PRO) model of ENSO is employed to determine the effects of the seasonally varying background state instability on a variety of ENSO synchronization metrics. Analytical solutions for the model’s seasonal variance, amplitude modulation, and phase synchronization are obtained, and are shown to match well with observations and the variety of ENSO behavior identified in state-of-the-art coupled general circulation models (CGCMs). The parametric model is also shown to explain spectral peaks at combination tone frequencies that are present in both the observed and modeled ENSO spectra. Additionally, the ENSO recharge oscillator model is extended to include external periodic forcing by the annual cycle and a cubic damping term, in order to produce frequency-locked model solutions. The extended model corresponds to the well-known van der Pol oscillator (van der Pol 1927), which exhibits frequency locking through so-called Arnol’d tongues (Arnol’d et al. 1983), with global behavior similar to the ENSO models used in previous studies. The van der Pol oscillator is used to examine the relevance of the frequency locking scenario to ENSO seasonal synchronization, ultimately demonstrating that the observed ENSO lacks the characteristics of a frequency-locked oscillation.

The remainder of the paper is organized as follows: Section 2 discusses features of ENSO synchronization, including seasonal variance, amplitude modulation, phase synchronization, and secondary peaks in the ENSO spectrum. These features are examined both for the observations and within a range of state-of-the-art coupled general circulation models. Section 3 discusses the analytical solution of a simplified neutrally stable version of the parametric recharge oscillator following the perturbation expansion of An and Jin (2011). It is demonstrated that the observed features of ENSO–annual cycle interactions arise directly from this solution. Numerical experiments with the parametric recharge oscillator are presented in section 4, and the solutions within the relevant subset of the model’s parameter space are discussed. A derivation of the van der Pol oscillator from the extended oscillator model is shown in section 5, which is used to examine the likelihood that ENSO synchronization is due to subharmonic frequency locking. The observed ENSO is shown to be lacking the characteristics of a frequency-locked oscillation. A discussion of the implications and caveats concerning the work herein is presented in section 6, and the paper concludes with a summary of major results in section 7.

## 2. Features of ENSO seasonal synchronization

To test the two leading theories of ENSO synchronization, it is necessary to construct a detailed picture of annual cycle–ENSO interactions using a variety of different metrics. This study examines the synchronization theories within the framework of a simple oscillator model; to allow for direct comparison with model results, the synchronization metrics are calculated from the detrended Niño-3.4 index time series, defined as the area average of SSTA within 5°N–5°S, 120°–170°W. The metrics were also calculated using time series of ENSO derived from complex empirical orthogonal functions (CEOFs; not shown), and the results are very similar due to the fact that ENSO anomalies are large-scale and occur in phase across the central Pacific (Stein et al. 2011).

Figure 1 (top, black) shows the monthly Niño-3.4 SSTA index based on data from the National Oceanic and Atmospheric Administration’s (NOAA’s) Extended Reconstructed Sea Surface Temperature version 3b dataset (ERSST.v3b; Smith et al. 2008) from the year 1950 to 2011, with the monthly mean climatology (bottom left) and long-term trend removed. The Niño-3.4 SSTA index *T*(*t*), can be expressed as a cyclostationary process *T*(*y*, *m*), where *y* is the year, *m* the month, and *T*(*y*, *m* + 12) = *T*(*y* + 1, *m*). The monthly variance is then determined as , where *E* indicates the expected value and we have made use of the fact that the monthly means have been removed from the time series. Figure 1 (bottom middle) shows the monthly variance of the Niño-3.4 SSTA index, which has minimum variance in March–April and highest variance in December, reflecting the tendency for ENSO events to peak in boreal winter.

The seasonally modulated variance is an expression of either a phase synchronization of ENSO with the annual cycle (Stein et al. 2011) or a seasonal modulation of ENSO’s amplitude, or a combination of the two. To separate the processes, it is necessary to construct a state space that allows for the definition of ENSO magnitude and phase. A well-known method for defining amplitude and phase from a dataset, and the one most naturally applied to climate data, is to construct the analytical signal (Gabor 1946) of the data using the Hilbert transform (Pikovksy et al. 2000). The analytical signal of the Niño-3.4 SSTA index is defined as

where is the Hilbert transform of *T*, and . The top panel of Fig. 1 shows the time series of the Hilbert transform of the Niño-3.4 SSTA index (gray dashed). The amplitude and phase of the index can be calculated from the complex analytical signal , based on a Cartesian to polar coordinate transform

where Arg is the principal value of the arg function of complex numbers, defined such that . For consistency, all phase values considered here will be calculated modulo 2*π*, and therefore on the interval [0, 2*π*).

The monthly mean amplitude of the analytical signal of the Niño-3.4 SSTA index, *α*_{m} = *E*[*α*(*y*, *m*)], is compared to the monthly variance of *T*(*t*) in Fig. 1 (bottom center). There is an indication of amplitude modulation of the complex signal by the seasonal cycle, with maximum mean amplitude of the complex ENSO signal occurring in boreal winter, but with a minimum in September, which does not correspond to the minimum in monthly variance of the Niño-3.4 SSTA index. Moreover, the overall seasonal amplitude modulation of analytical signal is insufficient to reproduce the observed seasonal variance of *T*(*t*).

To investigate the phase synchronization of ENSO with the annual cycle, we define the generalized phase difference

where *ϕ*_{e} is the ENSO phase^{1} [Eq. (3)], *ϕ*_{a} is the annual cycle phase, and *k*, *l* ∈ □^{+}. If the phase difference were constant the annual cycle and ENSO would be perfectly synchronized, and if the values of *δϕ*_{k,l}(*t*) were equally distributed throughout the 2*π* range then the two modes would have no phase relationship. Phase locking is defined as a bounded *δϕ*_{k,l} difference, that is, if |*δϕ*_{k,l}(*t*) − *s*| < *c*, where *c* < 2*π* is a constant and *s* is the average phase shift between the two time series (Pikovksy et al. 2000). Because the monthly mean climatology was removed from the Niño-3.4 SSTA index, the annual cycle is fixed, and we will define the phase of the annual cycle to be *ϕ*_{a}(*t*) = (*ω*_{a}*t*)mod(2*π*), where is the annual cycle angular frequency.

Figure 1 (bottom right) shows a histogram of *δϕ*_{2,1}, indicating the strength of the 2:1 phase synchronization of ENSO to the annual cycle, where the notation 2:1 indicates that the average phase progression of the annual cycle proceeds at twice the rate of ENSO. The 2:1 ratio was the only one found to show significant phase synchronization between ENSO and the annual cycle (see section 5). As can be seen in Fig. 1, the phase differences between the two signals are not bounded, indicating that ENSO is not strictly phase-locked to the annual cycle throughout the observable record. However, certain phase differences are 3 times more likely than others, which is evidence of partial 2:1 phase synchronization of ENSO to the annual cycle.

In terms of existing ENSO theory, such behavior could be explained by intermittent periods of frequency locking of ENSO to the annual cycle. However, the 2:1 phase synchronization of ENSO to the annual cycle is not associated with a distinct peak in the Niño-3.4 spectrum at 0.5 yr^{−1} (Fig. 2), as would be expected if the phase synchronization behavior was due to frequency locking^{2} (see section 5). Instead, the primary ENSO spectral peak occurs at 0.23 yr^{−1}, along with a second secondary peaks at 0.77 and 1.23 yr^{−1}, which occur at combination tones frequencies *ω*_{a} ± *ω*_{e}, although the higher-frequency peak is not statistically significant based on the spectrum from a first-order autoregressive fit of the Niño-3.4 time series (Fig. 2).^{3} The combination tones are indicative of nonlinear interaction between ENSO the annual cycle (Stuecker et al. 2013) and arise directly from the modulation of ENSO growth rate within the parametric recharge oscillator model (An and Jin 2011; see section 3).

To demonstrate that the above features of ENSO’s seasonal synchronization are robust, we examine the synchronization of a variety of different model-based representations of ENSO and the annual cycle, as simulated by state-of-the-art coupled general circulation models participating in the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al. 2012). The detrended Niño-3.4 SST index was calculated from CMIP5 historical run output for each model, covering model years 1901–2000. The annual cycles of the Niño-3.4 SST index for each model are compared in Fig. 3, the monthly variance of the Niño-3.4 SSTA indices are compared in Fig. 4, and the *δϕ*_{2,1} phase differences with the annual cycle are compared in Fig. 5. The models have been ordered according to the strength of the seasonal synchronization of ENSO to the annual cycle (as measured by the *υ* index, defined below).

As can be seen in Fig. 3, the CMIP5 CGCMs show a variety of differing amplitudes for the annual cycle of the Niño-3.4 SSTA index. The annual cycles for all the models are more symmetric than the observations (Fig. 1, bottom left), showing a warming for six months out of the year rather than four, but the timing of the warm to cold phase transitions for the models is in agreement with the observations. In particular, there are no semiannual seasonal cycles, allowing for the use of a single sinusoidal function at the annual frequency as the proxy for the annual cycle in the models, as with the observations. The models also display a variety of behavior in terms of overall ENSO variance and the strength of the seasonal modulation of ENSO variance (Fig. 4). For most models, total ENSO variance and the magnitude of the seasonal modulation of ENSO variance are both weaker than observations, which is reflected in the flatter distributions of phase differences between ENSO and the annual cycle (Fig. 5). The seasonal variance of each model ENSO (Fig. 4) shows a peak in December/January, and so the most likely *δϕ*_{1,2} phase differences for each model are all near the same value (Fig. 5).

To examine any systematic relationship between synchronization characteristics across all the CMIP5 CGCMs, we define three indices relating to ENSO’s seasonal variance, 2:1 phase synchronization, and amplitude modulation of the corresponding complex analytical signal. The first index simply captures the range of the monthly ENSO variance,

and will be referred to as the seasonal variance index. The index ranges from 0 to 1, where 0 indicates that each month has exactly the same variance throughout the year and 1 indicates that the variance of a particular month drops to zero.

An index for the strength of the 2:1 phase synchronization between ENSO and the annual cycle can be defined as

where *ϕ*_{a} is the annual cycle phase, *ϕ*_{e} is the ENSO phase, and 〈⋅〉_{t} indicates temporal averaging (Kralemann et al. 2008). The index is a measure of the length of the vector in the complex plane that results from the temporal averaging of unit vectors with angles equal to the phase difference *δϕ*_{2,1}. The index varies from 0 to 1, with 0 indicating that the phase of the two time series are completely independent and 1 indicating perfect phase synchronization. Similarly, one can define an index that captures the strength of the complex amplitude modulation of ENSO by the annual cycle as

where *α* is the amplitude of the complex ENSO signal [Eq. (2)]. The index is the temporal average of the vectors in the complex plane defined by the ENSO amplitude and annual cycle phase, which is then normalized by the mean ENSO amplitude. The index varies from 0 to 1, with 0 indicating that the complex ENSO amplitude is equal across all phases of the annual cycle, and 1 indicating that ENSO’s analytical signal only has a finite amplitude at a single time of the year, and has zero amplitude at all other times. As such, the values of the Ψ index would be expected to be much smaller than the values of the other two indices, and these indices should not be directly compared, such that a value of *χ* that is larger than Ψ is interpreted as stronger phase synchronization than amplitude modulation. Rather, the indices measure the changing strength of these processes across the CMIP5 models.

Figure 6 shows scatterplots of the seasonal variance index (*υ*) versus the complex amplitude modulation index (Ψ, top) and the 2:1 phase synchronization index (*χ*, bottom) for each of the CMIP5 models and the observations. Most models simulate ENSO variability that is more weakly synchronized with the annual cycle than in the observations. The seasonal variance of ENSO (*υ*) is linearly related to the strength of the phase synchronization of ENSO to the annual cycle (*χ*), while the strength of the modulation of the complex ENSO signal (Ψ) is not as closely related, though the Ψ index does tend to be larger for models with a larger range of seasonal ENSO variance. The models that display the strongest synchronization of ENSO to the annual cycle also display peaks at one or both of the *ω*_{a} ± *ω*_{e} combination tone frequencies, as seen in Fig. 7. Five of the first six models (CNRM-CM5, FGOALS-g2, CanESM2, BCC-CSM1.1, GISS-E2-R; see expansions of model names in the appendix) show a significant peak at the *ω*_{a} − *ω*_{e} combination tone, and all but three of the models show a spectral peak at at least one combination tone, although the peaks are statistically significant in only 9 of the 16 models (CNRM-CM5, FGOALS-g2, CanESM2, BCC-CSM1.1, Giss-E2-R, NorESM1-ME, Giss-E2-H, GFDL-ESM2M, CCSM4). Four of the 16 models (CNRM-CM5, FGOALS-g2, CanESM2, GISS-E2-R) also display a spectral peak at the 2-yr period, but each of these models also shows a peak at a combination tone frequency. It is unclear at this point whether the spectral peak at periods of 2 yr in these models is due to frequency locking, parametric resonance, or another mechanism. However, it is apparent that overall ENSO seasonal synchronization is more often associated with peaks in the ENSO spectrum at combination tone frequencies, as opposed to a peak at ½ yr^{−1}. In the next section, we will show how the *ω*_{a} ± *ω*_{e} spectral peaks, along with ENSO’s seasonal variance, amplitude modulation, and partial 2:1 phase synchronization, all arise directly from the modulation of ENSO’s growth rate by the annual cycle within the parametric ENSO model.

## 3. Amplitude and phase modulation of the parametric recharge oscillator

The parametric recharge oscillator model is based upon the recharge oscillator model of ENSO (Jin 1997) derived with seasonal varying coefficients (Stein et al. 2010), taking the form of a stochastic parametric oscillator

where *T* represents eastern equatorial Pacific SST anomalies, *H* represents the zonal mean equatorial Pacific thermocline depth anomalies, *λ*(*t*) and *ω*(*t*) are the seasonally varying growth rate and angular frequency parameters of the oscillator, the constant *R* relates to the time scale of the geostrophic adjustment of the thermocline to wind stress anomalies, and *ξ*(*t*) is Gaussian white noise representing forcing by the atmosphere. The parameters *λ* and *ω* can be derived from a statistical-dynamical estimation of the linearized upper ocean heat budget based on the Bjerknes index (Jin et al. 2006), and model runs utilizing the estimated growth parameter reproduce the observed seasonal cycle of ENSO variance (Stein et al. 2010).

Modulation of the angular frequency parameter (*ω*) was shown to have little effect on the seasonal variance of ENSO (Stein et al. 2010), so the model can be further reduced to

An approximate analytical solution was obtained for the neutrally stable, unforced case of the PRO model, where *ξ*(*t*) = 0, *λ*(*t*) = *λ*_{0} cos(*ω*_{a}*t*), and *ω*_{a} is the annual frequency in Eq. (10). Setting *λ*_{0} = *ϵ*, 0 < *ϵ* ≪ *ω*_{a}, and using the perturbation method, the first-order approximation to the solution for *T* can be obtained (An and Jin 2011) as

where *C* is a constant dependent on the initial conditions, , and . Without loss of generality, we can set *C* = 1. Note that due to the denominators of *A* and *B*, *T* → ∞ as , which reflects parametric resonance in the system. For all values of , the parameter *B* is larger than *A*, with the ratio of for realistic ENSO periods.

An example of the *T* and [*T*] time series, seasonal variance, and amplitude and phase for Eq. (12) is shown in Fig. 8. For this example, the ENSO period was set to 3.75 yr and the amplitude of the seasonal modulation was set to *ϵ* = 2 yr^{−1} in order to match the range of the observed monthly Niño-3.4 SSTA index variance. Although there is an arbitrary phase shift in terms of the time of year that the modulation occurs, the equation can reproduce the salient aspects of ENSO seasonal synchronization. The ranges of the seasonal variance, amplitude modulation, and the PDF of the phase difference with the annual cycle agree well with the observations (Fig. 1), especially considering the model’s simplicity. Utilizing Eq. (12), analytical solutions can be obtained for the Hilbert transform of *T*, as well as the seasonal variance and amplitude modulation of the *T* time series. Additionally, it can be shown that the ENSO signal expressed in Eq. (12) will have a preferred phase difference with the annual cycle. This analysis will demonstrate that the various aspects of ENSO synchronization can result directly from the modulation of ENSO’s growth rate by the annual cycle.

First, to calculate the seasonal variance of *T*, write the time series as the cyclostationary process *T*(*y*, *m*), with a monthly time series chosen for consistency with the observational data and coupled model output used in this study. The seasonal variance is then

where *E* indicates the expected value and we have made use of the fact that *E*(*T*) = 0. For the relevant ENSO time scale, the *A*^{2} term in Eq. (13) is negligible. With a long enough sampling time, the expected value of *ω*_{e} terms will be constant with respect to the annual cycle and we can set

resulting in

Thus, the variance of ENSO will be periodic with respect to the annual cycle, with this periodicity arising from the second cross term between the annual cycle and ENSO in Eq. (12). The third term in Eq. (14) represents a semiannual contribution to the seasonally modulated variance.

To calculate the Hilbert transform of the *T* time series, first note that Eq. (12) is equivalent to

This form of the analytical solution shows explicitly that modulation of ENSO by the annual cycle will result in energy at the combination tone frequencies *ω*_{a} ± *ω*_{e}, which account for the secondary peaks in the observed and modeled ENSO spectra. In particular, the approximate solution [Eq. (15)] explains the dependence of the frequencies of the secondary peaks on the mean ENSO frequency, and shows that *ω*_{a} − *ω*_{e} component of the time series will be larger than the *ω*_{a} + *ω*_{e} component. Both of these characteristics may be seen in the observed Niño-3.4 SSTA spectrum (Fig. 2).

We then utilize the fact that the Hilbert transform is a linear operator, and that the Hilbert transform of the cosine and sine functions are known, namely [cos(*x*)] = sin(*x*) and [sin(*x*)] = −cos(*x*). The Hilbert transform of *T* is then

The seasonal variance of the Hilbert transform of *T* can also be calculated as it was for the *T* time series itself,

This periodic variance of the *T* and [*T*] time series are related to the seasonal amplitude modulation. The amplitude of the analytical signal of *T* is defined as

and the seasonal amplitude modulation is defined by taking the cyclostationary mean of the amplitude, namely

which can be solved readily to first order using Eqs. (12) and (16). Because the complex amplitude modulation is proportional to the difference of the *A* and *B* parameters, the signal is relatively small.

Additionally, the modulation of the growth parameter in Eq. (10) results in a preferred phase difference between ENSO and the annual cycle. The result may seem counterintuitive, but the reason that modulation of the growth rate parameter in the PRO model can lead to phase synchronization is that the modulation of the growth parameter occurs on a shorter time scale than the intrinsic ENSO period. This is in contrast to a standard amplitude modulated (AM) signal, in which the modulation occurs at a much lower frequency than the carrier signal, and therefore no phase synchronization can result.

The phase of ENSO is defined as

where the tangent half-angle formula is used to calculate the principal value of the argument in a uniform matter. Neglecting the *A* terms yields

which is form of the solution that was used for calculating the analytical phase difference *δϕ*_{2,1} in Figs. 8–10. However, a simplified version of the equation is more helpful in demonstrating conceptually the relationship between ENSO seasonal variance and phase synchronization with the annual cycle. As shown before, the seasonal amplitude modulation is small, so the value of *α*(*t*) can be fixed at 1. Neglecting the *ϵ* term in the numerator then yields

Solving Eq. (14) for *ϵ**B* and inserting the solution into Eq. (23) results in a solution for ENSO phase modulation in terms of the seasonal variance:

where we have made use of the fact that sin(*ω*_{a}*t*) = sin(*ω*_{a}*m*) by definition.

With no modulation, the phase of ENSO would proceed monotonically as expected, *ϕ*_{e}(*t*) = *ω*_{e}*t*. The stronger the modulation of the growth rate parameter (*ϵ*), the stronger the seasonal modulation of ENSO variance [Eq. (14)], and therefore the stronger the modulation of the ENSO phase progression [Eq. (24)], which leads to the phase synchronization with the annual cycle. Equation (24) thus explains the linear relationship between the magnitude of the seasonal cycle of ENSO variance and the strength of the phase synchronization with the annual cycle that occurs across CMIP5 models and observations (Fig. 6).

The analysis of the PRO model shows that the annual modulation of the growth rate parameter can account for all the observed features of ENSO synchronization: seasonal variance, amplitude modulation, phase synchronization, and combination tones in the ENSO spectrum. The secondary spectral peaks at combination tones between ENSO and the annual cycle result directly from the modulation of ENSO’s growth rate by the annual cycle, and are seen explicitly in the first-order solution to the system [Eq. (15)]. Analytical solutions for the seasonal variance [Eq. (14)], amplitude modulation [Eq. (19)], and phase synchronization [Eq. (22)] of ENSO are also obtainable, with low-order forms of the solutions agreeing well with the statistical analysis (Fig. 8). The simplified Eqs. (23) and (24) for ENSO phase capture the essential mechanism that produces the preferred 2:1 phase difference with the annual cycle and explicitly relates ENSO’s phase synchronization to its seasonal variance. The counterintuitive result that modulation of the growth parameter can result in a preferred phase difference is due to the fact that the modulation of ENSO occurs at a higher frequency than ENSO itself. In such a system, amplitude modulation and phase synchronization are not independent, but result directly from the annual cycle modulation, and are dependent on the two free parameters of the system: the strength of the modulation (*ϵ*) and ENSO’s intrinsic frequency (*ω*_{e}). In the next section, the simplified analytical equations are shown to apply to more realistic numerical runs of the PRO model, supporting the proposition that the equations capture the salient features of ENSO synchronization in the natural system.

## 4. Numerical confirmation of the analytical solutions

The preceding analysis of ENSO synchronization is based on a first-order perturbation expansion to a particular case of the PRO model, namely the neutrally stable system with no stochastic forcing. To examine the relevance of these solutions, numerical experiments of the full PRO model [Eqs. (10) and (11)] were performed. For these experiments, the ENSO period was set to 3.75 yr, the amplitude of the modulation was set to *ϵ* = 2 yr^{−1}, a daily time step was used with stochastic forcing applied at each time step, and monthly averages of the *T* time series were saved as output. For simplicity, the model year consists of 360 days, 12 months of 30 days each.

First, the analytical solution was compared to a numerical integration of the neutral, undamped PRO model. The model was initialized with values of *T* = 0, *H* = 1.35, resulting in a variance for the *T* time series of 1.029°C. Figure 9 shows the *T* and [*T*] time series, the seasonal variance of the *T* time series, the amplitude modulation of the analytical signal, and a histogram of the *δϕ*_{1,2} phase difference between the analytical signal and the annual cycle, from a 60-yr integration of Eqs. (10) and (11) with *ξ*(*t*) = 0. Analytical solutions for the seasonal variance and amplitude modulation are shown with dashed lines for comparison. The numerical and analytical solutions for the time series agree well, and the analytical solutions of the amplitude modulation, seasonal variance, and phase difference are very good approximations to the values computed statistically from the analytical signal of the *T* time series, justifying the simplifications made in obtaining those solutions.

The analytical solution was then compared to an integration of the damped, stochastically forced PRO model. The damping parameter was set to , with the mean damping rate set to based upon a statistical-dynamical fit of the PRO model to output from a high-resolution ocean GCM reanalysis (Stein et al. 2010). Gaussian stochastic noise with a variance of 0.04°C was applied at each time step of one model day. Figure 10 shows the synchronization statistics from a 100-member ensemble of 150-yr runs of the stochastically forced model, with an example of the time series output above. For each ensemble member, the model was run for 200 yr, with the last 150 yr of output used for analysis. The ensemble mean monthly variance, monthly amplitude, and *δϕ*_{1,2} phase differences are shown, along with the 90% confidence intervals for each statistic. Although the time series of model output for the two cases are quite different, the ensemble mean statistics of the stochastic case are very similar to the damped, unforced case and agree with the analytical solutions obtained in the previous section. Thus, the analytical solutions obtained for the neutral, unforced case of the PRO model are valid for the long-term behavior of the damped, stochastically forced case. The stochastic forcing adds a large degree of uncertainty about the statistics of a particular ensemble member. Notably, the seasonal variance and preferred phase difference that result from the annual cycle modulation will be apparent for each ensemble member, while the same cannot be said for the amplitude modulation, which is characterized by a very small signal-to-noise ratio for the stochastically forced case. For time series of ≤150 yr, the stochastic forcing may mask the underlying amplitude modulation, but the seasonal variance and preferred phase difference that result from the annual modulation of ENSO will still be apparent. These results explain the weak correspondence between seasonal variance and amplitude modulation across the CMIP5 models (Fig. 6), as reflected in Eq. (19).

The analytical solutions presented in the previous section thus capture the essential aspects of ENSO synchronization in the PRO model for both the neutral case and for the forced, damped case where the mean damping rate is based on a fit to reanalysis data (Stein et al. 2010). We can therefore examine the dependence of ENSO synchronization statistics on the degree of annual cycle modulation and intrinsic ENSO frequency using Eq. (12). Figure 11 shows contours of the three synchronization indices across a section of the PRO model parameter space. For ENSO modes with periods longer than three years, the index contours are nearly horizontal, indicating that the strength of all the synchronization indices is determined almost entirely by the strength of the annual cycle modulation. In this regime, ENSO seasonal variance, amplitude modulation, and phase synchronization increase nearly monotonically with increasing modulation of the growth rate parameter. Model ENSOs with a period of <3 yr move toward a regime where increased modulation of the growth rate parameter results in increased seasonal ENSO variance and phase synchronization, but not an increase in amplitude modulation of the complex analytical signal. In this regime, the model becomes parametrically destabilized, and nonlinear saturation is necessary to limit overall ENSO variance.

Using Fig. 11, one can predict the strength of phase synchronization of ENSO to the annual cycle, given the mean ENSO period and seasonal variance, which provide an estimate of the strength of the annual cycle modulation [i.e., the range of *λ*(*t*)]. This calculation was done for the observations and the CMIP5 model ENSOs, using the *υ* index to derive the strength of the seasonal variance and defining the mean ENSO frequency using the phase time series:

An estimate of the strength of the modulation by the annual cycle based on the analytical solution to the PRO model can be obtained using *υ* and *ω*_{e} (Fig. 11, top), and the values of *ω*_{e} and *υ* then correspond to a predicted phase synchronization strength *χ* (Fig. 11, middle). The predicted and measured strengths of the 2:1 phase synchronization are compared in Fig. 12. The two values agree well, with the value predicted by the analytical PRO model accounting for 75% of the variance in the strength of the phase synchronization in the CMIP5 models. The agreement indicates that the measured 2:1 phase synchronization arises from the modulation of ENSO’s growth rate by the annual cycle, and provides further evidence that the analytical solutions of the PRO model capture the essential features of ENSO synchronization in the observations and CMIP5 models.

## 5. ENSO–annual cycle synchronization through frequency locking

To compare the above results to the frequency locking scenario for ENSO, Eqs. (8) and (9) are extended to include periodic forcing and a nonlinear damping term. The damping rate *λ* is held constant with respect to the annual cycle, and the behavior of the model is examined in the unstable regime above the Hopf bifurcation at *λ* = 0. The extended model is not offered as a particularly realistic representation of ENSO, but rather as a system that has similar global behavior to the simple (Tziperman et al. 1994; Wang and Fang 1996; Liu 2002) and intermediate (Chang et al. 1994; Tziperman et al. 1995; Jin et al. 1996) models of ENSO on which the frequency locking scenario is based, and which admits the same suite of possible model solutions.

The extended model reads as follows:

where the sinusoidal forcing term represents external forcing by the annual cycle that is independent of eastern Pacific SST anomalies (Liu 2002) and the cubic damping term represents nonlinear saturation of SST anomalies due to, for example, the nonlinear dependence of subsurface temperature on thermocline depth (Jin 1997).

Equations (25) and (26) can be reduced to the well-known forced van der Pol oscillator (van der Pol 1927) with two simplifications. First, rescale the equations by *τ* = *ω*_{e}*t*, and second, fix the strength of the cubic damping to , resulting in

where , , , and . By setting , the growth rate and cubic damping terms have been combined into the nonlinear damping term , which also controls the intrinsic frequency of the unforced oscillation. The periodically forced van der Pol oscillator [Eq. (27)] has been thoroughly studied [see Mettin et al. (1993) for a complete description of the bifurcation structure], so we will limit our discussion here to the types of solutions that are allowed by the system and their applicability to the seasonal synchronization of ENSO. We begin by providing examples of each type of model solution, using a value of because of the known chaotic solutions at that driving frequency (Mettin et al. 1993), and then proceed to examine the solutions that occur within the parameter space relevant to ENSO.

In the absence of external forcing, the solution of the self-sustained oscillator follows a stable limit cycle with an amplitude of 2°C (Fig. 13a), which is comparable to the amplitude of the Niño-3.4 index (Fig. 1, top), so the choice to fix the value of the cubic damping parameter is reasonable. As the nonlinear damping parameter becomes larger, the period of oscillation increases and the system moves toward a relaxation oscillation (Fig. 13b). With the inclusion of periodic forcing, quasiperiodic solutions exist for relatively small values of the nonlinear damping (Fig. 13c). Further increasing the nonlinearity, results in the oscillator becoming frequency-locked to a rational multiple of the driving frequency (Fig. 13d). With strong enough forcing, chaotic oscillations can occur when several frequency-locked solutions overlap, and the oscillation jumps irregularly between the possible resonant frequencies (Fig. 13e). The chaotic solutions observed in the driven van der Pol oscillator result from period doubling bifurcations (Mettin et al. 1993), a subcase of the quasiperiodic route to chaos (Jin et al. 1996) followed by ENSO models used in previous studies.

To examine the model solutions within a realistic parameter range, we fix the value of the forcing *F* in Eq. (25) and then examine model solutions with various values of the growth rate *λ* and neutral ENSO period . The range of the annual cycle of the Niño-3.4 index provides a reasonable upper limit of 1.2°C for the annual cycle forcing *F* in Eq. (25) (Liu 2002). The neutral ENSO period *T*_{e} is varied from 2 to 5 yr, and the *λ* parameter is varied from 0 to 6, corresponding to a maximum growth rate of 3°C yr^{−1}. The annual cycle frequency is fixed at *ω*_{a} = 2*π*, and hence the driving frequency is equal to the mean ENSO period *T*_{e}.

Figure 14 displays contours of the ratio of the frequency of the *T* time series output from Eq. (27) to the driving frequency , based on runs with rescaled parameters , , and calculated from *F*, *λ*, and *ω*_{e}. Regions of quasiperiodic solutions are separated by frequency-locked solutions that occur within Arnol’d tongues (Arnol’d et al. 1983), which become larger with increasing growth rate, and hence nonlinear damping . The Arnol’d tongues slope up and to the left because the shorter the neutral ENSO period *T*_{e}, the weaker the effective forcing for a fixed *F* and *ω*_{a}. As can be seen, the driven van der Pol oscillator favors frequency-locked solutions at odd multiples of the driving frequency; frequency-locked solutions at even multiples of the driving frequency are very unstable (Mettin et al. 1993). No chaotic solutions were found in the neighborhood of realistic parameter ranges, reflecting the rarity of chaotic solutions within more complex ENSO models (Jin et al. 2006). We will therefore only compare the observed ENSO synchronization characteristics to the probable model solutions: quasiperiodic oscillations, frequency-locked oscillations, and frequency-locked oscillations with added stochastic forcing.

Figure 15 shows an example of a quasiperiodic oscillation (top) obtained by setting *T*_{e} = 2.76 and *λ* = 0.5, along with the corresponding spectrum (bottom left). The strength of the *k*: *l* phase synchronization between the *T* time series and the periodic forcing for *k*, *l* ∈ [1, 10] is shown at the bottom right, defined as

For a quasiperiodic solution, the ENSO frequency is still largely determined by the nonlinear damping rate , independent of the forcing frequency, so the oscillation will not have a seasonal cycle of variance or display any phase synchronization with the periodic forcing (Fig. 15, bottom right). Quasiperiodic solutions therefore produce no synchronization between ENSO and the annual cycle.

Increasing the growth rate to *λ* = 2 moves the system into the 3:1 Arnol’d tongue (Fig. 14), resulting in a regular 3-yr ENSO cycle as shown in the spectrum (Fig. 16, bottom left). Such subharmonic frequency locking will by definition produce phase synchronization as well, which can be seen most strongly at 3:1, but also at the other 3*n*: *n* ratios (Fig. 16, bottom right). The weaker phase synchronization at other ratios results from the fact that the oscillation is not perfectly sinusoidal. The frequency-locked oscillation is far more regular than the observed ENSO, so stochastic forcing is necessary to make the solution more realistic (Fig. 17). The inclusion of the stochastic forcing broadens the spectrum of the *T* output, and weakens the phase synchronization between the *T* time series and the periodic forcing, particular for larger values of *k*, *l*. The only remaining phase synchronization occurs at the ratio 3:1, associated with the peak in the spectrum at the 3-yr period.

To examine the possibility that the observed ENSO results from a stochastically perturbed frequency-locked oscillation, we compared the above synchronization characteristics to those of the observations. Figure 18 shows the phase synchronization indices for *k*, *l* ∈ [1, 10] as calculated from the observed Niño-3.4 index. As mentioned before, the 2:1 phase synchronization is by far the strongest in the observations, and there is no indication of particularly strong synchronization at any other ratio. The analysis of the CMIP5 CGCMs yields the same result (not shown). However, the 2:1 phase synchronization is not accompanied by a spectral peak at the biennial frequency in either the observations (Fig. 2) or the majority of CMIP5 GCMs (Fig. 7), as would be the case if the phase synchronization were due to frequency locking. This observation, combined with the ability of the PRO model to reproduce seasonal ENSO variance, amplitude modulation, and phase synchronization, suggests that the seasonal synchronization of ENSO is due to the periodic modulation of ENSO’s growth rate by the annual cycle, rather than frequency locking of ENSO to periodic forcing by the annual cycle.

## 6. Discussion

In addition to the PRO model’s ability to capture the various features of ENSO synchronization, the results presented here complement other recent observational studies on annual cycle–ENSO interactions. A newly identified climate mode, corresponding to the second principal component of an EOF decomposition of western Pacific winds, has been shown to have energy at the *ω*_{a} − *ω*_{e} combination tone frequency (Stuecker et al. 2013). The spatial pattern of the mode captures the development of the Philippine Sea anticyclone (Wang et al. 1999), as well as the southward shift of westerly wind anomalies in the western Pacific (McGregor et al. 2012) during El Niño events. These two features have been proposed as possible mechanisms for ENSO phase transition via the modulation of the atmospheric response to SST anomalies and thereby ENSO’s growth rate, captured heuristically in the PRO model via the modulation of *λ*. In turn, the modulation of the growth rate (*λ*) in the PRO model results directly in energy at combination tone frequencies [Eq. (15)] due to the nonlinear interaction of ENSO and the annual cycle. It is important to note that the combination tone is proposed to be due to the interaction of El Niño events with the annual cycle, so there is an apparent asymmetry in the observed system not accounted for in the PRO model.

The combination of the analytical PRO model results and the observation of the combination climate mode suggests that the most important ENSO–annual cycle interaction occurs in the atmosphere over the western Pacific, which helps to narrow down the myriad possible physical interactions that could result in the modulation of ENSO’s growth rate by the annual cycle. This idea could be tested by examining characteristics of the various CMIP5 CGCM simulations of the western Pacific and relating them systematically to the strength of the synchronization of ENSO to the annual cycle within each model (Fig. 6). The analysis of the CMIP5 output suggests that there is a universal scaling of the models according to the strength of ENSO synchronization that might be related to, for example, biases in the location and strength of western Pacific atmospheric convection, SST, surface winds, etc. If such a connection were made, it would useful as a way to evaluate and improve coupled general circulation models and provide insight into the physics of ENSO–annual cycle interaction.

Although the PRO model appears to capture the salient features of ENSO synchronization as presented here, and is further supported by the observed ENSO–annual cycle combination climate mode, the model has essential limitations that should be discussed. Throughout this work we have considered only the two leading theories for ENSO synchronization to the annual cycle: that the synchronization results from the seasonal modulation of ENSO’s growth rate or from frequency locking of ENSO to the periodic forcing by the annual cycle. These two theories correspond to limiting cases of the following model:

which includes the seasonal variation of ENSO’s growth rate, nonlinear saturation, periodic forcing by the annual cycle, and stochastic atmospheric forcing at fast time scales. The effect of the seasonal modulation of ENSO’s growth rate was represented by the parametric recharge oscillator model of ENSO, which corresponds to the limiting cases of *c* = *F* = 0. Frequency locking of ENSO to the annual cycle was captured by the limiting case .

Besides these two limiting cases, there are other possible scenarios that are outside the scope of this study but which may warrant further investigation. For example, one could consider the interaction of the seasonally varying growth rate [*λ*(*t*)] and the nonlinear saturation term (−*cT*^{3}) in the system, or the effect of seasonally varying (Hendon et al. 2007) or state-dependent stochastic noise forcing (Kessler et al. 1995; Kessler and Kleeman 2000; Jin et al. 2007; Levine and Jin 2010). The importance of the various interactions within the full system [Eqs. (29) and (30)] would be dependent on the balance between the mean damping/growth rate , the strength of the nonlinearity (*c*), the strength of the periodic forcing (*F*), and the variance of the stochastic forcing [*ξ*(*t*)]. Clearly, the complexity increases rapidly even within the simplest model representation of ENSO, and the increased complexity must be justified by an increased understanding of ENSO or ability to explain features of ENSO beyond the stochastically forced model with seasonally varying growth rate. Of the possible avenues of investigation, the inclusion of state-dependent (multiplicative) noise forcing (Levine and Jin 2010) in the system seems the most fruitful, as it can reproduce ENSO skewness (Burgers and Stephenson 1999; An and Jin 2004), a feature that is not present in the symmetric models presented here.

## 7. Summary and conclusions

The two leading theories of ENSO seasonal synchronization are that ENSO is frequency-locked to periodic forcing by the annual cycle (Jin et al. 1994; Tziperman et al. 1994) or that the seasonal modulation of ENSO’s growth rate leads to the synchronization of ENSO events (Philander et al. 1984; Hirst 1986). These two theories were tested here within a consistent model framework based on the recharge oscillator model of ENSO (Jin 1997). The observed ENSO synchronization, as well as the synchronization simulated by state-of-the-art CGCMs, was described in terms of ENSO’s seasonal variance, amplitude modulation, phase synchronization, and secondary peaks in the spectra at “combination tones” with the annual cycle. These synchronization metrics were then compared to those of the parametric recharge oscillator (PRO) model of ENSO, which captures in the simplest form the seasonal modulation of ENSO’s growth by the annual cycle. Analytical results of the neutral, unforced case of the PRO model show that the annual modulation of the growth rate parameter results directly in all of the synchronization features, and the strength of each of the synchronization features agree well with those of the observations. The analytical solutions were shown to be very good approximations of the behavior of PRO model numerical simulations for both the neutral case, and for the long-term behavior of the damped PRO model excited by stochastic forcing. Further, analytical and numerical results from the PRO model explain the relationship between synchronization features across the range of the CGCMs. The results are further supported by the recent identification of a climate mode with energy at the combination tone frequency that results from the nonlinear interaction of ENSO and the annual cycle (Stuecker et al. 2013).

Next, the scenario that ENSO synchronization results from a manifestation of subharmonic frequency locking to the annual cycle was examined within the same recharge oscillator model framework. To do so, the behavior of the oscillator model was examined with the inclusion of periodic forcing by the annual cycle and a cubic damping term representing the nonlinear saturation of SST anomalies. The system was shown to be equivalent to the periodically forced van der Pol oscillator (van der Pol 1927), which has similar global behavior to the ENSO models on which the frequency locking scenario is based, and which allows the same suite of model solutions: quasiperiodic, frequency-locked, and chaotic oscillations. The characteristics of each type of solution were then compared to the observed ENSO. Quasiperiodic oscillations will not produce synchronization with the annual cycle, and so cannot account for the observed ENSO characteristics. Chaotic oscillations were not found within the parameter space of the van der Pol oscillator relevant to the observed ENSO, and occur only rarely in the parameter space of other ENSO models (Jin et al. 1996), so it is unlikely that the observed ENSO characteristics are due to chaotic interaction with the annual cycle. Moreover, chaotic solutions retain the subharmonic peaks associated with frequency-locked solutions (Jin et al. 1996), so the chaotic case can be considered along with the frequency-locked case. Frequency-locked solutions could account for the observed ENSO synchronization, but are far more regular than the observed ENSO, so a stochastically perturbed frequency-locked oscillation remains the only realistic model solution. In such a system, the frequency locking will be accompanied by phase synchronization at the same rational multiple of the annual cycle. In the observations, the only evidence of frequency locking occurs at a ratio of 2:1, but this phase synchronization is not associated with a peak in the spectrum at the biennial frequency. The results indicate that the annual modulation of the coupled stability of the equatorial Pacific ocean–atmosphere system is by far the more likely mechanism generating the synchronization of ENSO events to the annual cycle.

## Acknowledgments

This research was supported by the Office of Science (BER), U.S. Department of Energy, Grants DE-FG02-04ER63862 and DE- FG02-07ER64469. We gratefully acknowledge additional support through the sponsorship of research at the International Pacific Research Center by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC), by NASA through Grant NNX07AG53G, and by NOAA through Grant NA17RJ1230. We thank Prof. Soon-Il An for insightful discussions, and reviewers Prof. Paul Schopf and Prof. Eli Tziperman for their time and effort.

### APPENDIX

#### Expansions of CMIP5 Model Names

ACCESS1.0 Australian Community Climate and Earth-System Simulator, version 1.0

BCC_CSM1.1 Beijing Climate Center, Climate System Model, version 1.1

CanESM2 Second Generation Canadian Earth System Model

CCSM4 Community Climate System Model, version 4

CNRM-CM5 Centre National de Recherches Météorologiques Coupled Global Climate Model, version 5

CSIRO Mk3.6.0 Commonwealth Scientific and Industrial Research Organisation Mark, version 3.6.0

FGOALS-g2 Flexible Global Ocean–Atmosphere–Land System Model gridpoint, version 2.0

GFDL CM3 Geophysical Fluid Dynamics Laboratory (GFDL) Climate Model, version 3

GFDL-ESM2G GFDL Earth System Model with Generalized Ocean Layer Dynamics (GOLD) component

GFDL-ESM2M GFDL Earth System Model with Modular Ocean Model, version 4 (MOM4) component

GISS-E2-H Goddard Institute for Space Studies Model E2, coupled with the Hybrid Coordinate Ocean Model (HYCOM)

GISS-E2-R Goddard Institute for Space Studies Model E2, coupled with the Russell ocean model

HadCM3 Hadley Centre Coupled Model, version 3

HadGEM2-CC Hadley Centre Global Environment Model, version 2–Carbon Cycle

HadGEM2-ES Hadley Centre Global Environment Model, version 2–Earth System

INM-CM4.0 Institute of Numerical Mathematics Coupled Model, version 4.0

IPSL-CM5A-LR L’Institut Pierre-Simon Laplace Coupled Model, version 5A, low resolution

IPSL-CM5A-MR L’Institut Pierre-Simon Laplace Coupled Model, version 5A, mid resolution

IPSL-CM5B-LR L’Institut Pierre-Simon Laplace Coupled Model, version 5B, low resolution

MIROC-ESM Model for Interdisciplinary Research on Climate, Earth System Model

MIROC-ESM-CHEM Model for Interdisciplinary Research on Climate, Earth System Model, Chemistry Coupled

MIROC5 Model for Interdisciplinary Research on Climate, version 5

MRI-CGCM3 Meteorological Research Institute Coupled Atmosphere–Ocean General Circulation Model, version 3

NorESM1-M Norwegian Earth System Model, version 1 (intermediate resolution)

NorESM1-ME Norwegian Earth System Model (intermediate resolution) with carbon cycle

## REFERENCES

*Geometrical Methods in the Theory of Ordinary Differential Equations.*2nd ed. Springer-Verlag, 351 pp.

*Geophys. Res. Lett.,*

**33,**L23708, doi:.

*Phys. Rev. E,*

**77,**066205, doi:.

*Phys. Rev. Lett.,*

**107,**128501, doi:.

*London Edinburgh Dublin Philos. Mag.,*

## Footnotes

^{1}

The *ϕ*(*t*) time series as calculated from the analytical signal are in fact the protophase of the time series, which may differ from the true phase (Kralemann et al. 2008). However, the differences between the protophases and phases calculated in this study are negligible and have no effects on the results, so a discussion of this complexity is omitted.

^{2}

Note that the definition of frequency locking, , where is the annual cycle frequency, is the ENSO frequency, and *k*, *l* ∈ □^{+} is stronger than that for phase locking.

^{3}

There is an indication of another peak in the ENSO spectrum at ~0.35 yr^{−1} and the corresponding combination tones. The spectral estimation methods are not in close agreement in those frequency regions, so these peaks are not considered robust.