Does ENSO Regularity Increase in a Warming Climate?

The impact of a warming climate on El Niño–Southern Oscillation (ENSO) is investigated in large-ensemble simulations of the Community Earth System Model (CESM1). These simulations are forced by historical emissions for the past and the RCP8.5-scenario emissions for future projections. The simulated variance of the Niño-3.4 ENSO index increases from 1.48C in 1921–80 to 1.98C in 1981–2040 and 2.28C in 2041–2100. The autocorrelation time scale of the index also increases, consistent with a narrowing of its spectral peak in the 3–7-yr ENSO band, raising the possibility of greater seasonal to interannual predictability in the future. Low-order linear inverse models (LIMs) fitted separately to the three 60-yr periods capture the CESM1 increase in ENSO variance and regularity. Remarkably, most of the increase can be attributed to the increase in the 23-month damping time scale of a single damped oscillatory ENSO eigenmode of these LIMs by 5 months in 1981–2040 and 6 months in 2041–2100. These apparently robust projected increases may, however, be compromised by CESM1 biases in ENSO amplitude and damping time scale. An LIM fitted to the 1921–80 observations has an ENSO eigenmode with a much shorter 8-month damping time scale, similar to that of several other eigenmodes. When the mode’s damping time scale is increased by 5 and 6 months in this observational LIM, a much smaller increase of ENSO variance is obtained than in the CESM1 projections. This may be because ENSO is not as dominated by a single ENSO eigenmode in reality as it is in the CESM1.


