Monte Sierra Nevada stands at 2,527 m above sea level (a.s.l.), with ice capping the highest peaks and extending partially down several east-facing cirque valleys (Figure 2). The site is located to the west of the Cordillera Principal, which is made up of the granitoid North Patagonia Batholith (The Cura-Mallin Formation) of Cretaceous and Miocene age (Glodny et al., 2006). Pleistocene and Holocene volcanoes (e.g., Llaima, Lonquimay, Sierra Nevada, and Tolhuaca) make up the active volcanic arc of the Cordillera Principal that formed on the eroded, and relatively flat-lying Miocene-aged rocks (Muñoz et al., 2011). Regionally, this southern volcanic zone is predominantly composed of basalts and basaltic andesites (Hickey et al., 1986; Ferguson et al., 1992; Muñoz et al., 2011). Additionally, these active Holocene volcanos might have deposited volcanic ash on glaciated terrain, potentially affecting glaciers’ response to climate changes by two competing effects associated with debris thickness and temperature: insulation from incoming energy or enhancement of ablation (Davies and Glasser, 2012).

In our study area, there are a total of five morainal loops along the valley floor over a distance of ∼350 m (Figure 3). Mapped moraines exhibit variable effects of weathering and erosion suggesting relative age differences. We number each moraine, starting with the oldest and most distal from the current glacier front (M1), progressing westward to the youngest, most ice-proximal moraine (M5) (Figure 3). Visually, moraine M2 is the largest and most intact of the moraine complex, while moraine M1 is the sole moraine in the valley with observed soil development, and exhibits a subtle dual-crest to the south of the field site. On the ice-proximal side of moraine M5 chaotic hummocks are observed, with glacial scree deposits further up valley (Figure 4). Additionally, several striated boulders and roches moutonnées are observed, both within and outside of the moraine complex.

The geomorphology of Monte Sierra Nevada yields information on the relative timing of glacial landform development due to the stratigraphy observed in the valley. For example, the distal end of the valley is truncated by a steep cliff at the contact between the relatively flat-lying Pleistocene volcanic layers and the underlying Miocene plutonic rocks (Figure 4). The steep topography below this contact likely precludes the establishment of moraines from glacial advances that were more extensive than the presumed Late Holocene moraines that end just before a steep drop in topography to the valley below. On Monte Sierra Nevada, moraines are generally limited to elevations between ∼1,703 and 1,722 m a.s.l (Table 2) and the moraines have since been bifurcated by Holocene stream development.

Directly south of the valley with the mapped moraine complex is a ridge separating it from the adjacent, topographically higher valley (Figure 3). On the ridge separating the two valleys is a well-defined saddle and chute, both of which are currently ice free. Along this same ridge is a trimline indicating the most recent highest extent of ice thickness in the valley. This trimline extends down valley intersecting moraine M1, as well as up valley past the saddle and chute towards the mountain’s peak. The chute exhibits a cross-cutting relationship with the trimline, having eroded any recognizable lateral features. Notable talus deposits are identifiable along the northern and southern valley walls, particularly below the trimline and along the chute. The modern glacial ice cap on Monte Sierra Nevada extends into the uppermost reaches of both valleys.

Central Chile has a Mediterranean-type climate with mild, wet winters and dry, hot summers. The modern climate is largely modulated by Pacific and Southern Ocean processes, resulting in a strong seasonality of precipitation south of ∼35°. winter storms dominate here, bringing ∼2,000 mm of annual rainfall to the region (Sinclair, 1996; Garreaud et al., 2009). However, Monte Sierra Nevada (38°S) is located near what is considered to be a transition zone between Mediterranean-subtropical and temperate-humid climatic zones, identified at ∼37.5°S (Figure 1; Muñoz et al., 2016). Precipitation at this location is therefore impacted by shifts in the subtropical South Pacific Subtropical High (SPSH) and extratropical Southern Westerly Wind Belt (SWWB) year-round (Miller, 1976; Espinoza et al., 2020).

Glaciers in CSCh predominantly receive precipitation sourced from the Pacific Ocean due to the moisture barrier produced by the Andes mountains (Montecinos and Aceituno, 2003). Isotopic analyses of precipitation suggest that Atlantic-derived moisture rarely penetrates through to the western flanks of the Andes (Hoke et al., 2013). While regionally both the height and width of the CSCh Andes are smaller in comparison to the Northern Cordillera, orographic precipitation leads to continental precipitation values that are up to three times larger than corresponding oceanic amounts (Garreaud et al., 2009).

In Central Chile the SPSH typically fluctuates between ∼25 and 35°S and delivers moist air masses to the western slope of the Andes (Ancapichún and Garcés-Vargas, 2015). Active year-round, the position of the SPSH is modulated by the Bolivian High, which migrates northward during the austral winter (May-August) as the jet stream moves equatorward. This increase in winter storm tracks drives midlatitude precipitation along the central Chilean coast (Sinclair, 1996; Lamy et al., 2001; Garreaud et al., 2009; Espinoza, et al., 2020). From February-April the SPSH shifts south, resulting in increased precipitation off the coast of CSCh (∼37°S) (Ancapichún and Garcés-Vargas, 2015). This hydrologic regime results in an enhanced probability of aridity in CSCh, as streamflow peaks in winter due to rainfall and again in late spring as seasonal snow melts in the mountains.

The regional climate in the CSCh Andes is also affected by the Antarctic Circumpolar Current (ACC). The ACC splits as it approaches the western coast of South America between ∼40 and 45°S (Boltovskoy, 1976; Villalba et al., 2012), with the Peru-Chile Current flowing north and the Cape Horn Current flowing south (Lamy et al., 2007). The Peru-Chile Current transports cold, low-saline waters equatorward (Strub, 1998), with seasonal changes in upwelling driven by the position and intensity of the SPSH (Garreaud et al., 2009; Ancapichún and Garcés-Vargas, 2015). The SPSH is thus one of the dominant controls on sea surface temperatures (SSTs) associated with the Peru-Chile Current, with the ability to increase upwelling and offshore Ekman transport. The position of the ACC also results in strong winds directed at the coast, in the form of the SWWB. Situated between the subtropical anticyclone and subpolar cyclonic air masses, these westerlies bring strong winds and abundant precipitation to CSCh and the adjacent northern Patagonian Andes, intersecting the western coast of South America between 40°S-55°S (Bertrand et al., 2014). During Austral summer the SWWB strengthens, contracts, and shifts to the southern extent of its range in response to southern-ocean processes, while during the winter months the SWWB weakens and expands both poleward and towards the equator (Bentley et al., 2009). Comparatively, the SWWB is stronger than the storm tracks associated with the SPSH at the study location, as it is confined to a smaller latitudinal range and more directly affects regional precipitation.

