Agreement between reconstructed and modeled boreal precipitation of the Last Interglacial

See allHide authors and affiliations

Science Advances  20 Nov 2019:
Vol. 5, no. 11, eaax7047
DOI: 10.1126/sciadv.aax7047


The last extended time period when climate may have been warmer than today was during the Last Interglacial (LIG; ca. 129 to 120 thousand years ago). However, a global view of LIG precipitation is lacking. Here, seven new LIG climate models are compared to the first global database of proxies for LIG precipitation. In this way, models are assessed in their ability to capture important hydroclimatic processes during a different climate. The models can reproduce the proxy-based positive precipitation anomalies from the preindustrial period over much of the boreal continents. Over the Southern Hemisphere, proxy-model agreement is partial. In models, LIG boreal monsoons have 42% wider area than in the preindustrial and produce 55% more precipitation and 50% more extreme precipitation. Austral monsoons are weaker. The mechanisms behind these changes are consistent with stronger summer radiative forcing over boreal high latitudes and with the associated higher temperatures during the LIG.


The Last Interglacial (LIG) is a primary target period for climate research because its climate holds one of the closest, although partial, analogies to possible climates of the coming decades and centuries (1). Compared to today, the LIG was slightly warmer. Atlantic Ocean circulation was broadly similar (2, 3), although reconstructions of its northern branch conflict (4). Sea levels were several meters higher, implying that polar ice sheets were smaller (5). Yet, LIG greenhouse gas concentration was comparable to the preindustrial (PI) period, with carbon dioxide around 275 parts per million (ppm) and methane around 700 parts per billion (ppb) (6). Much research has used geological evidence to estimate the temperature of the LIG, in an effort to understand the thermal response to different orbital forcing (Fig. 1 and fig. S1) and constrain the internal feedbacks of the climate system (7). Compilations of climate proxies estimate that, during the peak of the LIG [ca. 127 to 125 thousand years (ka) ago], climate may have been warmer than PI by 2.0° ± 0.1°C over the whole Earth surface (8), although this assumes geographically synchronous peak LIG warming, which might be unrealistic. More recently, proxies have shown a global anomaly of only 0.5° ± 0.3°C in the ocean surface (9), with substantial regional differences (10). In the mid and high northern latitudes, land masses were considerably warmer by about 2° to 5°C (8) [see (11) and references therein], a pattern that also emerges in projections of future global warming (12). However, certain thermal features of LIG limit its validity as analog of future climates (13): Higher temperatures were probably rather a seasonal signal (11), reflecting the seasonality of the orbitally induced insolation anomaly, and land temperatures of the tropics and the Southern Hemisphere are hardly constrained. Previous climate model simulations of the LIG capture these large-scale patterns but tend to yield temperatures lower than those derived from proxies (14). Aside from this, the LIG offers a unique case to study the response of the hydrological system to differently distributed radiative forcing and associated higher-surface temperature. This has a wide societal relevance, since climate models project that global precipitation may increase by 3% per degree of global warming (15) because of the higher concentration of atmospheric water vapor. Furthermore, extreme precipitation has already increased over the majority of land areas, scaling with global warming (16), and models consistently simulate more extreme future precipitation outside of the tropics (12).

Fig. 1 Zonal mean anomalies between LIG and PI climates for each climate model and the ensemble mean.

Anomalies of near-surface air temperature over land and ocean (top) and precipitation over land only (bottom) are shown for the annual mean and each season (defined after the angular calendar). Annual and seasonal insolation anomalies [between 127 ka and 1850 CE; (6)] are plotted as dashed lines.

On the basis of proxy evidence, it has been suggested that the LIG was wetter than the Holocene in the Arctic region (17) and the boreal midlatitudes (18, 19). In addition, the Asian (20) and North African (21) monsoons were probably stronger, a finding that was replicated in LIG experiments with individual global climate models (22, 23) and ensembles of models over the Arabian Peninsula (24) and Southeast Asia (25). To gain confidence in LIG simulations, precipitation from models needs to be systematically compared with the best available empirical evidence. Beyond the interest for paleoclimate research, such comparison between proxies and models in the context of the LIG offers a rare opportunity to test the capabilities of global climate models to correctly simulate hydroclimate in a different climate state. Namely, it can inform us on the confidence that we can place in future precipitation projections, which are carried out with the same family of models. Yet, because proxy evidence of LIG precipitation is much less abundant than for more recent climates (26), a large-scale comparison was never attempted.