Introduction
As the dominant mode of tropical interannual variability with global teleconnections, El Niño-Southern Oscillation (ENSO) not only is the leading source of forecast skill on seasonal and interannual time scales but also plays an important role in the global dynamics of climate change. ENSO is associated with an irregular oscillation of sea surface temperatures (SSTs) in the tropical Indo-Pacific Ocean with periods in the 3-7-yr range. Given its prominent role in climate variability, it is of great interest to determine how ENSO may change in a warming climate. Capotondi and Sardeshmukh (2017) and Aiken et al. (2013) studied the change of ENSO in the observational record and found a general increase of ENSO variability, as well as of SST spectral power in the 3-7-yr band. Their studies were limited by the length of their observational record of about 50 years. They partly addressed this limitation, when assessing the statistical significance of their conclusions, by generating large synthetic dynamically consistent samples using linear inverse models (LIMs) fitted to different segments of the observational record.
To address the problem of limited sampling, the climate response is usually studied in ensembles of climate model simulations with prescribed radiative forcing. With regard to ENSO, while models are powerful tools Denotes content that is immediately available upon publication as open access.
for representing the complex interactions among trade winds, atmospheric convection, and ocean dynamics-all contributing to and modulating ENSO dynamics-the delicate balance between damping and amplifying feedback processes differs across models. There has consequently been no clear consensus on how ENSO amplitude and frequency will change under global warming (Cai et al. 2015;Collins et al. 2010). Out of the 17 models of phase 3 of the Coupled Model Intercomparison Project (CMIP3), for instance, about half showed an increase and half a decrease in ENSO variability (Collins et al. 2010). Kim et al. (2014) argued that in models with a good representation of ENSO feedbacks, ENSO variability increases in the period before 2040, in which SSTs warm faster in the eastern Pacific Ocean than over the Maritime Continent region, but decreases thereafter. An increase in variability was also claimed by one of the first studies of this topic (Timmermann et al. 1999).
The emphasis in previous studies has generally been more on changes in ENSO amplitude than on ENSO period (its peak in the 3-7-yr band) and ENSO regularity (the sharpness of that peak). The causes of changes in specific ENSO properties have usually been sought in terms of changes in the slowly evolving background state, such as the weakening of the mean Walker circulation or faster warming in the eastern equatorial Pacific and the Maritime Continent than in the central Pacific (Cai et al. 2015, and references therein).
Changes in ENSO variability need not always be associated with changes in the background state. Christensen et al. (2017), for example, reported a marked improvement in simulated ENSO variability by adding a stochastic parameterization of unrepresented subgrid variability, including that of atmospheric diabatic heating, to the coupled NCAR-CCSM4 model. The stochastic parameterization affected the model variability, but not the background state. This was explained by Berner et al. (2018) using a simple conceptual damped linear oscillator model to show how perturbations to its parameters could alter ENSO variability without impacting the mean. Such a damped linear oscillator arguably provides the simplest description of ENSO, and may be identified with the damped oscillatory eigenmodes of a linear inverse model fitted to observational data or climate model output.
The influential work of Penland and Sardeshmukh (1995) showed that the evolution of observed tropical SSTs is well captured by a low-order LIM (see also Penland 1989;Penland and Magorian 1993). In particular, LIMs can capture the temporary growth of ENSO through modal interference between nonorthogonal eigenmodes (e.g., Alexander et al. 2008;Zanna 2012;Gehne et al. 2014;Capotondi and Sardeshmukh 2017;Newman and Sardeshmukh 2017). Newman et al. (2009) showed that on seasonal to interannual scales the power spectra of synthetic ENSO time series obtained from LIM integrations generally compare better with the observed spectrum than do spectra from the CMIP3 coupled models.
In this paper, we analyze the large-ensemble simulations (LENS; Kay et al. 2015) of the 1921-2100 period generated using NCAR's Community Earth System Model, version 1 (CESM1). These coupled simulations use historical forcing for the past and the RCP8.5-scenario forcing for future projections. We analyze the modeled long-term changes in ENSO properties in the framework of the simple damped linear oscillator considered by Berner et al. (2018). To this end we fit low-order LIMs to the past and future simulated periods [period 1 (P1): 1921-80; period 2 (P2): 1980-2040; and period 3 (P3): 2040-2100] and quantify the changes not only in the variance but also in the temporal autocorrelation. The latter is directly linked to the width of the spectral peak, or equivalently to the regularity of the oscillation. This has implications for altering the predictability of interannual variations in a warmer climate.
Specifically, our goal here is to address the following questions: We recognize that a single-model study cannot resolve the ambiguity of the ENSO response to greenhouse warming reported in the literature. Nonetheless, we believe that using a large ensemble of simulations with a single model can establish the response in that model with greater statistical confidence, and using low-order LIMs can help us understand it better. The manuscript is organized as follows. Section 2 summarizes the datasets and their preprocessing. Section 3 provides a brief overview of linear inverse modeling and the simple damped linear oscillator used to interpret the changes of ENSO in the CESM1 simulations (section 4). A discussion and conclusions follow in section 5.

Data
The CESM1 LENS utilized here comprise an ensemble of 33 separate simulations of the period 1920-2100 differing only in initial conditions. The simulations are forced by historical emissions of carbon dioxide, ozone, and aerosols in 1920-2005 and by the RCP8.5-scenario emissions from 2005 onward. Figure 1 shows the simulated annually averaged global mean temperature. The global mean temperature varies around 3.58C with no pronounced trend until the 1970s or so, and rises markedly thereafter. We additionally utilize a 900-yr control simulation (CNTL) generated using the same model, but forced with constant preindustrial emissions.
ENSO variability is often characterized by the time series of a scalar index. The Niño-3.4 index is defined as the equatorial SST anomaly averaged over the region 58-58N, 1708E-1208W (e.g., Trenberth 1997). We remove the climate trend and seasonal cycle by removing the ensemble mean and applying a 3-month runningmean filter before calculating the index. Alternatively, an ENSO-index can be defined as the first principal component (PC) of detrended SST anomalies in the tropical belt, 208S and 208N (e.g., Penland and Sardeshmukh 1995). To ensure that the empirical orthogonal function (EOF) basis has no climate change signal, we project the LENS simulations onto the EOFs of the control simulation. The time series of these two indices are shown in Fig. 2 for five arbitrarily picked ensemble members. Both indices pick up the same ENSO variations; indeed they are correlated at 0.95. We use the PC-based index from here on.
We characterize the observed ENSO variations using monthly SSTs from the HadISST2 dataset, which spans the period 1900-2010. These SSTs were preprocessed analogously, except that a linear trend was removed and the EOFs obtained from observed SSTs. The observational ENSO indices in the period 1920-2010 are also shown in Fig. 2.

