Research ArticleGEOLOGY

Pacific warm pool subsurface heat sequestration modulated Walker circulation and ENSO activity during the Holocene

See allHide authors and affiliations

Science Advances  14 Oct 2020:
Vol. 6, no. 42, eabc0402
DOI: 10.1126/sciadv.abc0402


Dynamics driving the El Niño–Southern Oscillation (ENSO) over longer-than-interannual time scales are poorly understood. Here, we compile thermocline temperature records of the Indo-Pacific warm pool over the past 25,000 years, which reveal a major warming in the Early Holocene and a secondary warming in the Middle Holocene. We suggest that the first thermocline warming corresponds to heat transport of southern Pacific shallow overturning circulation driven by June (austral winter) insolation maximum. The second thermocline warming follows equatorial September insolation maximum, which may have caused a steeper west-east upper-ocean thermal gradient and an intensified Walker circulation in the equatorial Pacific. We propose that the warm pool thermocline warming ultimately reduced the interannual ENSO activity in the Early to Middle Holocene. Thus, a substantially increased oceanic heat content of the warm pool, acting as a negative feedback for ENSO in the past, may play its role in the ongoing global warming.


The equatorial eastern Indian and western Pacific oceans with a persistent sea surface temperature (SST) above 28°C [termed the Indo-Pacific warm pool (IPWP)] represent a major oceanic heat source for the atmosphere, characterized by deep atmospheric convection accompanied with heavy rainfall (1). The heat storage in the IPWP is essentially formed by the accumulation of warm surface waters driven by the equatorial trade winds (2) and modulated by the convergence of subsurface ocean heat anomalies from North (W-E) and South subtropical and eastern Pacific (Fig. 1A) (3). The west-east thermal asymmetry across the tropical Pacific and associated Walker circulation play a key role in both the interannual variability of El Niño–Southern Oscillation (ENSO) (2, 4) and the decadal to multidecadal Pacific climate changes (3, 4, 5). In association with the greenhouse warming over the 20th century, the equatorial Pacific W-E thermal gradient was possibly reduced, the Walker circulation slowed down (6), and the ENSO variability increased (4, 7). More recently, the slowdown of surface air warming between ~2000 and 2014 AD (aka “global warming hiatus”) was featured by substantial cooling of equatorial Pacific SST and strengthening of the zonal thermal gradient and the Walker circulation (3, 5), which has been attributed to enhanced heat storage in the equatorial Pacific thermocline (8, 9). It is not fully resolved, however, over decadal and even longer time scales, how the upper-ocean heat anomalies will vary with changes in the Walker circulation, ENSO activity, and the shallow overturning circulation from the subtropics. Understanding these relationships can greatly improve the prediction of the future climate change.

Fig. 1 Time series of thermocline and SST anomalies in the IPWP compared to global climate indices during the past 25,000 years.

(A) Site locations of paired SST and thermocline water temperature (TWT) records (white circles) and SST-only records (blue triangles) (table S1). Shadings indicate temperatures at 120-m water depth. (B) Mean TWTA (red) of the IPWP records. Solid black arrows mark the two major warming phases of TWTA between 22 to 19 ka and 13 to 11 ka, respectively. (C) Precession (dashed purple) and obliquity (orange) parameters (47). (D) Atmospheric pCO2 derived from West Antarctic Ice Sheet Divide ice core [gray dots; (26)]. ppmv, parts per million by volume. (E) Mean SSTA (blue) of the IPWP records and the global mean SST anomaly [ΔT, dark gray line; (48)]. (F) Mean IPWP G. ruber δ18OG anomaly (δ18OG-A, green) and LR04 benthic δ18O stack [gray line and symbols; (49)]. Shadings of proxy records show the 1σ SD. Vertical dashed lines denote the timing of the deglacial onset of SST (~19 ka, blue) and TWT (~22 ka, red), the onset of the second deglacial warming step (gray), and the Early Holocene peak of TWT (EH-peak, ~10.8 ka, red). Dotted red arrow denotes the Middle Holocene peak of TWT (MH-peak, 7 ka). B.P., before the present.

