Research ArticleCLIMATOLOGY

Glacial changes in tropical climate amplified by the Indian Ocean

See allHide authors and affiliations

Science Advances  12 Dec 2018:
Vol. 4, no. 12, eaat9658
DOI: 10.1126/sciadv.aat9658


The mechanisms driving glacial-interglacial changes in the climate of the Indo-Pacific warm pool are poorly understood. Here, we address this question by combining paleoclimate proxies with model simulations of the Last Glacial Maximum climate. We find evidence of two mechanisms explaining key patterns of ocean cooling and rainfall change interpreted from proxy data. Exposure of the Sahul shelf excites a positive ocean-atmosphere feedback involving a stronger surface temperature gradient along the equatorial Indian Ocean and a weaker Walker circulation—a response explaining the drier/wetter dipole across the basin. Northern Hemisphere cooling by ice sheet albedo drives a monsoonal retreat across Africa and the Arabian Peninsula—a response that triggers a weakening of the Indian monsoon via cooling of the Arabian Sea and associated reductions in moisture supply. These results demonstrate the importance of air-sea interactions in the Indian Ocean, amplifying externally forced climate changes over a large part of the tropics.


The Indo-Pacific warm pool (IPWP)—the region of warm waters and strong rainfall straddling the Indian Ocean (IO) and Pacific Ocean—experienced large climate changes during the Pleistocene glacial periods (13); yet the drivers and mechanisms behind these changes remain unclear. This is an important question because a better understanding of these past intervals could help evaluate models and theories used to predict future changes in this key component of the global hydrological cycle. These glacial changes are expected to be driven by changes in the Pacific Ocean (1, 4), a basin that dominates IPWP climate variability on year-to-year time scales associated with the El Niño/Southern Oscillation phenomenon. The patterns in the paleodata, however, are spatially heterogeneous and could reflect multiple processes, such as the influence of seasonal insolation forcing (5), the exposure of the Sunda and Sahul landmasses due to lowered sea level (2, 610), and/or ice sheet topography (3). Models simulate a large diversity of responses during the most recent glacial period, many of which do not match inferred signals from paleoclimate proxies (10), limiting our ability to attribute the drivers of IPWP climate on glacial-interglacial time scales.

Here, we explore the drivers and mechanisms of IPWP climate change using a suite of model simulations and two multiproxy reconstructions of rainfall and sea surface temperature (SST) of the Last Glacial Maximum (LGM; ca. 21,000 years ago). We focus our study on the LGM because it is firmly dated and well sampled, allowing us to infer robust patterns of hydroclimate change and ocean surface cooling. To identify key atmospheric and oceanic processes, we use a synthesis of rainfall-sensitive proxies and a reconstruction of LGM SST anomalies derived from sedimentary alkenone Embedded Image, TEX86, Mg/Ca, and δ18O data. Each reconstruction provides an independent target over land and ocean, allowing us to make robust inferences on coupled ocean-atmosphere dynamics. To diagnose the mechanisms driving glacial changes in IPWP climate, we compare these reconstructions to simulations performed with the Community Earth System Model version 1.2 (CESM1). Our modeling approach includes a simulation of the LGM climate and a series of “single-forcing” simulations with different combinations of greenhouse gas (GHG) concentrations, ice sheet topography, ice cover, and sea level (table S3). We also perform a hierarchy of simulations with different degrees of air-sea coupling to isolate the importance of atmospheric versus ocean processes. See the Supplementary Materials for further details on the proxy syntheses, the model simulations, and the coupling hierarchy.

Different methodologies are used to quantify model-proxy agreement. The hydroclimate reconstruction is categorical, and therefore, we estimate the model-proxy agreement using the weighted Cohen’s κ statistic, following previous work (10). Cohen’s κ values are computed as a function of thresholds used to categorize the model changes as a drier, unchanged, or wetter LGM (relative to preindustrial conditions). Our main metric of proxy-model agreement is the maximum κ value with 2σ significance (κmax) resulting from this sensitivity analysis. The ocean temperature reconstruction is quantitative, and therefore, we simply calculate the pattern correlation between proxy and simulated values at the core sites. Further details on the proxy syntheses, implementation of the glacial boundary conditions, the single-forcing simulations, the coupling hierarchy, and the proxy-model agreement quantification may be found in the Supplementary Materials.


