A Review of Paleo El Niño-Southern Oscillation

: The Earth has seen El Niño-Southern Oscillation (ENSO)—the leading mode of interannual climate variability—for at least millennia and likely over millions of years. This paper reviews previous studies from perspectives of both paleoclimate proxy data (from traditional sediment records to the latest high-resolution oxygen isotope records) and model simulations (including earlier intermediate models to the latest isotope-enabled coupled models). It summarizes current understanding of ENSO’s past evolution during both interglacial and glacial periods and its response to external climatic forcings such as volcanic, orbital, ice-sheet and greenhouse gas forcings. Due to the intrinsic irregularity of ENSO and its complicated relationship with other climate phenomena, reconstructions and model simulations of ENSO variability are subject to inherent difﬁculties in interpretations and biases. Resolving these challenges through new data syntheses, new statistical methods, more complex climate model simulations as well as direct model-data comparisons can potentially better constrain uncertainty regarding ENSO’s response to future global warming.


Introduction
El Niño-Southern Oscillation (ENSO) is the leading mode of natural climate variability that fluctuates irregularly on interannual time scales. It is known that ENSO is associated with tropical ocean-atmosphere interactions [1][2][3][4][5]: an initial warming in the eastern equatorial Pacific sea surface temperature (SST) tends to reduce the surface pressure, trade wind strength, and, in turn, the oceanic equatorial upwelling, further amplifying the initial warming and eventually leading to an El Niño event. The oceanic waves (westward-propagating Rossby waves and eastward-propagating Kelvin waves) accompanying the initial warming tend to generate a delayed cooling in the subsurface ocean, which eventually switches the surface warming to a cooling, or a La Niña event. This episodic warming and cooling of SST can cause unusual weather patterns such as droughts and floods across the globe [6] through atmospheric teleconnections [7], leaving a widespread impact on human society [4].
Studies on ENSO dynamics and their impact have progressed dramatically in the last 40 years, in which ENSO is either characterized as a self-sustained oscillator (various oscillator models are thus proposed) or a stable mode punctuated by stochastic forcing (see a review by Wang et al. [8]). Together with advanced understanding of ENSO dynamics, operational ENSO predictions using both dynamical and statistical models have shown substantial skills at seasonal time scale, as exemplified by the prediction of the last 2015-16 El Niño event [9]. However, the model prediction skill of ENSO is

Reconstructing Paleo ENSO Using Proxies
The paleoclimate proxy reconstructions that are closely related to ENSO variability are either direct proxies in the core ENSO region (most of them are coral and mollusk records and some are marine sediments) or indirect ones from remote regions that are affected significantly by ENSO through atmospheric and oceanic teleconnections (lake sediment, tree ring, speleothem and ice core) ( Figure 1; Table 1). Circles/triangles denote proxies that cover a period older/newer than 1 ka. Larger green triangles are for tree-ring synthesis clusters. Black boxes show Nino 3.4 region and Nino 1 + 2 region.  [23] Lake sediment 10 ka Bainbridge Crater Lake, Galápagos Islands Conroy et al., 2008 [24] Lake sediment 10 ka El Junco Crater Lake, Galápagos Islands Z. Zhang et al., 2014 [25] Lake sediment 10 ka El Junco Crater Lake, Galápagos Islands Thompson et al., 2017 [26] Lake sediment 10 ka Bainbridge Crater Lake, Galápagos Islands Rein et al., 2005 [27] Marine sediment 21 ka 78 • W, 12 • S Leduc et al., 2009 [28] Marine sediment 50 ka 84 • W, 8 • N Koutavas and Joanides, 2012 [29] Marine sediment 21 ka 90 • W, 1 • S Sadekov et al., 2013 [30] Marine sediment 21  A direct assessment of past SST variability in the tropical Pacific spanning several decades can be retrieved from oxygen isotope measurements on fossil corals, collected from the tropical central and eastern Pacific [33,35,[39][40][41], or from mollusk shells well preserved in the archeological sites on the coast of Peru [42]. The carbonate oxygen isotopes (δ 18 O) in such marine records reflect variations in both SST and the hydrological cycle at as high as monthly resolution [43]. Most of these high-resolution records cover the past several thousand years. Some marine sediment records also have been measured and also been linked to the past ENSO variability, going back to tens of thousands of years. They include continuous time series of carotenoid pigments and fine-grained lithic [27] and scattered records of oxygen isotope (δ 18 O) and magnesium/calcium (Mg/Ca) ratios in sedimentary foraminifera shells [28][29][30][31]44].
South America coastal regions and nearby islands in the tropical eastern Pacific, located in the ENSO center of action, are ideal regions to search for information on past ENSO variability. Early archaeological records [45] and later high-resolution lake sediment records showing variations in sediment composition [22][23][24]26,46] have been collected. These radiocarbon-dated records are usually continuous in time and can date back thousands of years in the past. To reconstruct ENSO history over the last few centuries, tree ring data has been used (see a synthesis by Li et al. [36]). Strong ENSO signals have been detected in tree rings from the tropics and the mid-latitudes. Recently, tree-ring cellulose δ 18 O, instead of the conventional width, was also measured to infer the SST change over the tropical central Pacific [37,47]. On rare occasions, oxygen isotope ratios from ice cores at low latitudes can provide a high-resolution record of climate variability in the tropics during the last millennium [38]. All these circum-Pacific proxy reconstructions are interpreted as anomalous regional weather or climate patterns, which are linked to the SST variability in the tropical Pacific.

Simulating Paleo ENSO Using Climate Models
Since the earlier studies, paleo ENSO simulations have been carried out in both intermediate models that require prescribed mean climate state [48,49] and coupled general circulation models (CGCM) [50]. Current CGCMs can reproduce ENSO's amplitude, spatial pattern, asymmetry, phase locking and irregularity that are generally within the range of its natural variation but are heavily model dependent [51][52][53]. Recently, water isotope-enabled CGCMs allowing direct model-data comparison of oxygen isotopes have been applied to paleo-ENSO [54][55][56]. Paleo ENSO modeling has benefitted significantly from the Paleoclimate Modelling Intercomparison Project (PMIP) [19] (https:// pmip4.lsce.ipsl.fr/doku.php/exp_design:index) for two periods: the mid-Holocene (MH,~6000 years ago) [57] and the Last Glacial Maximum (LGM,~21,000 years ago) [58]. The MH simulations provide an opportunity to explore ENSO response to orbital forcing, while the LGM simulations focus on ENSO sensitivity to the presence of vast continental ice sheets, reduced greenhouse gas (GHG) concentrations and a generally cold mean climate. Additionally, the Pliocene Model Intercomparison Project (PlioMIP) [59,60] has been created to investigate the mid-Pliocene warm period (~3.3-3 Ma). The PlioMIP ensemble mainly aims to advance our knowledge about the ENSO response to increased GHG concentrations (400 ppm of CO 2 ) and a different land-sea distribution. These CGCM experiments are all snapshot experiments that simulate the climate at a particular period of a few hundred years. Recently, transient CGCM simulations have been used to study the evolution and abrupt change of the ENSO response during the entire last deglaciation, such as in the TRACE21 project that simulates the transient climate evolution since the LGM in CCSM3 [61,62] and the transient simulation of the last 142,000 years in ECHO-G model forced by the orbital forcing accelerated by a factor of 100 [63]. Ensemble transient simulations are also widely performed to study the ENSO variability and forcing mechanisms in the last millennium [64][65][66][67][68], in which even more external forcings are included (solar activities, volcanic eruptions, land use, GHG, orbital parameters and ozone-aerosol).
Next, we will discuss the specific ENSO reconstructions and possible reasons for inconsistencies among the data and review the present status of paleoclimatic numerical experiments and model-data comparison of ENSO. Then, we examine recent developments in understanding the forcing mechanisms of paleo ENSO in terms of the mean climate and, in turn, coupled instability, the nonlinear interaction of ENSO with the annual cycle and the impact of stochastic forcing. We proceed in time order of different periods (and the related external forcing factors), with particular reference to the last millennium, the Holocene, the glacial-deglacial and warm periods beyond the last glacial.

