High- and low-latitude forcings drive Atacama Desert rainfall variations over the past 16,000 years.

Late Quaternary precipitation dynamics in the central Andes have been linked to both high- and low-latitude atmospheric teleconnections. We use present-day relationships between fecal pellet diameters from ashy chinchilla rats (Abrocoma cinerea) and mean annual rainfall to reconstruct the timing and magnitude of pluvials (wet episodes) spanning the past 16,000 years in the Atacama Desert based on 81 14C-dated A. cinerea paleomiddens. A transient climate simulation shows that pluvials identified at 15.9 to 14.8, 13.0 to 8.6, and 8.1 to 7.6 ka B.P. can be linked to North Atlantic (high-latitude) forcing (e.g., Heinrich Stadial 1, Younger Dryas, and Bond cold events). Holocene pluvials at 5.0 to 4.6, 3.2 to 2.1, and 1.4 to 0.7 ka B.P. are not simulated, implying low-latitude internal variability forcing (i.e., ENSO regime shifts). These results help constrain future central Andean hydroclimatic variability and hold promise for reconstructing past climates from rodent middens in desert ecosystems worldwide.

INTRODUCTION

The dry central Andes in South America (15° to 27°S), encompassing the semiarid Altiplano and high Andes >~3600 m and the adjacent hyperarid Atacama Desert, are potentially vulnerable to future warming. Regional climate simulations under extreme greenhouse warming scenarios indicate a reduction in summer rainfall by 30% over the coming century (, ), bringing about severe impacts on ecosystems and human societies (, ). The central Andes, however, are remote and topographically complex with only sparse weather observations, making it difficult to discriminate different climate outcomes from independent lines of evidence. Precisely identifying episodes of past rainfall variability and their myriad drivers could improve forecasts of future hydroclimatic variability, as well as aid the understanding of past human-environment interactions in this region. Paleoclimatic records spanning the past 50 ka B.P. (thousands of years before present) have shown that major hydroclimatic variations occurred during the late Quaternary (–). The timing of these changes, however, suggests a complex set of forcing mechanisms that could be responding to regional (, ) and/or interhemispheric drivers ().

The Atacama Desert and Pacific slope of the Andes have been particularly sensitive to past rainfall variability, which left an indelible footprint on the coupled environmental and cultural history of the Atacama (). Permanent arid to hyperarid conditions in this region have been interrupted in the past by the incursion of wet events that are widely replicated in multiple late Quaternary paleoecological (, ) and paleohydrological (, , ) records. The largest of these events over the past 100 ka B.P., in terms of both magnitude and spatial extent, was the Central Andean Pluvial Event (CAPE) originally constrained from >15.9 to 9.7 ka B.P. () and later expanded to 18 to 9 ka B.P. based on additional records (, ). As currently defined, CAPE consists of two extensive pluvials (CAPE I and CAPE II) interrupted by a dry phase between ca. 14 and 13 ka B.P. (). Increased precipitation during the two CAPE pluvials has been attributed either to increased local summer insolation () or to permanent or semipermanent “La Niña–like” states (, ) driven by anomalously cool tropical Pacific sea surface temperatures (SSTs) and are referred to here as low-latitude forcing. Conversely, prolonged dry periods are thought to have resulted from El Niño–like conditions. More recent evidence suggests that changes in North Atlantic SSTs (i.e., high-latitude forcing) during Heinrich Stadial 1 (HS1) could have played a much larger role during CAPE, likely through southward displacements and intensification of the South American summer monsoon (SASM) (, ). Hence, periods of decreased SST in the North Atlantic throughout the Holocene, such as during Bond events (), could have affected precipitation regimes in the central Andes.

Hydroclimatic variability throughout the Holocene following CAPE has been largely attributed to low-latitude forcing. During the Middle Holocene (~8 to 4 ka B.P.), for instance, most records indicate prevailing hyperarid conditions that had major impacts on human settlements (, ). El Niño–Southern Oscillation (ENSO) activity was severely limited between 5 and 4 ka B.P., most likely driven by orbital changes (–). Major ENSO regime shifts have also been linked to hydroclimatic variability during the medieval climate anomaly (MCA) and the Little Ice Age (LIA) between 1.0 and 0.8 ka B.P. and 0.75 and 0.55 ka B.P., respectively (). Discrepancies exist, however, on whether El Niño or La Niña prevailed during MCA with the opposite ENSO stage following afterward during LIA (–). Changes in SASM intensity have also been suggested for these same periods in the Bolivian and Peruvian Andes (, ). More recent “historic” pluvials during the early 19th and mid-20th centuries described in records from the Altiplano and the Atacama have been attributed to dominant La Niña conditions ().