Our simulation using all LGM boundary conditions combined shows drying of the central IPWP, northern Australia, India, and southeastern Asia, and wetter conditions over the western IO (Fig. 1A). These patterns closely resemble the rainfall-sensitive paleorecords, which also show a dipole of a drier IPWP center and a wetter coastal equatorial Africa, in addition to drying over Asia and inland equatorial Africa (Fig. 1A, symbols). Our model also simulates an altered zonal SST gradient along the IO, with enhanced cooling in the east and subdued cooling in the west (Fig. 1B), a pattern evident in the SST proxy data (Fig. 1B, circles). Likewise, both the model simulation and the proxy data show enhanced cooling over the Arabian Sea, where several records show cooling in excess of 4 K. The simulated patterns show statistically significant (P < 0.05) agreement with the proxy data, with a maximum Cohen’s κ = 0.28 for hydroclimate and a pattern correlation r = 0.34 for SST, suggesting that CESM1 provides a realistic representation of the glacial climate state of the Indo-Pacific. Hereafter, we unravel the mechanisms driving the distinct spatial patterns seen in the proxy data, assessing proxy-model agreement with select single-forcing simulations.

Fig. 1 Simulated and reconstructed climate of the IPWP during the LGM.

Change in (A) rainfall and (B) SST (shading) and surface wind stress (vectors) simulated by CESM1 under LGM boundary conditions relative to preindustrial conditions. The color schemes for SST and rainfall change are not linear. Proxy-inferred changes in hydroclimate (symbols) indicate drier (brown), wetter (green), or unchanged (white) conditions at the LGM relative to the Late Holocene (top). Triangles represent merged sites within a radius of 100 km. Dark gray circles or triangles indicate disagreement between different proxies at a given site or merged sites. SST changes (bottom) are derived from geochemical proxies (see the Supplementary Materials) and are plotted as circles on the same color scale as the simulated changes. Proxy-model agreement for hydroclimate is measured by a Cohen’s κ value of 0.28. Proxy-model agreement for SST is measured by a pattern correlation of 0.34. Coastlines correspond to the LGM 120-m drop in sea level.

Response to lower sea level

To isolate the impact of lowered sea level on IPWP glacial climate, we focus on the SLMC response, which captures the effect of lower LGM sea level over the Maritime Continent (exposure of the Sunda and Sahul shelves, altered Indonesian throughflow passages, and tidal mixing) but holds GHGs and ice sheet extent at preindustrial conditions (see the Supplementary Materials for details). These sea level effects explain many patterns of hydroclimate change at the LGM inferred from our proxy network, including the large-scale drying over the Maritime Continent and northern Australia, as well as the wetter conditions over equatorial East Africa (Fig. 2A). These rainfall changes are accompanied by colder SSTs over the eastern IO and warmer SSTs over the western IO (Fig. 2B, shading), also evident in the proxy data after removing the mean IO cooling (Fig. 2B, circles). The simulated dipole of drier/wetter conditions across the equatorial IO reflects a greatly weakened Walker circulation in the IO, with increased ascending motion and convection in the western IO and reduced ascending motion and convection over the eastern IO and the Maritime Continent (Fig. 3A).

Fig. 2 Glacial drivers of IPWP hydroclimate change.

Annual mean climate response to (A to C) effect of lowered sea level on the geography of the Maritime Continent, (D to F) glacial concentrations of GHGs, and (G to I) ice sheet albedo. Maps show, from left to right, annual mean changes in rainfall, SST and surface wind stress (vectors), and upwelling (shading) and thermocline depth (contours). Overlaid circles show proxy-inferred changes in hydroclimate (left) and SST (center). The simulated and proxy SST change averaged over the domain was removed to emphasize the changes in horizontal gradients. The value of the simulated SST cooling averaged over the tropics (30°S to 30°N) is shown in blue. Upwelling is defined as the ocean vertical velocity at a 50-m depth. Thermocline depth is computed as the depth of the maximum vertical temperature gradient. Thermocline depth change contours are at 5-m intervals. Note the different proportions of the right panels focusing on the equatorial IO. The color schemes for relative SST and rainfall change are not linear.

