Lake Tauca highstand (Heinrich Stadial 1a) driven by a southward shift of the Bolivian High.

Heinrich events are characterized by worldwide climate modifications. Over the Altiplano endorheic basin (high tropical Andes), the second half of Heinrich Stadial 1 (HS1a) was coeval with the highstand of the giant paleolake Tauca. However, the atmospheric mechanisms underlying this wet event are still unknown at the regional to global scale. We use cosmic-ray exposure ages of glacial landforms to reconstruct the spatial variability in the equilibrium line altitude of the HS1a Altiplano glaciers. By combining glacier and lake modeling, we reconstruct a precipitation map for the HS1a period. Our results show that paleoprecipitation mainly increased along the Eastern Cordillera, whereas the southwestern region of the basin remained relatively dry. This pattern indicates a southward expansion of the easterlies, which is interpreted as being a consequence of a southward shift of the Bolivian High. The results provide a new understanding of atmospheric teleconnections during HS1 and of rainfall redistribution in a changing climate.

INTRODUCTION

A major uncertainty in our understanding of 21st century climate concerns the global redistribution of rainfall and modification of monsoonal systems in response to global warming (). Past global climatic oscillations are of particular interest for understanding the teleconnections that lead to rainfall modifications during periods of global scale climatic change. These oscillations constitute a baseline against which various scenarios can be tested (, ) and provide case studies of abrupt redistributions of rainfall and major regional hydrologic changes. The Heinrich stadials that characterized late Pleistocene climate were marked by massive iceberg discharges accompanied by global reorganizations of the atmosphere-ocean system, involving major modifications of monsoon dynamics and the positions of the westerly wind belts in both hemispheres (–). Understanding the climatic processes at work and their interconnections during these periods is therefore crucial. In South America, Heinrich Stadial 1 [HS1; 18.5 to 14.5 thousand years (ka) before present (BP) ()] presents one of the most notable examples of abrupt rainfall redistribution and major hydroclimatic change: the transgression, highstand, and regression of the giant paleolake Tauca in the endorheic Altiplano basin (Bolivian Andes). Lake Tauca reached a depth of 120 m and covered an extent of 52,000 km2 during the second half of HS1 (HS1a; 16.5 to 14.5 ka BP). It is the largest paleolake expansion in the Altiplano in the last 130 ka (–).

The distribution of precipitation responsible for such marked hydroclimatic changes remains unknown. Most previous studies of South American paleorainfall were based on oxygen isotope records from speleothems in the northeast (NE) and southwest (SW) Amazonian lowlands. These studies show that atmospheric circulation over subtropical Brazil changed in pace with the Heinrich stadials (, –) and numerous sites record wetter conditions during these episodes [see compilations (, )]. Nevertheless, these records do not provide spatially distributed information on the regional precipitation changes that induced the Tauca highstand. Over the Altiplano, two studies that couple lake and glacier mass balance models propose a range of overall temperature shifts of between −5° and −7°C and an overall annual rainfall increase of between 70 and 150% (, ) during the Tauca highstand. However, these studies were unable to determine the spatial variability of the precipitation during this period, which precludes identification of the moisture sources and mechanisms associated with the substantial redistribution of rainfall concurrent with the worldwide climatic modifications of HS1. The persistent La Niña–like conditions that are assumed to have occurred in the central Pacific during the Tauca cycle have been suggested as a possible forcing for the precipitation anomaly (), but none of paleoclimatic data available to date have been able to confirm this theory.

Today, two modes of moisture transport, involving different atmospheric features of the South American summer monsoon (SASM), have been identified for the Altiplano basin (). In the northern Altiplano basin, rainfall anomalies result from advection by the easterlies of tropical moisture driven south by the low-level jets. Precipitation anomalies located further south on the Altiplano result from southward deflection of the Bolivian High (BH), an upper troposphere high-pressure cell centered over Amazonian Bolivia that develops during the austral summer (, ) and promotes moisture transport from the Gran Chaco. The interplay between these two rainfall modes creates interannual variability, which has already been discussed in the context of central Andes paleoclimate (, ). However, the relative contribution of each component during the hydroclimate changes responsible for the Lake Tauca highstand remains unknown. Constraining the spatial distribution of precipitation during this period is therefore crucial to (i) track the source of humidity, (ii) understand the connections between the Altiplano and the South American climate dynamics, and (iii) identify possible interhemispheric forcings that link the SASM to North Atlantic cooling events (, ).