Located in the transition zone between Mediterranean and temperate climates, both ENSO and SAM play important roles in driving changes in temperature and moisture availability to Monte Sierra Nevada. ENSO is characterized by shifting oceanic-atmospheric conditions in the tropical Pacific, the effects of which are superimposed on the climate system of CSCh. El Niño (La Niña) years are associated with a weakening (strengthening) of the Peru-Chile Current, due to a reduction in the east-west pressure gradient (Montecinos and Aceituno, 2003), forcing a weaker SPSH into more northerly position (Villalba, 1994). This northward displacement of the SPSH reduces moisture availability in CSCh. When SSTs are lower the high-pressure system strengthens, and westerly trade winds intensify in response (Bertrand et al., 2014). In the northern portion of the transition zone (35–40°S), ENSO is characterized as having a stronger influence on the SPSH (Rutllant and Fuenzalida, 1991), whereas south of ∼40°S interannual precipitation driven by the SWWB is heavily modulated by the strength/phase of the SAM (Aravena and Luckman, 2009; Garreaud et al., 2009; Villalba et al., 2012).

The SAM is measured by changes in sea-surface pressure between the mid- and high-latitudes (Villalba et al., 2012; Ancapichún and Garcés-Vargas, 2015). Both the strength and position of the SWWB is dictated by SSTs, and sea surface pressure gradients associated with the ACC and the Southern Ocean, specifically the SAM (Gillett et al., 2006; Reynhout et al., 2019). Enhancement of precipitation in Patagonia south of ∼55°S is predominantly driven by a southern displacement of the SWWB, subsequently resulting in decreased precipitation around 38–43°S. During negative (positive) phases of the SAM the SWWB shifts equatorward (poleward), resulting in more humid (arid) conditions in CSCh (Garreaud et al., 2009; Bertrand et al., 2014; Muñoz et al., 2016).

Proxy paleo-climate evidence suggests that Late Holocene glacial variability in the northern Patagonia Ice Field was driven by both precipitation and temperature (Bertrand et al., 2012). A tree-ring reconstruction of summer temperatures from northern Patagonia suggests there was a cold interval from ∼1.1 to 0.95 ka, followed by warming until ∼0.77 ka and another phase of cool and wet conditions peaking between ∼0.68 and 0.38 ka (Villalba, 1994). Proxy paleoclimate records also suggest that from ∼0.65 ka until the last century was the wettest interval of the last ∼5,000 years in the southern Chilean Andes (Moreno et al., 2009). Moreno et al. (2018) strengthened this analysis, and documented that regional glacial advances in the southern Patagonia Ice Field were forced from negative SAM phases (cold and wet), that were modulated on regional scales by the position and intensity of the SWWB. Recently, Davies et al. (2020) compiled geomorphic evidence and recalibrated chronological data from dozens of studies assessing fluctuations of the Patagonian Ice Sheet from 35 ka to the present. Records detailing the Late Holocene from 41°S to 52°S noted climatic conditions alternated between cold/wet and warm/dry on centennial time scales (Elbert et al., 2013; Álvarez et al., 2015; Moreno and Videla, 2016; Moreno et al., 2018; Davies et al., 2020). Well-defined Late Holocene moraines are apparent regionally (Glasser et al., 2008; Glasser and Jansson, 2008; Davies and Glasser, 2012; Davies et al., 2020), but these advances overall do not seem to surpass moraines emplaced at ∼5 ka. While these Late Holocene landforms exhibit clear evidence of later formation, such as vegetation differences and degradation (Davies and Glasser, 2012; Davies et al., 2020), it is possible that positions of earlier advances were reoccupied (Sagredo et al., 2016; Davies et al., 2020).

Additional records of Neoglaciation in CSCh and to the north are also coarsely resolved. For example, glacier advances have been dated in the nearby Rio Grande basin (35°S) on the eastern (Argentinian) slope of the Andes: radiocarbon dating from peat bogs in contact with terminal moraines, which are located at elevations around 2,500 m, yielded ages between ∼0.5 and 0.2 ka (Espizua and Pitte, 2009). On Monte Tronador (41°S), located to the south, tree-ring dating of frontal and lateral moraines (elevations around 1,000 m) of glaciers flowing toward the eastern side of the Andes indicated ages between ∼1 and 0.4 ka, with some smaller ridges suggesting advances between the late 19th and early 20th centuries (Masiokas et al., 2009). On the Chilean slope, ¹⁰Be exposure dating of boulders located atop moraines ridges (elevations between 1,700 and 2,000 m) in the foreground of the Cipreses glacier (34°S) yielded older ages compared to the other sites, clustered between ∼1 and 0.6 ka, but with some samples older than 8.0 ka (Sagredo et al., 2016). Thus, while regionally distributed records suggest ice advanced at or around the time of the LIA, comprehensive studies that combine multiple independent dating methods, such as tree-ring records, with direct moraine ages in the same valley are lacking (Fernández et al., 2022).

Moraines were initially mapped and targeted for sampling using Google Earth Pro imagery. The location and extent of these moraines was documented in greater detail in the field utilizing handheld GPS methods during our field campaign of November 2018, with additional geomorphic observations and mapping of the valley conducted in January 2022, when directional measurements were taken from striated rock surfaces to assess paleo ice-flow directions. Sample locations were measured with a single-frequency GPS receiver (Emlid reach). The receiver was placed on a rod and positioned on top of or to the side of each boulder, and recorded continuously for 10–15 min in RINEX format. Subsequently, these data were post-processed using the closest GPS base station available, located in Temuco, about 100 km from Monte Sierra Nevada (Figure 1). This differential correction delivered uncertainties in the order of ±0.8 m.

Basaltic boulders located atop moraine ridges were sampled for surface exposure dating, with a total of 10 processed for ³⁶Cl measurement. These boulders heights ranged from ∼0.42 to ∼2.3 m, and included basaltic conglomerates and fine-grained crystalline basalts (Figure 5). We preferentially selected the largest boulders perched on or near moraine crests, which increases the likelihood that the boulders have been stable and experienced minimal snow cover since deposition (Supplementary Figure S1.; Ivy-Ochs and Kober, 2008; Balco, 2011). Additionally, boulders were assessed for weathering characteristics to distinguish moraine boulders from rock fall (Zreda and Phillips, 2000; Briner et al., 2001). No suitable sample material for exposure dating was found on the expression of the proximal (M5) or near proximal (M4) moraines to the north of the stream, as an immature drainage network had eroded large sections of these moraines (Figure 3). Samples (∼1 kg) were chiseled from the upper few centimeters of these basaltic boulders. Topographic shielding was measured by collecting azimuth and zenith angles every 15° at each sample location using a Suunto tandem compass/clinometer.