Proxy Records
The existence of many long, high-resolution paleoclimate records from the circum-Pacific regions paints a relatively detailed picture of ENSO variability in the last millennium. In particular, paleo ENSO studies for the past millennium have been focusing on resolving the dynamical connections between ENSO processes and centennial Northern Hemisphere (NH) climate anomalies [69]-notably the warm Medieval Climate Anomaly (MCA,~900-1450 CE), the cold Little Ice Age (LIA,~1450-1850 CE) and then an industrial global warming. Volcanic eruptions inject aerosols into the atmosphere and impact the global climate through radiative forcing and complicated interactions with clouds [70,71]. The last millennium is an ideal period to investigate the ENSO sensitivity to radiative forcing induced by large tropical volcanic eruptions that are well documented [72] (Figure 2a, black bar).
ENSO variability is the dominant source of δ 18 O variability during the last millennium in the Palmyra fossil-coral δ 18 O records [39] (Figure 2b). The mean climate in the central tropical Pacific shifted from relatively cool and dry during the tenth century to increasingly warmer and wetter climate in the 20th century. ENSO exhibits considerable natural variability but its behavior is poorly correlated with the estimates of mean climate. The most intense ENSO activity occurred during the mid-seventeenth century. It is also found that ENSO exhibited considerable natural variability at decadal time scale. These results imply that the majority of ENSO variability over the last millennium might have arisen from internal ENSO dynamics.
Carotenoid pigments and fine-grained lithic rate in the marine sediment off Peruvian coast [27] (Figure 2c; see also Figure 12c in [27]) shows that at the beginning of the last millennium El Niño was not active and was persistently suppressed since its maximum at~3-2 thousand years before present (ka). Later, sediment flux regained pre-anomaly values for a short period (~1300-1500 CE) after the MCA. The last five hundred years of the sediment is rather uniform, implying relatively weak El Niño. However, their conclusion is contradictory to that from another multidecadally-resolved marine sediments on Galápagos in the eastern equatorial Pacific [32]. In Rustic et al. [32] sample variance of δ 18 O in individual Globigerinoides ruber (G. ruber) shells was calculated and compared with modern sample benchmark (Figure 2d), as well as with the occurrence of warm and cold outliers beyond high-probability ENSO events over the continuous period~1000-2009 CE. Both the frequency of outliers and δ 18 O variance indicate a mid-millennium (~1500 CE) shift from a less volatile state to higher ENSO variability. Also, the averaged modern δ 18 O variance is found to be intermediate between the millennium maxima and minima but the exceptionally intense 1997-1998 El Niño has no match in the earlier part of the millennium.  [39], modified from Otto-Bliesner et al. [66]; (c) Lithic concentration in a marine sediment from off Peru (% of maximum) [27]; (d) δ 18 O variability from marine sediments from the tropical eastern Pacific [32]; (e) Tree-ring synthesis from circum-Pacific sites (modified from Li et al. [36]); (f) δ 18 O variability in tree-ring from Taiwan (modified from Liu et al. [37]). Vertical bars denote the Medieval Climate Anomaly (MCA), the Little Ice Age (LIA) and the industrial warming, respectively.
Li et al. [36] synthesized 2222 tree-ring records from both the tropics and mid-latitudes (Asia, New Zealand and North and South America) to reconstruct ENSO variability for the past seven centuries (Figure 2e). This comprehensive tree-ring network includes several individual datasets in previous studies on this topic [73,74]. Strong ENSO signals and substantial interdecadal-centennial variability is resolved in these records. The reconstructed ENSO variance is coherent with the NH mean climate. ENSO weakened in the early LIA of the 1300s-1550s, then strengthened during the late LIA of the 1550s-1880s. Moreover, ENSO activity was anomalously high in the late 20th century, exceeding that exhibited over the past seven centuries. The tree-ring compilation also suggests a robust ENSO response to large tropical eruptions, with anomalous cooling in the east-central tropical Pacific in the year of eruption, followed by anomalous warming one year after. The new δ 18 O reconstruction from tree cellulose in Taiwan [37] (Figure 2f) also reveals a late 20th century enhancement of interannual variability compared to any other intervals of the 818-year-long reconstruction, associated with relatively warmer Nino 4 sea-surface mean temperature. Before the 20th century, the largest interannual excursions occurred during the early to mid-seventeenth century. The highest single tree-ring δ 18 O value of the entire reconstruction was in 1651 AD and might correspond to an exceptionally large El Niño event.
A compilation of the last millennium records shows some degree of similarity in the highly volatile ENSO variance on decadal to centennial time scales. Although ENSO variability was consistently above the average level over the past several centuries, there still seems to be not enough evidence to conclude that the late 20th-century ENSO intensity is exceptionally high, which is only seen in the tree-ring synthesis [36]. The sensitivity of ENSO to volcanic forcing (the cooling effect) is more evident in tree-ring synthesis [36] than in the fossil-coral records [39] and marine sediment records [32], probably because of the compensation of increased solar irradiance (the warming effect) [32]. Changes in the ENSO variability during the MCA and LIA remain unclear in current reconstructions: despite a mid-millennium shift captured in most continuous records [27,32,36], the sign of change is inconsistent. Moreover, most records do not show a significant difference between the end of LIA and the beginning of the industrial warming period. In short, substantial discrepancies still exist in current reconstructions, probably because of the large internal/spatial variation of ENSO and the relative weak forcing during the last millennium, illustrating the difficulty of quantitatively assessing ENSO sensitivity to the mean climate anomalies for this period.