Here, we gather a large set of new cosmogenic exposure ages and previously published ages to identify a set of nine valleys distributed over the Bolivian Altiplano that display moraines left by glacial stillstands or readvances synchronous with the Lake Tauca highstand (Fig. 1). By combining a precipitation-temperature–dependent glacier model at these sites with lake budget modeling of Lake Tauca, we then use an inverse method to derive a paleoclimate reconstruction that incorporates a regional cooling and a quantitative precipitation field over the Altiplano during the Lake Tauca highstand. The spatial pattern of precipitation in the reconstruction enables identification of the atmospheric processes involved in this unparalleled redistribution of rainfall linked to global climate change.

Fig. 1

Paleoglacial extent of glacial cover during the Lake Tauca highstand (16.5 to 14.5 ka BP) at the nine paleoglaciated mountains where a moraine synchronous with the Lake Tauca highstand was identified.

(A) Nevado Sajama (data from this study). (B) Cerro Pikacho (data from this study). (C) Cerro Luxar (data from this study). (D) Cerro Uturuncu (). (E) Cerro Tunupa (, ). (F) Cerro Tambo (data from this study). (G) Cerro Azanaques (). The additional feature mapped here is the Challapata fan delta and its boulder field (light orange; see section S1.2.1). (H) Zongo Valley. The white arrows indicate the former ice-flow direction (). (I) Laguna Aricoma (). Aerial photographs are from Google Earth. Paleo-ELAs were reconstructed using the AAR method (see the Supplementary Materials for more details).

Settings and approach

We can explore past climate conditions by simultaneously reconstructing the fluctuations of lakes and glaciers because these variations both depend on precipitation (rainfall and snowfall) and temperature (soil and lake evaporation and ice/snow melting). The abundant shoreline records and ubiquity of glacial morphologies over the Altiplano make this location particularly suited to this approach. The nine glacial sites we rely on for this study span from 14.3°S to 22.3°S and 66.5°W to 69.8°W, ensuring good coverage of the Altiplano basin (Fig. 1). The Laguna Aricoma, Zongo Valley, Cerro Azanaques, and Cerro Tambo sites spread from north to south in different subranges of the Eastern Cordillera and present various petrologies, from siliciclastic metasedimentary rocks to granite and granodiorite. The other sites correspond to isolated andesitic and dacitic volcanoes located close to the western (Nevado Sajama, Cerro Pikacho, and Cerro Luxar) and southern (Cerro Uturuncu) edges of the basin. Cerro Tunupa occupies a central position close to the center of paleolake Tauca. All of these sites contain a detailed sequence of moraines, probably formed during the Last Termination (from the Last Glacial Maximum to the Holocene), that can be used to produce a precise reconstruction of the former successive glacial extents (see the Supplementary Materials).

Given the available lithologies along the basin, the total exposure age data set includes both 10Be ages in quartz and 3He ages in pyroxenes. For both isotopic systems, recent independent and convergent studies reported robust and locally calibrated in situ production rates (–). For consistency, all ages were systematically recalculated using the online cosmic ray exposure program (CREp) calculator (), following a scaling procedure shown to be appropriate for the high tropical Andes (see the Supplementary Materials) ().

For all sites, determination of the extent of ice coeval with the Lake Tauca highstand required a comprehensive view of the broad deglaciation sequence. Stratigraphic relationships between successive or synchronous glacial landforms provided additional relative time constraints that allowed one moraine coeval with the highstand to be identified at each site (Fig. 2 and see the Supplementary Materials). The equilibrium line altitude (ELA) was then calculated for each paleo-ice extent using the accumulation area ratio (AAR) method (). The AAR used for this calculation, 0.63 to 0.73, was derived from monitoring of modern glaciers in the high tropical Andes (see the Supplementary Materials for the method) (, ).

Fig. 2