On longer time scales, from the Last Glacial Maximum [LGM; 19 to 25 ka (kilo annum, thousand years before 1950 AD)] to the Early Holocene (9 to 12 ka), a warming of up to 4°C has been observed in the thermocline temperature records over the IPWP (10, 11, 12), possibly leading to a reduction of ENSO activity in the Early to Middle Holocene (13, 14). In addition, stalagmite records from Borneo suggest that Walker circulation was relatively weak during the last deglacial and strengthened in the Middle Holocene (4 to 7 ka) (15) when ENSO activity was suppressed (16). However, while the Borneo stalagmite record suggests similar-to-modern ENSO activity during the Early Holocene (16), other proxy records of ENSO activity from eastern equatorial Pacific indicate little or no ENSO activity during this period (17, 18). Such a discrepancy is largely caused by the sparsity of sedimentary archives and thus leads to an incomplete understanding of past changes in the upper-ocean heat structure of the IPWP.

Here, we present a compilation of more than 30 sedimentary proxy records (3 from this work) from the equatorial Pacific (Fig. 1A and table S1) to comprehensively examine changes in the upper-water (thermocline and above) temperature of the IPWP over the past ~25,000 years. We excluded five records [e.g., those from southwest Sumatra (19)] that are predominately influenced by local processes such as upwelling and do not reflect the general characteristics of IPWP’s subsurface (table S1). The records with an average temporal resolution of ~150 years constrained by a total of 217 radiocarbon dates (see Materials and Methods and table S2) covering the course of the LGM through the Holocene. The thermal structure of IPWP is examined by reconstructing temperatures of the thermocline water (TWT) and sea surface (SST) using shell Mg/Ca of two planktic foraminifera, the upper-thermocline dweller Pulleniatina obliquiloculata (20, 21) and the mixed-layer dweller Globigerinoides ruber (see Materials and Methods) (21). To minimize the possible interlaboratory and intercalibration biases and the effect of different cleaning protocols (see the Supplementary Materials), we calculated the SST and TWT anomalies relative to the average value of each temperature record over 6 to 10 ka (denoted as SSTA and TWTA; Fig. 1, B and D). The records of G. ruber δ18O from 22 cores are processed in the same way to acquire the mean δ18OG. ruber anomaly (δ18OG-A; Fig. 1F) and are calculated for seawater δ18O by subtracting the amount related to changes in local temperature and global ice volume [see the Supplementary Materials (22)]. Three deep thermocline temperature records of the eastern equatorial Pacific, estimated by the shell Mg/Ca of Neogloboquadrina dutertrei, are also analyzed (table S1).


Mean TWTA and SSTA variations

The mean SSTA continuously warms by ~2.8° ± 0.6°C since ~19 to ~10 ka and cools by ~0.3° ± 0.5°C from ~9 to 0 ka (Fig. 1E), consistent with previous estimates (23, 24). The deglacial onset of positive SSTA occurs at ~18 ka with site-specific differences (25), which is generally synchronous with the onsets of the atmospheric pCO2 rise (Fig. 1D), global mean SST warming (Fig. 1E), and decreases in the IPWP mean G. ruber δ18O and the global benthic δ18O stack (Fig. 1F).

The mean TWTA warms by 3.0° ± 0.6°C from ~22 to 11 to 9 ka and cools by 1.0° ± 0.7°C from 9 to 0 ka (Fig. 1B). The TWTA at 11 to 9 ka is highest and synchronous with the precession minimum and obliquity maximum (Fig. 1C). The onset of deglacial warming occurs at ~22 ka in TWTA, in phase with the turning point of precession parameter (Fig. 1C) and precedes the deglacial pCO2 rise (26) by ~4000 years (Fig. 1D). The deglacial mean TWTA warming mainly occurs in two phases: a first warming between 22 and 19 ka, coeval with the initial decrease of precession parameter, and a second warming between 13 and 11 ka, coinciding with the minima of precession parameter (Fig. 1B). The later warming phase is also synchronous with the final phase in the deglacial rise of the atmospheric pCO2 (Fig. 1D). The overall trend and the timing of the deglacial warming illustrate that orbital-driven insolation forcing controls the TWT change in the IPWP.