Model Simulations and Model-Data Comparison
The volcanic impacts on the tropical climate and ENSO variability have first been studied using some early models (see a review by Timmreck [71]). It is proposed that the cooling induced by volcanic eruptions can force the climate system to a state where El Niño-like conditions are favored in the year following a large volcanic eruption, through the 'thermostat' mechanism [75,76]. This finding is supported by intermediate model [75] and GCM [77] simulations but with some exceptions [78]. According to the thermostat mechanism, a sudden cooling caused by a volcanic eruption induces a relative warming in the eastern equatorial Pacific, because of the dilution of cooling in the eastern equatorial Pacific by the upwelling and diverging surface currents. This relative warming induces a reduction of the trade wind and in turn further anomalous warming through the Bjerknes feedback [1].
In the recent efforts of the Community Earth System Model-Last Millennium Ensemble project (CESM-LME), however, the ensemble mean indicates different ENSO phase response to individual volcanic eruptions [66,68]. It is found that eruptions in the Northern Hemisphere lead to a subsequent weakening but distinct El Niño-like pattern; eruptions in the Southern Hemisphere enhances the probability of a La Niña; eruptions in the tropics has little systematic impact on the generation of either El Niño or La Niña events during the winter after the eruption. Furthermore, probability of the El Niño events is found to increase one boreal winter following Tropical and Northern eruptions and two boreal winters following all eruption types. The mechanism for the different ENSO responses to volcanic eruptions in different hemispheres in CESM-LME seems to be more complex than the 'thermostat' mechanism and related to the different responses of the Intertropical Convergence Zone (ITCZ) [68]. In the CESM that the strongest dynamic response to Northern eruptions is an equatorward migration of the ITCZ and associated trade wind weakening, which favors an El Niño initiation; the reverse may explain the enhanced La Niña probability after Southern events. Nevertheless, this interpretation is consistent with a study using another CGCM [79]. They suggest that after a high-latitude eruption in the Northern Hemisphere, the ITCZ shifts southward and implies a weakening of the easterly trade winds (weakest in the proximity of the ITCZ) along the equator in the central and eastern equatorial Pacific, thus leading to a reduction in the east-west temperature contrast across the tropical Pacific and an El Niño-like anomaly.
The CESM-LME also provide a more comprehensive view on ENSO variability under realistic transient forcings (solar activities, volcanic eruptions, land use, GHG, orbital parameters and ozone-aerosol) during the last millennium [66]. The ENSO variability in the CESM-LME is the dominant internal variability, with little sign of a systematic change and inconsistent linear trends in variance among all realizations. Some ensemble members exhibit an increasing trend in ENSO variance during the 20th century but the signal is not robust until the very end of the simulation. ENSO variance responses to varying solar irradiance and volcanic eruptions are obvious. CESM-LME can also detect substantial impacts of forced changes on ENSO diversity identified as the relative proportions of the Eastern Pacific (EP) and Central Pacific (CP) El Niño events [67]. It is found that anthropogenic changes in GHG and ozone/aerosol emissions [80] modify the persistence of both EP and CP El Niño events, although the effect is opposite; land use changes [81,82] most strongly amplify EP El Niño events while orbital forcing most strongly damps CP El Niño events [67].
Lewis and LeGrande and Hope et al. [64,65] also investigated the ENSO variability during the last millennium based on the Last Millennium experiments of Coupled Model Intercomparison Project-Phase 5 and Paleoclimate Modelling Intercomparison Project-Phase 3 (CMIP5-PMIP3) [83] with the focus on ENSO's spatial pattern and spectral characteristics, respectively. Evident decadalto centennial-scale El Niño-and La Niña-like episodes during the last millennium simulations were identified. ENSO amplitude in the historical period  was found to be greater than that of the last millennium, due to the increased atmospheric GHG concentrations. The gridwise correlation of Tropical Pacific (East, Central, West) temperature and precipitation was found temporally and regionally dependent.
The forced change in ENSO diversity in the CESM-LME was quantitatively investigated [67] using heat budget analysis adopted by Capotondi [84]. They suggested that EP El Niño events are most strongly forced by Land use/land cover changes that result in strongly amplified EP El Niño events, due to stronger westerly current anomalies and enhanced zonal advective feedback during an El Niño development. EP El Niño events are preferentially terminated under increased GHG, due to more efficient upwelling-induced cooling in a more strongly stratified equatorial eastern Pacific; the reverse effect of ozone/aerosol forcing enhances EP El Niño persistence. On the other hand, CP El Niño events also exhibit compensating influences from GHG and ozone/aerosol forcing, with increased GHG inducing a slight and marginally significant enhancement in zonal advective feedback that prolongs an El Niño; the opposite effect is shown from the latter. But CP El Niño events are most strongly forced by orbital changes, which preferentially cool the eastern Pacific during boreal winter and spring, enhance the zonal SST gradient during El Niño termination and lead to much more strongly negative zonal advective feedback [67].
Large uncertainty of ENSO variability in the CESM-LME and Last Millennium experiments of CMIP5-PMIP3 along with strong internal variability of ENSO in proxy records limits the confidence of a vigorous model-data comparison. The only qualitatively consistent conclusion among different simulations seems to be that the 20th century ENSO was intensified [37,[64][65][66]. Additionally, the hydroclimate responses show agreement in volcano-induced ENSO anomaly patterns between CESM-LME and tree-ring data following the eruption year but they somehow contradict during the eruption year [68]. Model simulations suggest that the interpretations of proxy records should be considered in conjunction with paleo-reconstructions from within the central Pacific to deconvolve changes occurring in the Nino 3.4 region from those related to local climatic variables [65].
In summary, relatively intensified ENSO after the industrial warming compared to the last millennium is recognized in both proxy records and model simulations while its magnitude is still uncertain. Consistent sensitivity of the ENSO variability to centennial mean climate anomalies such as the MCA and LIA is yet to be found from reconstructions and models.

Proxy Records
The tropical climate exhibits significant variations on orbital time scales (see a review by Schneider et al. [85]) due to changes in latitudinal and seasonal distribution of solar insolation [86]. The slowly modulated ENSO evolution in the Holocene has been studied most intensively. Early studies of archaeological records from the coast of Peru suggest that active ENSO variability similar to today can be traced back to the MH [45]. Later some continuous lake sediment archives offer the first prospect of understanding past ENSO dynamics and its responses to altered external forcing throughout the Holocene. Rodbell et al. and Moy et al. [22,46] presented a continuous record of ENSO variability based on measuring different reflectances between light-colored inorganic clastic and dark-colored organic-rich silt in sediment cores in the center of Laguna Pallcacocha in Ecuador. The change of intensity of red (light) color is associated with ENSO-induced precipitation and suggests an overall increase in ENSO frequency and magnitude that began about 6.5 ka (Figure 3b).
Another record was collected from a sediment core from hypersaline Bainbridge Crater Lake, Galápagos Islands [23]. Using a similar method to scan the optical features in short intervals, millennial-scale trend in El Niño frequency and intensity was analyzed. It is estimated that the MH El Niño activity was infrequent, with 23 strong to very strong and 56 moderate events occurring between 6-4 ka. El Niño frequency then increased between 4-3 ka (7 strong/very strong and 83 moderate events). The last 3 ka before the 1990s are characterized by an increase in the frequency (266) and intensity (80 strong/very strong) of events. Recently, Thompson et al. [26] presented long-term monitoring results of the local climate and limnology in the same region. In their revised interpretation of this record, brown-green, organic-rich, siliciclastic laminae reflect warm and wet conditions typical of El Niño events (previously interpreted as strong El Niño events), whereas carbonate and gypsum precipitate during cool and dry La Niña events and persistent dry periods (previously interpreted as moderate El Niño events). They found that ENSO events of both phases were generally less frequent between 6-4 ka relative to the last 1.5 ka and were intensified from 2-1.7 ka until 0.7 ka (Figure 3c).
Photospectrometer-measured changes in carotenoid pigments and fine-grained lithic reflectance in a marine sediment record spanning the last 20,000 years on the continental shelf off the coast of Peru are hypothesized to reflect El Niño-induced wind/upwelling and precipitation anomalies [27]. This result also suggests reduced El Niño activity during 8-5.6 ka (Figure 3d). Then the sediment flux rates recovered from the minimum to gradually reached a maximum around 3-2 ka; afterwards they were again persistently weakened. During the early Holocene, strong El Niño floods were abundant as inferred from strong lithic peaks recurring with ENSO frequencies.
Based on grain size data from a more recent continuous lake sediment record from El Junco Crater Lake, Galápagos Islands, Conroy et al. [24] produced a new record of past lake level variability and likely seasonal precipitation and ENSO frequency history through the Holocene. The grain size is interpreted as the amount of terrestrial versus aquatic organics in the lake sediments, corroborating lake levels indicated by carbon/nitrogen (C/N) ratios (Figure 3e,f). It implies increased precipitation intensity prior to~9 ka and after~4 ka, as well as a two-step increase in precipitation at 3 ka and 2 ka. Maximum Holocene precipitation and ENSO variability occurred between 2-1.5 ka. However, changes in the concentration and hydrogen isotope ratio of alga lipids in a sediment record retrieved from the same location [25] show only 5.6-3.5 ka was characterized by minimum El Nino activities and moderate ENSO events persisted during the early Holocene (Figure 3g).  [22]; (c) El Nino (thin red) and La Nina (thin blue) events [26] in Bainbridge Crater Lake [23]; (d) Lithic flux in a marine sediment from off Peru (% of maximum) [27]; (e-g) El Junco sand record and C/N ratio [24] and botryococcene concentration [25]; (h-j) Coral and mollusk synthesis of δ 18 O variability from Western Pacific (WP), Central Pacific (CP) and Eastern Pacific (EP) clusters [33]; (k) δ 18 O variability from speleothem [34]; (l-n) δ 18 O variability in marine sediment records from the tropical eastern Pacific [28][29][30]44].
Apart from these sediment records, a more direct proxy for ENSO is fossil coral records from the Northern Line Islands in the tropical central Pacific spanning the past 7 ka [40]. The interannual variability in fossil coral δ 18 O was compared against the late 20th century counterpart to assess the ENSO variance in the past. The coral records suggest highly variable ENSO activity in the last 7 ka, with no evidence for a systematic trend in ENSO variance. The average fossil coral ENSO variance is lower but with some rare cases exceeding or comparable to that in the 20th century. A reconstruction of ENSO in the eastern tropical Pacific Nino 1 + 2 region (90 • -80 • W, 0 • -10 • S) of the past 10,000 years derived from oxygen isotopes in mollusk aragonitic skeletons from Peru was presented by Carré et al. [42]. A composite Holocene record of ENSO variance and skewness from 180 mollusk shells and seven archaeological sites was analyzed. The results suggest that ENSO variability was volatile throughout the Holocene, with highest variance during the Holocene in the late 20th century, even if the influence of the 1982-1983 and 1997-1998 extreme El Niño events were excluded, and lowest variance between 5-4 ka. ENSO variance between 10-6 ka was variable but was not statistically different from that in the late Holocene.
Emile-Geay et al. [33] synthesized high-resolution, well-dated oxygen isotope records of paleo ENSO variance across the tropical Pacific spanning the Holocene (Figure 3h-j). Records from 65 sites in the core ENSO region in the Pacific were compiled, including the above coral and mollusk reconstructions [40,42]. Most of the individual records are comparatively short in coverage (50 years or less) and together they constitute combined 2162 sample years out of the past 10,000 years. The records were clustered into three groups according to their regions: the tropical western Pacific, the tropical central Pacific and the tropical eastern Pacific. Within each region, changes in ENSO variance were quantified by computing the ratio of fossil to modern variance in the 2-7-year band. Most records suggest paleo ENSO variance was lower than modern ENSO variance but there are no coherent millennial-scale changes among ENSO reconstructions from the three regions, possibly because paleo ENSO was highly variable. However, this dataset does show consistently reduced MH ENSO variance, for example, during the early Holocene and 6-2 ka in the Western Pacific (WP), during the past 7 ka and during 3-5 ka (a reduction of 64%) in the Central Pacific (CP) and during~4.6 ka in the Eastern Pacific (EP). A newly discovered seasonally-resolved stalagmite δ 18 O variance record [34] (Figure 3k) consistently shows a significant reduction in interannual variability during the MH (~3 ka and~5 ka) relative to both the late Holocene (~2 ka) and early Holocene (~7 ka but not applied to~8 ka).
So far, it is still difficult to draw a clear conclusion on the general trend of ENSO evolution during the Holocene from the diverse lines of evidence. High-resolution ENSO records as direct estimate of the tropical monthly SST (Figure 3h-j) and salinity and convection show no clear trend of the ENSO amplitude throughout the Holocene and recognize minimum ENSO strength during the MH [33,34]. In comparison, conservative sediment records suggest an increased frequency of ENSO events towards the late Holocene, although the timing of the onset of active ENSO somewhat differs (Figure 3b-f). Moreover, all sediment records suggest the period of 2-1 ka experienced extremely high frequency of ENSO activities, which is almost absent in high-resolution records. It should also be noted that other scarce glacial ENSO reconstructions which happen to cover some periods of the Holocene (see more details in Section 5.1) generally support an intensified ENSO during the Holocene (Figure 3l,m) but it should be treated with caution because of the small sample size. In short, the intensification of ENSO from the mid-to late Holocene seems to be qualitatively consistent among ENSO proxies but the ENSO evolution from the early to mid-Holocene remains controversial.
The discrepancies among proxies may result from the high irregularity of ENSO [87], the sparse data sampling and the caveats attached to data interpretation. The short length and small number of high-resolution oxygen isotope records denote a greater possibility that they may only capture some irregular signals in random ENSO events rather than its long-term trend [61]. For the sediment records, they are often not sensitive to La Niña (or drought in the sampling regions) conditions [22,24,25,27] and lack the resolution to resolve ENSO explicitly. In addition, the representation of physical processes that is the basis of paleo-ENSO proxy data interpretation could be altered in the past due to changes in ENSO centers of action [33,42] or ENSO seasonal phase locking [88]. For example, the spatial patterns of ENSO-induced atmospheric teleconnection for precipitation-sensitive proxy records, or the salinity-isotope relation for temperature-sensitive proxy records may have been different from what they are today. Ultimately, high-resolution paleoclimate data from ENSO centers of action (the central-eastern equatorial Pacific) must be expanded significantly to resolve these issues and provide particularly effective benchmarks for model simulations of ENSO variability.