Whole rock samples were sawed to a measured thickness of between 0.4 and 2.0 cm before crushing and sieving to collect grains between 250 and 125 µm. Chemical preparation of samples was conducted at the University of New Hampshire Cosmogenic Isotope Lab using a modified version of the protocols presented in Stone et al. (1996) and Licciardi et al. (2008). A key component of the modified procedures includes measuring ³⁵Cl/³⁷Cl on a ∼1 g aliquot of rock prior to ³⁶Cl extraction and measurement on the full sample. The ∼1 g aliquot was spiked with a ³⁷Cl-enriched solution (“LLNL Spike A”; ³⁵Cl/³⁷Cl = 0.93) and ∼4,000 µg of Br, which served to increase the size of the final precipitate. Total sample Cl was determined through isotope dilution methods. ³⁶Cl was extracted as Ag(Cl,Br) from full rock samples in two analytical batches following standard procedures. All samples received a carrier solution containing ∼4,000 μg of Br and a natural ratio Cl carrier (³⁵Cl/³⁷Cl = 3.127). This technique was used to better control the total Cl within the analytical batch and reduce potential AMS memory effects (Arnold et al., 2010; Finkel et al., 2013).

Measurements of ³⁵Cl/³⁷Cl and ³⁶Cl/³⁷Cl were conducted using the 10 MV tandem accelerator at the Lawrence Livermore National Laboratory Center for Accelerator Mass Spectrometry (LLNL-CAMS). Analytical uncertainty on ³⁵Cl/³⁷Cl measurements ranged from 0.42 to 0.45%. Analytical uncertainty on ³⁶Cl/³⁷Cl AMS measurements was between 4.3 and 13%. Sample ³⁶Cl concentrations are corrected for lab background using an average blank value for the two analytical batches (4.12 × 10⁴ atoms ³⁶Cl, n=2). Major and trace element concentrations in the samples were determined using x-ray fluorescence and inductively coupled plasma mass spectrometry at SGS Mineral Services in Burnaby, British Columbia, Canada. Analytical data used to determine surface exposure ages are provided in Table 3.