Besides the warm TWTA peak at ~11 ka, a second peak is found around 7 to 6 ka (Fig. 1B). TWT features observed in sites from open ocean differ from those within the Maritime Continent waters. The TWTA from open-ocean sites gradually warms from 22 ka and peaks at 11 to 10 ka (Fig. 2A), which we define as the Early Holocene peak (EH-peak) type. The near-equator TWTA records from the Maritime Continent waters are characterized by a rise after 15 ka and a Middle Holocene warm peak around 7 ka (Fig. 2B), defined as the Middle Holocene peak (MH-peak) type. The principal components analysis confirms the distinction, with a first principal component (Fig. 2C) yielding positive loadings for all records (Fig. 2D) and a second principal component (Fig. 2C) yielding different signs of loadings among the sites (Fig. 2E). The linear combinations of the first and second principal components resemble the two types of TWTA change (Fig. 2F), suggesting that “PC1 − PC2” represents the feature of EH-peak type TWTA and “PC1 + PC2” represents the MH-peak type.

Fig. 2 The two types of thermocline temperature anomaly (TWTA) records in the IPWP since the LGM.

(A) Average TWTA (brown) and the original TWTA records of the open-ocean sites. (B) Same as (A) but for the near-equator sites in the Maritime Continent waters. The TWT records in (A) and (B) are defined as (EH-peak and MH-peak types), respectively. (C) First (blue) and second (red) principal components (PCs) of all the TWTA records. PC1 and PC2 explain 62 and 16% of the total variance, respectively. (D and E) Loadings of PC1 (D) and PC2 (E) for each site. (F) Linear combinations of PC1 and PC2 that resemble the Early Holocene peak type (PC1 − PC2) and Middle Holocene peak type (PC1 + PC2), respectively.

The Early Holocene TWTA peak

The EH-peak type, consistent with the mean TWTA trend (Fig. 1B), is in phase with the changes in Earth’s orbital configuration of precession and obliquity (Fig. 3A). The relationship between precession/obliquity and IPWP’s thermocline change has been found before (1012, 27, 28) and was explained by some regional oceanographic processes. Such a common Early Holocene peak, however, implies a common driving mechanism over the entire IPWP thermocline. The western equatorial Pacific thermocline water originates from the basin-wide shallow overturning circulation of the Pacific Ocean (1, 29), which is fed by the subduction of relatively salty, warm surface waters in the subtropical North and South Pacific (1), and is primarily regulated by the surface wind stress curls determined by the meridional SST gradients (30). Over the past 25 ka, the gradient between the southwestern Pacific SSTA (from 45.5°S, 174.9°E) (31, 32) and the IPWP mean SSTA (Fig. 3B) resembles the EH-peak type TWTA, possibly reflecting relatively warmer midlatitude and thereby enhanced warm water transport of shallow overturning circulation in the Early Holocene. In addition, foraminifera δ13C records of the equatorial Pacific also suggest an Early Holocene intensification in the advection of southern-sourced subsurface waters (33, 34). An Early Holocene peak also appears in the southern Pacific and Antarctic temperature records because of the June (austral winter) insolation maximum at precession minimum (fig. S3). In modern observations (35), the northern Pacific shallow overturning circulation also contributes to the IPWP thermocline water, but the northwestern Pacific SSTA record shows no direct linkage to the Early Holocene peak in IPWP thermocline (fig. S3). Therefore, we can only propose that the overall trend of IPWP’s thermocline evolution over the LGM-Holocene may be dominated by the southern Pacific shallow overturning circulation, under the control of changing meridional insolation gradient induced by orbital forcing [i.e., precession and obliquity (36)].

Fig. 3 Precession-forced Early Holocene TWTA peak.

(A) Mean TWTA of the Early Holocene peak type (brown), precession (red dashed line), and obliquity [gray dashed line; (47)]. (B) Meridional SSTA gradient between southwest Pacific [SWP; site MD97-2120, from 45.5°S, 174.9°E; (3132)] and the IPWP. Holocene, last deglaciation, and LGM are separated by dashed vertical lines. CESM-simulated responses of Pacific subsurface temperature to June insolation in the Early Holocene are shown in (C to F): horizontal temperature anomaly distributions of upper-thermocline [at 120 m in (C)] and deeper thermocline [160- to 180-m water depth average in (D)] and meridional upper-water temperature anomaly profiles in the open Pacific [140°E to 140°W in (E)] and the Maritime Continent waters [100°E to 140°E in (F)]. Temperature anomalies in (C) to (F) are shown as regression coefficients against the standardized time series of the June insolation at precessional band in experiment CESM_GHG. White shadings mask insignificant results below 95% confidence level (t test). EQ, equator.