Map of the Altiplano basin showing locations of the nine paleoglaciated sites and associated geochronological constraints (3He and 10Be dating in thousand years BP) and ELA (in meters above sea level).

CRE ages (presented in age probability density function plots) at each site (yellow points) constrain the extent of glacier cover coeval with the Tauca highstand (16.5 to 14.5 ka BP, above 3760 m, dark blue contour on the map, and vertical blue bands on the age plots). The horizontal axis of each graph shows the CRE age (in thousand years before 2010). The blue number in each age plot is the associated ELA in meters above sea level (see the Supplementary Materials)

The determination of paleoclimatic conditions was achieved by exploring incremental values of uniform cooling (ΔT) over the Altiplano basin. For each incremental value, we derived local precipitation values at each glacial site using a robust empirical relation that links the ELA position to local precipitation and temperature (see fig. S16) (). By extrapolating from these nine precipitation values, we built a precipitation grid over the Altiplano. We then used a distributed lake budget model (, ) with a quarterly time resolution and finally retained the optimal values of the uniform cooling (ΔT) and the precipitation grid over the Altiplano to allow us to balance the annual hydrologic budget of Lake Tauca at its highstand. Present-day climatic observations over the Altiplano basin support the assumption of uniform cooling required by this algorithm (see the Supplementary Materials). In addition, because the precipitation grid was interpolated from precipitation values from nine glacial sites, it cannot capture small-scale precipitation anomalies such as that currently observed over Lake Titicaca (Fig. 3A). To compare present-day and Tauca precipitation regimes, we therefore downsampled the present-day rainfall map at a similar resolution (see the Supplementary Materials). Finally, for estimating the uncertainties in the map reconstruction, we applied a Monte Carlo method to propagate both those associated with the present mean temperature and annual rainfall at each site and that associated with the ELAs.

Fig. 3

Annual precipitation maps for the Altiplano—modern and Tauca highstand.

(A) Current annual rainfall over the Altiplano (). White dots, sample sites; blue contours, Lakes Titicaca and Tauca. (B) Tauca highstand annual rainfall derived from our reconstruction. (A) and (B) share the same color scale. (C) Difference between the Tauca and modern annual rainfall. (D) Rainfall amplification during the Lake Tauca highstand compared to Present. Precip., precipitation. (E) Relative uncertainty (uncert.) in the precipitation results (uncertainties in the present-day climate and ELA were propagated using Monte Carlo methods; see the Supplementary Materials). (F) Hydrological and glacial map.

RESULTS AND DISCUSSION

Our reconstruction provides quantitative constraints on the climatic conditions associated with the Tauca highstand. We derive a uniform cooling of 2.9° ± 0.2°C from modern for the whole endorheic catchment and an average value of 900 ± 200 mm for the annual rainfall, corresponding to a 130% increase compared to the present. The annual paleorainfall grid (Fig. 3B) shows a clear and strong N-S to NE-SW gradient characterized by abundant precipitation not only over the Titicaca watershed in the northern Altiplano but also further south along the Eastern Cordillera, defining a band of abundant precipitation (1200 to 1300 mm year−1) spreading from 14°S to 19°S. In contrast, the southwestern part of the Altiplano is the driest zone with a mean annual precipitation of 100 to 500 mm year−1 (Fig. 3B). The relative uncertainty in the annual precipitation value ranges from 9 to 23%, with a mean value of 14 ± 3% (Fig. 3E), confirming the robustness of the overall spatial pattern that we report for the basin. Given the strong contrasts that we calculate, even an overall relative uncertainty of 20% would not significantly alter this spatial distribution or the related discussion. Additional sensitivity tests on the AAR and the ELA-P-T relation were performed to establish the robustness of the results and can be found in the Supplementary Materials.

A map showing the precipitation difference (Fig. 3C) exhibits several important features: (i) During the Tauca phase, the whole central Altiplano received on average significantly more precipitation than today (+500 mm; that is, a 130% increase, or a 2.3× amplification); (ii) the maximum increase is observed along the central-eastern part of the basin, whereas the southwestern part of the basin was almost as dry as today; and (iii) the maximum rainfall along the Eastern Cordillera was clearly further south than it is in the present spatial distribution; however, this southward shift is limited in extent to 4° of latitude as it does not affect the extreme south of the basin (Fig. 3B). The asymmetric precipitation anomaly during the Tauca phase arises from a regional modification of moisture transport—in this case, an enhanced eastern input of moisture from the southwestern Amazonian basin, across the Eastern Cordillera.