Modern precipitation in the central Andes occurs predominantly (>80%) during the austral summer, with two different upper air forcing mechanisms causing an anti-phased pattern in the northern and southern extents of this region (). Precipitation in the north (15° to 23°S) results from advection of tropical moisture from the Amazon basin and the tropical Atlantic (), and its interannual variability depends on both the intensity of the SASM () and SST changes in the tropical Pacific (). Further south (24° to 27°S), convective activity is modulated by the strength and location of the Bolivian high (), which relies on the moisture content of lowlands east of the Andes (i.e., the Gran Chaco) (). Moisture in the Gran Chaco depends not only on the intensity of the SASM and the transport of warm, moist air from the Amazon basin by the South American low-level jet but also on extratropical cold air incursions (). Intensification of the SASM is therefore crucial for rainfall episodes in the central Andes, yet its contribution to past pluvials has received little attention.

Here, we present a new proxy record of rainfall variability from the central Atacama Desert (22° to 23°S; Fig. 1A) that spans the past 16,000 years. Our record is based on a new quantitative method that estimates past changes in rainfall using a correlation between mean annual rainfall (MAR) and fecal pellet average diameters of present-day ashy chinchilla rats (Abrocoma cinerea) in northern Chile. We then apply the inverse correlation as a transfer function to A. cinerea pellets obtained from a series of radiocarbon-dated fossil middens (or paleomiddens) to develop a record of the magnitude and temporal extension of past rainfall episodes. Paleomiddens are amalgamations of plant and animal debris encased in crystallized urine (amberat), which enhances their preservation for tens of thousands of years in caves and rock shelters of arid environments and have been studied extensively to infer past climates in arid regions of western North and South America (). To establish the extent of either high- or low-latitude forcing of past hydroclimatic variability, we compare our results to other regional and extraregional paleoclimate records and to a transient climate simulation of rainfall in the Atacama Desert that spans the past 22,000 years, which is mostly driven by high-latitude forcing ().

Fig. 1.

Maps of the study area showing precipitation patterns, sites of modern middens and paleomiddens, and the linear regression used to estimate past annual rainfall.

(A) Map of study area showing major trends in mean annual rainfall (MAR), which positively correlates with elevation and negatively correlates with latitude. Modern midden sites (red stars: QLH, Quebrada La Higuera; SHO, Salar del Huasco; CDT, Cordón de Tuina) and DGA meteorological stations (white triangles) are also shown, as well as the current distributional range of A. cinerea (hatched shading). Map was created using QGIS 3.16 software using the SRTM30 elevation model (source: U.S. Geological Survey) and monthly precipitation data spanning 50 years (1950–2000) [source: ()]. (B) Satellite image of the eastern slope of the Calama basin in the central Atacama Desert showing paleomidden sites (black circles) along with isohyets. CDM, Cerros de Minta; CDA, Cerros de Aiquina; CDT, Cordón de Tuina; ACO, Arroyo Coya; EHT, El Hotel; BSA, Barros Arana; PVI, Pampa Vizchachilla (table S4). (C) Scatterplot of modern A. cinerea fecal pellet diameters and MAR showing the linear model used for estimating MAR in paleomiddens (MAR = 67.62*[Paleomidden N20C] – 148.63; r2 = 0.80; P < 0.001; prediction accuracy = 0.97). Shaded area corresponds to the 95% CI of the model, and error bars represent the 95% standard error CI.

Fossil fecal pellets and past climate change

Why would fecal pellet size be related to climate? Changes in body size profoundly affect the life history, ecology, and evolution of organisms and are a heritable and easily measured morphological trait in animals that scales allometrically with fecundity, energetic requirements, diet, territory, and home range size (). Identification of the environmental factors that drive long-term (i.e., at millennial scales) body size changes is therefore of particular interest for the ecology and evolution of species. Fecal pellet size (diameter) is positively correlated with body size in modern pack rats (Neotoma cinerea), and negatively with ambient temperature, and has been used to estimate body size variations in pack rats in response to temperature variations during the late Quaternary as a test of Bergmann’s rule (, ). In other arid lands, however, modern studies indicate that the relationship between mammal body size and environmental factors is variable and that precipitation (through net primary production) rather than temperature may be the main limiting factor ().