Fig. 3 Walker circulation response to LGM boundary conditions.

Changes in vertical velocity (ω) over the equatorial Indo-Pacific simulated by CESM1 in response to the following changes in boundary conditions at the LGM: (A) sea level effects over the Maritime Continent, (B) lowered GHGs, and (C) full LGM effects combined. Contours are annual mean ω simulated in the preindustrial control experiment. Vertical velocity values are averaged over the 10°S to 10°N latitude band. Contour intervals are 5 hPa day−1. Note that the color scale is not linear.

These changes represent a marked departure from the modern Walker circulation, albeit not a complete reversal of the large-scale overturning over the IO. The changes in SST are a critical part of this response, as the warmer western IO becomes more favorable for ascending motion—and deep convection—while the colder eastern IO favors subsiding motion and drying. This SST dipole drives anomalous easterly winds along the equator (Fig. 2B, vectors), which drive the ocean into a state that also represents a major departure from modern conditions. The equatorial thermocline, which is uniformly deep across the modern equatorial IO, responds to the easterly winds, establishing a large-scale east-west tilt with shoaling of more than 15 m over the eastern side of the basin (Fig. 2C, contours). In addition, the easterly winds drive equatorial upwelling (Fig. 2C, shading), also in contrast to the present-day climate, in which upwelling is restricted off the coast of Sumatra. These two processes act to cool the eastern IO, reinforcing the SST gradient and the easterly winds, consistent with a positive feedback loop akin to the Bjerknes feedback (11). This feedback is excited by an atmospheric response to exposure of the Sahul shelf, with drying and anomalous subsidence over the shelf and anomalous easterly winds over the eastern IO. These wind changes shoal the thermocline in the east, initiating the cooling (12). In addition, anomalous westward currents advect the climatological temperature gradient (warmer in the east, less warm in the west), effectively warming the western IO and establishing the characteristic SST dipole seen in the simulations and proxy data.

There is significant proxy-model agreement between SLMC and both hydroclimate and SSTs (κmax = 0.26 and r = 0.38), indicating that lowered sea level during the LGM relative to modern conditions is responsible for key features of the proxy patterns (Fig. 4). Among the sea level boundary conditions, exposure of both the Sahul and Sunda shelves is the main driver of these responses, as shown by the superior agreement with both the hydroclimate (κmax = 0.26; Fig. 4A) and the SST (r = 0.27, P < 0.1; Fig. 4B) proxies. Exposure of the Sunda shelf explains drying over Borneo, Sumatra, and Thailand, while the Sahul shelf explains the drying over New Guinea and northern Australia and wetter equatorial East Africa (fig. S2). Statistically significant κ values occur for a limited range of thresholds used to categorize the model changes into drier, unchanged, or wetter LGM (fig. S3), suggesting that each shelf in isolation cannot explain data. Proxy-inferred SST changes are best explained by exposure of both shelves combined (Fig. 4B), with Sahul dominating the agreement, a result that demonstrates that this shelf is most effective at exciting the Bjerknes feedback (12). Further inclusion of GHG and orbital forcing does not change the agreement for hydroclimate (Fig. 4A) and degrades the agreement with SSTs (Fig. 4B). Together, these results demonstrate that shelf exposure and Bjerknes amplification best explain the proxy-inferred hydrological changes over the IO and Maritime Continent.

Fig. 4 Model-proxy agreement.

Maximum Cohen’s κ (left) and pattern correlation coefficient (right) for each climate response simulated by CESM1. Cohen’s κ quantifies agreement between simulated patterns of rainfall change and a multiproxy synthesis of hydroclimate change at the LGM. The maximum κ value obtained from the threshold sensitivity analysis is reported. The correlation coefficient (r) quantifies agreement between simulated patterns of SST change and our reconstruction of LGM temperature changes derived from Mg/Ca and Embedded Image records. A κ or r value of one indicates perfect model-proxy agreement. A value of zero indicates lack of agreement. Solid circles indicate statistically significant (P < 0.05) κ and r values. Error bars indicate the 95% confidence intervals for the κ values. Uncertainty due to calibration of SST proxies is reported as the minimum-maximum range of the r values computed using each reconstruction in our ensemble (blue error bars). The 95% confidence interval of the r values computed using the ensemble mean reconstruction is also shown (gray error bars). Further details on the computation of these proxy-model agreement estimators, as well as the definition of the climate responses, are provided in the Supplementary Materials.