This study aims to (i) provide the first empirical global view of LIG precipitation; (ii) use this to quantify the ability of state-of-the-art climate models to reproduce precipitation in a different climate; and, with these models, (iii) study the seasonal response of LIG precipitation and its extremes. We compiled the existing empirical evidence from studies covering decades of geological fieldwork into the first LIG precipitation database. Parallel to this, we simulated the LIG within the framework of lig127k experiment of the Paleoclimate Modelling Intercomparison Project phase 4 (PMIP4) (27), which is part of the Coupled Model Intercomparison Project phase 6 (CMIP6) (28), and of the piControl experiment of CMIP6. Boundary conditions for the LIG and the PI experiments differ in the seasonal and geographical distribution of solar insolation (Fig. 1 and fig. S1), resulting from different orbital parameters, and secondarily, in greenhouse gas concentration (see Methods). The configuration of polar ice sheets is prescribed and equal in the two simulations. The common forcing protocol adopted across all models enables treating simulations as a coherent ensemble (6).


Annual global precipitation from proxies and models

Our new database includes 138 published records, with near-global coverage (see Methods and Fig. 2). We selected proxy signals at ca. 127 ka to ensure comparability with the model simulations. From these records, we extracted semiquantitative anomalies between the LIG and the present/PI. Some of these records reflect precipitation amount [e.g., loess deposits; speleothem δ18O according to some interpretations (29)], while others reflect the net balance between precipitation and evapotranspiration (e.g., pollen and lake levels). As shown in our models, there is a tight linear relationship between anomalies in precipitation and precipitation minus evapotranspiration (fig. S2); therefore, in the remainder of this paper, we refer to “precipitation proxies” for simplicity. A clear picture emerges from the proxy database of a wetter LIG over most parts of the boreal continents (Fig. 2): over Central and Northern Africa, the Arabian Peninsula and the Levant (where also some drier LIG anomalies are present), south and northeastern Asia, southern and eastern Europe, the western midlatitudes and the northwestern part of North America, and the northernmost part of South America. Over the austral continents, about 60% of the records indicate positive anomalies (i.e., LIG wetter than PI), and 30% indicate negative anomalies (i.e., LIG drier than PI). The south and northeast of Australia are wetter during the LIG; austral South America is wetter at its western coast over regions corresponding to central Peru and northern Chile and mostly drier elsewhere; and southern Africa has a spatially varying signal with both wetter and drier climates.

Fig. 2 Annual precipitation anomaly between LIG and PI from models (ensemble average in contoured colors) and proxies (filled markers)

Blue colors indicate higher precipitation during the LIG from models or proxies and conversely for red colors. Proxy anomalies are on a semiquantitative scale (see Methods): dark blue (much wetter LIG), light blue (wetter), white (no noticeable anomaly), light red (drier), and dark red (much drier). Different markers represent different types of proxy records. Hatched areas indicate regions where one (single hatching) or at least two (double hatching) of seven models disagree with the sign of the ensemble average anomaly.

Global mean precipitation of our model ensemble (see Methods) is almost the same for the LIG and PI simulations and so is global mean near-surface air temperature (presented in Supplementary text S1) and atmospheric water vapor. Integrating over land only, precipitation is higher in the LIG by 2.7% (Fig. 3), which is in agreement with a similar increase in atmospheric water vapor. The highest precipitation anomaly between the LIG and PI simulations is in the Northern Hemisphere (NH) tropics, with zonal maximum around +500 mm/year around 10° to 15°N (Fig. 1). LIG precipitation is higher over Central and Northern Africa, the Arabian Peninsula, the Indian subcontinent, Northeast Asia, Alaska, and Central America (Fig. 2). The only boreal continental region that is consistently modeled drier in the LIG is approximately the Mississippi River basin. The Southern Hemisphere tropics are drier in the LIG, with a maximum negative zonal anomaly of ca. −300 mm/year around 20°S. All austral land receives less precipitation in the LIG, except the western part of Central Africa, parts of Borneo and New Guinea, and New Zealand. For all these regional anomalies, all seven models agree on the sign of change.

Fig. 3 Percentage anomalies between LIG and PI.

Left: Anomalies in area-weighted precipitation. Right: Anomalies in monsoon domain area and monsoon precipitation for each monsoon and each hemisphere. Monsoon domains were calculated after Wang and Ding (43) and are reported in fig. S5. Monsoon precipitation is calculated within the monsoon domains. Only grid cells over land are considered.