Model Simulations and Model-Data Comparison
One of the most robust features of paleo ENSO simulation is the reduction of ENSO during the early and mid-Holocene. This feature has been simulated by most climate models, including the earlier simulations in intermediate coupled models [49,89] and a CGCM [50] and later across most CGCM simulations [90][91][92][93][94][95][96] and, systematically, in the PMIP2/3 simulations [19,20,97,98]. However, there are exceptions. For example, 1 out of the 5 PMIP2 experiments [20] and 8 out of the 28 PMIP2/3 experiments [98] do not show a significant reduction in ENSO amplitude during the MH [97].
In the transient simulation TRACE21, ENSO variance is modulated by the precessional forcing and it increases gradually by~15% from the mid-to late Holocene [61] (Figure 4e, red). The magnitude of the ENSO strengthening simulated by TRACE21 lies within the ensemble spread of PMIP2/PMIP3 simulations with a median of 10-15% [19] (Figure 4e, green bars). The prominent precessional modulation on ENSO evolution (the overall high-low-high trend) is also consistent with another transient simulation [63].
The ENSO intensification over the Holocene simulated by the transient and snapshot simulations is qualitatively consistent with several key ENSO proxies, for example, the lake sediment records [22][23][24]26] and the marine sediment record [29]. However, most of the simulations fail to capture the magnitude of the ENSO reduction inferred in the proxy, which could amount to 30-60% in the MH relative to the late Holocene. Furthermore, none of the models is able to simulate a minimum ENSO in the MH (i.e., the MH ENSO is weaker compared to both the early Holocene and the late Holocene) as suggested by recent high-resolution ENSO reconstructions [33,34,40,42]. This model-data discrepancy may result, partly, from the small sample size of the high-resolution ENSO proxies, relative to the high level of intrinsic variability of ENSO [61]. It can also be caused by model deficiencies, either in or outside the coupled ocean-atmosphere system. The coarse resolution in both the atmosphere and ocean models and the persistence tropical biases in current CGCMs [99] may contribute to the deficient ENSO simulation. A recent study further suggests that a vegetated and less dusty North Africa during the MH could also suppress ENSO significantly through its impact on the atmosphere and in turn atmospheric teleconnection into the tropical Pacific [95]. Moreover, the shift of ENSO centers of action [33,61] and the change in timing of ENSO events [88,93] further complicate this problem for both more effective data interpretations and more reliable model simulations. To resolve the model-data inconsistency regarding the evolution of the ENSO variability during the Holocene, much more high-resolution paleoclimate data from ENSO centers of action, model simulations with more complete physical processes and better model-data comparison/synthesis are still needed.