Methodology
a. Linear inverse model Penland and Sardeshmukh (1995) showed that observed tropical SST anomaly dynamics may be approximated by a linear stochastically forced system of the form where x(t) is the time-evolving state variable, L the linear feedback matrix, and S the amplitude of the additive (i.e., state-independent) white-noise forcing e. When L and S are obtained empirically, (1) is called a ''linear inverse model'' (LIM). Note that while (1) is formally linear, nonlinear processes are not discarded, but approximated as linear terms plus noise. The statistics of any such system are fully represented in the time-lag covariance matrix C t (Penland 1989): which can be used to estimate of feedback matrix as where t is the time lag, log is the matrix logarithm, and the hat (Á) denotes that the operator is estimated from observational data or climate mode output. In practice, C t is estimated for a particular time lag t 5 t 0 , so that the validity of the underlying model assumptions can then be tested for other lags via the so-called t test (Penland and Sardeshmukh 1995). An estimate of the noise covariance Q 5 SS T can be obtained from the Lyapunov equation (Gardiner 1983;Farrell and Ioannou 1993;Penland and Matrosova 1994) by inserting the estimatedL andĈ 0 and solving forQ. Because C 0 and Q 5 SS T are by definition positive-definite, (4) can be satisfied only if L is a stable operator (i.e., if its eigenvalues are either real and negative or come in complex conjugate pairs with negative real parts). Each oscillatory eigenmode is complex and describes an evolving structure that can be represented as a combination of two patterns varying in quadrature. In other words, each eigenmode has the form of a damped linear oscillator (see below). Of special interest is the least damped oscillatory mode with period and spatial patterns close to those associated with ENSO (Penland and Sardeshmukh 1995), which we henceforth refer to as the ''ENSO eigenmode.'' While fitting a linear inverse model enables a modeby-mode comparison of different datasets, it does not follow the system can be decomposed into a set of independent modes. Indeed, the power of the LIM lies in its ability to capture the transient energy growth through interference between nonorthogonal eigenmodes. The ratio between energy at time t and initial time is given by (Penland and Sardeshmukh 1995)  possible growth through modal interference and the associated left and right eigenvectors denote the optimal initial and evolved structures.

b. The simple model of a linear damped harmonic oscillator
We are interested here in relating the statistical properties of a damped linear oscillator to the parameters of the oscillator. The evolution of such an oscillator forced by noise is governed by (1): The state vector x comprises the amplitude and rate of displacement of an oscillator with frequency v and damping rate n. In order for L to be stable, the damping rate has to be positive, n . 0. We assume for simplicity that the noise forcing of the two oscillator components is uncorrelated and has the same amplitude, S 5 [s « , 0j0, s « ].
As discussed by Berner et al. (2018; see also references therein) the zero-lag covariance of the two-component state vector x in such a system may be obtained using (4) as Hence an increase in the noise forcing amplitude increases the variance of x. Interestingly, the covariance C 0 does not depend upon the oscillator frequency. This means that a change in frequency does not change the variance. On the other hand, the temporal autocorrelation is independent of the noise amplitude s « . Thus a change in the strength of the forcing changes the autocovariance but not the autocorrelation of the damped oscillator. Of particular interest is the decorrelation (or e-folding) time t d at which the autocorrelation decays to 1/e. For our simple model the decorrelation time scale is simply which is a function of the damping rate but not on the frequency. According to the Wiener-Khinchin theorem (e.g., Gardiner 1983), the power spectrum [or power spectral density (PSD)] normalized by the variance is the Fourier transform of the autocorrelation function. This theorem holds for a wide range of so-called wide-sense-stationary random processes, including the physical system investigated here. Consequently, a change in the temporal autocorrelation translates directly to a change in the (normalized) spectrum, and vice versa.
The spectral matrix of a system given by (1) is given as (Gardiner 1983) where f is the frequency variable. For our simple system, (1) and (5), the spectrum of each component simplifies to From (10) we see that the spectrum has a maximum at f 5 v and the spectral width is a function of the damping rate n. The larger the damping rate, the broader and lesspeaked is the spectrum. These simple considerations motivate us to use not only the variance but also the autocorrelation (or alternatively power spectral density) to characterize and interpret the changes of ENSO variability in response to global warming.