Response to GHGs

We followed the same approach to determine whether lower GHGs could explain any part of the patterns seen in the proxy data. Here, we use our single-forcing simulation that reduces GHG concentrations to LGM levels but holds ice extent and sea level effects at preindustrial values. The associated response shows wetter conditions over the Maritime Continent and drying over the western IO and equatorial East Africa (Fig. 2D, shading). GHG forcing also drives drying in inland East Africa, explaining the proxy-inferred patterns over that region. The simulated wetter conditions over the Maritime Continent are not supported by the proxies (Fig. 2D, circles), resulting in low model-proxy agreement (Fig. 4). Reduced GHGs drive a tropical mean cooling of −1.9 K and produce an SST pattern that is largely the opposite of the sea level response, with a warmer eastern IO (Fig. 2E, shading) and anomalous westerlies along the IO (Fig. 2E, vectors). Accordingly, the pattern correlation with the SST proxy data is negative but not significant at the 1σ level (r = −0.13, P > 0.33; Fig. 4B). Together, these changes are consistent with stronger ascending motion over the IPWP center (Fig. 3B), anomalous westerly winds along the equatorial IO (Fig. 2E, vectors), suppression of equatorial upwelling (Fig. 2F, shading), and a deeper thermocline in the east (Fig. 2F, contours). These responses also indicate Bjerknes amplification, as in the sea level case, but with circulation and SST gradient changes in the opposite direction, hence the lack of proxy-model agreement (Fig. 4).

Response to ice sheets

We isolate the response to ice sheets (topography and albedo) using single-forcing simulations that change only these aspects. Additional simulations allowed us to isolate the responses to ice sheet albedo (ISalbedo) and topography (IStopo). The ISalbedo response shows that the higher albedo of the LGM ice sheets causes dry conditions over the northern IPWP, India, and southeastern Asia, a key feature of the proxy-inferred patterns that cannot be explained by shelf exposure (Fig. 2G). Ice sheet albedo drives SST changes with enhanced cooling in the Arabian Sea, Bay of Bengal, and South China Sea (Fig. 2H). These changes lead to a significant pattern correlation with the SST proxies (r = 0.34, P < 0.05; Fig. 4B), which show pronounced glacial cooling in these areas (Fig. 2H, circles). The changes in response to ice sheet topography are generally opposite and weaker over the IPWP domain, confirming the preeminence of ice sheet albedo as a driver of glacial changes in IPWP climate (see section S7.1 for full details on the responses to albedo and topography and their nonlinear interactions).

Glacial response of the Indian monsoon

To further explore the mechanisms by which the Northern Hemisphere (NH) ice sheets influence IPWP climate, we performed a hierarchy of simulations with varying degrees of ocean-atmosphere coupling. The hierarchy consists of uncoupled (atmosphere only), thermally coupled (atmosphere coupled through heat fluxes to a mixed layer ocean), and fully coupled simulations. All uncoupled simulations were run with SSTs from the preindustrial control (see section S7.1 for further details). All three types of simulations show a large-scale weakening of the Afro-Asian summer monsoon system with drying stretching from the Sahel to northwest India (Fig. 5, top, and also fig. S9). The fact that widespread continental drying occurs in the uncoupled case (Fig. 5A and also fig. S9C) indicates that SST changes do not play a role in linking ice sheet cooling with the large-scale monsoon weakening. We characterize this large-scale response as a monsoonal retreat, with the exception of south India, which shows increased rainfall in the uncoupled case (Fig. 5A). Inclusion of air-sea coupling causes drying over south India (Fig. 5, B and C), suggesting that SST changes play a role on this regional aspect of the response. The large-scale pattern of drying from Africa to central Asia can be attributed to enhanced “ventilation” of mid-latitude air as follows. Cold and dry air generated over the Laurentide ice sheet (LIS) is advected by the westerly winds around the NH, particularly over northern Africa and central Asia. This air is mixed equatorward, making the environment less conducive to deep convection. This limits the summertime advance of the monsoon rain belts seen as reductions in summertime rainfall (see section S7.1 for further details).