A further observation that can be established from the map showing the ratio of the Tauca to the present-day precipitation field (Fig. 3D) is that the rainfall amplification is much greater in the center of the paleolake Tauca (where the maximum anomaly reaches a factor of 4) than it is in the surrounding regional field (where the rainfall amplification is smaller and relatively homogeneous). As it is centered on the paleolake, this anomaly probably results from local recycling of moisture in a daily cycle of evaporation of water during the day and rainfall at night. These anomalies are also observed today over large lakes, such as Lake Titicaca or Lake Victoria (). The occurrence of this type of anomaly was previously postulated for Lake Tauca (), and our study confirms this hypothesis. While this local rainfall anomaly contributed to the Lake Tauca highstand, the spatial mismatch between the maximum precipitation ratios (Fig. 3D) and the maximum absolute increase in rainfall along the eastern margin (Fig. 3C) highlights the dominant contribution of moisture from sources located outside the Altiplano. Overall, our study quantifies the inner South American paleoprecipitation at a regional scale with unprecedented spatial resolution.

Today, the SASM is the main driver of precipitation over the Altiplano basin, and our results should therefore be discussed in the framework of the present-day SASM and its past variations. The Altiplano today receives between 100 and 900 mm of annual rainfall from the SW to the NE, mostly during the austral summer [December, January, and February (DJF)], which corresponds to the mature phase of the SASM (Fig. 4) (). During this period, Atlantic moisture advected by the trade winds penetrates the tropical latitudes of the Amazonian basin and is channeled southward along the Andes by low-level jets. South of this, strong convective activity over central and southern Brazil is linked with the development of the South Atlantic convergence zone (SACZ), a southeast-oriented band of precipitation that extends from southern Amazonia toward southeastern Brazil ().

Fig. 4

Paleoclimatic implications for South America during HS1—Tauca highstand.

Left: Modern features of South American climate with particular focus on the summer monsoon. SWW, southern westerly winds; AMOC, Atlantic meridional overturning circulation; ITCZ, intertropical convergence zone; SALLJ, South American low-level jets. The colored background shows the mean annual rainfall. Asterisk indicates the color scale truncated at 3500 mm. Blue contours show the DJF/annual precipitation ratio. Precipitation data are from ERA-Interim (), with mean values between 1979 and 2016. Right: Modification of South American climate during the second half of HS1, as suggested by our data. Compared to present, the BH is intensified and located further south. The SACZ is intensified, and the SWW are also displaced southward. These features all enhanced the transport of moisture onto the Altiplano. Concurrently, the ITCZ shifted southward, and AMOC intensity was reduced. Blue and red symbols correspond to climatic records that report wetter and drier conditions, respectively, during the Tauca highstand. Circles, δ18O results from speleothems; rectangles, δ18O from lacustrine or marine cores. Points 8, 10, and 11 are locations of the Jaragua, Paixão, and Lapa Sem Fim speleothems, respectively (, ). White star indicates the location of a core that implies a southward shift of the SWW (). See the Supplementary Materials for all site names and references.

Over the Altiplano, during the same period (DJF), the weakening of the dry westerly winds in favor of wet easterlies promotes moisture transport from tropical and subtropical Brazil (, , , ). This westward transport is not uniform over the Altiplano basin; however, in the north of the Altiplano rainfall is fed by tropical moisture brought to the basin from the northeast (Fig. 3A). The El Niño–Southern Oscillation (ENSO) conditions modulate the present-day interannual variability of this NE rainfall contribution over the Altiplano. La Niña periods are characterized by strengthened easterlies and precipitation anomalies from the NE precipitation source. In the center and south of the basin, easterly moisture transport is linked to the formation of the BH. The anticlockwise circulation associated with this high-pressure cell generates upper-level easterly flows toward the Altiplano basin, which advect moisture from the Gran Chaco and generate precipitation (Fig. 4). This second mode of moisture advection, modulated by the position of the BH, is independent of ENSO conditions (, ).