Most proxies in our database reflect terrestrial precipitation in its annual mean, so we compare them with the annual mean of the model ensemble over land. Overall, modeled precipitation anomalies agree in sign with 64% of the proxies. This hit rate is significantly higher than in a comparison with randomized data with the same probability distribution (see Methods). It also far exceeds the hit rate calculated for the proxy dataset versus anomalies between the PI and an idealized simulation in which wet areas simply become wetter and dry areas become drier. This suggests that proxies are closer to precipitation patterns from the LIG ensemble mean than to a situation in which anomalies are only driven by thermodynamic effects. Agreement is even higher over vast regions of the NH: Namely, models and proxies agree on more LIG precipitation over Central and Northern Africa, the Middle East, South Asia, northern South America, western and northwestern North America, northeastern Asia (see the map in Fig. 2 and the regional “hit rate” between the models and proxies in Fig. 4), and, secondarily, the Mediterranean. Signals diverge over northeastern Australia (with drier LIG model anomaly and wetter LIG proxy anomaly) and over the southwestern part of western Africa (with wetter models and drier proxies). This is likely related to deficiencies in the simulation of the meridional position of the intertropical convergence zone and associated rain belt at those longitudes. Over the east Atlantic, this may imply inaccurate representation of high-latitude/Arctic warming (14). Model and proxy signals also diverge over parts of Europe: The models do not replicate proxy signals of higher LIG precipitation over eastern Europe and lower LIG precipitation over Scandinavia. Over central Asia, some wetter LIG proxies are not mirrored by the models. Over southeastern Africa, the models do not replicate the higher LIG precipitation from most of the proxies. Mismatches between LIG precipitation, as documented by proxies and as simulated by models, have multiple potential sources that cannot be assessed with the present knowledge: deficiencies in models’ representation of reality, inaccuracy in proxy reconstructions and interpretation, and inaccuracies in the experimental design of the LIG simulations due to insufficient constrains on key boundary conditions such as the Greenland ice sheet (see Supplementary text S3 for a discussion of the limits of a proxy-model comparison) (30).

Fig. 4 Precipitation anomaly between LIG and PI, comparison of models and proxies.

The hit rate is the percentage of matches between the sign of the anomaly in the model ensemble average and in the proxy records in the database (see Methods). n is the number of comparisons per region.

Opposite monsoon changes in the Northern and Southern Hemispheres

Only few proxies report seasonal anomalies of LIG precipitation. For most of these locations, the sign of the seasonal anomaly in our model ensemble agrees (Supplementary text S2). Therefore, the discussion of seasonal precipitation is mainly based on the model results. Months and seasons are defined according to the angular calendar (31), i.e., they cover the same angle in Earth’s orbit, regardless of orbital configuration. In contrast to the classic, fixed-date calendar, the angular calendar enables the comparison, across geological periods, of climate along the same portion (e.g., season) of Earth’s orbit. There is a strong seasonal contrast in precipitation anomalies, with much higher (lower) boreal (austral) summer and autumn rainfall during the LIG than during the PI (Figs. 1, 3, and 5). Zonal seasonal precipitation anomalies are more accentuated over land than over the ocean and are robust across all models. With few exceptions, models disagree only where anomalies are weak or absent. The seasonal contrast in LIG precipitation is strongest in the tropics, with maximum zonal June-July-August (JJA) anomaly of +1300 mm/year at 10°N and December-January-February (DJF) peak drying of −800 mm/year at 20°S. In JJA, LIG upflow (i.e., negative omega; fig. S3) is stronger than PI mostly over the Sahel and Sahara, the Arabian Peninsula, northern India, and northeast South America, vastly coinciding with the areas of increased precipitation. The months of September-October-November (SON) also show strong tropical anomalies (Fig. 5): wet over Central America, across a band from Central and Northern Africa to South Asia and until East Asia and dry over South America and southern Africa. Anomalies of March-April-May (MAM) are more moderate. Strong SON anomalies from our analysis are accentuated compared to previous studies (22, 23), which is partially explained by the effect of the angular definition of seasons that we apply (see Methods and fig. S4).

Fig. 5 Model ensemble average of seasonal precipitation anomalies between LIG and PI.

Seasons are defined after the angular calendar. Hatched areas are as in Fig. 2.

Terrestrial domains of boreal monsoons are 42% larger in the LIG simulations (Fig. 3 and fig. S5). Conversely, austral monsoons have 11% smaller domains, with very narrow intermodel spread. Total precipitation and extreme precipitation of monsoons, expressed with the RX5day index (i.e., maximum precipitation in five consecutive days, calculated from daily time series; see Methods), have the same sign of anomaly as the monsoon domains but vary to a larger extent (Figs. 3 and 6). Results for each monsoon system are discussed in the following paragraphs.