Forcing Mechanisms
For the modulation of ENSO variance in the early and mid-Holocene, the orbital forcing is a likely cause because ice volume and GHG concentrations changed little [80,100,101]. Starting from the early Holocene, precession induced a weakened (enhanced) annual cycle of insolation in the Northern (Southern) Hemisphere following approximately half of a precession cycle (Figure 4a, grey). In addition, when obliquity was reduced, leading to less (more) annual mean insolation over the tropics (mid and high latitudes) and a usually secondary but globally weakened annual cycle in insolation.
Multiple hypotheses have been put forward to explain the reduction of ENSO during the early to mid-Holocene by orbital modulation. Some studies have attributed the weakened ENSO to an enhanced seasonal forcing through the nonlinear mechanism of frequency entrainment [102]: when a strong equatorial annual cycle forcing is present, the amplitude of ENSO oscillation tends to be reduced due to its nonlinear interaction with the annual cycle. The enhanced seasonal forcing, however, could be attributed to either the seasonal insolation forcing locally over the Pacific [49], the seasonal cycle of cloudiness in the off-equatorial regions and the resulting meridional SST gradient in the eastern equatorial Pacific [63], or the remote forcing from the enhanced NH monsoons [50,98], all under precession modulation. The latter also induces enhanced trade winds over the equatorial central Pacific and suppresses the development of warm ENSO events through upwelling cooling [20,50,90,91,94]. However, it should be pointed out that the inverse relationship between the amplitude of the annual cycle and ENSO strength during the Holocene is not found in high-resolution reconstructions [33] or in some other model simulations [33,61,90,96].  Orbital forcing may modulate the tropical mean climate and, in turn, the strength of ocean-atmosphere feedback and change ENSO intensity in the Holocene. A recent study in a transient CGCM simulation [61] suggests that the precessional modulation of the coupled ocean-atmosphere feedback dominates the ENSO evolution in the Holocene. This mechanism of the reduced ENSO in the MH involves the subduction of South Pacific subtropical water, the weakening of the equatorial thermocline [103] and in turn positive ocean-atmosphere feedbacks [61]. Relative to the late Holocene, the boreal winter (or austral summer) insolation in the early and mid-Holocene is reduced locally in the Southern Hemisphere (SH). The reduced seasonal cycle of insolation in winter in the SH warms the SST in the SH and, in turn, the subducted water in the permanent thermocline [104] in the South Pacific, which is eventually transported to the equatorial thermocline, forming a warmer subsurface water in the tropics, weakening the oceanic stratification and the strength of positive ocean-atmosphere feedbacks (dominated by Ekman upwelling feedback). This mechanism is now supported by a quantitative examination of ENSO feedbacks employing Bjerknes stability (BJ) analysis [105][106][107], as well as by comparing full forcing and individual forcing sensitivity experiments. Quantitative analyses further illustrate that the precessional forcing can reduce ocean-atmosphere feedbacks through other mechanisms, such as the ocean's SST response to wind stress forcing, the increased thermal damping and weakened atmospheric response to SST [61,98]. For example, the annual mean cooling in the equatorial Pacific in the MH reduces the saturation moisture and in turn atmospheric response to a given SST anomaly, leading to a reduced Bjerknes feedback [89].
It has been proposed that ENSO can be reduced in the MH by the reduced intensity of extratropical atmospheric synoptic variability [92], in addition to local mechanisms in the coupled ocean-atmosphere in the tropical Pacific. The extratropical atmospheric variability acts as a stochastic forcing to ENSO through equatorward propagation of the coupled wind-evaporation-SST feedback [108] or the so-called Pacific Meridional Mode (PMM) [109,110]. The PMM represents the coupled teleconnection process through which the North Pacific atmospheric variations may trigger ENSO events in the present day observation [111].
ENSO can be further weakened by mechanisms involved outside the ocean-atmosphere system. A recent study using a set of GCGM experiments with different vegetation covers and dust concentrations [95] shows that ENSO in the MH can be weakened by an increased vegetation cover and a reduced dust aerosol level over the North Africa. They highlight a link between the intensity of the climatological West Africa Monsoon (WAM), Walker circulation strength and position and their indirect effects on ENSO activity. The increased vegetation cover strengthens the WAM, which then shifts the ascending branch of the Walker circulation in the Atlantic westward in boreal summer, leading to an increased divergence in the central equatorial Pacific through the westward propagation of atmospheric equatorial Rossby waves. The associated wind field deepens the equatorial thermocline in the east and causes a shoaling of the thermocline in the central Pacific, both reducing the Bjerknes feedback. Given that active vegetation cover and the associated dust emission reduction have not been included in previous modeling studies on ENSO in the Holocene, this study suggests that ENSO may be suppressed further by feedback processes that are missing in current models and the weaker suppression of MH ENSO activity in current CCGMs (~10%) than in proxy records (30-60%) may be caused partly by this deficiency of current CGCMs [95].
Recent modeling studies further suggest that, in addition to the reduction in the variance, ENSO is also affected in its phase-locking during the MH such that its peak season is delayed towards spring by approximately 2 months, resulting from changes in the seasonality of ENSO feedbacks [93]. This feature of ENSO phase shift is also implied by paleoclimate proxies from the central Pacific [88]. Under the MH insolation forcing, a delayed ENSO phase is remotely forced by the West Pacific, where a weakening of the trade winds in early boreal spring in the MH initiates an anomalous downwelling annual Kelvin wave, which reaches the eastern Pacific and deepens the thermocline during the boreal fall [112], weakens the upper ocean stratification and results in reduced upwelling feedback during ENSO development season [93].
Finally, the transient simulation of Liu et al. [61] indicates that the orbital time-scale ENSO variability is the most sensitive to the precessional forcing, while the obliquity effect does not play a significant role, qualitatively consistent with a previous transient simulation of the last 140 ka [63].
In summary, ENSO suppression in the early to mid-Holocene seems to be robust across models. The mechanism of the suppression, however, is less certain, ranging from reduced coupled feedbacks, interaction with the annual cycle, the impact from the extratropics, to factors even outside the ocean-atmosphere system. The reduced ENSO in the MH in the models is consistent with proxy records. However, there is little evidence in models so far that ENSO intensity reached minimum in the MH, as implied from some high-resolution proxy records.