In the context of the Tauca highstand, these present-day atmospheric controls can aid interpretation of our rainfall reconstruction. Although tropical moisture inputs are responsible for precipitation anomalies restricted to the north of the basin, a southward shift of the BH during the austral summer creates significant precipitation anomalies further south in the basin (). Therefore, the rainfall pattern that we report for the Tauca period strongly suggests that a southward shift of the seasonal BH was the major driving force behind the Tauca highstand precipitation.

Today, precipitation anomalies linked to a southward shift of the BH are also associated with an intensification of the atmospheric circulation (, ). Our paleoprecipitation estimates also show a large overall increase in precipitation (130% average increase), and we suggest that the southerly shift of the BH was accompanied by an intensification of the BH allowing greater transport of moisture to the Altiplano during the Tauca highstand (Fig. 4).

Intensification and southward movement of the BH have already been predicted from atmospheric modeling to be the main contributor to rainfall anomalies in the central Andes during the Little Ice Age (). The possible role of the BH in glacial dynamics during the late glacial of the central Andes has also been discussed (, ). Our results provide independent constraints and additional validation of the crucial role of the southward shift of the BH. Furthermore, our interpretation is consistent with several modifications of the SASM during HS1. The seasonal onset and withdrawal of the BH presently proceeds in pace with the southward and northward migration of the dry westerlies, which prevents moisture transport from the Amazonian basin to the Altiplano during the dry austral winter (June, July, and August) (). In a paleoclimatic perspective, a long-term southward shift of the BH is thus consistent with the prolonged southward shift of the westerlies () observed during periods of AMOC weakening, such as HS1a (). This common southward movement underscores the consistent responses of two synoptic climatic features of South American climate to the Heinrich stadials.

The BH today develops in response to condensational heating over the Amazon basin (, ) and its north-south position is related to the position and intensity of the SACZ (). The modern climatology, which ties the BH to the SACZ dynamic, highlights the probable influence of the SACZ on the Tauca highstand. The SACZ has experienced periods of intensified activity over the lateglacial period [mega-SACZ events ()], one of which (16.1 to 14.7 ka BP) was coeval with the highstand. These events are periods of major convective activity in the SACZ and are responsible for large amounts of precipitation over central and southern Brazil. Similarly, the δ18O signal of the Jaragua Cave speleothem (Fig. 4) (), located between the SACZ influence zone and the Altiplano basin, also records a wet period between 16.0 and 14.8 ka, also concurrent with the highstand (Fig. 4). Hence, the second half of HS1 seems to be a distinctive period of continental-scale modification of SASM activity, during which the atmospheric conditions over central and southern Brazil were favorable to a modulation of the BH dynamics and provided significant levels of moisture in the source region.

Observations suggest that the precipitation transport outlined here is independent of ENSO modulations (, ) and is, instead, controlled by humidity levels in the source area located east of the Altiplano (, ). Our results do not exclude potential superimposition of a Pacific forcing on the influence of the BH, given that La Niña conditions bring increased precipitation to the north of the basin. However, the locations of the precipitation maxima (1250 to 1300 mm year−1) at between 16°S and 19°S rather than in the north of the basin (1000 mm year−1 from 14°S to 15°S) do not support this hypothesis. As our results highlight the driving role of the BH in the precipitation regime, this suggests that atmospheric changes that led to the Tauca highstand were dominantly controlled by the moisture levels observed over the central Brazil–Gran Chaco region, under the influence of the SACZ and, in the end, an Atlantic forcing ().