In the hyperarid Atacama Desert (21° to 24°S), paleomiddens built by ashy chinchilla rats (A. cinerea) are identified by the size and shape of their copiously abundant fecal pellets (fig. S1). At present, this species occurs at elevations >3500 m, where MAR and mean annual temperature (MAT) are ca. 60 to 70 mm/year and 5° to 6°C, respectively (Fig. 1A) (). Several late Quaternary radiocarbon-dated middens record downslope displacements of this species, relative to its modern range into what is currently known as “absolute desert” between 2400 and 3100 m (Fig. 1B) (). As with Neotoma middens, fecal pellets from Abrocoma middens also exhibit large variations in size and have been linked to past changes in precipitation (). Contrary to North American pack rats, however, estimating body size changes in A. cinerea from fecal pellet diameters is hindered because individuals of this species are very shy, are difficult to capture and observe in the field, and do poorly in laboratory experiments that control for diet and ambient conditions (temperature and relative humidity). On the basis of similar studies with Neotoma, we assume that a relation exists between fecal pellet size and body size in A. cinerea and use this theoretical link to test for modern correlations with temperature and/or precipitation along an environmental gradient in northern Chile between 18° and 24°S (Fig. 1A and see the “Modern calibration and past estimations of MAR” section in Materials and Methods).

RESULTS

Past rainfall anomalies over the past 16,000 years in the central Atacama

Our record (Fig. 2) shows that persistently very wet and almost uniform conditions prevailed from 15.9 to 14.8 ka B.P. (e.g., during CAPE I) followed by a conspicuous absence of paleomiddens between 14.8 and 13.1 ka B.P. This conspicuous lack of middens could be due to several causes, including a very dry period, which would have forced midden building rodents upslope from our midden sites, but we remain cautious with our interpretations until more evidence is available (see also the “Criteria for defining the timing of wet and arid phases” section in Materials and Methods). In contrast, rainfall was much more variable from 13.0 to 8.6 ka B.P. during CAPE II (Fig. 2), characterized by large centennial-scale fluctuations of very wet (from 100 to 60 mm/year, the largest anomalies in our record) to mildly wet or dry periods (60 to 10 mm/year) until 8.6 ka B.P. Aridity increased gradually from a very wet period at 9.3 ka B.P. to near-modern conditions at 8.4 ka B.P. A wet phase ensued at 8.1 ka B.P., reaching a maximum of ~88 ± 11 mm/year around 7.5 ka B.P. This was followed by a ~2500-year gap in our record bounded by similar-to-present and positive rainfall anomalies at 7.4 and 5.0 ka B.P., respectively. As with the previous gap, one possible explanation for the lack of middens from 7.4 to 5.0 ka B.P. is the prevalence of extremely arid conditions.

Fig. 2.

Temporal variation of A. cinerea pellet diameters (top) and estimated rainfall anomalies (bottom) in the central Atacama.