a. Variance and autocorrelation of the ENSO index
Even a casual inspection of the ENSO time series in Fig. 2 suggests that it is nonstationary. Both frequency and amplitude appear to increase with time. To quantify these changes we divide the simulations into three time periods of 60 years each. P1 covers the years 1921-80, P2 the years 1981-2040 (gray shading), and P3 the years 2041-2100. This choice was made so that in the first period the global mean temperature has no pronounced trend (Fig. 1), and also so that all three periods have the same number of years.
Guided by the simple oscillator model, the variance and autocorrelation of the ENSO index were computed. The temporal variance averaged across all ensemble members increases from 1.48C 2 in the first 60 years to 1.9 and 2.28C 2 in the subsequent 60-yr periods (Fig. 3, Table 1). According to an F test, these differences are significant at the 95% confidence level, even if the samples are made independent and identically distributed (iid) by only using every 48th sample to account for temporal correlations. While the shift in the variance is pronounced, the width of the probability distribution of the variance obtained in different ensemble members ( Fig. 3b) points to large uncertainty. Note that the observed variance of 0.78C 2 in 1921-80 is smaller than that of any LENS member in that period. We will return to this issue in section 4d, which contains a detailed comparison of the model results with observations over the historic period.
The autocorrelation functions for P1-P3 (Fig. 3c) are indicative of oscillations with periods of 45-49 months. The autocorrelation decays less rapidly in the later periods, P2 and P3, indicating less damping and more ENSO memory in P2 and P3 than in P1 (Fig. 3c). When fitting an exponential envelope to the amplitude of the autocorrelation function (shading), we can determine the decorrelation time as time where this envelope intersects the 1/e 5 const curve (dashed line). From the graph, we determine the decorrelation times to 18 (P1), 27 (P2), and 34 (P3) months (Table 1). The LENS simulations thus indicate that in a warming climate ENSO events will have both greater amplitude and greater memory, with implications for potentially greater predictability.

b. Spectra
To facilitate comparison with other studies, we also present spectra in addition to autocorrelation functions.  To this end we computed the power spectral density (PSD) not as the Fourier-transformed autocorrelation function but directly from the squared Fourier coefficient of the time series obtained using Welch's periodogram method (which is included the Matlab software package). The spectra are characterized by a distinct peak in the 3-7-yr ENSO band (Fig. 4). Note that our preprocessing of the data filters out variability on the longest time scales (through detrending) as well as the annual cycle. The mean spectrum averaged over all ensemble members shows sharper peaks in P2 and P3, with spectral densities of 110 and 1508C 2 cpm compared to 658C 2 cpm in P1 (where cpm is cycles per month).
The uncertainty in the power spectral density is indicated by the standard deviation across all ensemble members for each frequency (dark shading) and by the maximum and minimum density range (light shading). The uncertainty is large and roughly proportional to the power (i.e., largest in the ENSO band). Intercomparisons across P1-P3 show that the mean power spectrum in each period is within the sampling uncertainties in the other periods. Note that a sample spectrum tracking the lower bound at all frequencies has a much smaller total variance (given by the weighted area under the PSD curve) than the total variance associated with the mean spectrum.
To focus on the uncertainty of ENSO, we plot histograms of PSD across the entire frequency band between 3 and 7 years. The PSD distributions are all positively skewed toward high spectral densities (Fig. 5a). The means for P2 and P3 are higher than for P1 and the right tail is much heavier. We note that the change from P1 to P2 is more pronounced than from P2 to P3, pointing to a possible saturation effect. The LENS results suggest that the ENSO amplitudes in the future will be on average larger than in the past.