Fig. 6 Model ensemble average of seasonal precipitation anomalies between LIG and PI.

Top four panels: Precipitation extremes (RX5day index, i.e., maximum precipitation during five consecutive days, in millimeters). Hatched areas are as in Fig. 2. Bottom four panels: Synthesis of climatological average precipitation and precipitation extremes. Green indicates more average precipitation during the LIG. Blue indicates more precipitation extremes during the LIG. Purple indicates both more average precipitation and more precipitation extremes during the LIG.

Stronger LIG North African monsoon is the clearest from both proxies and models. While this has been suggested for a long time (32), we provide here the most robust and geographically extensive evidence yet. The modeled domain of the LIG monsoon expands by 88% (Fig. 3), to the north and to the east, and encompasses the Arabian Peninsula, while the summer tropical rain belt shifts north over the continent west of 30°E (fig. S6). Average summer precipitation and extreme summer precipitation are much stronger by 123 and 138%, respectively. Abundant proxy evidence in our database supports the north- and westward expansion of an area of higher precipitation. As in previous simulations (24), wind anomalies of our models indicate that the moisture for monsoon precipitation over the Arabian Peninsula comes from North Africa.

The Asian monsoon is also unambiguously stronger from our proxy database and models. Proxies document increased LIG annual precipitation over eastern China, also specifically for summer (Supplementary text S2) [e.g., (20)], and over the Indian subcontinent. Average and extreme simulated monsoon precipitation increase by 40 and 29%, respectively, with very narrow intermodel spread. Increase is greatest over the Himalayas and Bangladesh/Myanmar and less prominent over the south of India. This is potentially linked to a more northern LIG position of the prevailing JJA westerlies over the northwest Indian Ocean and an attendant northward shift of the intertropical convergence zone (fig. S6). The modeled monsoon domain is 26% larger; its northward expansion is supported by a wet phase of the Kesang cave record in western China (see reference in external database S1).

Our ensemble shows higher average precipitation (+34%) and extreme precipitation (+26%) of the LIG North American monsoon. The spatial extent and the intensity of the precipitation anomaly are larger than those found in a recent model study (22) and are particularly clear over northern South America. There, half of the proxies in our database document positive annual anomaly, and marine sediments from Cariaco Basin suggest a more northern intertropical convergence zone (see reference in external database S1). However, the anomaly for this monsoon is the smallest among NH monsoons, with a LIG domain only 14% larger and with considerable intermodel spread.

The Australian monsoon is much weaker in our LIG simulations, with reductions in domain area (−36%) and in average (−49%) and extreme (−38%) monsoon precipitation. This may be linked to a more northerly tropical rain belt, visible in several models (fig. S6). One proxy record shows decreased annual precipitation over the west flank, and two records show increased precipitation over the northeast of the monsoon region. The LIG domain of the South African monsoon is 9% smaller and retreats from the southeast of the continent, which is partially supported by reduced annual precipitation documented in proxies only at the south of the continent. Pollen records across the southwest Africa region have shown expansion of dry biomes during Marine Isotopic Stage 5e (33). Average monsoon precipitation and extreme monsoon precipitation decrease by 26 and 17%, respectively. The LIG South American monsoon is simulated only 3% smaller but receives less average (−22%) and extreme (−18%) precipitation, in agreement with several annual proxies. Drying is already under way before the monsoon season during SON. Over the western monsoon flank, intermodel agreement in annual mean precipitation is lower but does not contradict annual proxies at the Peruvian coast, Atacama Desert, and Salar de Uyuni (see references in external database S1).

Drivers of precipitation anomalies