During the first half of HS1 (18.5 to 16.5 ka BP), the Altiplano basin was not as dry as today, as attested to by shoreline deposits. The period from 18.5 to 17 ka BP exhibits shorelines around 3660 m, indicating a shallow and small extent lake (maximal depth, <20 m; area, <10,000 km2). The 17 to 16.5 ka correspond to the transgression of the lake (). This progressive increase in the lake level indicates that the atmospheric mechanisms at work during the highstand have been evolving during the first half of HS1 to reach a stable state that permitted the highstand of the Tauca in the Altiplano basin. The midwestern and central-eastern Brazilian speleothem records show convergent behaviors during the first half of HS1 (, ), indicating a consistent dynamic of the SASM for this period as well. As it is designed, our approach only allows us to produce results for the Tauca highstand. For this reason, our discussions and interpretations are restricted to this 2-ka period (16.5 to 14.5 ka BP). Considering the short-time response of the lake to a rainfall change [0.1 to 0.2 ka; see ()] compared to the durations of the highstand and transgression, precipitation was necessarily much higher during the highstand than during the transgression. The rapid (<1 ka) regression of Lake Tauca after 14.5 ka BP () implies a sudden return of the BH to a more northerly location and possibly weaker circulation.

These abrupt changes in regional-scale atmospheric dynamics are supported by other paleoclimatic studies [for example, (, , –)]. Considering the global climatic changes reported during HS1, a southward shift of the BH is consistent with the concurrent southward shifts of (i) the ITCZ, as observed in offshore and continental records for same period (, ); (ii) the North Hemisphere (NH) westerlies as reproduced in modeling (, ) under the influence of NH ice sheets; and (iii) the southern westerlies (, ). Thus, our new results support the hypothesis of a global southward shift of atmospheric circulation features during the Heinrich events. They highlight the propensity of the atmosphere to propagate and modulate climatic perturbations across both hemispheres and to drive major redistributions of rainfall over very short time scales. Ultimately, our results emphasize the versatility and interdependence of Earth’s hydroclimatic systems under a changing climate.

MATERIALS AND METHODS

For a detailed presentation of the geological settings, methods, results, and sensitivity tests, see the Supplementary Materials.

Cosmogenic nuclide measurements

Samples for 10Be dating were prepared and processed at Centre de Recherches Pétrographiques et Géochimiques (CRPG; all steps from quartz extraction through to chemical insulation of the in situ produced cosmogenic 10Be component) and then analyzed at the ASTER French National AMS (accelerator mass spectrometry) Facility [Centre Européen de Recherche et d’Enseignement en Géosciences de l’Environnement (CEREGE)]. For each analytical session, the blank correction was taken as the averages of 2 to 5 10Be/9Be full-process blanks realized during the same analytical session. This average value was then subtracted to the 10Be/9Be ratio of each sample of the session to calculate the 10Be concentration. It induced a correction of 1.4 ± 1.5% on the values (maximum of 5.5%). Uncertainty on 10Be concentration of a sample includes the uncertainty of the 10Be/9Be ratio of the sample and the uncertainty over the blanks (taken as the standard deviation of the blank values of the session).

Samples for 3He dating were processed and analyzed at CRPG. Black and green pyroxenes were first examined and identified using a scanning electron microscope and then handpicked under a binocular microscope. Total 3He concentrations were measured on an split flight tube mass spectrometer, as described in (). Pyroxene aliquots were fused in a vacuum furnace at 1400° to 1500°C for 15 min. The extracted gas was purified using activated charcoals, getters, and a cryogenic pump, before measurement of the 3He and 4He abundances in the mass spectrometer. The furnace blanks induced a mean correction of 4 ± 3% (maximum of 12%).

Cosmic ray exposure ages calculation

For the sea-level high-latitude production rate used to calculate 10Be CRE ages, we used the uncertainty-weighted mean of the production rates of Blard et al. (), Martin et al. (), and Kelly et al. (). All of these studies, particularly, Blard et al. () and Martin et al. (), provide robust and local production rates for the high tropical Andes over time spans that are similar to the age of Lake Tauca (ca. 15 ka). The production rate value used here can therefore be considered a synthetic local value that is both robust and able to limit age dependence on scaling procedures ().

The same method was used for calculation of 3He ages, using the weighted mean of the local 3He production rates of Delunel et al. () and Blard et al. (). Age calculations were based on the Lal-modified time-dependent scaling scheme (, ), along with the ERA-40–spatialized atmosphere () and the geomagnetic database of Muscheler et al. (), and were performed on the CREp calculator ().

ELA calculations