Proxy Records
Following the LGM [113] the deglacial climate change was accompanied by a retreat of ice sheets over the globe, notably over North America and Northwestern Eurasia [101], a rise of atmospheric CO 2 by~100 ppm [114] and changes of orography, surface conditions (e.g., albedo) and sea level. These climatic forcing factors have been shown to have major influences on the atmospheric and oceanic circulations [115][116][117][118][119], including the coupled ocean-atmosphere system in the tropics [120][121][122]. Moreover, meltwater fluxes from the retreating ice sheets were released to the North and South Atlantic (see a review by Carlson and Clark [123]), causing abrupt climate changes, for example, Heinrich event 1 (~19-14.5 ka), meltwater pulse 1A (~14.6 ka), the Bølling-Allerød (~14.7-12.7 ka), the Younger Dryas (~12.9-11.7 ka) and the 8.2 ka event (~8.4-8.2 ka). All these make the ENSO evolution in the last deglacial period more complex.
The paleoclimate records of ENSO variability before the Holocene are scarce. Some rarely preserved corals from Papua New Guinea suggest that the ENSO variability has existed for the past 130,000 years of glacial and interglacial times [41], which is consistent with recently retrieved giant clams (Tridacna sp.) from the Western Pacific warm pool [124] over 60 ka old. For the last deglaciation period (since~21 ka), there is only one continuous high-resolution ENSO record, a marine sediment core from off Peru [27]. Nevertheless, non-continuous marine sediment records have been used to infer ENSO variability across a wide range of times, albeit covering only some intervals separated by gaps [28][29][30][31]. A common method adopted by these marine records is individual foraminifera analysis, which compares the variance of samples of repeatedly measured δ 18 O in foraminifera shells between a past time period and the modern time or the late Holocene.
The continuous ENSO record from Peru [27] (Figure 3d) suggests that El Niño activity was weaker during the LGM than the present-day one and it started to strengthen around 17 ka, coinciding with the retreat of the Laurentide Ice Sheet during Heinrich event 1. ENSO experienced centennial-scale variability during the deglaciation, as Peru climate was characterized by centennial-scale drier (less El Niño activity) and wetter (more El Niño activity) periods; then ENSO strengthened sharply throughout and after the Younger Dryas before it weakened again towards the MH.
Leduc et al. [28] performed repeated isotopic analysis on δ 18 O in shells of thermocline-dwelling Neogloboquadrina dutertrei (N. dutertrei) at eight intervals of the last 50 ka (Figure 3l). Their record suggests marginally reduced ENSO activity during the LGM compared to the Holocene, while it increased during the last deglacial period, with the lowest ENSO reached during the late LGM. The highest ENSO variance was at the beginning of Marine Isotope Stage 3 (MIS3) during Dansgaard-Oeschger (DO) interstadial 14 (~52.1 ka). Furthermore, ENSO activity did not respond to rapid climate shifts such as the Heinrich events and DO interstadials, especially during Heinrich event 4 (~39.6 ka) and DO8 (~38.9 ka).
Sadekov et al. [30] also investigated total variance in specimens of subsurface-dwelling N. dutertrei in the eastern equatorial Pacific at intervals of the modern (1.6 ka), the early Holocene (9.1 ka), deglaciation (12.5, 15.1 and 17.9 ka) and the Last Glacial Maximum (20.7 ka) (Figure 3m). ENSO variability was found to be~20% higher than the modern-day during the LGM and increased significantly during early deglaciation (for example, at 17.9 and 15.1 ka interval). It then declined, reaching levels close to modern at 12.5 ka, lowest in the early Holocene (~27% lower than modern level at 9.1 ka).
Using individual foraminifera analysis of the surface-dwelling Globigerinoides ruber (G. ruber) in the eastern equatorial Pacific over thirty-nine narrow intervals in the last 25 ka (Figure 3n), Koutavas et al. [29,44] suggested that LGM ENSO was 27% greater than the late Holocene mean (0-4 ka) and 7% greater than the late 20th century. The highest LGM variance at 20.7 ka exceeded the lowest MH value at 5.3 ka by four-fold.
It should be noted that the variance of the foraminifera temperature, if treated as the strength of ENSO, can be used to infer the broad trend of ENSO variability through different time periods [28][29][30]44]. However, since these snapshots are distributed most likely randomly in time, their variance most likely represents the total variability across all the time scales including the seasonal cycle, interannual variability, decadal and longer variabilities. A test using a last millennial sediment shows that there is only a slight influence of annual cycle on the total variance [32]. However, this representation is likely to depend on location, time period and the inferred variable. In fact, statistical and modeling studies suggest that changes in the total variance during the LGM can be very different from changes at the interannual time scale, especially over the eastern equatorial Pacific [56,125] where the annual cycle is prominent [126].
Ford et al. [31] used both surface-(Globigerinoides sacculifer) and subsurface-dwelling (Globorotalia tumida) foraminifera in cores in the equatorial Pacific to reconstruct tropical SST variability during the late Holocene (<6000 years ago) and the LGM. They claimed that the quantile-quantile (Q-Q) plots of populations of ocean temperature from individual foraminifera analysis were able to separate changes in ENSO variability from seasonality during the LGM. The LGM records in the eastern equatorial Pacific suggest a reduced ENSO variability and an enhanced seasonal cycle compare to present day. They also suggest although thermocline temperature proxy shows enhanced variability during the LGM, it might not be the best indicator of ENSO strength because there could be substantial changes in thermocline depth.
In summary, although evident ENSO variability existed within the glacial climates [41,124], proxy records so far show contradictions on the strength of ENSO during the glacial and deglacial periods. For example, surface [29] and subsurface [30] proxies have been interpreted as strengthened ENSO variability at the LGM, while other proxies [28,31] show weakened ENSO. These contradictions suggest that ENSO strength during the glacial period thus remains unclear.