At the large scale, atmospheric circulation follows changes in temperature gradients, in turn affecting precipitation patterns. Precipitation is also strongly affected by thermodynamic processes. Our new proxy database indicates widespread, although not ubiquitous, increase in boreal precipitation, concurrent with warming of all land masses north of ca. 45°N (8). The most robust empirical knowledge about large-scale LIG temperature patterns is that surface air was considerably warmer in the boreal high latitudes. This is also seen in our model ensemble (Fig. 1 and fig. S7). The causes of this are high-latitude albedo and water vapor feedbacks triggered by higher boreal insolation from April to August (34), the months of perihelion during the LIG (fig. S1). In particular, sea ice cover remains lower than the PI year-round, which acts to increase temperatures in the other seasons as well. The warming in high latitudes and associated increase in evaporation rates and in air’s water-holding capacity, confirmed in our models, can explain the moderate increase in LIG precipitation over high latitudes. The patterns of anomalies of seasonal extreme precipitation, as per the RX5day index, closely follow those of mean precipitation (Fig. 6), suggesting that they are driven by the same physical processes described above. The magnitude of the extreme precipitation changes is, however, much larger by almost a factor of 3 (fig. S8). This is roughly consistent with the concept that extreme precipitation changes after the Clausius-Clapeyron equation, whereas mean precipitation is limited by evaporation and therefore changes at a much reduced rate (16). The much amplified boreal monsoonal precipitation is related to a decreased latitudinal temperature gradient during the monsoon season (May to September), which constitutes a major driver of monsoon circulation (35). Idealized simulations with minimum precession and maximum obliquity (as during 127 ka) have yielded enhanced North African (36) and Asian (25) monsoons. The Pleistocene speleothem record of the Asian monsoon seems in fact to be modulated by NH insolation changes (29), although recently different proxies have suggested that Asian monsoon intensity responds rather to NH climate and ice sheet variations (37) (note that the latter are not represented in our models). Our ensemble of simulations suggests a link between hemispheric thermodynamics and precipitation of the North African monsoon: The latter tends to be larger for models that simulate higher LIG temperatures in the NH and in the Arctic, and lower LIG temperature gradients between NH low latitudes and the Arctic (fig. S9). The model EC-EARTH has the strongest NH and Arctic summer warming and displays the greatest North African monsoon expansion.

In our simulations, boreal monsoons last longer and persist until September, when boreal LIG insolation is only slightly higher than PI’s (fig. S1), and marginally even until October. Arctic amplification reaches its maximum during SON (Fig. 1 and fig. S7) and may therefore drive the persistence of boreal monsoon anomalies. An additional driver, as noted before for the North African monsoon (22), can be the increased local water availability in an altered water cycle. Namely, (seasonal) precipitation anomalies can be enhanced in models equipped with dynamic vegetation, which, in our ensemble, are the models MPI-ESM and NUIST-CSM (see Methods). In MPI-ESM, plant functional types transition from grass to forest over the Sahel and southern Sahara, where JJA and SON cooling is the strongest of the ensemble. This may concur to explain positive anomalies in the NH monsoon domain and precipitation in this model, which are the largest of the ensemble. NUIST-CSM shows widespread LIG increase in total leaf area index but not specifically over boreal monsoon regions, and its monsoon anomalies are slightly above the ensemble average. The model IPSL-CM has semidynamic vegetation and shows only moderate and sparse LIG increase in grass and forest types, which do not seem to affect remarkable monsoon anomalies.

Last, models show robust annual and summer drying over the midlatitudes of North America and across the Atlantic Ocean and summer drying over Europe. These anomalies are only partly supported by proxies on land. Across the midlatitude of the North Atlantic, a LIG reduction in upflow is visible (i.e., positive omega; fig. S3), indicating weakening activity of storm tracks that contribute precipitation to the adjacent continents (Figs. 4 and 6). This might also be the result of the enhanced high-latitude warming, which reduces the poleward temperature gradient and therefore baroclinicity.


There is an urgent need to evaluate the ability of climate models to reproduce actual precipitation (12). Using the LIG as a testbed, we have shown that models of the latest generation can reproduce precipitation of a regionally warmer world over large parts of the continents. Proxies and models point to stronger boreal precipitation, over the high latitudes and especially over monsoon areas. Boreal LIG monsoons in models are larger, more intense, and longer. For the austral hemisphere, proxies and models do not coincide in a clear message on LIG precipitation.

Radiative forcing of LIG climate is orbital and has a pronounced seasonal cycle with strongest forcing in boreal summer. Greenhouse gas radiative forcing of future climate is expected to be more uniform seasonally and regionally. Nevertheless, projections of future climate change bear essential similarities to anomalies of the LIG: in temperature and, as we have shown here, in boreal precipitation. As reconstructions of the LIG, future projections for low-emission scenarios feature warmer NH temperatures, especially in the high latitudes (13), and enhanced precipitation and extreme precipitation over high-latitude continents and over boreal monsoons (12). The match between proxies and models over those regions provides a degree of confidence in models’ capability of faithfully reproducing precipitation changes. In turn, this informs the reliability of projections of future precipitation changes over those regions.


Model ensemble