ELAs were determined using the AAR method. The AAR parameter was derived from glacier monitoring observations conducted by the GLACIOCLIM–IRD (Les GLACIers, un Observatoire du CLIMat–Institut de Recherche pour le Développement) National Observation Service at (i) Huayna Potosi, summit of the Zongo Valley, (ii) Antizana (Ecuador, 0.5°S to 78.1°W), and (iii) Artesonraju (Perú, 9.0°S to 77.6°W), over the periods 1991–2010, 1995–2010, and 2004–2011, respectively (, ). To derive the AAR values, we performed a linear regression on the AAR = f(mass balance) relationship and then computed the AAR0 value (AAR for mass balance = 0). This AAR0 value represents the AAR for a steady-state glacier with a null mass balance, as is required to study moraine stillstands during the Tauca highstand. We calculated AAR0 values of 0.65 for the Zongo glacier (r2 = 0.73), 0.72 for the Antizana (r2 = 0.83), and 0.68 for the Artesonraju (r2 = 0.77; fig. S12). We used an AAR range of 0.63 to 0.73 to take into account the uncertainty associated with the paleo-ELA determination. The hypsometry data required for the AAR method were derived from the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global digital elevation models from NASA–U.S. Geological Survey (USGS).

The glaciers used for the calibration have various orientations, and the AAR range that we used only resulted in 30 to 50 m uncertainty on the ELA values, making the spatial ELA gradient highly significant. As AAR values are sensitive to the climate settings, we moreover propose an iterative use of our model that optimizes spatialized values of AARs over the basin in section S4.3.7. This alternative way to reconstruct the Tauca highstand precipitation yields similar results to those presented here, which support the robustness of the calculation.

Glacier modeling

We combined glacier and lake modeling in an inversion algorithm to identify a precipitation grid for the Altiplano basin that satisfied both glaciers and lakes extents during the Lake Tauca highstand. For the glacier model, we used the ELA = f(P,T) model of Condom et al. ()where ELA is the elevation of the ELA in meters above sea level, P is the annual rainfall in millimeters, T is the mean annual temperature in degrees Celsius, LR is the atmospheric lapse rate in degrees Celsius per meter, and z is the altitude of the temperature measurement in meters above sea level. Present-day climate variables for the nine glacial sites were derived from meteorological station data (). We derived identical results (precipitation map and uniform cooling) when replacing the model of Condom et al. () by the one of Greene et al. () in our algorithm. Those kinds of empirical relationships are more efficient to link the ELA and the climate than the PDD (positive degree day) for the Altiplano where sublimation is important and where precipitation can be low (section S2.3.1).

Lake modeling

The hydrological budget of Lake Tauca was modeled using a distributed model derived from that of Condom et al. (). The model is run at a quarterly time resolution and at a spatial resolution of 5 km. Temperature and precipitation inputs were derived from the New et al. () data set. Elevation was based on the SRTM 1 Arc-Second Global digital elevation model from NASA-USGS. The model parameters were calibrated on the Titicaca watershed for the present period.

Paleoclimatic inversion algorithm

To determine a precipitation field for the Tauca highstand, we applied the following algorithm (fig. S16). First, we assumed a uniform ΔT cooling over the whole basin and applied this ΔT cooling at each glacial site. Second, using the present-day annual rainfall, mean annual temperature, the Tauca ELA, and the chosen cooling, we were able to compute the paleomean annual rainfall at each site using Eq. 1. From these nine rainfall values, we then interpolated a precipitation grid for the whole Altiplano. Third, using this interpolated precipitation grid and the actual mean annual temperature () modified by the uniform ΔT, we computed the hydrological budget of the basin with the modified version of the model by Condom et al. (). Finally, depending on the final value obtained, we modified the initial ΔT cooling and then repeated the same calculations until the hydrologic lake balance reached a null value (characteristic of the Tauca highstand). This condition allows identification of a precipitation field that is representative of the Lake Tauca highstand and that satisfies the concurrent Tauca glacial extents at each site. We based the spatial interpolation of precipitation on the square inverse of the distance. We used a Monte Carlo method to propagate the uncertainties attached to the paleo-ELAs and to present-day climatic conditions (precipitation and temperature) in the paleoprecipitation grid.