Model Simulations and Model-Data Comparison
As with the paleo reconstructions, there is no consensus on the intensity of ENSO at the LGM versus the late Holocene among CGCMs [19,20,56,94,127,128]. The PMIP2 simulations show contrasting ENSO strength at the LGM compared to the present day: there is one model showing an ENSO weakening of more than 30%, one model showing a strengthening of more than 10% and two models showing slightly insignificant weakening (<5%) [20]. The median of 10 LGM simulations from the PMIP2/PMIP3 shows a weakening of the ENSO variability around 10% but the ensemble spread is twice as large as that of the MH simulations [19] (Figure 4e, green bars). Despite the large uncertainty between different climate models, there might be some consistency between different generations of climate models from a single group, highlighting the model dependence of ENSO sensitivity. For instance, the simulated weakened ENSO at the LGM in a recent isotope-enabled CGCM [56] is in qualitative agreement with results from its previous versions [20,127].
TRACE21 also simulates ENSO evolution during the entire last deglaciation [61] (Figure 4e, red), which is characterized largely by a precessional modulation, with an initial weakening at the early Holocene and a subsequent strengthening towards the late Holocene. There is pronounced millennial variability of ENSO intensity during the early deglaciation, associated with meltwater pulses in the North Atlantic. The meltwater forcing tends to intensify the ENSO amplitude, agreeing with a study using multiple models [129]. Furthermore, the presence of large continental ice sheets during the glacial period tends to suppress the ENSO variability, while the lower concentrations of GHG also at glacial time seem to strengthen it, resulting in a modestly reduced ENSO at the LGM relative to the late Holocene.
The model-data comparison regarding glacial ENSO strength is rather challenging. First, as discussed above, information from reconstructions using individual foraminifera analysis is inconclusive and these reconstructions are usually subject to complicated uncertainties. For example, changes in total variance from individual foraminifera analysis could reflect changes in the annual cycle [56,125], vertical migration of the foraminifera [56], variations in the isotope composition of seawater [54], intershell variability and changes in seasonality of sedimentation rate [130]. Second, all climate models have to employ lots of parameterizations and they all have systematic biases in simulating the mean climate in the tropics and the extratropics, for example, the double ITCZ problem and the bias in simulating stratiform clouds in the subtropics [99]. Last but not least, variables from proxy records are usually not directly simulated in climate models, making it hard to address the inconsistencies between models and reconstructions and among themselves.
Using a state-of-the-art water isotope-enabled climate model, Zhu et al. [56] made the first direct model-data comparison regarding ENSO strength at the LGM. By realizing the possible uncertainties in reconstructions and combining evidence from both model and data, they argued that the ENSO strength at the LGM was most likely weaker than that at the preindustrial. Their water-isotope enabled modeling suggests that changes in the total variance from records in the eastern equatorial Pacific using the individual foraminifera analysis [29] could merely reflect changes in the annual cycle instead of ENSO variability. Furthermore, their modeling study suggests that the subsurface records using the individual foraminifera analysis [30] could be subject to substantial complications from the possible migration of the habitat depth of the thermocline-dwelling foraminifera. Consequently, the weaker ENSO at the LGM suggested by the isotope-enabled simulation is in a qualitative agreement with the ensemble mean of recent generation of climate models [19,20,127], the transient simulation TRACE21 [61] and other proxy reconstructions [28,31,41,131].

Forcing Mechanisms
Compared to the Holocene, glacial and deglacial climate changes are associated with more complicated forcings, for example, the GHG, ice sheets, freshwater forcing, surface boundary conditions and ocean bathymetry. The ENSO strength at the glacial period and its deglacial evolution should be a result of interactions with these climatic forcings. The evolution of ENSO during the last deglaciation is explained by response to GHG, ice-sheet and melt water forcings through the nonlinear frequency entrainment mechanism in a CGCM [61] (Figure 4g,h). During the deglaciation, the increased GHG concentrations led to an asymmetric annual mean warming in the tropical Pacific with a stronger warming to the north of the equator [132], which enhanced the equatorial asymmetry and in turn the annual cycle, favoring weakened ENSO amplitude [133]. In contrast, a retreat of continental ice sheets in the NH led to a northward displacement of the NH jet stream, increased sea-ice coverage over the North Atlantic and North Pacific and a subsequent cooling over the equatorial northeast Pacific and in turn a weaker annual cycle, favoring a stronger ENSO [121]. These competing effects on ENSO by the retreat of ice sheets and increased GHG concentrations could explain the large uncertainty of the ENSO strength at the LGM across multiple models. This interpretation is consistent with PMIP2 simulations which show less constraint of the seasonal cycle on the El Niño characteristics at the LGM than in the MH, probably because of competing heat fluxes and local responses and the resultant zonal SST gradient change within the Pacific Ocean [20].
So far, reconstructions on the ENSO strength during abrupt climate changes (e.g., the Heinrich Event 1, the Bølling-Allerød warming and the Younger Dryas) are still lacking. These abrupt climate changes are thought to be related to variations in the Atlantic Meridional Overturning Circulation (AMOC) [134]. However, an enhanced ENSO variability in response to a weakening of the AMOC seems to be a robust feature in almost all the climate models [61,129] (Figure 4e, red; Figure 4g, blue). The models suggest that a meltwater pulse in the North Atlantic weakens the AMOC, shifts the ITCZ southwards and creates a more symmetric annual mean climate in the eastern equatorial Pacific. The more symmetric mean climate is associated with a weaker north-south difference in the SST and a weaker equatorial annual cycle, resulting an amplified ENSO through the nonlinear mechanism of frequency entrainment. The influence of north-south asymmetry and the shift of ITCZ may not be as sufficient during the LGM due to either only subtle southward-shift of the ITCZ [135], or the remote forcing of ice-sheets [121].
In contrast to these studies that rely on frequency entrainment mechanism, a recent modeling study [56] attributes the weakened ENSO at the LGM to weaker coupled ocean-atmosphere feedbacks. Feedback analysis shows a decrease in all three positive feedbacks, namely the Ekman pumping feedback, thermocline feedback and zonal advection feedback. A further quantification of these feedbacks reveals that the weakened feedbacks are caused by the reduced mean thermal stratification, the deepened thermocline and the generally reduced sensitivity of the ocean dynamic response to surface wind anomalies (probably due to the increased thermodynamic and dynamic inertia). Additional sensitivity experiments with single climatic factors suggest that both lowered GHG and the presence of the glacial ice sheets contribute to the weakened ENSO at the LGM but the role of GHG is in contrast to that in the transient simulation in Liu et al. [61].
In summary, the response of ENSO to deglacial climate changes remains uncertain in both models and reconstructions. Also, mechanisms responsible for the glacial and deglacial ENSO changes range from the nonlinear frequency entrainment to different details of coupled ocean-atmosphere feedbacks. Model simulations seem to suggest that the presence of glacial ice sheets tend to suppress ENSO variability, nonetheless only two models have done this exercise so far. The role of GHG forcing is still debatable in models. The response of ENSO to meltwater forcing is more robust in models, with a meltwater forcing inducing stronger ENSO. This response, however, has yet to be confirmed by observations.