Fig. 5 Boreal summer climate response to ice sheet albedo as a function of ocean-atmosphere coupling.

Rainfall (A to C) and surface temperature (D to F) change during July-August-September (JAS) simulated by CESM1 in response to ice sheet albedo. Surface wind stress changes (vectors) (D to F) are June-July-August (JJA) to highlight their role initiating the cooling of the Arabian Sea. Climate responses under varying degrees of ocean-atmosphere coupling are shown. Left to right: CESM1 with fixed SSTs, with a mixed-layer ocean model, and with a fully interactive ocean model.

The retreat of the monsoon rainbelts leads to reductions in cloud cover and relative humidity, causing massive changes in the energy balance of the land surface (fig. S10). Increased radiation and decreased evaporative cooling cause pronounced warming over a swath of land stretching from the Sahel to northwest India (Fig. 5, bottom, and fig. S9, top). This warming results in a large land-sea temperature contrast between the Arabian Peninsula and the Arabian Sea that drives anomalous southwesterly winds along the shore (Fig. 5, bottom, vectors). These wind changes are present in the uncoupled case (Fig. 5D), showing that they are directly driven by the large-scale monsoonal retreat. Thermal air-sea coupling plays a key role in enhancing the land-sea temperature contrast, as the increased wind speed cools the ocean surface over the Arabian Sea (Fig. 5E). Interactive ocean dynamics do not change the pattern and amplitude of this response (Fig. 5F), demonstrating that evaporative cooling is the essential mechanism cooling the Arabian Sea. Furthermore, this response explains the drying over the core Indian monsoon region (Fig. 5, B and C), especially in southern India, confirming that cooling of the Arabian Sea is the key feedback contributing to the glacial drying over India. Note that this thermodynamic effect overwhelms the effect of stronger winds, which would otherwise increase rainfall, as seen in the uncoupled response (Fig. 5A).

Drivers and mechanisms of LGM climate

Our simulations and quantitative proxy-model agreement show that the full LGM response is dominated by the combined effect of sea level and ice sheet albedo. These two drivers explain most of the large-scale features seen in the proxies, including the wetter/drier dipole along the equatorial IO and the drier conditions in the central IPWP, northern Australia, India, and southeastern Asia. The combined response (ISalbedo + SLMC) also shows the highest pattern correlation (r = 0.50, P < 0.05) with the SST proxies (Fig. 4B), which is higher than the correlation with the response to sea level alone, attesting to the enhanced improvement by inclusion of the ice sheet albedo response. Over the deep tropics, the hydroclimate changes are dominated by changes in the Walker circulation, which shows a massive weakening of its ascending branch over the Maritime Continent accompanied by increased ascending motion over the western IO (Fig. 3C). The spatial patterns of the ISalbedo + SLMC response are nearly identical to that of the full LGM response; however, the full LGM conditions result in lower κ and pattern correlation values (Fig. 4) because the exclusion of the GHG forcing artificially enhances the agreement of the model with the proxies. While statistically significant (P < 0.05), our r and κ metrics of proxy-model agreement show rather low values, particularly for the full LGM response. This indicates that only a small fraction of the proxy pattern variance is explained by our simulations. However, we note multiple cases where nearby proxy sites show lack of agreement, for example, in the eastern IO, where the magnitude of the cooling estimated by Mg/Ca and Embedded Image proxies does not necessarily agree with each other. Therefore, the low r values could also be indicative of proxy uncertainties, in addition to model uncertainties. Two robust patterns in the hydrological proxies, the muted changes in Sumatra and New Guinea, are not captured by our model, suggesting obvious model deficiencies there. Nevertheless, our conclusion that shelf exposure and ice sheet albedo are the main drivers of the LGM climate is not sensitive to these issues, stressing the importance of multiproxy approaches to studying past climate changes.