Our hypothesized mechanism for the Pacific subsurface temperature change is verified by a transient simulation of the Community Earth System Model [CESM1.0.4 (37)], forced by the orbital insolation and greenhouse gas (GHG) changes of the past 300,000 years (detailed in Materials and Methods)(38). The responses of the Pacific upper-ocean thermal structure point out the key role played by precession in forming the Early Holocene peak of TWT in the IPWP. At precession minimum during the Early Holocene, an intrusion of southern Pacific warm waters resulted in a drastic thermocline warming in the 30- to 200-m water depth of the open-ocean equatorial Pacific (140°E to 140°W; Fig. 3, C to E) and in the relatively deeper (below 120-m water depth) Maritime Continent waters (100°E to 140°E; Fig. 3F) (39). Noteworthy, the precession minimum also induces considerable cooling anomalies in the shallower Maritime Continent waters above 120-m depth in our simulation (Fig. 3F), in contrast to the paleo-proxy–based Early Holocene TWTA warming off the Philippines and in the Timor Sea. Therefore, the Early Holocene TWT warming in the IPWP may be a result of the precession-forced warming of deeper thermocline waters (Fig. 3D), which cannot be explained by the obliquity maximum with a cooling effect instead (fig. S4). Of course, the Early Holocene warming may also be induced by the influence of increased atmosphere pCO2, which results in a universal warming at all latitudes (fig. S4). In addition, the deglacial sea level rise can also deepen the thermocline and result in thermocline warming in the Maritime Continent waters (40). Thus, the Early Holocene warming of the IPWP subsurface water could have been caused by the combined effect of precession minimum, atmosphere pCO2 maximum, and sea level high stand.

The Middle Holocene TWTA peak

For the near-equator sites within the Maritime Continent waters, the most substantial thermocline warming occurs in the Middle Holocene. These records share two main features: (i) a cooling spell in 18 to 15 ka and (ii) a warming peak at ~7 ka (Fig. 4A). These sites are apparently less directly influenced by the southern Pacific sourced signal of the Early Holocene TWTA peak. The MH-peak type TWTA varies in phase with the equatorial September insolation change (Fig. 4A) that dominates ENSO-related activities in the tropical Pacific (41). For example, the overall pattern of the MH-peak type TWTA is consistent with the W-E zonal temperature difference in both the sea surface (42) and subsurface across the equatorial Pacific, which shows maxima in the Middle Holocene and minima in the early stage of the deglaciation (Fig. 4B). Likewise, the strength of the ascending limb of the Walker circulation, as indicated by the Borneo stalagmite δ18O records, shows a minimum (more positive δ18O) around 17 to 16 ka and a maximum (more negative δ18O) around 7 ka (Fig. 4C), suggesting an enhanced atmospheric convection over Borneo in the Middle Holocene.

Fig. 4 Time series and simulated temperature and rainfall anomalies in the IPWP since the LGM.

(A) Mean TWTA of Middle Holocene peak type (green) and the September 21st insolation at the equator (dashed red line) and obliquity [gray dotted line; (47)]. (B) Zonal temperature gradients at subsurface (Δsub-TA, in gray) and at sea surface [ΔSSTA, in violet; (42)] shown as the difference between the western equatorial Pacific and eastern equatorial Pacific. (C) δ18O records of northern Borneo stalagmites (15, 16) (dark and light green) and simulated annual mean rainfall (millimeters per day) over Borneo (dark gray, this study). PDB, Pee Dee belemnite. (D) Mean anomaly of seawater δ18O of IPWP (Δδ18Osw, dark gray, shading shows the 1σ error of the records). Shadings, vertical bars, and dashed lines are as in Fig. 3. Simulated response of the Pacific subsurface temperature and atmospheric variables to September insolation maximum are shown in (E to H): (E) Annual mean TWTA at 120-m water depth. (F) Depth profile of the annual mean temperature anomaly across the Pacific between 5°S and 5°N. (G) Late-autumn (October to December) anomalies of mean rainfall (colors, in millimeters per day) and horizontal winds at 850 hPa (arrows, in meters per second, reference arrow on top right). (H) Late-autumn mean Walker circulation anomalies between 5°S and 5°N across the Pacific, as indicated by anomalies in wind (arrows, in meters per second, reference arrow on top right) and in vertical velocity (colors, in pascal per second). Positive values in red indicate upward motion, and negative values in blue indicate downward motion. These anomalies in (E) to (H) are shown as regression coefficients against the standardized time series of the September insolation at precessional band in experiment CESM_GHG. SMOW, standard mean ocean water.