We used equilibrium simulations of climate at 127 ka ago. In PMIP4, part of the sixth round of the CMIP6 (28), the community has recently agreed that 127 ka should capture peak NH summer insolation and potentially, therefore, also peak warmth of the LIG in most areas (6). The forcing of the climate models follows the PMIP4 prescriptions for the lig127k experiment, i.e., eccentricity = 0.039378; obliquity = 24.04°; perihelion-180° = 275.41°; [CO2] = 275 ppm; [CH4] = 685 ppb; [N2O] = 255 ppb; and the rest as for the piControl experiment of CMIP6. Note that the configuration of ice sheets is the same in the two experiments, because proxy reconstructions of Greenland and Arctic ice sheets do not allow conclusive constrains for their representation in models for 127 ka [see discussion in (6)]. Our ensemble includes seven Earth system models (table S1): CESM1.2, EC-EARTH3.2, HadGEM3-GC3.1 (38), IPSL-CM6-LR, MPI-ESM1.2.01p1-LR (39), NorESM1-F (40), and NUIST-CSM (41). Within the model ensemble, we used five types of atmospheric modules (among these, two versions of CAM and two versions of ECHAM6), four types of ocean modules (among them, four versions of NEMO), and five types of vegetation/land modules (among them, two versions of CLM and two versions of JSBACH). Main relevant differences are in the resolution of the atmospheric modules (from T159 to T63 spectral truncation and from 26 to 79 vertical levels) and in the representation of vegetation dynamics: dynamic, semidynamic, or prescribed. Vegetation in MPI-ESM and NUIST-CSM is dynamic, and in IPSL-CM, it is semidynamic, meaning that the maximum fraction of each plant functional type is prescribed, but the leaf area index for each type responds to changes in radiation and moisture and interacts with the carbon cycle. In addition, the hydrology in IPSL-CM is based on 11 vertical layers to more realistically simulate changes in evaporation and moisture availability. We calculated climatological statistics (LIG and PI climatological averages are presented in fig. S10) and extremes using the time series of daily output from the models. These time series reflect climate close to equilibrium, i.e., after a spin-up model period of several hundred years when the climate is allowed to adjust to the altered forcing. Spin-up time was considered sufficient when certain metrics, such as top-of-the-atmosphere radiation imbalance, were close to zero. The trend stationarity of the time series used here was statistically verified with the augmented Dickey-Fuller test. The length of the time series is ca. 300 years for the models CESM, NorESM, and MPI-ESM and 100 years for NUIST-CSM, IPSL-CM, EC-EARTH, and HadGEM. For the analysis, months and seasons were defined according to the angular calendar, also known as celestial calendar (31). The angular calendar is preferable to the classic, fixed-date calendar when comparing seasons across geological periods, such as LIG and PI, since it enables to compare climate along the same portion (e.g., a season) of Earth’s orbit. Seasonal precipitation and temperature anomalies between the LIG and the PI, obtained by applying either the classic or the angular calendar to our time series, show non-negligible differences and even significantly different geographic patterns for SON (fig. S4). We calculated total precipitation as the sum of the model fields corresponding to convective and large-scale precipitation, including liquid and snow (this is simply called “precipitation” in the main article). Near-surface air temperature is at 1.5 or 2 m of elevation. To estimate changes in atmospheric water vapor content, we used the model fields corresponding to either total vertically integrated precipitable water or total vertically integrated water vapor, depending on the model. With either of these model fields, all models show positive anomaly of the global annual mean in the LIG simulation. Indices of extreme precipitation were calculated with the Python package ICCLIM (Indice Calculation CLIMate), developed at CERFACS, France. We calculated the standard indices: RX5day (maximum precipitation in five consecutive days), RX1day (maximum daily precipitation), R10mm (number of days with more than 10 mm of precipitation), and R20mm (number of days with more than 20 mm of precipitation). From these, we calculated extreme values for the entire time series at the return periods of 20, 50, and 100 years. Although values of these different indices differ, the geographic pattern of anomalies in extreme precipitation of the LIG with respect to the PI remains largely the same. We therefore only showed here extremes with the RX5day index (maximum rainfall in five consecutive days). Although the skill of models in reproducing extremes is limited, because of known challenges in representing convection and to general biases (42), expressing precipitation extremes as anomalies of LIG from PI climate may partially compensate the bias. In addition, to avoid overinterpreting fine-scale phenomena where precipitation and extremes are vulnerable to model bias and low skill, we limited our discussion to the large-scale patterns. Values of the Monsoon Precipitation Index and monsoon domain areas were calculated according to the method proposed in (43), whereby the threshold of the Monsoon Precipitation Index and the threshold of the seasonal precipitation range used to calculate the monsoon domain areas have been modified to 0.8 and 600 mm, respectively, to obtain reasonable areas for the PI simulation. We applied the angular calendar to define monsoon months. Ensemble means were calculated after regridding the results of each model by bilinear interpolation. Ensemble median values are negligibly different from the mean. For the ensemble, the bias of the PI simulation was assessed via comparison of its ensemble mean precipitation with the climatological average of the NOAA-CIRES 20th Century Reanalysis V2c dataset for the period of 1851–1880 (closer in time to the PI) and for the period of 1900–1929 (for which the confidence on the product is higher) from and with the climatological average of the ERA-20C Reanalysis for the period of 1900–1929 (fig. S11). The bias is very similar when calculated with the two periods of the NOAA-CIRES product. However, there are differences between the bias with the NOAA-CIRES and the ERA products. The following are most notable: over central North America, over central regions of South America, over the southern part of Western Africa and over parts of Central Africa, over northern Europe and Siberia, and over Southeast Asia. In general, the ensemble PI mainly seems to overestimate PI precipitation over northern Australia, the northeast of South America, northwestern America, southern Africa, and central Asia and to underestimate precipitation mainly over small parts of South Asia and South America.