Moraine boulder exposure ages were calculated using a developmental version of the CRONUS-Earth ³⁶Cl calculator (Balco et al., 2008; http://stoneage.ice-d.org/math/Cl36/v3/v3_Cl36_age_in.html) and LSDn scaling (Table 2; Lifton et al., 2014). The ³⁶Cl calibration datasets applied by the calculator can be found in the ICE-D database at http://calibration.ice-d.org/allcds. Based on field observations, no corrections were made for rock surface erosion, vegetation cover, or snow coverage. A bulk density value of 2.7 was used for all samples. All cosmogenic ages here are presented as years prior to sampling (AD 2018).

We made corrections for nucleogenic ³⁶Cl production using a rock formation age of 130 ka (Muñoz et al., 2011). Nucleogenic ³⁶Cl is produced in moderate amounts when ³⁵Cl absorbs neutrons released by the decay of U and Th isotopes (Gosse and Phillips, 2001). In rocks with high Cl, U, or Th concentrations, nucleogenic ³⁶Cl can thus constitute a substantial proportion of the total ³⁶Cl. In typical exposure age calculations, the amount of nucleogenic ³⁶Cl in a sample is calculated assuming steady state production, and that value is subtracted from the total ³⁶Cl to obtain the cosmogenic component (Marrero et al., 2016b). However, in young rocks (<1 Ma formation age), production and decay of nucleogenic ³⁶Cl may not be in equilibrium, and erroneously assuming steady state could lead to uncertain and/or inaccurate surface exposure age determinations (e.g., Sarikaya et al., 2019; Anjar et al., 2021). To assess the sensitivity of our ³⁶Cl ages to the time of rock formation, we also calculated ³⁶Cl ages assuming a rock formation age of 30 ka. Calculated surface exposure ages presented in the results and discussion of this study use a rock formation age of 130 ka to represent the last stage of volcanism on Monte Sierra Nevada based on lack of visual evidence of younger vents or deposits to support a formation age of 30 ka (Muñoz et al., 2011). Surface exposure ages representative of a younger date of rock formation of 30 ka are presented in Supplementary Table S1, as a possible late Pleistocene formation age has not been dismissed by available evidence.

Dendrochronology was used to obtain an approximate date of establishment for Araucaria araucana (Pewén) trees growing on the most distal moraine in the valley, to provide an independent minimum age of moraine formation and deglacial processes. To assess the minimum age of the moraine, increment borers were used to collect samples from mature A. araucana trees atop the outermost moraine (M1) of the proglacial area of the study site (Figure 3). These trees were selected with consideration to morphological attributes such as bark, height, diameter, as well as size and form of branches, all visual characteristics associated with old A. araucana trees. Two cores per tree from A. araucana (n=35) were extracted from living trees growing on the site using the aforementioned criteria, forming the “NEV” chronology. Samples were taken at 1.3 m from the soil, to assess a minimum age of the trees. Tree cores were then sanded for visual identification of growth rings (Stokes and Smiley, 1996), in order to be counted and measured in the Laboratory of Dendrochonology and Environmental Studies at Pontificia Universidad Católica de Valparaíso, Chile (https://www.pucv.cl/uuaa/site/edic/base/port/dendrolab.html). Ring width was measured to ±1 µm resolution using a Velmex system (https://www.velmex.com/), after which the COFECHA software (Holmes, 1983) was used to assess and corroborate the dating of each tree using the ring width measurement patterns across the suite of tree cores. Of the 35 sampled trees, 30 were able to be properly cross-dated. From these 30 trees, 20 exhibited curvatures close to the pith. We selected the five oldest trees that exhibited this curvature, suggesting a minimal age close to the total age of the tree, and using the Duncan method (Duncan, 1989) estimated the number of missing rings. The age of the oldest tree was used to provide an independent and cross-dated minimum age of the glacial landform from which the trees had colonized.

We applied a modified version of DETIM by Hock (1999), which is available at http://regine.github.io/meltmodel/. As inputs, we used a climate file from the Global Historical Climatology Network (GHCN) (https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network-daily) with daily data from AD 2004–2013 (Temuco station, Figure 1), as well as the gridded NASADEM surface. The model also requires gridded surfaces for aspect and slope values that were calculated using the spatial analyst tools in ArcGIS 10.4.1. The model ran on the Cygwin64 C compiler following the approach of Stansell et al. (2022). The model ran with degree-day factors (DDFs) for ice and snow of 7 and 3.5 mmw. e. K⁻¹ d⁻¹, values typically used for Chilean glaciers (Möller and Schneider, 2010).

The climate input file requires some pre-processing of climate data, including the calculation of a precipitation gradient and temperature lapse rate. The precipitation amount was set at a value of ×1.5 the values recorded from the Temuco station, which is similar to what is observed at other GHCN stations of similar elevation offsets in the region. No lapse rate was applied to the precipitation amounts along the glacier profile.

The model also requires an estimated temperature offset (Tdiff) corresponding to the location of the oldest mapped moraine (M1), which we independently determined by finding the difference between present (pr) and past (pa) Equilibrium Line Altitudes (ELAs) and applying a simple linear relationship using a lapse rate of 6.5°C/100 m: T d i f f = ( E L A p r − E L A p a ) × 0.0065

We used the Area Altitude Balance-Ratio (AABR) method (Osmaston, 2005) to determine ELAs. The AABR method estimates the ELA using the hypsometry of the glacier or the proglacial area within the limits of former ice extents, usually demarked by moraines. We used the 12.5 pixel size DEM available from ALOS PALSAR imagery as a source of topographical information to build the hypsometric curve (ASF DAAC, 2015), and a 3-m resolution true-color composite from the Planet Constellation to interpret the glacier margin in 2013 and to derive ELA _pr . For the past hypsometry, we trimmed the ALOS DEM extending the current margin by following the pattern of the trimline to the moraine crest closest to the cliff (M1). For the calculation, the method iteratively searches the elevation of the zero-balance, assuming that the mass balance gradient is linear, with a steeper slope in the ablation zone (Kaser and Osmaston, 2002). A key component of this approach is to select a suitable ratio between mass balance gradients, which can be computed from mass balance observations and assumed to hold in the past. In our case we face the challenge that almost no mass balance data exists for the region. The few direct measurements available from the Mocho Glacier (Rivera et al., 2005; Schaefer et al., 2017), situated about 115 km to the south, suggest a balance ratio of 3.2, larger than figures commonly used, such as 1.75 or 2 (Osmaston, 2005). A recent study found that a value of 1.56 represents glaciers at a global scale, although with important spread (Oien et al., 2021). Guided by these studies, we tested several ratios, opting for 2.38 as a best approximation to represent this region. We notice, however, that the difference in ELAs between ratios of 1.56 and 3.2 is 47 and 77 m for the present and the past, respectively, which translates to an uncertainty of 0.2°C. Using the median AABR value of 2.38 results in an ELA _pr and an ELA _pa of ∼2,267 m a.s.l. and ∼2029 m a.s.l. This change in ELA (ΔELA) of -238 m translates to T _diff of ∼−1.5°C.

We then ran the DETIM with multiple climate scenarios. First, the model was run over the interval from A.D. 2004 to 2013 with no change in temperature and precipitation relative to today. Next, we ran the model with 1.5°C cooling and a 20% increase in precipitation as a first approximation to simulate Neoglacial conditions that are represented by the AABR method. We assumed wetter conditions were associated with ice advances because proxy and modeling studies suggest that this region had increased precipitation amounts at times when it was colder during the Late Holocene (Bertrand et al., 2014). In order to test glacier sensitivity, we also ran the model with temperature changes of +1°C with a doubling of precipitation amounts. A series of other tests were run by changing only one variable at a time using ±1°C change in temperature, ±20% change in precipitation, and lapse rates ranging from 0.55°C/100 to 0.75°C/100. Additionally, to test the influence of DDFs, we adjusted the ice and snow values by ±1, following the approach of Möller et al. (2007).

Finally, we applied the flow model of Plummer and Phillips (2003) to simulate glacial thickness using the input of the mass balance grid. While the mass balance model identifies areas of ablation and accumulation, the flow model predicts the shape, size, and extent that a glacier stabilizes at for those conditions at steady-state. The model applies the standard flow law for ice using the shallow ice approximation, and using multiple interpolation methods. The results of the combined energy mass balance and flow models are shown in Figure 6. For all model runs, we present the lowest altitude of ice limit, maximum and average ice thickness values (Table 4).

The surface exposure ages of moraines on Monte Sierra Nevada range from 4.2 ± 0.2 ka to 0.28 ± 0.05 ka, and represent multiple phases of glacial growth and decay (Figure 7; Table 2). Sample AFCL18-16B, located on the most distal moraine in the valley, M1, is dated to 4.2 ± 0.2 ka. A second boulder sampled from M1, AFCL18-15, yielded a younger age of 0.65 ± 0.12 ka. Two samples dating to 2.1 ± 0.1 and 0.79 ± 0.08 ka are located on moraine M2, near the cliff on the north end of the valley. On the center moraine of the complex, M3, samples dated to 0.39 ± 0.05 and 0.28 ± 0.05 ka. Three samples from M4 yielded exposure ages of 2.6 ± 0.1, 0.46 ± 0.07, and 0.42 ± 0.06 ka. On M5, the innermost moraine in the valley, a single sample dated to 0.82 ± 0.09 ka. No apparent age biases (e.g., too young or too old) were detected between the fine-grained crystalline vs the basaltic conglomerate source rocks.

The cross-dating process resulted in an inter-correlation value of 0.55, a high correlation for A. araucana ring width chronologies (Muñoz et al., 2014). The five oldest trees comprising the NEV chronology, which showed a complete curvature close to the pith, were established between five and six hundred years ago. The first rings in samples from this group of trees were developed around the years AD 1407, 1462, 1511, 1559, and 1609. With the addition of missing rings due to the exhibited curvature (Duncan, 1989), the first ring for each of these five trees results in an approximate year of AD 1318, 1378, 1482, 1535, and 1580, respectively. The oldest tree of the NEV chronology yielded a date of AD 1318, providing a minimum limiting age of moraine M1 at 0.70 ka. This tree cross-dated to 0.54 with the complete set of trees measured, a value close to the mean inter-correlation value at this site, corroborating the age of this tree. Finally, due to the height of the extraction (1.3 m from the ground), a number of rings were missing with respect to the complete age of the trees, leading to an underestimation of the age of each tree. As such, a number of years, not determinate, should be added on to the determined age after ring counting, necessary to obtain the growth height observed. Following Lusk and Le Quesne (2000), the mean height growth rates of A. araucana were 169 mm/year in sun-grown trees, and 111 mm/year in shade-grown individuals. Considering these estimations, a minimum period of around 10 years should be added to the age of each tree. However, depending on competition and site requirements, this period could be longer than 10 years.

The combined glacier mass balance and flow models produced simulated ice extents that mimic the modern and paleo-conditions. Using the GHCN temperature and precipitation values produced an ice margin at an elevation (∼1905 m) that is similar to today. Using a cooling of 1.5°C and a 20% increase in precipitation produced a modeled ice extent that is similar to the outer-most mapped and dated moraines on Monte Sierra Nevada (Figure 6). As the moraine complex represents an elevation change of only ∼19 m over a distance of ∼350 m, distinguishing between moraines on the gridded volume proves difficult. In general, the model produces the thinnest ice near moraine positions and thicker ice elsewhere. Glacial coverage on Monte Sierra Nevada also tends to be thicker on the eastern slopes that are less steep than the western slopes. For modeling purposes, we assumed that periods of colder temperatures were associated with wetter conditions over the last several millennia (see Discussion). Nevertheless, the sensitivity analysis indicates that varying temperature alone by ±0.5°C shifted the elevation of the lowest ice position by ∼±130 m, while adjusting precipitation by ±20% shifted the ice margin by ∼−40 and +125 m (Table 4). Shifting the lapse rate by ±1.0°C/100 m offset the lowest ice position by −230 and +330 m. Lastly, shifting the DDF values by ±1 offset the lowest ice limit by +120 and −110 m.

The exposure ages presented here are some of the youngest to be developed using ³⁶Cl methods on glacial boulders in the Andes. While our dataset includes only a modest number of ages, these results demonstrate the feasibility of using ³⁶Cl to date young glacial deposits of similar compositions elsewhere in the Chilean Andes. Due in part to the small number of boulders, definitive and discreet ages cannot be confidently assigned to each moraine position. However, the chronology established for this moraine complex provides clear documentation for multiple phases of retreat and readvance of Monte Sierra Nevada glacier during the Late Holocene.

Multiple sources of uncertainty must be considered when using ³⁶Cl to date Late Holocene moraine boulders. Processes such as post depositional erosion and boulder rotation can contribute to uncertainties in surface exposure dating (e.g., Balco, 2020). In addition, the presence of nucleogenic (i.e., non-cosmogenic) ³⁶Cl in rock surfaces can introduce considerable uncertainty to exposure age calculations. In a setting such as Monte Sierra Nevada, where the target landforms for dating and the lithologic units are likely relatively young, low-energy neutron capture can also produce non-negligible amounts of ³⁶Cl. In exposure age calculations, the nucleogenic ³⁶Cl must be subtracted from the total ³⁶Cl to obtain the concentration of cosmogenic ³⁶Cl. At this time, ³⁶Cl production rates from low-energy neutron capture are the least well-characterized production pathway (Marrero et al., 2016a), partially because low-energy neutron flux depends on difficult-to-constrain paleoenvironmental factors such as snow cover and presence of water (Zreda et al., 1993; Phillips et al., 2001; Schimmelpfennig et al., 2009; Zweck et al., 2013; Dunai and Lifton., 2014). Improved constraints on the timing of volcanic episodes of Monte Sierra Nevada will help to fully evaluate the nucleogenic production of ³⁶Cl in the Monte Sierra Nevada moraines, as these Pleistocene deposits have likely not reached steady-state production for nucleogenic ³⁶Cl (Sarikaya et al., 2019; Anjar et al., 2021). With these considerations in mind, we applied recent advances in ³⁶Cl methods to produce robust ages for the last several millennia.

The cosmogenic ³⁶Cl ages from Monte Sierra Nevada suggest glaciers in the CSCh Andes fluctuated throughout the latest Middle and Late Holocene. Ice margin retreat from the oldest and most distal moraine, M1, likely commenced by ∼4.2 ± 0.2 ka based on sample AFCL18-16B (Figure 3). A second boulder sampled from M1 to the north of the stream yielded a much younger age (AFCL18-15; 0.65 ± 0.12 ka), which aligns with the establishment of A. araucana trees on this glacial landform by 0.70 ka. Additional exposure samples are needed from moraine M1 to further assess the potential of nuclide inheritance due to unknown timing of boulder stabilization over melting debris covered ice, as well as to produce a statistically robust age determination of this landform. However, the morphology of the M1 moraine might aid in deciphering this complex glacial setting. The subtle duel-crested nature of the southern extent of M1 and observed pedogenesis on the moraine surface (Jenny, 1944; Mason and Jacobs, 2018) may suggest that M1 demarcates an expansion of the ice front prior to 4.2 ± 0.2 ka, and again before 0.65 ± 0.12 ka. Prior to 0.65 ± 0.12 ka the glacier likely expanded to the position of moraine M1, resulting in reoccupation of this landform, abutted the ice-proximal expression of M1 formed in the mid-Holocene, and shaped the dual-crest observed on the southernmost extent. The glacier then commenced its final retreat from the M1 position by ∼0.65 ± 0.12 ka.

The composite nature of moraine M1, with a potential formation in the mid-Holocene followed by a latest Holocene glacial advance reoccupying this landform, would allow for a period of several thousands of years for soil development to occur. While soil may develop rapidly on deglaciated terrain, this process is largely dependent on environmental factors such as humidity and precipitation, as well as weathering rates of the original parent material (Mason and Jacobs, 2018). The lack of soil development on moraines M2-M5, along with an initial colonization age of moraine M1 by A. araucana trees of 0.70 ka, suggests that attributing a younger moraine formation age of ∼0.65 ± 0.12 ka would preclude sufficient time for pedogenesis to support the colonization of A. araucana on that landform. When combined, these independent datasets suggest glaciers were at more advanced positions in the region just prior to ∼4.2 ka, and supports assigning M1 an initial mid-Holocene formation age.

The trimline abutting M1 along the southern ridge of the valley suggests these features formed concurrently, and while no samples were dated from up-valley of the moraine complex, we suggest the trimline formed contemporaneously with M1. The continuation of the trimline along the southern ridge up-valley of the chute suggests that ice from the cirque directly up-valley from the dated moraines contributed to the formation of M1 in the mid-Holocene. However, this up-valley source of ice was likely minimal during the Late Holocene, with ice predominantly flowing over the southern ridge from the adjacent valley, as interpreted from the cross-cutting stratigraphy of the chute eroding a portion of the trimline (dashed blue, and black lines respectively, Figure 3). The modeling results suggest that greatest ice thicknesses in the valley were obtained when the ice margin expanded past the mapped extent of the moraine complex, which based on both the exposure ages and dendrochronology, has not occurred in the Late Holocene. Thus, this expression of high ice thickness in the valley likely formed during a period of glacial expansion during the mid-Holocene, predating glacial retreat from the M1 position at ∼4.2 ± 0.2 ka. The glacier then retreated to an up-valley position before retreat halted. During this interval of glacial recession, ice from the cirque likely retreated up valley, separating from the glacial tongue flowing over the saddle on the southern ridge.

Four samples (AFCL18-01; 2.6 ± 0.1 ka, AFCL18-05: 0.81 ± 0.09 ka, AFCL18-12; 2.1 ± 0.1 ka, and AFCL18-13; 0.79 ± 0.08 ka) date from 2.6 ± 0.1 ka to 0.79 ± 0.08 ka, and are located at the northern and southernmost extents of the moraine complex (Figure 3). The exposure ages obtained from these samples support a fluctuating ice margin, and initial exposure from the subglacial environment predating a final glacial advance in the valley during the latest Holocene. The mapped locations of these samples at the outermost extents of the moraines on which they are perched suggests reworking and transport to their present locations by one or more readvances of the glacier. Additionally, the older exposure ages of these marginal samples compared to samples closer to the center of the moraine complex supports the likelihood of multiple stages of ice marginal expansion. Ice flow directions assessed from striated surfaces (Figure 3) indicate a radial advance of glacial ice in the valley, which would have obliterated previously formed moraines, and redeposited these boulders in more lateral positions of the moraine complex. These samples may have originally been deposited on the valley floor during recession of the ice margin, or were potentially exposed from subglacial conditions on the southern ridge. Additionally, these samples may have variable amounts of ³⁶Cl inheritance, which could explain some of the scatter seen in the available ages. Finally, these boulders may have experienced destabilization and subsequent rotation as underlying ice melted, an inherent pitfall to interpreting surface exposure ages on ice-cored moraines.

The youngest surface exposure ages (AFCL18-07, 0.42 ± 0.12 ka; AFCL18-08, 0.46 ± 0.07 ka; AFCL18-09, 0.40 ± 0.05 ka; AFCL18-11B, 0.28 ± 0.05 ka; and AFCL18-15, 0.65 ± 0.12 ka), and the development of moraines M2-M5 provide evidence for multiple phases of ice expansion in the valley over the latest Holocene. Two scenarios are plausible to explain the formation of moraines M2-M5 between this final glacial advance and the modern. The first would entail a glacier with minimal to no debris coverage commencing a continuous pulsed retreat of the ice marginal environment from 0.65 ± 0.12 ka through the modern, emplacing moraines M2-M5. However, it is more likely that during this interval, and perhaps the entirety of the Holocene, that the lower portion of the glacier was debris-mantled, as debris coverage on many regional glaciers has been increasing in modern times (Glasser et al., 2016). If the glacial snout was indeed debris-mantled, the response of the ice marginal environment to climate-driven changes could be more muted (Kirkbride and Warren, 1999; Naito et al., 2000; Davies and Glasser, 2012).

Critically thin ice margins along the frontal edge of the glacier likely resulted in instability of slope angles on moraine M1, and potentially along the northern and southernmost lateral positions of the other moraines in the complex, leading to partial or full collapse and deposition of supraglacial debris. Several exposure ages towards the center of moraines M3 (0.40 ± 0.05 ka) and M4 (0.42 ± 0.12 ka and 0.46 ± 0.07 ka) are in close agreement. The location of these samples may have been better insulated from more rapid fluctuations of the ice margins at the time collapse was initiated, as surface melt rates often vary spatially depending on the thickness of debris coverage (Östrem, 1959; Moore, 2018). The timing between the final abandonment of moraine M1 at 0.65 ± 0.12 ka and the exposure of boulders on the central portions of moraines M3 and M4 at ∼0.4 ka may represent an interval when climatic conditions were cooler and wetter, resulting in temporary stabilization of the ice margin at moraine M3. Furthermore, the cluster of ages on moraines M3 and M4 between 0.46 ± 0.07 ka and 0.40 ± 0.05 ka might suggest that latest Holocene warming and/or drying resulted in ice surface lowering, the melt out of supra- and englacial debris, and the continued recession of the ice margin up valley. However, as debris can either increase insulation or enhance ablation of the glacial ice depending on both the thickness and extent of coverage (Davies and Glasser, 2012; Bartlett et al., 2021), further data and modeling are needed to substantiate these tentative interpretations.

On the ice-proximal side of moraine M5 to the north of the stream, hummocks were surveyed (Figure 3, Figure 4). The area over which hummocks are observed has been substantially dissected by an immature drainage network that has also eroded large sections of moraines M4 and M5. Development of this chaotic terrain suggests a rapid retreat from the location of moraine M5 (Owen and England, 1998; Figure 3). The development of hummocks also supports the idea of a debris-covered glacier in the lower reaches of the valley, particularly along the glacial snout (Anderson, 2000; Bartlett et al., 2021). An abundance of observed glacial scree and talus deposits in the valley, particularly along the southern ridge of the valley, would provide sufficient source material to generate debris coverage of the snout of the paleo-glacier (Shakesby et al., 1987; Giardino and Vitek, 1988; Moran et al., 2016; Charton et al., 2021b). As the talus is most extensive beneath the trimline and covering the chute along the southern ridge of the valley, debris coverage of the lowest extents of the glacier most likely increased following the initial abatement from the M1 location and thinning of the glacial surface.

Moraines M4 and M5 to the north of the stream are also weakly preserved, potentially indicating that these moraines were ice-cored, where dead ice would possibly be entrained along with glacial debris (Evans, 2009). Deflation of these glacial landforms would occur as the ice core melted, resulting in the poor preservation observed, and the destabilization of perched boulders (Putkonen and Swanson, 2003; Briner et al., 2005). A debris-mantled glacier would present a possible explanation in the scatter of the cosmogenic ages from this study site, such as the younger exposure age of sample AFCL18-11B (0.28 ± 0.05 ka) on M3, as it may have experienced some rotation due to debris coverage limiting ablation and the enhanced likelihood of boulder rotation upon destabilization.

The glacier modeling results and geomorphic evidence provides additional support that the present-day ice cap system was more extensive at times during the Late Holocene. Our model indicates that as the ice cap expanded, it likely covered more ridges and had additional source-ice regions with more complex ice-flow dynamics down-valley than what is seen in present-day (Figure 6). These complex ice patterns likely explain part of the scatter of available ages as the glacier did not follow a standard down-valley glacier pathway in the latest Holocene. Additionally, an eroded trimline indicates ice was flowing from south-to-north in the valley from a ridge that is currently ice-free (Figure 3). Combined, this evidence supports the initial formation of moraine M1 via a combination of glacial ice flow from the cirque of the modern glacier, merged with ice sourced from the valley directly to the south via the saddle along the southern ridge, and possibly the north side of the valley as well (Figure 3). The eroded trimline feature also suggests that following the abandonment of the M1 position at ∼4.2 ± 0.2 ka, glacial ice was likely sourced solely from the southern ridge of the valley with little to no contribution directly down-valley from the glacial cirque. That ice then flowed outward in a radial piedmont-style pattern across the table-top-like valley floor to reoccupy M1, and to the positions of moraines M2-M5 as supported by the striated surfaces throughout the lower extents of the valley.

Our chronology of Late Holocene glacial fluctuations from Monte Sierra Nevada is further augmented with independent tree-ring data collected near the margins of the apparent trim-line (Figure 3). The ³⁶Cl ages suggest that a Late Holocene retreat of the frontal/lateral margin occurred at ∼0.65 ± 0.12 ka. While the chronology of the oldest tree provided by the NEV record only extends to 0.70 ka, it provides a minimum-limiting age of glacial retreat from M1. It is possible that this landform has been colonialized by successive generations of A. araucana trees, as very few specimens have been dated to older than 500 years (Aguilera-Betti et al., 2017). Additionally, the period necessary between moraine stabilization and seedling establishment depends in part on the soil evolution of the landform to provide conditions for germination and growth, as well as individual species capacities to grow under limited conditions. This period could take several decades in most cases. However, A. araucana is a recognized species that can grow on low developed soils (Veblen, 1982). This, combined with the time required for the tree to grow to the height at which cores were extracted, results in a likely formation age of moraine M1 at minimum a few decades prior to 0.70 ka, further suggesting the true age of the moraine is older than 0.70 ka.

The modeling results presented here provide a physical and quantitative basis for interpreting glacial variability at Monte Sierra Nevada, and despite not being fully capable of handling the complex dynamics of a debris-mantled glacier, the Temperature-Index approach provides an important first-order estimate at quantifying paleo-environmental conditions on Monte Sierra Nevada. Additionally, the interpretation of a debris-covered glacier for this valley concerns only the snout, not the entire ablation area, and even less of the whole glacier surface. On modern decadal timescales, the dynamics of debris-covered glaciers are observed to operate differently than clean glaciers, but very little research exists to determine whether debris cover exerts a significant effect on the temperature sensitivity at centennial to millennia scales. Furthermore, when climate trends are strong, debris-cover does not cause a significant insulation effect (Burger et al., 2019), suggesting that under certain scenarios, in particular considering relatively persistent climate changes, debris-coverage does not alter glacier sensitivity significantly.

This model was parameterized with a range of representative temperature and precipitation values from the AABR calculations, suggesting temperatures were ∼1.5°C colder and precipitation was ∼20% higher during peak Neoglacial ice advances over the last several millennia. It is presently unclear, based on the available chronology, whether the peak Neoglacial ice advance at Monte Sierra Nevada was associated with the LIA or an earlier advance. However, the modeling outcome is generally the same regardless of the chronology when recognizing that it was wetter when it was colder at multiple times during the Neoglacial (Bertrand et al., 2014). Following this assumption, we were able to achieve ice margin reconstructions that match the mapped moraine positions on the landscape (Figure 6). The sensitivity analysis also highlights that a 20% increase in precipitation offsets a temperature shift by ∼0.5°C, which is less than other modeling studies that suggest a doubling of precipitation is needed to offset a ∼1°C temperature change (Figure 6; Sagredo et al., 2014). The results also highlight that multiple combinations of climate parameters are capable of explaining glacial variability in this region, and these modeling scenarios will improve as independent proxy evidence of Late Holocene environmental conditions, such as tree-ring records, continue to be developed.

The results presented here are not of sufficient resolution to attribute confident alignment of this record of glacial fluctuations to either the “Aniya-type” or “Mercer-type” chronology of Neoglacial events (Figure 8), however, certain similarities exist. Event I of Aniya (2013) is proposed to have occurred at ∼5.0 to 4.3 ka, while Mercer (1976) designates a first Neoglacial advance between 4.6 and 4.1 ka. Thus, the timing from both chronologies is similar to the oldest, outermost moraine age of 4.2 ± 0.2 ka. No samples dated in this study fall within the event II range of the “Aniya-type” chronology (3.7–3.4 ka), suggesting that either the climatic conditions in CSCh did not support a glacial advance during this interval, or that evidence of this event locally was destroyed by subsequent Neoglacial ice fluctuations. As the two exposure ages from this study (2.6 ± 0.2 ka and 2.1 ± 0.1 ka) that overlap with Mercer’s second event and Aniya’s III event are from proposed reworked boulders, we cannot assess any overlap in glacial advances from these samples. Similarly, while two exposure ages developed from our study (0.81 ± 0.09 ka and 0.79 ± 0.08 ka) fit into the event IV range (1.4–0.7 ka) of the Aniya chronology, we are unsure if these samples were recycled from a previous exposure due in part to sample placements at the far northern and southern lateral extents of the moraine complex. The latest Holocene glacial readvance is assigned an approximate range of 0.57 to 0.27 ka in the Mercer (1976) chronology and 0.52 to 0.3 ka (event IV) in the Aniya (2013) chronology. Four of the developed exposure ages from this study fit into this interval (0.46 ± 0.07 ka, 0.42 ± 0.12 ka, 0.40 ± 0.05 ka, and 0.28 ± 0.05 ka) of the Southern Hemisphere expression of the LIA. However as previously discussed, these ages are distributed across two moraines, M3 and M4 (Figure 3), precluding a definitive assessment as to which moraine represents the LIA with the current data set.

Slightly different regions were assessed to build the afore mentioned chronologies. The Mercer-type chronology evaluated a latitudinally longer corridor of glacial records from 39°S to 53°S, while the Aniya-type chronology was slightly more restricted and was compiled from records sourced from between ∼46°S and 52°S. Additionally, the inclusion of records from tide-water glaciers as opposed to a restriction of the analysis to fully terrestrial based glacial systems could potentially result in more variability in age control in these records. As highlighted by Glasser et al. (2004), the differences between the Aniya-type and Mercer-type chronologies are also likely more related to dating uncertainties and the lack of bracketing age control rather than actual differences in interhemispheric glacier behavior. The “Aniya-type” chronology though widely recognized, was developed largely from ¹⁴C dating of organic material combined with tree-ring data and lichenometry (Aniya, 1995; Aniya, 1996). Since the establishment of these chronologies numerous advancements have been made in various dating methods (i.e., cosmogenic nuclide exposure age dating, optically stimulated luminescence, tephrochronology, etc.) providing researchers a lens through which age constraints on regional records can be reexamined. Nevertheless, additional multi-proxy approaches similar to what we present here for Monte Sierra Nevada are needed as researchers work to develop comprehensive chronologies of past glacial variability in the high-latitude Andes.

The data reanalysis on the Patagonian Ice Sheet published by Davies et al. (2020) includes recalibrated chronologies of glacial fluctuations, but notes that the Late Holocene is not yet well understood over large regions of southernmost South America. However, three intervals of Neoglacial advances were well supported from the studies reviewed within: in the mid-Holocene at 5 ka, a smaller Late Holocene advance between 2.0–1.0 ka, and a LIA advance between 0.5–0.2 ka (Davies et al., 2020). Our proposed advance prior to 4.2 ± 0.2 ka, and the cluster of ages that fall within the LIA (0.46 ± 0.07 ka, 0.42 ± 0.12 ka, 0.40 ± 0.05 ka, and 0.28 ± 0.05 ka), are supported by this large-scale assessment from the Patagonia Ice Sheet. Nevertheless, our ability to disentangle how much past glacial variability in the Chilean Andes was driven by climate-related processes versus hypsometric controls will improve as we continue to apply more comprehensive dating methods, including additional direct ages on moraine boulders.

Andean glaciers have varying sensitivities to temperature and hydrologic changes depending on location, geometry and other local and regional factors (e.g., Sagredo et al., 2014). CSCh glaciers are typical of mid-latitude glaciers, with mass balance sensitivities that are seasonal and seemingly most sensitive to summer temperatures and winter precipitation (Masiokas et al., 2008). These temperature and precipitation anomalies are tightly coupled, with cold intervals corresponding to wetter periods and vice versa (Masiokas et al., 2009). Modern ice fields and outlet glaciers in Chile have historical mass balances that are sensitive to both atmospheric temperature and precipitation changes (Bown and Rivera, 2007; Casassa et al., 2007). Since 1975, the regional freezing level (0°C isotherm) has risen by ∼120 m in winter and 200 m in summer (Carrasco et al., 2005), with an overall increase of 0.25°C/decade in mean annual temperature since 1979 (Falvey and Garreaud, 2007). There has also been a strong negative trend in winter precipitation amounts over the last century (Masiokas et al., 2008).

Even though the annual cycles of temperature and precipitation vary across the Andes, it is possible that Pacific oceanic-atmospheric forcing synchronized climate variability across the low and high southern latitudes during much of the Late Holocene. For example, it has been proposed that SAM and ENSO are teleconnected and possibly related on centennial time-scales, providing one possible explanation for Southern Andean glacial variability (Moreno et al., 2018). It is also possible that changes in the position of the SWWB drove much of the observed pattern of Neoglaciation as these winds deliver critical precipitation to Patagonian glaciers (Glasser et al., 2004; Kaplan et al., 2016). Millennial-scale SAM variability began after ∼5.8 ka (Moreno et al., 2018), and it is suggested that insolation combined with SAM drives millennial and centennial-scale glacier variability in southern South America through the Holocene (Reynhout et al., 2019). An increase in Southern Hemisphere summer insolation (assessed at 37°S) from the start of the Holocene through the modern (Berger and Loutre, 1991; van den Bos et al., 2018), would support the largest glacial expanse recorded in the valley, predating a retreat as early as ∼4.2 ± 0.2 ka. A warm and dry period during the Late Holocene (∼4.1–2.8 ka) has been associated with prolonged SAM-positive phases, resulting in less extensive glaciers (Reynhout et al., 2019), while after 2.8 ka conditions appear to reflect more SAM-negative conditions, becoming colder and wetter, similar to the proposed intervals of SAM phases by Moreno et al. (2018). These SAM negative phases roughly correspond to intervals of assessed glacial expansion on Monte Sierra Nevada (Figure 8). Additionally, ENSO variability might make SAM stronger due to the effect of atmospheric Rossby waves across the Pacific. These SST anomalies in the tropics drive pressure anomalies in the high latitudes, as positive phases of ENSO lead to negative SAM, and vice versa (Moreno et al., 2018).

There are multiple proxy records of wetter conditions and evidence of a northward shift of the SWWB from ∼0.82 to 0.22 ka (Figure 8; Lamy et al., 2001; Jenny et al., 2002; Bertrand et al., 2005; Fletcher and Moreno, 2012; Bertrand et al., 2014). Perren et al. (2020) suggests that the extratropical westerlies migrate northward during cold Neoglacial periods, and therefore, it seems likely that periods of increased precipitation were associated with colder temperatures. The tree-ring and glacier modeling results presented here further suggest that both increased precipitation and colder conditions drove the pattern of latest Holocene glacier variability.

Regardless of the degree to which temperature vs. hydroclimate were the dominant factors, Late Holocene glacial variability in the Chilean Andes was likely influenced by complex oceanic and atmospheric processes. For example, there appears to have been a close relationship between Pacific and Southern Ocean SSTs and atmospheric circulation during the Late Holocene. Oceanic SST records (44°S) demonstrate the connectivity between the Southern Ocean with increased ENSO variability, resulting in an ∼2.5°C cooling trend between ∼1.2–0.7 ka off the central Chilean coast (Figure 8; Collins et al., 2019). Additionally, precipitation proxy records suggest an equatorward or northward shift of the SWWB associated with colder SSTs and a southward shift associated with warmer SSTs (Figure 8; Bertrand et al., 2014). The SWWB delivers critical precipitation to temperate glaciers, especially on the windward side of the Andes (Bertrand et al., 2012) where a strong correlation exists between precipitation amounts and westerly wind anomalies (Garreaud, 2007; Moreno et al., 2009). Historical archives compared to meteorological data support this observation and suggest that some glaciers on the drier side of precipitation gradients in northern Patagonia have varied over the last several decades mostly in response to precipitation changes (Warren and Sugden, 1993). There is also a correlation between Pacific SSTs and freezing level height in the Andes (Bradley et al., 2009), suggesting that shifting ocean circulation and temperatures could drive both temperature and hydroclimate variability in CSCh.

Proxy and modeling results suggest that decreased radiative forcing due to reduced solar activity causes an equatorward shift in the mean position of the SWWB (Varma et al., 2011). During the LIA, reduced solar forcing combined with increased volcanic activity led to net global cooling (Mann et al., 2009; Miller et al., 2012). This cooling would likely have produced a similar equatorward shift in the SWWB and wetter conditions in CSCh Andes, driving the observed pattern of ice advances. The higher temperatures that followed the LIA likely led to a poleward shift in the SWWB, and drier conditions that combined, would have caused ice to retreat.