The hydroclimate changes revealed by Borneo stalagmite are supported by our CESM simulation of annual mean rainfall time series over Borneo forced solely by orbital insolation change (Fig. 4C). In addition to the Borneo stalagmite records, the surface seawater δ18O stack (δ18Osw) of the IPWP shows positive excursions in the last deglacial and a negative peak in the Middle Holocene (Fig. 4D), indicating a strengthened convective precipitation over evaporation in the Middle Holocene. Thus, we argue that the Middle Holocene thermocline warming of the near-equator IPWP is dynamically linked to the equatorial Pacific ENSO-like changes (e.g., enhanced Walker circulation and strengthened W-E zonal thermal contrast in the Middle Holocene). Our model simulations verify that September insolation maximum forces a warming in the IPWP thermocline (Fig. 4E) and a stronger zonal thermal difference across the equatorial Pacific (Fig. 4F). The atmospheric response to an increased zonal thermal gradient leads to increased rainfall over western equatorial Pacific (Fig. 4G) and a stronger Walker circulation (Fig. 4H).


The long-term evolution of the tropical Pacific mean state—including the IPWP’s thermocline temperature, the W-E temperature gradients, and the western equatorial Pacific hydroclimate—has the potential to shape shorter-term climate oscillations, i.e., interannual ENSO activity, as suggested by simple model simulations (41, 43). An Early to Middle Holocene depression of ENSO activity associated with strengthening of the Walker circulation relative to modern is evidenced by several proxy records and model simulations (13, 16, 44). Our findings suggest that the evolution of the equatorial Pacific climate in response to precession forcing could be understood in analogy to the modern seasonal development of the equatorial Pacific air-sea coupled system (2). That is, in the Early Holocene under the precession minimum, the thermocline of the open-ocean IPWP warmed widely, thereby likely suppressing ENSO activity. During the Middle Holocene, maximal September insolation may have caused an overall thermocline warming, increased precipitation, and decreased sea surface salinity in the IPWP and strengthening of the Walker circulation (Fig. 4D). A maximum in W-E upper-ocean thermal contrast (Fig. 4B) ultimately led to an extreme reduction of ENSO activity in the Middle Holocene.

The response of ENSO activity to future global warming and consequences to Earth’s climate evolution are not well constrained by either modern observations or model simulations (45), thus necessitating additional observations from paleoclimate records. Our study shows that warming of the western equatorial Pacific thermocline coupled with increased W-E thermal gradient and strengthened Walker circulation may have ultimately led to the reduction in ENSO activity during the Early and Middle Holocene, when climate was arguably slightly warmer than at present (39, 46). This inference raises the possibility that enhanced anthropogenic heat sequestration in the western equatorial Pacific subsurface waters, through the shallow overturning cell and equatorial Pacific air-sea coupled system, may further augment heat uptake in the eastern equatorial Pacific cold tongue due to reduced ENSO activity. In the near future, these may subsequently lead to an intermittent slowdown of surface warming, likely for short periods, in a pattern akin to the global warming hiatus between 2000 and 2014 (3, 5, 9).