Our findings have important implications for understanding mechanisms driving tropical climate change. First, the LGM test case demonstrates that IO dynamics amplify external perturbations through the Bjerknes feedback, resulting in marked large-scale hydroclimate changes extending over one-third of the global tropics, from eastern Africa, across the entire Maritime Continent, to the western Pacific. This positive feedback loop was first proposed by Bjerknes (11) to explain the growth of El Niño events. A similar mechanism has been hypothesized to control the response of the Pacific Ocean to external forcings (4, 13). However, proxies, observations, and models are inconclusive regarding Bjerknes amplification of past or future climate changes in the Pacific (1417). In contrast, here we present robust evidence that coupled ocean-atmosphere interactions can amplify externally forced climate changes in the IO—an ocean that is generally considered to play a passive role in tropical climate change. While there is no future analog for shelf exposure, the mechanisms uncovered here are directly relevant for projecting changes in the IPWP in response to greenhouse warming. Our simulations and proxy reconstructions show that small wind perturbations driven by shelf exposure are amplified by changes in the ocean through changes in the east-west SST gradient. The dynamics of this externally forced response are potentially unique to the IO and indicate that this ocean is capable of amplifying small asymmetries introduced by any external forcing, for instance, as hypothesized in response to greenhouse warming (18, 19). This result also revises the long-standing belief that IPWP climate is mainly forced by the tropical Pacific (4, 13, 2023) and supports the emerging view that the IO can behave autonomously on multidecadal and longer time scales (24), with important implications for regional hydroclimate prediction.

Second, our results reveal important new details on monsoon dynamics. Previous studies have argued that monsoon weakening during glacial periods is simply the terrestrial manifestation of global shifts in the Intertropical Convergence Zone (ITCZ) in response to NH cooling (2527). Climate models consistently show that the ITCZ shifts south when the NH cools, as a compensation for the alteration in energy balance (2831). Our simulations show that albedo cooling of the NH does shift the ITCZ in the Atlantic and eastern Pacific (fig. S9D), but this response requires thermal (or full) air-sea coupling. The monsoonal retreat, in contrast, occurs in the absence of air-sea coupling and without an ITCZ shift (fig. S9C), demonstrating that these responses obey distinct mechanisms. We posit that albedo cooling is also communicated to the tropics by land-atmosphere processes operating over Africa and Asia. Cold, dry air generated over the LIS is advected around the NH and mixed equatorward, inhibiting convection during the seasonal advance of the NH monsoons. This part of the response is consistent with the ventilation mechanism, proposed as a limit for the poleward extent of the monsoons (32).

The ventilation mechanism has been proposed to explain weakening of the African and North American monsoons in response to hemispheric cooling (33, 34), two monsoon systems that are not topographically insulated from mid-latitude air. The Indian monsoon, in contrast, is insulated by the Himalayas (35), explaining the more muted drying over India in the absence of air-sea coupling (Fig. 5A and fig. S9C). Enhanced glacial cooling of the Arabian Sea is thus the critical link between albedo cooling of the NH and the rainfall reductions over India. This cooling is tied to glacial ventilation, which causes warming of the Arabian Peninsula, establishing an anomalous thermal gradient relative to the Arabian Sea. This gradient drives stronger southwesterly winds increasing the evaporative cooling over the Arabian Sea. The colder SSTs reduce moisture transport into the core monsoon region—a thermodynamic effect that dominates over the dynamic effect of increased southwesterly flow. This mechanism allows glacial ventilation to “bypass” the insulating effect of the Himalayas and cause widespread drying over India. This mechanism could also explain drying over India during other intervals characterized by NH cooling. Previous studies have found a robust link between Arabian Sea cooling and India drying (36, 37), but the ultimate link to high-latitude cooling remains unclear. Monsoon ventilation and its effect on Arabian Sea surface temperatures could provide that link. Proxy evidence for monsoon ventilation could help establish this link; however, we anticipate that the evidence for land surface warming would be challenging to detect due to its seasonality and the cooling effect of other glacial forcings.