(Top) Plot of corrected pellet diameter from paleomiddens over the past 16,000 years. Error bars represent the 95% standard error CI, and horizontal bars represent calibrated age errors at two-sigma. (Bottom) Plot of MAR anomalies estimated from corrected pellet diameters from paleomiddens (top) using the linear model shown in Fig. 1C. Hydroclimatic phases for the Atacama are indicated in blue (wet) or red (dry) shading (see Results and Discussion). Major climate events relevant for the Atacama as established by previous research (8 and 13 for CAPE and 21 for MCA and LIA) are shown for reference. Quaternary epochs and stages follow the International Stratigraphic Chart v.2019 (http://stratigraphy.org).

A brief wet episode occurred at 5.0 ka B.P. followed by increasingly large wet anomalies during the Late Holocene. Sustained but decreasing rainfall from 2.5 to 1.7 ka B.P. was followed by increased rainfall variability from 1.4 to 0.76 ka B.P. with large wet anomalies (~83 ± 11 mm/year) at 0.86 ka B.P., coeval with MCA (between 1.0 and 0.7 ka B.P.) (). Some of the largest changes in rainfall anomalies occur in paleomiddens with almost identical calibrated 14C dates. For example, an 83 ± 11 mm/year anomaly occurs at 0.862 ka B.P. (960 ± 15 14C yr B.P.), followed by an ~11 ± 20 mm/year anomaly from another midden dated to 0.857 ka B.P. (1015 ± 15 14C yr B.P.). One possible explanation is that these changes could be due to large decadal-scale variability in precipitation, which would not be resolvable by direct accelerator mass spectrometry (AMS) dating of individual middens (table S3 and fig. S5). Four middens dated to between 0.5 and 0.1 ka B.P. [e.g., during LIA; (21)] have anomalies indistinguishable from present-day values, except during the past 200 years when a slight increase (45 to 50 mm/year anomaly) in rainfall estimated from two middens matches a regional recent historical pluvial event in the early 1800s (, ).

We used a global transient climate simulation () of the past 22,000 years to further analyze rainfall in the Atacama Desert (centered at 24°S) (Fig. 3I). This simulation shows two major pluvials with >20% increase in rainfall during the late Pleistocene, from 22 to 14.5 ka B.P. and from 13.8 to 11.4 ka B.P. The simulation further shows a precipitation minimum from 14.3 to 14.0 ka B.P. After 11.3 ka B.P., precipitation rises almost linearly throughout the Holocene, from an average deficit of −15% to modern values. A less intense pluvial interrupts this secular trend between 8.9 and 7.6 ka B.P., although average rainfall was still below the present-day average.

Fig. 3.

Comparison of regional and extra-regional paleoclimate records.

(A) Mid-month January (summer) insolation 15°S (). (B) δ18O from GISP2 Greenland ice core (). (C) CH4 ice core record from EPICA Dome C Antarctica (). (D) δ18O speleothem record from Botuverá cave in southeast Brazil (). (E) Lithics concentration in a marine sediment core offshore central Peru (). (F) Modeled precipitation of the southern Bolivian Altiplano (). (G) MAR anomalies in the central Atacama Desert (this study). (H) Transient simulation of precipitation (TraCE 21k) over the central Atacama (24°S) plotted as 100-year running average anomalies from present-day climate (this study). Bond events 0 to 8 and the 8.2k event are shown by vertical gray shading (). Timing of major regional climatic events follows 8 and 13 for CAPE and 21 for MCA and LIA. Quaternary epoch stages are after the International Stratigraphic Chart v.2019 (http://stratigraphy.org). p.p.b.v., parts per billion by volume; p.p.m.v., parts per million by volume.

DISCUSSION

Our evidence for extensive pluvials at 15.9 to 14.8 and 13.0 to 8.6 ka B.P. in the central Atacama (15° to 27°S) further constrains the temporal extension of CAPE. Our improved chronology also shows that pluvials during the Early Holocene were possibly absent between 7.5 and 5 ka B.P. and resumed during the Middle to Late Holocene. These younger pluvials where much shorter in duration (<1 ka) that gradually started at ~5 ka B.P. and intensified at 2.5 and 1.1 ka B.P. To what extent were variations in present-day climate mechanisms (e.g., ENSO) responsible for causing these regional changes versus changes in summer insolation? To what extent is high-latitude forcing responsible for these pluvials?

Local changes in summer insolation (the amount of solar radiation at the top of atmosphere, typically at 20°S for the central Andes) were among the first mechanisms proposed to explain the formation of paleolakes in the central Andes (, ). The timing of the pluvial events seen in our record, however, shows no support for local orbital insolation forcing as these occur during an insolation minimum (at 11 ka B.P.) (Fig. 3A). Late Holocene Atacama pluvials, however, increased in magnitude over time, similar to other paleo-ENSO records (Fig. 3, E and F) (, ), pointing to a possible role for increasing insolation through its impact on ENSO ().

Until recently, past changes in rainfall described in different paleorecords from the Atacama were attributed to permanent or semipermanent La Niña–like states to explain large increases in precipitation during the late Quaternary (, , ). Wet events (or “pluvials”) occurring during these periods of low-latitude forcing should thus be characterized by “cold phases” in the tropical Pacific and strengthened west-east SST gradients. Our timing of the 15.9 to 14.8 ka B.P. pluvial (CAPE I) does occur during La Niña–like conditions (). Yet, our record shows that this was not the case during the second part of CAPE II (11 to 9.6 ka B.P.) or at 7.5 ka B.P. Additional explanations are necessary apart from this low-latitude forcing.

Recent studies have examined how changes in North Atlantic SSTs affect rainfall variability over the central Andes (, , , ). A correlation between North Atlantic cooling and precipitation increases over the central Andes during the Holocene has been previously suggested (). Two major mechanisms are likely involved in regulating interhemispheric temperature contrasts in the Atlantic and its subsequent impacts on central Andean climate: (i) a southward displacement of the intertropical convergence zone (ITCZ) and (ii) enhanced moisture transport from the tropical Atlantic (). These changes lead to a southward displacement and intensification of the SASM, thus increasing available moisture in the Amazon basin, leading to increased easterly advection over the central Andes (). Several paleorecords of SASM intensity have been correlated with North Atlantic SSTs over the past decade (, ) and link massive ice discharges in the North Atlantic to corresponding increases in rainfall over the central Andes (, ). An intensified and southerly displaced Bolivian high resulting from a southerly shifted SASM could partly explain the intensity and magnitude of CAPE I ().

We support the hypothesis that CAPE likely resulted from an intensified SASM triggered by the cooling of the North Atlantic by massive ice discharges during HS1 (between 16 and 14 ka B.P.) and the Younger Dryas event (YD; between 12.8 and 11 ka B.P.). This, in turn, drove an intensification and southward shift of the SASM (Fig. 3, B and D) (, –). These two periods of increased and widespread SASM activity were partly coeval with the two CAPE phases: CAPE I during HS1 and the first half of CAPE II during the YD (Fig. 3G). Furthermore, records of equatorial Pacific SSTs suggest La Niña–like conditions (, ), which would have even further increased precipitation over the central Andes. The Atacama rainfall transient climate simulation also suggests cold North Atlantic forcing during this interval (Fig. 3I). One key difference, however, is that ice discharge was far less intense during the YD compared to HS1 (), which would have resulted in a less-intensified SASM during CAPE II (). This difference in SASM intensity appears to have affected most of the central Andes, with deep paleolakes in the Altiplano and increased, sustained rainfall over the Atacama during CAPE I and shallow paleolakes and increased but variable rainfall during CAPE II (Fig. 3, G to I).

Precipitation increases appear to have been more intense in the southern Altiplano during CAPE II, which has been attributed to changes in moisture source (, , , ). If CAPE II involved an intensification and southward displacement of the Bolivian high, as with CAPE I, then this would shift precipitation maxima southward with moisture transport from the central Brazil–Gran Chaco region (). This implies that the differences observed between the CAPE phases likely resulted from weakened and/or more variable monsoonal activity during CAPE II compared to CAPE I and not from changes in moisture source, which is consistent with the transient climate simulation (Fig. 3I) and our MAR anomaly record (Fig. 3H).

Our results imply that hydroclimatic variability in the Atacama appeared to switch from high- to low-latitude forcing during the Middle Holocene. The positive anomaly peaks in our record between 11 and 7.5 ka B.P. correlate with SASM intensification during Bond events (, ), suggesting that a high-latitude control remained even after the end of CAPE II (Fig. 3, D and H). This relationship breaks down as the Laurentide ice sheet disappeared around 7 to 7.5 ka B.P. along with possibly reduced ENSO activity between 10 and ka B.P. (Fig. 3, E and F) (, , , ). The Middle to Late Holocene pluvials were markedly similar in terms of intensity and temporal extension to changes in ENSO event frequency and intensity (Fig. 3, E and F), implying a switch to low-latitude forcing sometime between 8 and 5 ka B.P.

Increased orbital insolation since 5 ka B.P. resulted in enhanced ENSO activity through its impact on equatorial Pacific SSTs (Fig. 3, A and F) (, ). SASM intensification between 3 and 2.4 ka B.P., however, could have also been driven by high-latitude forcing () through southward displacement of the ITCZ (). No single driver, however, appears to explain all the rainfall variations over the past few millennia in the central Atacama, and these events are noticeably absent from the transient simulation. Because these modes of variability are internal and intrinsic to the climate system, their exact timing cannot be reproduced by a climate model, which can only respond to a given external forcing or driver. Nevertheless, a stronger low-latitude forcing is more clearly visible during the Common era (i.e., last 2 ka B.P.), as the timing seen in our record (with positive rainfall anomalies during MCA and negative ones at the onset of LIA) corresponds with major ENSO regime shifts. The nature of hydroclimatic variability over the central Andes during MCA and LIA is a topic of ongoing debate, much of it centered on which ENSO phase predominated during each period. Wet conditions during MCA have been attributed to continuous La Niña–like states [with subsequent El Niño dry conditions during LIA; (, , )], whereas the opposite [El Niño followed by La Niña–like conditions, i.e., (, )] has also been suggested.

In the Atacama, evidence points to a pluvial event during MCA (, ), although some recent studies also show evidence for pluvials during LIA (, ). As most records point to a weak SASM during MCA, which later strengthened during LIA (Fig. 3D) (, –), a wet Atacama during MCA would imply a low-latitude driver for hydroclimatic variability that favored moisture advection across the central Andes. Although a recent historic pluvial could also be driven by changes in ENSO regimes (), the timing of these changes indicates a more complex set of forcing mechanisms with differing regional impacts ().

Last, increased moisture availability from 13 to 8.6 ka B.P. in the Atacama Desert may have attracted the first human settlers who were expanding at the time across the different ecosystems of the central Andes (). Increased environmental stress likely drove human societies to abandon certain territories while dispersing to others from 8 to 5 ka B.P., where they developed different socioecological strategies including camelid domestication (). Although shorter in duration, Late Holocene pluvial events (<1 ka B.P.) had major impacts on human societies that thrived in the Atacama hyperarid core (, ). Major technological improvements for coastal and inland food procurement and processing, coupled with major sociocultural changes including investments in permanent hamlets with solid and multifunctional architectures, irrigated and fertilized desert farming, major population growth, and increased social complexity and inequalities, were all accelerated during these Late Holocene pluvials (, , ).

Hydroclimatic changes in the central Andes during the Holocene have been a matter of debate for decades, and regional differences between the western Andes and the Bolivian and Peruvian Altiplano have been suggested (, , ). Putting aside complexities in dating and interpretations of these different proxies, unraveling hydroclimatic variability in the region requires understanding the relative contributions of high- and low-latitude forcing mechanisms during the Holocene. Depending on time scales, we propose that changes in North Atlantic SSTs were as relevant as those occurring over the tropical Pacific in controlling rainfall variability over the central Andes (Fig. 3), and neither driver can independently explain our entire record. For instance, Holocene rainfall variability in the Atacama resulting from SASM modifications alone would generate pluvials that would only track austral summer insolation as shown by the transient climate simulation, instead of becoming less frequent and more intense as seen in our record (Fig. 3, D and I).

We link episodes of massive ice discharge during HS1 and YD, and brief discharges during Bond events until 7 ka B.P. to (i) increased Atlantic trade winds and southward displacement of the ITCZ, (ii) intensification and southward displacements of the SASM, and (iii) pluvials in the Atacama Desert. This implies that these pluvials likely had a North Atlantic driver from 16 to 7 ka B.P. instead of (or complimentary to) the La Niña–like states, which had been previously proposed for CAPE. The alteration of the Bolivian high as a response to these atmospheric changes during CAPE I likely caused the most intense pluvial (CAPE I) in the central Andes over the past 100 ka, providing the best conditions for inland late Pleistocene human societies. It is reasonable to assume that similar mechanisms operated during the YD and subsequent Bond events throughout the Holocene. Although changes in equatorial Pacific SSTs could have additive effects coupled with episodes of North Atlantic cooling, our record suggests that ENSO regime drivers became the dominant control of precipitation only since the Middle Holocene, which is likely linked to increased solar insolation and its influence on the tropical Pacific.

Paleomiddens are used widely to reconstruct past climate and ecological change throughout the world’s arid regions, but the discontinuous nature of the midden record has hampered comparisons with other continuous proxies. By pooling several paleomidden series from nearby localities (such as throughout a basin) covering a large geographical area and time scale, our study overcomes this restriction by providing an increased number of samples available for replication. We show that rodent fecal pellets from ancient middens can be used to obtain quantitative estimates of past climate change. Here, we identify the timing and magnitude of multiple wet and dry episodes in the Atacama Desert over the past 16,000 years and propose different mechanisms of climatic forcing. We note that the methods used could be readily exported to other arid regions of South America where chinchilla rat middens have been collected, which could provide further examination of the timing of events observed here.

Our results here highlight the complexity of the controls of climate variability over the central Andes. Although our interpretations are simplified, they are essential for assessing whether models of future climate change can correctly estimate the full range of future regional hydroclimatic responses. These will be critical for regional governance and development of sustainable policies for the decades to come.

MATERIALS AND METHODS

Midden collection and processing

We selected and analyzed a total of 81 previously collected chinchilla rat (A. cinerea) middens (, , ) in northern Chile (Fig. 1, A and B, and table S1). Midden agents (i.e., the species of rodents that build middens) were identified on the basis of fecal pellet size and shape (fig. S1) (). Following standard procedures (), indurated middens were placed in 10-liter buckets of water for 5 to 7 days to remove the crystallized urine (amberat) and then wet-sieved and dried at 50° to 60°C in a drying oven. Fecal pellets and other organic debris (plant macrofossils, vertebrate bones, and insect remains) were sorted and stored for various research purposes. Radiocarbon dating was performed on all middens by either conventional methods or AMS using between 0.1 and 0.3 g for AMS or 3 and 10 g for conventional dating of fecal pellets (table S1). Midden 14C ages were calibrated at two-sigma with the Southern Hemisphere curve () using the CALIB 7.0.2 program () to thousands of calendar years before 1950 (ka B.P.) to facilitate comparisons across records. Out of the total set of available middens, 13 were used to establish correlations with modern climate variables (table S2), and 68 middens were used for our MAR anomaly reconstruction (table S3). To ensure that local rainfall in the vicinity of the middens was the only source of water, the middens used in this study (both modern and fossil) were collected away from water sources (wetlands, point springs, etc.).

Pellet measurements

After separation from other materials, the width of the largest ~200 pellets (or ~20% of those available, whichever was greater) was measured to the nearest 0.1 mm using a digital caliper (Mitutoyo, USA) following protocols previously described for wood rat (N. cinerea) middens in North America (, ). Statistics on the largest pellets were used instead of total measurement statistics as this allows for discrimination of the confounding effect of juveniles (). Measuring 10 to 20% of the sample is sufficient to characterize the largest pellets; above this maximum threshold, pellet size does not vary with sample size (). When fewer than 200 pellets were present, all were measured. The mean, SD, SE, and 95% confidence interval (CI) of all the measurements (NTOT) and the largest 20% of the pellets measured (N20) were calculated.

Modern calibration and past estimations of MAR

Modern middens were defined here as those directly dated with 14C resulting in ages <150 cal (calibrated) yr B.P. or containing “bomb radiocarbon.” These middens were mostly collected and archived during fossil midden surveys carried out between 1996 and 2012 covering an extensive latitudinal range (18° to 26°S). During these collection campaigns, approximately 1600 middens have been recovered and are currently stored at the Laboratory of Paleoecology and Paleoenvironments of the Pontificia Universidad Católica de Chile. Application of criteria, such as sufficient number of measurements (>200 largest pellets or ~20% of those available), discarding middens collected in proximity to nearby water sources other than rainfall, and limiting the dataset to those collected within the latitudinal range covered by weather stations (see below), results in only 13 middens that could be used to establish the relationship between modern MAR and pellet diameter values.

Mean daily temperatures and monthly accumulated rainfall data in northern Chile (17 to 30°S) were retrieved from weather stations of the DGA (Dirección General de Aguas) available at the Center for Climate and Resilience Research website (MAT: http://cr2.cl/datos-de-temperatura/; MAR: http://cr2.cl/datos-de-precipitacion/). We included only those above the inversion layer boundary at ca. 800 m above sea level, spanning the longest continuous period possible without any gaps (table S4). Because of differences in collection of rainfall and temperature data between stations, we applied different filters for each dataset. For temperature, we used stations spanning a 13-year period (1976–1989) between 17° and 28°S with more than 76% of available data, except for stations located where gaps in elevation or latitude existed. Stations with less than 300 days of data information were also discarded. In the case of rainfall data, we used stations covering a 30-year period (1987–2016), including stations with more than 98% of data between 17° and 24°S. We considered water years from 1 November to 31 October as this groups both summer and the following winter rainfall peaks into the same year. After filtering both datasets, we performed analyses of variance (ANOVAs) to identify whether elevation or latitude better explained changes in temperature or rainfall.

We found that MAT varies inversely with elevation but very little with latitude between 17° and 25°S (fig. S2A), in agreement with previous research (). MAR presents a more complex relationship varying in both latitude and elevation (fig. S2B), decreasing toward the Pacific coast and southward until ~24°S, the southernmost region where summer rainfall regularly occurs (). Linear regressions to estimate temperature were then established on a yearly basis as a function of elevation (MAT = a – b*Elevation) and of latitude and elevation to estimate rainfall (MAR = a*Latitude + b*Elevation + c). These annual equations captured between 80 and 94% and 57 and 96% of the variance for MAT and MAR, respectively. We thus derived average MAT and MAR values for each modern midden site using latitude and elevation information (fig. S2, C and D). Estimates of means and error values of MAT and MAR for each modern midden site are included in table S2. N20 was measured in each A. cinerea modern midden and then plotted against MAR and MAT estimates for each site (fig. S3), which then enables us to link changes in pellet diameters as a function of these climatic variables (Eqs. 1 and 2).

Simple linear regression models show that MAR (Eq. 1) is a better predictor of changes in N20 pellet diameters than MAT (Eq. 2) in the Atacama Desert (fig. S3). We then use the inverse correlation as a transfer function to calculate past quantitative MAR estimates from fecal pellet measurements obtained from our paleomiddens (Fig. 1C)

Estimation of past rainfall anomalies

Fecal pellet measurements from 68 radiocarbon-dated paleomiddens collected from seven localities in the central Atacama Desert were used to estimate past rainfall anomalies (estimated past MAR minus present-day MAR at the same site; Fig. 1A and tables S1 and S3). Pellet diameters slightly decrease linearly with age, and a correction was applied to all N20 average measurements from paleomiddens (fig. S4). Accuracy of MAR estimates was assessed by calculating the difference in anomaly estimates between same-age midden pairs (within <10 14C yr difference of each other) collected from multiple sites and/or elevations (fig. S5). As most paleomiddens were collected along the hyperarid margin of the Atacama where no active A. cinerea populations exist today, a positive bias exists in our record implying that it may be better at tracking wet phases rather than dry ones (see the next section).

For estimation of rainfall anomalies, diameters from 16,071 fecal pellets were measured from 68 individually 14C-dated rodent middens (i.e., those older than 0.15 ka B.P.) collected from seven localities in the central Atacama Desert (Fig. 1B and table S3). Because pellet diameters slightly decrease linearly with age (Eq. 3 and fig. S4), a correction was first applied to all N20 average values for obtaining quantitative estimates of MAR in each midden. This decrease appears to be linear at least for the past 16 ka B.P. but may level off for older middens. This will need to be further explored with a larger dataset of much older paleomiddens

After correcting for changes in diameter with time, we then established a simple linear regression model based on the inverse relationship between modern fecal pellet diameters and MAR, which permits the prediction of MAR estimates using corrected paleomidden N20 mean values as a predictor variable (Eq. 4 and Fig. 1C)

Estimates of MAR values and 95% confidence and prediction intervals were thus obtained from each midden, and rainfall anomalies were calculated using modern MAR estimates for each paleomidden site (table S3). Replicability of MAR estimates was assessed by calculating the difference in anomaly estimates between same-age midden pairs (within <10 14C yr difference of each other) collected from multiple sites and/or elevations (fig. S5).

The use of a linear regression model was confirmed after checking for parametric assumptions through analyses of residuals, and additional validation was obtained in R () using the gvlma package (). Diagnostic evaluation of the regression model was performed by splitting the modern MAR and pellet diameter dataset into an 80:20 (training:test) subset for comparison of observed and predicted MAR values. These analyses suggested a min-max accuracy of 0.97 and a mean absolute percentage deviation of 0.029. Although a log function provides a slightly better fit [Akaike information criterion (AIC) = 109.2314] compared to a linear function (AIC = 110.9961), we opted for the linear function as (i) visual inspection of the log curve data remains linear within the range of values covered by the dataset with no indications of reaching a plateau, and (ii) a linear function provides a simpler way to interpret MAR estimates based on pellet diameter without data transformation. In addition, MAR estimates using both functions do not differ significantly, and therefore, the use of a linear model is appropriate under the dataset used. We do call for caution, however, as pellet diameters (and therefore body size) are restricted to increase indefinitely with increasing rainfall, and thus, a log model would be more biologically reasonable. Testing for this prediction should be explored in further studies using a wider range of fecal pellet diameters along a wider latitudinal gradient and/or from other regions of arid South America where A. cinerea occurs.

Criteria for defining the timing of wet and arid phases

Wet episodes are defined when MAR anomalies remained above our level of uncertainty (conservatively established at >40 mm/year; see fig. S5). These wet events differ significantly from the modern baseline (see previous section). Arid to hyperarid phases were established on the basis of either (i) low MAR anomaly estimates (<40 mm/year) or (ii) tentatively by the absence of any middens for prolonged periods of time (i.e., >1000 years; see Fig. 2). The reasoning for this is that many of our paleomiddens were originally collected with the purpose of establishing major vegetation changes along the hyperarid upper margin of the absolute desert (Fig. 1B), an area sparsely vegetated today with very low bioproductivity and no active A. cinerea populations (, ). As these areas would have been too dry for middens to form in the past during extended droughts, this could likely create the “gaps” in our record. Hence, a lack of middens for prolonged periods of time could be a response to very arid conditions. This also implies that our record is “positively biased” in that it does a much better job at tracking wet phases than dry ones. A possible test would be to develop a similar record at higher elevations (>3500 m) where A. cinerea populations remain active today.

Transient simulation of Atacama Desert rainfall

We used the TraCE 21k simulation () to analyze rainfall in the Atacama Desert over the past 16,000 years to enable comparison of our results to other regional paleoclimatic records. TraCE is a simulation of the transient (i.e., continuous) climate of the past 21,000 years. The transient climate simulation performed was centered at 24°S using a fully coupled model (National Center for Atmospheric Research Community Climate System Model version 3) that is forced by time-varying orbital parameters, greenhouse gases, time-varying ice-sheet extent and topography, and freshwater forcing to the oceans from the retreating ice sheets (). The model thus simulates local impacts of major high-latitude drivers such as the HS1 and YD stadials as well as the Antarctic Cold Reversal.