We analyzed the Mg/Ca and δ18O of G. ruber (250 to 350 μm) and P. obliquiloculata (350 to 440 μm) at the State Key Laboratory of Marine Geology, Tongji University, Shanghai, China. Mg/Ca measurements were conducted on an inductively coupled plasma mass spectrometry (Thermo VG-X7) with a measurement reproductivity of 2.2% for G. ruber (n = 311) and 4.8% for P. obliquiloculata (n = 302), estimated by replicate samples (n, total replicates of the three cores of this study; for details, see table S2). The shell δ18O of the two species was measured with a Finnigan-MAT253 mass spectrometer. Conversion to the international Pee Dee belemnite scale was performed using NBS19 standard, and the long-term variability of δ18O is better than 0.07 per mil. Details of pretreatments and procedures are described elsewhere (20).

The age models for the IPWP cores were all established mainly by linear relationships of radiocarbon dates, first corrected for the 14C reservoir ages by the Marine Reservoir Correction and then calibrated to calendar age using CALIB7.1 software ( (tables S1 and S3). The time series of proxies (SST, TWT, and δ18O) were then averaged at 150-year nonoverlapping bins using the staircase integration resampling method. The temperature gradients of IPWP relative to the eastern equatorial Pacific or extratropical seas are calculated by the differences between the respective temperature anomaly records and on temporal steps determined by the average temporal resolution of the corresponding records (150 years for W-E subsurface temperature gradient, 500 years for W-E SST gradient, and 600 years for South-Equatorial Pacific SST gradient).

Here, we use the CESM1.0.4 with T31_gx3v7 resolution [3.75° × 3.75° for atmosphere and nominal 3° resolution for ocean (37)] to simulate the response of Pacific upper-ocean thermal structure to the forcing of orbital configuration (obliquity and precession) and change in atmospheric GHG content (38). As a spin-up, the CESM was first run for 200 model years under orbital parameters and GHG of 300 ka and other boundary conditions in 1950 AD. Then, the model was integrated for 3000 model years with the transient orbital insolation forcing and GHG changes of the past 300,000 years, in which orbital parameters and GHG were advanced by 100 years at the end of each model year (experiment CESM_GHG). A similar transient accelerated experiment (CESM_ORB) was only forced by orbital insolation changes since 300 ka (38). The outputs in the last 3000 model years of these two experiments were both analyzed, and they exhibit similar responses to orbital insolation forcing. Thus, only the results from experiment CESM_GHG are shown. At first, ocean temperature, salinity, atmospheric circulation, and precipitation are extracted from original outputs along multiple profiles [i.e., the latitude-longitude profile at 120-m water depth, the longitude-vertical profile along the equator, and the latitude-vertical profile zonally averaged over the open Pacific (140°E to 140°W) or the western Pacific (100°E to 140°E)]. Then, these oceanic and atmospheric variables were linearly regressed onto the normalized time series of specific orbital forcing [i.e., obliquity parameter changes, GHG changes, and the June or September insolation changes defined by the solstice or equinox precessional mode, respectively (38)]. Associated regression coefficients represent the Pacific air-sea coupled responses between the maxima and minima of each orbital forcing. Statistical significance is assessed by the 95% confidence level of t test.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We appreciate the sea-board assistance of cruises MD122, MD181, SO185, and KX08-973 and the laboratory assistance by P. Qiao, W. Fan, and Z. Chu. Funding: This study is supported by the National Natural Science Foundation of China (grants 91958208, 41976047, 91428310, 41630965, and 91858106), the National Key Research and Development Program of China (grant 2018YFE0202401), and the Ministry of Natural Resources of China (grant GASI-GEOGE-04). Financial supports by the BMBF [grants 03G0806A (CARIMA) and 03G0864F (CAHOL) to M.M. and grants 03G0185A (VITAL) and 03G0217A (MAJA) to W.K.] are acknowledged. This study is implemented under the France-China Framework of LIA-MONOCL. Author contributions: H.D. wrote the draft of the paper, and all authors contributed to subsequent revisions. Z.J. designed the framework, and H.D. proposed the preliminary idea of this research. Z.J. took lead in the experiments, assisted by H.D. and L.Y. Y.W. carried out the numerical simulation. M.M. and Y.R. helped the interpretations. F.B. and W.K. provided part of the sediment materials and helped the accelerator mass spectrometry 14C measurements. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. All the data can also be obtained from NOAA ( and PANGAEA ( or from Z.J. (email: jian{at}

Stay Connected to Science Advances

Navigate This Article