Our study firmly establishes that ocean-atmosphere interactions play a central role in shaping externally forced responses in IPWP hydroclimate. We also show the importance of shelf exposure and ice sheet albedo relative to other glacial boundary conditions. Ice sheets should therefore be considered alongside orbital (precessional) forcing as a major driver of Late Quaternary climate change in the tropics. However, as we focused on a time period with a similar orbital configuration to today (the LGM), we cannot assess the relative importance of orbital versus glacial (ice sheet–related) forcings during the Pleistocene. Our results demonstrate that patterns of cooling (and arguably patterns of warming) play a first-order role in the response of tropical rainfall to external forcings. Furthermore, our finding of an active climatological Bjerknes feedback in the IO deserves further attention in the context of climate prediction, since this feedback could be activated by anthropogenic forcings. Likewise, the emerging importance of ventilation as a control knob on monsoon intensity is directly relevant to how we understand and predict future rainfall in seasonally arid lands. The LGM still therefore offers us important lessons about Earth’s climate system, despite the fact that we are heading toward a warmer, and not cooler, world.


Model-proxy agreement

The role of each glacial forcing and associated mechanisms was assessed by quantifying the agreement between the simulated changes with two multiproxy reconstructions. We assessed agreement between simulated rainfall responses (Table 1) against the updated version of our synthesis of LGM hydroclimate described in section S1. As this reconstruction is categorical, we estimated the model-proxy agreement using the weighted Cohen’s κ statistic, following previous work (10). Cohen’s κ ranged from κ = 1, if a simulated pattern is in complete agreement with the proxies, to κ = 0, if the agreement could be expected entirely by chance. The weighted Cohen’s κ uses a weight matrix that penalizes models for a total miss (e.g., drier when it should be wetter) more than a near miss or partial agreement (e.g., drier when it should be no change). A weight of 0.5 is given to these cases. We explored the sensitivity of the κ values by varying the thresholds of rainfall, over which we placed the model output into the same categories of change assigned to the proxies. We consider that the model-proxy agreement is robust when 2σ statistically significant κ values are obtained for large range of dry-wet thresholds used to categorize the model simulated changes. The maximum Cohen’s κ values for each climate response are reported in Fig. 4A along with their 2σ confidence intervals.

Table 1 Climate responses.

Approach used to compute the climate response to different LGM boundary conditions. First column identifiers are used throughout the study to identify each climate response. The third column is the operation between climatologies used to compute each climate response. See tables S3 and S4 for a list of simulations. PI, preindustrial; GHG, greenhouse gases; ITF, Indonesia throughflow; SL, sea level; MC, Maritime Continent.

View this table:

We also assessed the agreement between simulated SST changes and a new compilation of LGM temperature anomalies across the IPWP derived from sedimentary alkenone, Mg/Ca, and δ18O data. As this reconstruction is quantitative, we simply calculated the pattern correlation between proxy and simulated values at the core sites. We applied different temperature calibrations to the alkenone and Mg/Ca data to explore the sensitivity of the model-proxy agreement to calibration choice. Most calibrations show the highest significant correlations with the simulated SST changes in response to either full sea level boundary conditions or sea level plus ice sheet albedo (fig. S5). The minimum-maximum range of r values computed from this distribution of SST reconstruction is reported as uncertainty of the pattern correlation. Further details are available in sections S2 and S6.

Climate model simulations

We used CESM1 (38), a model that simulates realistic IPWP climate and is sensitive to changes in the configuration of land masses over the Maritime Continent (12)—a response that is important for simulating glacial climates (10). CESM1 consists of coupled general circulation models of the atmosphere and ocean, as well as sea ice and land models. Other components of the Earth system, such as the carbon cycle and marine ecosystems, can also be simulated using CESM1; however, we kept them inactive because our focus here is on the climate response to glacial boundary conditions. The climate at the LGM was simulated by prescribing the following boundary conditions: (i) reduced GHG concentrations, (ii) insolation changes due to the orbital configuration at 21 ka, (iii) orography changes due to ice sheets and corresponding roughness length, (iv) surface albedo changes due to ice sheets, and (v) changes in the land-sea distribution and altitude due to lower sea level. A simulation of preindustrial climate was used as control. A series of simulations forced with individual glacial boundary conditions, also known as “single forcing” runs, were used to isolate the climate responses to different glacial drivers. The climate responses and single forcing simulations used to compute them are listed in Table 1. Full details on the implementation of the LGM boundary conditions and the single forcing simulations are available in the Supplementary Materials.