Proxy Records
Sustained periods of global warmth beyond the last glacial period can provide extraordinary insights into the ENSO projections under future global warming. So far, there are very limited records allowing us glimpses of ENSO changes during these periods, such as the most recent warmer climate interval at the early/mid Pliocene (~4.5-3 Ma) [136,137] when the global temperature was~3 • C higher than that of the pre-industrial and the early Eocene (56-50 Ma) [138] when the global temperature increase was~13 • C warmer than modern temperatures [139].
Two monthly-resolved fossil porites coral records around 3.5-3.8 Ma in the Philippines reveals variability that is similar to the present ENSO variation [35]. Paleotemperature proxies from the early Eocene that can be used to assess ENSO variability have not been gathered from the tropical Pacific. Only annually-resolved varved sediments from midaltitude lakes show a strong interannual signal in varve thickness spectrum [140]. These pieces of evidence seem to support the presence of active ENSO in warm periods.
One highly debatable topic is the paradox of a 'permanent El Niño' at the Pliocene (see a review by Fedorov et al. [141]) but it should be regarded as an issue of the mean climate rather than climate variability. Wara et al. [142] first showed that during the warm early Pliocene (~4.5-3 Ma) the mean zonal SST contrast across the tropical Pacific was like a continual modern El Niño, using marine sediments from the western and eastern equatorial Pacific of the last 5 Ma. The reduced zonal SST gradient (less than 2 • C) during the early Pliocene is also supported by some other marine sediment records [141]. Extratropical land reconstructions in North America also show the distinctive climate signatures during the Pliocene warmth that resemble the teleconnection pattern of an El Niño [143]. However, a SST reconstruction using biomarkers suggests a zonal SST difference in the equatorial Pacific as large as 3-4 • C [144], in contrast to the 'permanent El Niño' hypothesis. Although this record showing a strong SST gradient was challenged later by Ravelo et al. [145], who argued that their results could have been misinterpreted because of the different between Quaternary and modern baselines and the coarse temporal resolution.
Currently observational evidence renders a 'permanent El Niño' (which also implies no ENSO variability) unlikely during these past warm periods in deep time [8]. However, the fundamental ENSO characters in response to past global warming and possibly weakened instability (weaker east-west SST contrast and thus weaker Walker circulation) remains inconclusive due to the lack of high resolution proxies.

Model Simulations and Model-Data Comparison
Modeling efforts have been made to answer the question of whether a 'permanent El Niño' existed and whether ENSO fluctuations were inactive in the past warm periods. So far models have difficulties in producing a 'permanent El Niño' mean state. Prescribed mid-Pliocene boundary condition fails to drive substantial reduction in zonal SST gradient without modifications of model physics [146]. Although earlier climate models were able to produce an active ENSO in past warm periods in deep time [140,147,148], both non-weaker [149,150] and weaker [151][152][153] ENSO variations compared to the pre-industrial were reported. The PlioMIP ensemble [60] shows overall reduced ENSO variability (7 models out of 8) but no robust changes in the ENSO flavor [146]. Song et al. [154] presented some model sensitivity experiments with seaway changes during the Pliocene. They found a reduced ENSO amplitude (15-20%) due to the closing of the Panama Seaway in the Kiel Climate Model (KCM), while prescribing changes in the Indonesian Passages has less impact. The ENSO sensitivity to the specific geometry changes highlights the challenge in model configurations to provide a reasonable representation of climatic features during these earlier warm periods.
The ENSO forcing mechanisms during the Pliocene warm period, which is often seen as a potential testbed for future ENSO projection, have not been studied thoroughly. Although several hypotheses [142,151,155] have been proposed to explain the mean SST gradient along the equator in the Pacific that is regarded very weak or absent (a 'permanent El Niño'), when it comes to the climate variability, however, there is a general lack of diagnostics due to the disagreement in ENSO variability in the Pliocene warm period [146,[149][150][151][152][153]. Particularly puzzling is why a mean climate with a 'permanent El Niño' condition can still maintain a robust ENSO cycle [35,140]. Some clues can be found in one special case which is missing in other modeling studies: closing Panama Seaway closing and/or narrowing Indonesian Passages during the Pliocene [154]. The closing of the Panama Seaway in the KCM can lead to an enhanced meridional SST gradient thus strongly strengthening southerly wind stress in the eastern Pacific. Other responses include shoaling of the thermocline across the tropical Pacific, slightly decreased zonal equatorial SST gradient and zonal wind stress. The annual cycle is also strengthened mainly due to the stronger southerly wind stress and can suppress ENSO through frequency entrainment [102]. In addition, BJ analysis shows relatively small change between the modern and Pliocene Panama Seaway experiments. It underlines some external forcing factors that can have considerable impacts on the tropical Pacific mean climate and ENSO during the Pliocene as documented in the proxy reconstructions, yet to be more systematically studied.
In summary, the paleo ENSO variability forced by natural variation of GHG in warm periods beyond the last glacial such as the early Eocene and the early/mid Pliocene likely exist, which provides possible analogues of its future behavior.

Discussion and Conclusions
Significant advances have been made in the paleoclimate ENSO studies in the last few decades. From the perspective of reconstructions, an increasing number of high-resolution proxy records have been recovered and synthesized across the tropical Pacific to provide a better picture of the past ENSO variability; new statistical methods have been introduced to target paleo ENSO more effectively. From the modeling perspective, numerical experiments using a newer generation of CGCMs with more complete physics have been widely coordinated and analyzed; continuous simulations of past climate for tens of thousands of years in a CGCM is now a possibility; isotope-enabled climate models have offered new opportunities to make direct model-data comparison and improve interpretations of proxy records. Held together through more intensive collaborations between communities, paleoclimate data scientists and modelers are now more likely to create conditions for unbiased interpretations of the results and portray a more confident outcome.
Nevertheless, we should keep in mind not to overestimate how much we actually understand ENSO under climate changes from a paleoclimate point of view. The first and perhaps the simplest case is the MH ENSO (Section 4). Sediment records disagree with recent high-resolution water isotope reconstructions on the magnitude of the ENSO variance reduction during the MH and it is not yet known very well why the latter suggest a 'mid-Holocene minimum.' CGCM simulations cannot quantitatively produce the reduced ENSO activity during the MH suggested by proxy data, as some models even generate a stronger MH ENSO. There are many proposed mechanisms but the underlying details are conflicting. For example, a simulated inverse relationship between the amplitude of the seasonal cycle and ENSO-related variance in sea surface temperatures in the eastern equatorial Pacific contradicts that of the reconstructions. Some new but important dynamic processes have not been thoroughly studied yet, such as the response of spatial shift (diversity) and seasonality of ENSO under external forcing and only after such study could we generate more convincing MH ENSO scenarios. The understanding of the last millennial ENSO and a clear detection of forcing signals are obscured by ENSO's highly variable features and somewhat minor forcings and whether the 20th century is a period of the strongest ENSO during the past millennium remains disputed. There is an additional need for rigorous examinations of more challenging cases such as the glacial and deglacial ENSO, or further back in time and isotope-enabled CGCM simulations could be the next application of reproducing these underlying comprehensive physical processes.
Past experience suggests that the next steps regarding paleo ENSO studies are the synthesis of much more high-resolution paleoclimate data from ENSO centers of action and new statistical methods to avoid any misinterpretation, the development of the next generation of Earth system models with more complete physical processes (e.g., reducing model biases such as the cold tongue bias, the use of isotope-enabled model, model with predictive vegetation module and producing ensembles, etc.) and above all, the continuous interaction between paleoclimate data scientists and modelers.