c. Realizations from linear inverse models
Previous studies have argued that a linear inverse model is able to capture observed interannual variability in the tropics and can be used to make predictions on this time scale (e.g., Penland and Sardeshmukh 1995;Newman et al. 2009;Capotondi and Sardeshmukh 2017 The LIM feedback matrix L is derived from the covariance matrix fitted to a particular lag t 0 . If the linear approximation is valid, we expect there to be a range of lags for which the LIM should capture the system's behavior. This range is bounded on the lower end by the time scale for which nonlinearities become important and on the upper end by the lag for which sampling errors dominate. We fitted LIMs for a range of time lags between 1 and 12 months and confirmed that the results were qualitatively similar. The following results are obtained for a time lag of t 0 5 3 months. LIMs for each period were fitted and used to generate 5000-member ensembles of 60-yr integrations. The agreement with the mean LENS spectra is overall excellent (Figs. 4d-f), confirming that the LIMs fitted to different periods are indeed able to capture the changes in tropical variability. On decadal time scales, the LIM appears to overestimate the uncertainty, although this could also be an artifact of the limited sample size of the LENS. In the ENSO band, the agreement between the LIM and LENS is excellent (Figs. 5b-d) and captures not only the shift in the mean and variance, but also the change in the skewness of the variance distribution. We conclude that the LIM is able to capture ENSO variability in a warmer climate and can be used to increase the signal-to-noise ratio in the power spectral density in the ENSO band (Fig. 5e). Next, we will investigate the eigenmodes of the linear feedback matrix L fitted separately to P1-P3 to gain insights into the mechanisms by which the ENSO dynamics change. The real part of the least damped oscillatory eigenmode (Fig. 6) is highly reminiscent of the ENSO peak pattern (Penland and Sardeshmukh 1995;Gehne et al. 2014). The spatial characteristics of this eigenmode (both real and imaginary parts) in all three 60-yr periods closely resemble that obtained from the observations (Figs. 6g,h).
The real and imaginary parts of the complex conjugate pair of eigenvalues l 1,2 5 2n 6 iv indicate the e-folding time t d 5 21/n and frequency v of each eigenmode. Scatterplots of e-folding times versus frequency for all eigenvalues (Fig. 7a) show that besides the ENSO mode, there are three other modes in the ENSO band, but with much shorter (;5-month) e-folding times. The ENSO mode has periods of 48, 50, and 51 months in P1, P2, and P3 respectively. The associated e-folding time increases from 23 months in P1 to 28 and 29 months in P2 and P3. For the linear oscillator, this increase in the e-folding time signifies an increase in variance, (6), as well as a narrowing of the spectrum, (10). This is in qualitative agreement with the behavior in the full climate model.
The observed ENSO mode has a period of 46 months and e-folding time of 8 months (Fig. 7b). While the modeled period is similar to the observed one, the decorrelation time in P1 is with 23 months almost 4 times as large as in the observations. This will be discussed in further detail in the next section.
To see if the changes in the eigenvalue of the ENSO mode are sufficient to explain the changes in ENSO dynamics quantitatively, we use the LIM fitted to P1 (1920-81) but replace the eigenvalue for the ENSO mode with that from the LIM fitted to P2 and P3, respectively. The modified feedback matrix are recomputed as L mod 5 VL mod V 21 , where V is the matrix containing the right eigenvectors of L and L mod is a modified diagonal matrix of eigenvalues with modification only of the ENSO complex conjugate eigenvalue pair.
As expected, the spectrum of the modified LIMs (dash-dotted lines) have more power in the ENSO band, although the peaks are slightly lower than those for LIMs directly fitted to P2 and P3 (Figs. 8a,b). These results suggest that an increase from 23 to 28 months (from P1 to P2) and from 23 to 29 months (from P1 to P3) in the damping time scale of the single ENSO mode can by itself explain almost all of the changes in the CESM1 spectra.
The temporal energy growth through modalinterference of the nonorthogonal eigenmodes is depicted by the ''maximal amplification (MA) curve,'' which is given by largest eigenvalue of G T G as function of the lag t (see section 3). The curves indicate maximal growth for a lag of 14-15 months in the climate simulations, but only 7 months for the observations (Fig. 9). Figure 9 also depicts the average error curves estimated as the trace of the error covariance matrix he(t)e T (t)i 5 C(0) 2 G(t)C(0) G T (t), normalized by the trace of the covariance (Penland and Sardeshmukh 1995).
The time lag t c at which these MA curve falls below the error energy curve gives an upper bound for the predictability limit in the absence of noise. We see that for observations this critical lag is 16 months, consistent with the literature. For LENS in the period 1921-80, this lag is with 32 months, roughly twice as long, which amounts to an overestimation of ENSO predictability in the climate model. For the future periods 1981-2040 and 2041-2100, t c increases to 43 and 46 months, which amounts to an increase of the predictability horizon by about 30%. This finding is consistent with our earlier result showing an increase in the decorrelation time for the future periods (Fig. 3c).

d. Comparison with observed ENSO-variability
This study focuses on the ENSO response to global warming as projected by the CESM1 under the RCP8.5 emissions scenario. The model projects a clear increase in ENSO variance and regularity, but this may be compromised by its biases in ENSO amplitude and damping time scale. To assess these biases, we computed an observed ENSO index using SSTs from the HadISST2 observations analogously to the model (Fig. 2). The observed variance of 0.78C 2 is smaller than that of any LENS ensemble member and the autocorrelation function also decays faster to zero (Table 1, Figs. 3b,c).
When compared with the LENS for the historical period, the observed spectrum has less power and a broader spectral peak (Fig. 4a). As already discussed, we cannot conclude that the observed spectrum is inconsistent with the LENS spectrum at any specific frequency, but the fact that the observed spectrum tracks the lower bound of the sample spectra at all frequencies further demonstrates that the total observed variance is outside the range of the sample LENS variances (Figs. 3b,e).
An LIM fitted to the HadISST2 observations in period P1 has an ENSO eigenmode very similar to that of the model (Figs. 6g,h), but with a much shorter 8-month e-folding time scale than the 23-month time scale of the model's ENSO mode in P1 (Fig. 7b). This is consistent with the smaller power and a broader peak in the observed spectrum (Fig. 10a). The differences between the LIMs fitted to observations and LENS for the period 1920-81 is especially evident in the ENSO band, where the LENS LIM is biased toward higher power spectral densities (Fig. 5f).
While there is compelling evidence for inconsistencies between the observations and the LENS for the historical period, we cannot overcome the inherent problem of the limited sample size in observations and our results have to be interpreted carefully. To illustrate this, we subjectively picked the LENS ensemble member whose spectrum in P1 is closest to the observed spectrum. An LIM fitted to this single member still has too much power in the ENSO band, although the discrepancy is much smaller than when using all members (Figs. 10b,c).
The sensitivity of the CESM1 projected changes in ENSO to its biases in ENSO dynamics and variability may be gauged using our observational LIM. As already noted, this LIM has an ENSO eigenmode with a much shorter 8-month damping time scale than the 23-month time scale of CESM1's ENSO eigenmode. When this 8month damping time scale is increased by 5 and 6 months, reflecting CESM1's projected 5-month increase from 23 to 28 months in P2 and 6-month increase from 23 to 29 months in P3, respectively, this observational LIM also ''projects'' an increase in ENSO variance and a narrower spectrum. However, these changes are much weaker than obtained in the CESM1 projections. This may be because ENSO is not as dominated by a single ENSO eigenmode in reality as it is in the CESM1, consistent with the fact that its observed decay time scale of 8 months is much closer to that of several other observational eigenmodes, whereas in the CESM1 the 23-month time scale of the eigenmode is much longer than that of other eigenmodes.

Conclusions
Given its dominant role in tropical interannual variability and predictability, it is of great interest to determine how ENSO responds to a warming climate. We addressed this question here by analyzing a large ensemble of coupled simulations of the 1921-2100 period generated using NCAR's Community Earth System Model (CESM1). To quantify the changes in ENSO over this long 180-yr simulation period, we focused on the changes in the ENSO index for three subperiods of 60 years each: P1 (1921-80), P2 (1981-2040), and P3 (2041-2100). According to this model, tropical variability will increase as global temperatures rise. Specifically, the model projects an increase of ENSO amplitude from 0.68C for P1 to 1.48C and 1.58C in P2 and P3 under the RCP8.5 emissions scenario. An increase of ENSO amplitude with warming is consistent with some previous studies (Timmermann et al. 1999;Kim et al. 2014;Cai et al. 2015;Capotondi and Sardeshmukh 2017), but there remains ambiguity (Collins et al. 2010). The temporal autocorrelation time scale of ENSO also increases with warming in this model, indicating increased system memory. This is associated with a narrowing of the spectral peak in the ENSO band, which points to an increase in ENSO regularity. The CESM1 simulations thus suggest potentially greater predictability of ENSO events in a warmer climate. Considering its widespread impact on characteristic patterns of rainfall and temperature with implications for flooding and drought especially in vulnerable regions, increased ENSO predictability is associated with considerable social and economic value.
Previous work has shown that linear inverse models (LIMs) are good at capturing observed interannual variability and predictability in the tropical belt (Penland and Sardeshmukh 1995;Newman et al. 2009). Here, we verify that LIMs are also able to capture changes in tropical variability under global warming. LIMs fitted separately to the 60-yr periods of P1-P3 capture the LENS spectra in those periods very well, especially in the ENSO band.
The changes in the CESM1 spectra and memory can be understood in terms of the changes in the eigenvalue of the least damped oscillatory eigenmode (the ENSO eigenmode) of the fitted LIMs. The increase in the decay time of this ENSO eigenmode from 23 months in P1 to 28 and 29 months in P2 and P3 accounts for almost all of the changes in the CESM1 spectra from P1 to P3. Since the spectra are closely related to the autocorrelation function, this also explains the increase in the memory of the autocorrelation function. This is consistent with the statistics of a stochastically forced simple damped linear oscillator, for which an increase the e-folding decay time can be directly linked to an increase in variance and regularity.
Our study was performed with one particular climate model, CESM1. It would be interesting to determine to what degree our results carry over to other models, such as those contributing to CMIP6. One can imagine that different models will project different ENSO responses to warming due to their different model biases. With regard to CESM1, we showed here that while CESM1 captures the main characteristics of ENSO variability, it overestimates the magnitude as well as the correlation time (i.e., its ENSO spectral peak is too high and too narrow). A separate LIM fitted to the 1921-80 observations has an ENSO eigenmode with a much shorter 8-month damping time scale, comparable to that of several other eigenmodes. When the mode's damping time scale in this observational LIM is increased by 5 and 6 months, corresponding to its increase in P2 and P3 in the CESM1 projections, a much smaller increase of ENSO variance is obtained than in the CESM1 projections. This may be because a single ENSO eigenmode is not as dominant in reality as it is in the CESM1. This is a subtle but important model error, whose existence in other climate models may similarly compromise changes of ENSO projected by those models. Reducing such model biases in ENSO variability, perhaps by implementing stochastic parameterizations of atmospheric diabatic processes Berner et al. 2017), will be necessary to build confidence in model projections of changes in ENSO with global warming.