The ensemble was augmented by simulations in which the CESM1 atmosphere and land models (CAM5 and CLM4) were coupled to ocean models of varying complexity. These simulations allowed us to explore the importance of ocean-atmosphere coupling in response to ice sheets. First, we replaced the fully dynamical ocean (POP2) with a model of the ocean mixed layer. In this model, the effect of vertical mixing, entrainment, and horizontal currents is prescribed as a seasonally varying heat source or Q-flux. Changes in SSTs computed by this “slab” ocean model can only be influenced by energy exchanges with the atmosphere, such as changes in evaporation or clouds. In the second case, the ocean model consists of climatological SSTs and sea ice extent from the preindustrial control. In this configuration, the air-sea heat fluxes are computed but cannot change the prescribed climatological SST and sea ice extent. Therefore, climate changes simulated in this configuration are due to “uncoupled” atmosphere or land changes. “Thermally coupled” and uncoupled responses were computed from each configuration as in the full-coupled cases by differencing the simulations as specified in Table 1.


Supplementary material for this article is available at

Section S1. Hydroclimate synthesis

Section S2. SST synthesis

Section S3. Climate model

Section S4. Simulations

Section S5. Initialization and climate equilibration

Section S6. Quantitative proxy-model agreement

Section S7. Climate response to ice sheet boundary conditions

Fig. S1. Climate response to exposure of the Sahul shelf as a function of land surface properties.

Fig. S2. Simulated patterns of hydroclimate change at the LGM.

Fig. S3. Proxy-model agreement for patterns of warm pool hydroclimate change at the LGM.

Fig. S4. Simulated patterns of IO cooling during the LGM.

Fig. S5. Impact of Embedded Image and Mg/Ca-δ18O calibrations on proxy-model agreement.

Fig. S6. Simulated rainfall response to LGM ice sheets.

Fig. S7. Simulated surface temperature response to LGM ice sheets.

Fig. S8. Simulated boreal summer response to LGM ice sheets.

Fig. S9. Boreal summer response to ice sheet albedo as a function of ocean-atmosphere coupling.

Fig. S10. Surface energy changes in response to ice sheet albedo as a function of ocean-atmosphere coupling.

Table S1. Updates to the synthesis of LGM hydroclimate.

Table S2. Compilation of SST proxies.

Table S3. Sea level simulations.

Table S4. LGM and single-forcing simulations.

References (39127)

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 thank T. Shanahan, C. Deser, and K. Thirumalai for comments on this work. We acknowledge members of NCAR’s Climate Modeling Section, CESM Software Engineering Group (CSEG), and Computation and Information Systems Laboratory (CISL) for their contributions to the development of CESM. Computing resources (ark:/85065/d7wd3xhc) were provided by the Climate Simulation Laboratory at NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation and other agencies. P.D.N. gratefully acknowledges CGD for supporting a long-term visit at NCAR. Funding: Funding for this work was provided by NSF (grants AGS-1204011 and OCN-1304910 to P.D.N.), the David and Lucile Packard Science and Engineering Fellowship to J.E.T., and the National Center for Atmospheric Research (NCAR) sponsored by the U.S. NSF (to B.O.-B., N.R., and E.B.). This study was supported by the Institute for Basic Science, Daejeon, South Korea, under the project “Earth System Dynamics” code IBS-R028-D1. Author contributions: P.D.N. designed, performed, and analyzed the climate model experiments. J.E.T. produced the new SST synthesis. J.E.T., A.T., and T.B. helped with the analysis of the climate model experiments. B.O.-B., N.R., and E.B. helped with the configuration of the model. P.D.N. and J.E.T. prepared the manuscript. All authors discussed the results and contributed to manuscript preparation. 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.

Stay Connected to Science Advances

Navigate This Article