Proxy database and model/proxy comparison

We compiled the evidence for LIG precipitation at 127 ka from the existing peer-reviewed literature. The database comprises 138 records of proxies based on the following types of paleoclimatic archives: 38 based on pollen, 28 based on lacustrine or marine sediment composition, 21 based on speleothems, 20 are multiproxy, 10 based on landscape features, 7 based on fossils other than pollen, 5 based on loess deposits, and 9 based on other types of archives. In the database (Annex 1), we included precipitation anomalies (between the LIG and the present/PI climates), as explicitly stated in the original publication and as interpreted by the authors’ expert judgment. Unlike for temperature, most proxies for precipitation in the literature are not quantitative. We assigned to each record a value on a semiquantitative scale to reflect the anomaly between the LIG and the present/PI as follows: +2 (much wetter LIG anomaly; dark blue in Fig. 2), +1 (wetter; light blue), 0 (no noticeable anomaly; white), −1 (drier; light red), and −2 (much drier; dark red). This was done on the basis of the authors’ expert judgment on the magnitude of the deviation between the signal of the LIG and of the present/PI/core top.

To calculate the hit rate between proxies and models, as shown in Fig. 4, we bilinearly interpolated the model ensemble average of precipitation anomalies for the coordinates of the proxy record; we considered as “hit” the cases in which the sign of ensemble average and of the proxy correspond; for proxy value equal to 0 (no noticeable anomaly), we considered as hit the cases in which the model anomaly comprised between −40 and + 40 mm/year. The significance of overall hit rate of 64% between proxies and models was tested as follows. (i) On the basis of the frequency distributions of anomaly values in the proxy database, we calculated the probability with which values fell within the ranges established by the above definition of the hit rate; we calculated probabilities in the same for values from the model ensemble average map (considering only grid cells on land). (ii) We generated one random value to represent the proxy anomaly and one random value to represent the model anomaly using the probabilities estimated in step 1 as weights. (iii) We considered the comparison between the randomly generated proxy value to the randomly generated model value a hit if the conditions described above were met. (iv) We repeated steps 2 and 3 138 times, corresponding to the actual number of proxy-to-model comparisons in this study and thus obtained a synthetic hit rate. (v) We repeated step 4 100,000 times. The average synthetic hit rate was 50.4%; a hit rate of 64% was obtained in only 0.06% of the iterations (P = 0.0006) and, as such, was very unlikely to be a result of chance. We did not account for spatial correlation in the proxy database because proxy locations, in general, are not clustered and are nearly randomly distributed in space with considerable distances from each other. We also calculated the hit rate between the proxies and a hypothetic result of a model simulation in which, based on PI climate, “wet gets wetter and dry gets drier,” i.e., areas that received below (above) a given threshold of precipitation in the PI were assigned a positive (negative) anomaly value. For plausible thresholds of 500, 800, and 1000 mm/year, the proxy/model hit rate was 42.8, 27.5, and 23.5%, respectively. The general results of the hit rate between models and proxies are robust with respect to different ways of calculating the hit rate, for example, by using percentage anomalies instead of absolute anomalies from models and using different thresholds to classify anomalies as null. Additional considerations regarding the limits of the proxy database and the proxy-model comparison are included in Supplementary text S3.


Supplementary material for this article is available at

Supplementary Text S1. Temperature anomalies from models.

Supplementary Text S2. Comparison of seasonal precipitation between proxies and models.

Supplementary Text S3. Limits of proxy dataset and of proxy-model comparisons.

Fig. S1. Anomaly in insolation between the LIG (127 ka) and the PI (year 1850 CE).

Fig. S2. Comparison of annual anomalies between LIG and PI.

Fig. S3. Ensemble average anomaly between LIG and PI for variables diagnostic of atmospheric circulation.

Fig. S4. Difference between seasonal ensemble anomalies (between LIG and PI) after the angular and the classic calendar.

Fig. S5. Monsoon Precipitation Index, in contoured colors, and monsoon domains, in contoured black lines.

Fig. S6. Approximate position of the ensemble mean seasonal intertropical convergence zone (over the oceans) and tropical rain belt (over the continents).

Fig. S7. Model ensemble near-surface air temperature anomalies between LIG and PI.

Fig. S8. Correlations between LIG-PI anomaly in average precipitation and in extreme precipitation, per season, for the ensemble average.

Fig. S9. Selected sets of correlations between anomalies in temperature and North African monsoon for each model in the ensemble.

Fig. S10. Precipitation from the ensemble mean of the PI and LIG simulations.

Fig. S11. Bias of the PI ensemble annual average precipitation with respect to the closest comparable reanalysis products.

Table S1. Details of the models in the ensemble, with main information regarding the modules relevant to this study.

External database S1. Proxy database.

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 acknowledge the World Climate Research Programme and its Working Group on Coupled Modeling that coordinated CMIP6. We thank the climate modeling groups for producing and providing model output and the Earth System Grid Federation (ESGF) for archiving the data and providing access. CESM1.2 simulations were carried out on the supercomputer of the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN). EC-EARTH simulations were performed at Swedish National Supercomputer Centre provided by the SNIC. The IPSL-CM6experiments were performed using the HPC resources of TGCC under the allocations 2016-A0030107732, 2017-R0040110492, and 2018-R0040110492 (project gencmip6) provided by GENCI (Grand Equipement National de Calcul Intensif). The IPSL-CM6 team of the IPSL Climate Modelling Centre ( is acknowledged for having developed, tested, evaluated, and tuned the IPSL climate model, as well as having performed and published the CMIP6 experiments. MPI-ESM simulations were conducted at the Deutsches Klimarechenzentrum (DKRZ). NorESM1-F simulations were performed on resources provided by the UNINETT Sigma2–the National Infrastructure for High Performance Computing and Data Storage in Norway (nn4659k and ns4659k). We thank P. Bartlein for providing angular calendar information. M. Sicard is acknowledged for support with the IPSL-CM datasets. We are grateful to the editor and to four anonymous reviewers for allowing us to improve an earlier version of the manuscript. Funding: P.S. acknowledges funding from the NWO (Nederlandse Organisatie voor Wetenschappelijk Onderzoek) under grant ALWOP.164. Q.Z. acknowledges funding from Swedish Research Council VR (2013-06476 and 2017-04232). C.S. acknowledges funding by the Helmholtz Climate Initiative REKLIM and the Alfred Wegener Institute’s research program Marine, Coastal and Polar Systems and provision and support for the MPI-ESM suite by the Max Planck Institute for Meteorology in Hamburg. P.Ba. and M.P. acknowledge support from the German climate modeling initiative PalMod (BMBF). M.-V.G. acknowledges the financial support of the NERC research grant NE/P013279/1. NorESM1-F simulations received funding from the European Research Council under the European Community’s Seventh Framework Program (FP7/2007–2013)/ERC grant agreement 610055, ice2ice project. B.O.-B. acknowledges the support of the National Science Foundation; this material is based on work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement no. 1852977. We thank the multiple funding agencies that support CMIP6 and ESGF. This study benefited from the ESPRI computing and data center (, which is supported by the CNRS, Sorbonne Université, Ecole Polytechnique, and CNES, as well as through national and international grants, particularly Labex L-IPSL (ANR 10-LABX-0018). Author contributions: P.S., P.Ba., H.R., P.J.W., and J.C.J.H.A. designed the research. P.Ba., C.G., C.S., Q.Z., P.Br., J.C., M.-V.G. performed the climate model experiments. P.S. compiled the proxy database, analyzed the model datasets, and prepared the figures. P.S. wrote the manuscript with the contributions from P.Ba., C.G., C.S., D.C., P.Br., M.P., P.J.W., H.R., M.K., and J.C.J.H.A. All authors revised the manuscript. 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. The proxy database is available in external database S1. Additional data related to this paper and datasets corresponding to the analysis of the climate model results may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article