Evolution of Glacial Lake Cochrane During the Last Glacial Termination, Central Chilean Patagonia (∼47°S)

1 Introduction

The Last Glacial Termination (T1: ∼18–11 ka, ka = 10³ calibrated years before present) is the most recent transition from glacial maximum to interglacial conditions and constitutes a key target for understanding the functioning of the global climate system through Quaternary ice ages. T1 featured a general trend of increase in atmospheric CO₂, along with abrupt temperature and atmospheric circulation swings at millennial timescales that led to the abrupt collapse of continental ice sheets from their LGM (∼26.5–19 ka, Clark et al., 2009) positions, and subsequent rapid sea-level rise (Denton et al., 2010). Recent studies have emphasized the importance of changes in the Southern Westerly Winds (SWW) as a key factor in the climatic transformations that terminated the last glaciation at regional, hemispheric, and global scale (Toggweiler et al., 2006; Lamy et al., 2007; Anderson et al., 2009; Moreno et al., 2018). The western portion of Patagonia (40°–55°S), in southern South America, is a key region to monitor T1 and associated SWW swings because it is the sole continuous landmass that extends south into the subantarctic realm, intersecting a large portion of the SWW belt. Furthermore, it harbored the largest ice mass in the Southern Hemisphere outside Antarctica during the LGM, the PIS, which allows the study of past changes in the cryosphere, atmosphere, hydrosphere, and biosphere, and their reciprocal relationships.

The PIS covered the Patagonian Andes and adjacent regions multiple times during Quaternary ice ages (Rabassa and Clapperton, 1990; Glasser et al., 2008). Large eastward-flowing outlet glacier lobes filled the basins of Lago General Carrera/Buenos Aires (LGCBA) and LCP in Central Patagonia (44°–48°S) during the LGM (Figure 1) (Turner et al., 2005; Bell, 2008; Hein et al., 2010). The Fénix and Río Blanco moraines define the LGM positions of these lobes, respectively (Douglass et al., 2006; Hein et al., 2010; Kaplan et al., 2011), and constrain the initial formation of ice-dammed proglacial lakes in both basins. These glacial lakes drained eastward to the Atlantic Ocean during the earliest stages of T1. Continuing ice retreat allowed the opening of new drainage routes and, eventually, both lakes merged through the upper Valle Río Baker (Turner et al., 2005; Hein et al., 2010; Glasser et al., 2016; Davies et al., 2018; Thorndycraft et al., 2019). The subsequent separation of the North and South Patagonian Icefields (NPI and SPI, respectively) and associated outlet glacier lobes, during the later stages of T1, allowed the drainage of this merged glacial lake via the lower Valle Río Baker towards the Pacific Ocean, establishing the drainage pattern that persists until the present (Turner et al., 2005). Deciphering the sequence of events through this process spurred a flurry of research over the last 17 years, leading to multiple interpretations (Turner et al., 2005; Hein et al., 2010; Bourgois et al., 2016; Glasser et al., 2016; Davies et al., 2018; Thorndycraft et al., 2019). We identify three main sources of divergence among these alternative models. First, although they agree on the identification of the principal steps that took place in the lake evolution, the timing of events varies considerably, owing mainly to differences in geomorphic/stratigraphic interpretations and the precision and accuracy of the available geochronological data. Second, most of the data come from the LGCBA basin and the upper Valle Río Baker, with relatively few field-based measurements and stratigraphic data from the LCP basin. And third, the majority of models do not consider the glacio-isostasy effects, which caused an important uplift on the western sector. The latter factor was discussed by Thorndycraft et al. (2019), who produced isostatically corrected lake surfaces. Their approach is based on digital elevation models (DEM) with coarse vertical (∼20 m) and horizontal (∼30 m) resolution along with selected geomorphologic, stratigraphic, and chronologic constraints.

FIGURE 1

FIGURE 1. Overview of the study region. Lake, river, glacier, and LGM limit shapes were obtained from Davies et al. (2020) for reference only. NPI, North Patagonian Icefield; SPI, South Patagonian Icefield; MSL, Monte San Lorenzo Ice Cap. LCP, Lago Cochrane/Pueyrredon; LGCBA, Lago General Carrera/Buenos Aires.

Here we report new data on shoreline-related landforms (raised deltas and paleoshorelines) in the LCP area, which is located ∼70 km east of the NPI (Figure 1), along with sediment cores from small closed-basin threshold lakes to decipher the evolution of GLC through T1. We provide Global Navigation Satellite System (GNSS) measurements with high precision for constraining the elevation at the brink points of the landforms, and the age for the transition from glaciolacustrine to organic-rich mud deposition in sediment cores from threshold lakes distributed at different elevations. We assess these results in the context of previously published geochronologies and interpretations, which we collated from the literature, recalculated, and recalibrated. We aim to 1) identify ancient lake levels in the LCP basin, 2) constrain their chronology, and 3) examine their response to glacial isostatic and climatic influences. These data allow assessment of the following questions: 1) When did the LCP glacier lobe abandon the LCP basin? 2) did recession of the LCP glacier lobe proceed steadily, or was it punctuated by stabilizations or reversals? 3) did variations in the extent and elevation of GLC represent responses to climate forcing through T1? and 4) what are the regional implications of the GLC evolution?

2 Previous Work

The first study of glacial landforms in the area was carried out by Caldenius (1932), who identified, named, and assigned tentative ages for several moraines east of LCP. A long hiatus in geomorphic studies followed until the publication of the Turner et al. (2005), Bell (2008), Glasser et al. (2008), and Hein et al. (2010) papers early during the 21st century. The latter study mapped multiple glacial limits around LCP and reported ¹⁰Be measurements from the Río Blanco moraines (Figure 2). These dates yielded a mean age of 20.7 ± 1.3 ka (recalculated) for the innermost crests (Hein et al.’s (2010) third limit). Subsequently, Boex et al. (2013) reported cosmogenic exposure ages from erratics at Cerro Oportus, with ages ranging between ∼19 and 17 ka (elevation range: 1,900–1,300 masl), Cerro Tamango, with ages ranging between ∼17–16 ka (elevation range: 1,500–500 masl), in the highlands that separate the LCP basin from Valle Chacabuco, and the María Elena moraines on the Valle Chacabuco floor (Figure 3), with a recalculated mean age of 16.2 ± 0.5 ka (elevation range: 600–450 masl). Davies et al. (2018) reported a recalculated mean age of 13.0 ± 0.4 ka for the Esmeralda lateral moraines and one sample from the Salto moraine (359 masl) (Figure 3) with an age of 12.5 ± 0.4 ka, both located in Valle Río Baker, southwest to LCP. They also dated two boulders below 350 masl on the Salto moraine with a recalculated mean age of 12.1 ± 0.4 ka, which they interpreted as anomalously young ages as a result of water shielding by the merged lake. Upstream in Valle Río Baker, the same authors reported a recalculated mean age of 14.0 ± 0.5 ka for the inner Lago Bertrand moraine (Figure 1), which they considered as a maximum limiting age for the merged lake phase, referred to in the literature as the “Lower United Lake” (Turner et al., 2005) and “Lago Chalenko” (Davies et al., 2018; Thorndycraft et al., 2019).

FIGURE 2

FIGURE 2. Distribution of shoreline-related landforms in the eastern sector of LCP, grouped according to their elevation in N5 (∼600–500 m), N4 (∼470–400 m), N3 (∼360–300 m), N2 (∼230–220 m), N1 (∼180–170 m), and N0, which correspond to active deltas. Lines represent paleo-shorelines (beaches, wave-cut scarps and benches, and delta brink-points) and polygons deltas. Geochronological data are shown in ka, without uncertainty (see details in Supplementary Tables S1–S3). Black lines represent ice limits.

FIGURE 3

FIGURE 3. Distribution of shoreline-related landforms in the western sector of LCP, grouped according to their elevation in N5 (∼600–500 m), N4 (∼470–400 m), N3 (∼360–300 m), N2 (∼230–220 m), N1 (∼180–170 m), and N0, which correspond to active deltas. Lines represent paleo-shorelines (beaches, wave-cut scarps and benches, and delta brink-points) and polygons deltas. Geochronological data are shown in ka, without uncertainty (see details in Supplementary Tables S1–S3). Black lines represent ice limits. Site 1 (Table 1) is shown as black circle. Red panels show sectors with field recognition (see Supplementary Figure S7 for more details).

Henríquez et al. (2017) studied sediment cores from Lago Edita (570 masl, GNSS), a small closed-basin lake located on the highlands that separate Valle Chacabuco from the LCP basin (Figure 3). They reported two identical radiocarbon dates from the basal organic sediments that yielded a median probability age of ∼19.4 ka. This constitutes a minimum limiting age estimate for local ice-free conditions and cessation of glaciolacustrine sedimentation at the Puesto Tejuela col (580 masl) (Figure 3). An identical age and interpretation apply to the Lago Augusta (440 masl, GNSS) (Villa-Martínez et al., 2012) record, a small closed-basin lake located ∼7 km northeast from Lago Edita in Valle Chacabuco (Figure 3). Turner et al. (2005) reported minimum limiting ages of kettle holes and shoreline features and concluded that the final drainage of the Lower United Lake occurred at ∼12.8 ka. Thorndycraft et al. (2019) combined radiocarbon dates reported by Turner et al. (2005) and cosmogenic ages reported by Glasser et al. (2012) to propose that the end of the Lower United Lake drainage occurred between ∼12.4–11.8 ka, using a Bayesian age model.

3 Materials and Methods

This study combines analyses of satellite data, field-based descriptions, and GNSS measurements of key sediment-landform assemblages, along with collection of sedimentary samples for documenting changes in depositional environments and obtaining absolute dates. We used high-resolution satellite imagery from Google Earth Pro and Esri™ World Imagery, ALOS PALSAR DEM (12.5 m horizontal resolution), and field data to construct a geomorphological map of selected sectors of the LCP basin (Figures 2, 3). We mapped shoreline-related landforms, including active deltas, raised deltas and their brink points, paleo-beaches, and wave-cut scarps and benches, using QGIS v3.10. We also described the sedimentology and stratigraphy of natural and artificial outcrops according to the depositional setting.

We acquired geodetic GNSS data (vertical datum WGS 84 ellipsoid; horizontal datum WGS 84) from the brink points of the landforms using a Trimble R1 GNSS receiver paired to a Juno 5 collector with Terrasync v5.85. The measurements were post-processed with the software Trimble Business Center, yielding an average vertical accuracy of ∼0.7 m. We extracted the elevation of geomorphic features that we were unable to visit on the field from ALOS PALSAR DEMs (vertical datum WGS 84 ellipsoid). Ten measurements were obtained from each linear shape. The DEM data accuracy was validated statistically by comparing it with our GNSS data and calculating their correlation coefficient, Root Mean Square Error (RMSE), and bias. The same process was applied to ASTER G-DEM to evaluate its vertical accuracy as well. We considered the NPI vertical elevation model proposed by Hubbard et al. (2005) and defined a simplified ice flow direction of 97° azimuth. We defined isobase zero orthogonally, passing through the Caracoles col (Figure 1). Finally, we plotted the elevation dataset in a scatter diagram against their distance to the isobase zero measured along the flowline of the LCP glacier lobe (we discarded landforms with strong evidence of erosion because of the difficulties in identifying their original brink point elevation). Using this plot together with field observations, we clustered landforms in six distinct groups according to their elevation (N5-N0, with N0 the lowest). We applied second-order polynomial regressions to N5, N4, and N3 using R (R Core Team, 2020) to obtain the glacio-isostatic gradients and paleolevels. Based on these results we adjusted the ALOS PALSAR DEMs in QGIS v3.10 to produce an isostatically corrected landscape reconstruction of the LCP area. We used these paleolake reconstructions, together with the inferred ice position, to estimate the water volume drained in each stage using raster analysis tools from QGIS. Some of the morphologies did not belong to the main levels and are therefore cataloged as sub-levels and not considered in the glacio-isostatic calculations. It is important to mention that although the aforementioned data processing was done using ellipsoidal elevations, all elevation data reported in the text are orthometric (EGM96).

We obtained multiple overlapping sediment cores from Lago Maldonado (−47.255, −72.512; 325 masl, GNSS) (Figure 3), a small closed-basin lake located on an ice-molded bedrock promontory ∼2 km west of LCP. Coring was carried out from an anchored raft equipped with a 7.5 cm diameter aluminum casing tube using a square-rod Livingstone corer. All cores were wrapped in plastic foil and stored at 4°C in a walk-in cooler. We obtained digital X radiographs to detect potential stratigraphic structures on the sediment cores, carried out loss-on-ignition analyses using the protocol described by Heiri et al. (2001) and obtained seven radiocarbon dates on bulk sediments, upon which we developed a Bayesian age model using the Bacon package of R (Blaauw and Christeny, 2011; R Core Team, 2020). We recalibrated the radiocarbon dates from Lago Augusta and Lago Edita (Figure 3) (Villa-Martínez et al., 2012; Henríquez et al., 2017) and updated the age models incorporating, in the case of Lago Augusta, the following modifications: 1) we excluded the radiocarbon date CAMS-144599 from the new age model considering it is anomalously young in the context of the over and underlying dates, and 2) we included the basal dates from Lago Edita (UCIAMS-133418, CAMS-144454) in the new Lago Augusta age model based on the Henríquez et al. (2017) discussion of the Lago Augusta record reported by Villa-Martínez et al. (2012), and their congruence with the basal date from Lago Edita (CAMS-144600).

We obtained two OSL dates from a raised delta front in Dos Arroyos (site 5) (Figure 3; Table 1). Measurements were performed on quartz following the single-aliquot regenerative-dose (SAR) protocols (Murray and Wintle, 2000, Murray and Wintle, 2003) at the Luminescence Laboratory, Indiana Geological and Water Survey, Indiana University. The SAR protocol was applied on 39 and 11 multi-grain aliquots (25–40 grains each). Equivalent dose and age distributions were estimated using the Minimum Age Model (Galbraith et al., 1999). The dose rate and age for each sample was calculated using the online DRAC v1.2 calculator tool (Durcan et al., 2015).

TABLE 1

TABLE 1. Summary of stratigraphic sections. See locations in Supplementary Figure S7.

We compiled 122 dates from previously published articles in the LGCBA and LCP basins, including cosmogenic, radiocarbon, and OSL dates (Supplementary Tables S1–S3). ¹⁰Be ages were recalculated with the CRONUS-Earth exposure age calculator v3, with version 2.2.1 of the constants files (Balco et al., 2008), and the high-resolution version of the geomagnetic framework (Lifton et al., 2008). All ages are presented using the time-dependent Lal/Stone scaling model (Lal, 1991; Stone, 2000) and the regional Patagonian production rate assuming zero erosion (Kaplan et al., 2011). We report cosmogenic dates with ±1σ analytical uncertainty. We decided not to use external uncertainties because we consider their use unnecessary given that all the ages presented here were calculated using a local production rate, which was empirically derived against radiocarbon dates. In this sense, for moraines, the arithmetic mean was calculated including the propagation of the analytical uncertainty and the local production rate uncertainties (3%) (Kaplan et al., 2011). All radiocarbon ages were recalibrated and converted to calendar years BP using the CALIB 8.2 (Stuiver et al., 2021) and the SHCal20 calibration dataset (Hogg et al., 2020) and are presented with ±2σ ranges. We report the calculated ages from OSL dates as reported in their original publications.

4 Results

4.1 Geomorphic Features and Sedimentology

Field mapping and satellite imagery revealed lake marginal features in LCP that include raised deltas, paleo-beaches, wave-cut scarps, and benches (Figures 2, 3). Raised deltas are distinguishable in the study area as flat surfaces composed of sedimentary materials, usually located upstream from active modern deltas (Figure 4). Their delta-top area range from 0.01 to 1 km², and their top slope varies between 1° and 10°. In general, they show well-defined brink points and steep fronts, and show well-developed topset and steeply inclined foreset beds (Figures 5, 6), which are the main features of Gilbert-type deltas, similar to those described in LGCBA by Bell (2008). Paleo-beaches show parallel lineations of gravel bars, gently concave towards the lake. Other paleo-shorelines are identified as narrow erosive or depositional terraces or scarps throughout the margins of LCP.

FIGURE 4

FIGURE 4. Field photographs from La Ponderosa (A), Frutillar (B,C), and Dos Arroyos (D–F) deltas. The location of these sites is expressed in Supplementary Figure S7.

FIGURE 5

FIGURE 5. Main stratigraphic and sedimentological features of deltas in Frutillar sector. (A) Outcrop showing the topset/foreset contact on site 2. (B) Trough cross-bedded gravel in topset. (C) Planar cross-bedding in topset. (D) Upward fining granules in foreset. (E) Topset/Foreset contact on site 3. Gm, clast-supported massive gravel; Gfo, deltaic foreset gravel; Gh, horizontally bedded gravel; Gt, Trough cross-bedded gravel; Gfu, Upward-finning gravel; GRm, massive granule; GRfu, upward-finning granule; GRp, Cross-bedded granule; GRfo, deltaic foreset granule; Sfo, deltaic foreset sand; Sm, massive sand.

FIGURE 6

FIGURE 6. Deltas from the Dos Arroyos sector. (A) Main stratigraphic features of site 4. (B) Detail showing main sedimentary facies of site 4. (C) Delta topset on site 5. (D) Topset and foreset on site 6. Sh, horizontally bedded sand; Gsi, matrix-supported imbricated gravel; Gms, matrix-supported massive gravel.

We mapped a total of 568 landforms (Figures 2, 3), which fall into six clusters of lake-marginal features with elevations between ∼600–500 (N5), ∼470–400 (N4), ∼360–300 (N3), ∼230–220 (N2), ∼180–170 masl (N1), and the currently active morphologies (N0), with their elevations close to the present-day water level (∼150 masl). The Río Blanco LGM moraines, located in the easternmost sector, feature shorelines at ∼500 (N5), 400 (N4), and 300 (N3) masl. They exhibit a steady rise in elevation toward the west, which explains the broad range of elevation within the clusters. In the following paragraphs we describe key sectors we visited during fieldwork or monitored from satellite imagery (Valle Río Baker).

Most of the landforms located in the Valle Río Baker (Figure 3) range between ∼360 and 350 masl, falling in the elevation range of the N3 cluster. The most conspicuous shoreline extends for ∼4 km, north of Cochrane township. The western side of Cerro Tamango contains shorelines with elevations of ∼460, ∼350, and ∼200 masl, associated to the N4, N3, and N2 clusters, respectively. We found glaciolacustrine rhythmites close to Río Cochrane at 173 masl (site 1 in Figure 3 and Table 1), composed of alternating millimetric clay and silt layers of dark and lighter sediments (Supplementary Figures S8A,B). Considering its elevation and the inferred ice position (discussed in Section 5.1), this deposit is most likely related to paleolevel N3 or N2.

We detect four raised deltas at La Ponderosa (Figure 3) which belong to clusters N4 to N1 (Figure 4A). N4 delta has a gently sloping upper surface and a steep delta-front slope, with a sharp and well-defined brink point at ∼470 masl. N3 presents a sharp brink point at 365 masl (GNSS) and erosional notches incised along its sloping front. The N2 delta extends over a larger area and has a diffuse brink point at 222 masl (GNSS). The N1 delta does not show a distinguishable brink point.

In the Frutillar sector (Figure 3) we detect two raised deltas (Figures 4B,C) which occupy small areas and have a moderately steep delta top slope (5–10°). The N4 delta has a somewhat diffuse and convex brink point at 455 masl (GNSS). The delta located within the N3 level does not show a clear brink point but its delta top has an elevation of 350 masl (GNSS). Outcrops in both deltas reveal the topset/foreset contacts for N3 and N4 landforms (sites 2 and 3, respectively) (Table 1). Topset beds at site 2 (Figures 5A–C) have a thickness of ∼3 m and correspond to alternating centimetric layers of darker matrix-supported pebbles in a matrix of coarse sand and layers of clast-supported sub-rounded well-sorted granules to pebbles, with normal gradation. We find scour and fill structures, planar and trough cross-bedding, and conglomerate lenses. The foreset beds are ∼2 m thick and exhibit steeply dipping (∼30°) decimetric layers of sub-rounded fine-grained to medium-grained massive sand with some isolated pebbles, intercalated with centimetric layers of coarse-grained sand to pebbles showing normal gradation (Figure 5D) and planar lamination. At site 3 (Figure 5E), the topset beds are ∼50 cm thick, and the contact between the topset and foreset is diffuse. Topset beds are composed mainly of subangular matrix-supported pebbles in a granule to sandy matrix, with a slight imbrication of clasts. Foreset beds have a thickness of ∼3 m and present a well-defined parallel lamination with an orientation of 25° NE. They correspond to alternating layers of massive matrix-supported granule clasts in a sandy matrix and massive matrix-supported pebbles in a granule sandy matrix.

The Dos Arroyos sector (Figure 3) contains the most extensive and well-preserved perched deltas in LCP (Figures 4D–F), which encompass elevations that correspond with the N5 to N1 clusters. The N5 delta top has a gentle slope and irregular surface, possibly due to fluvial erosion. The front is very eroded as well, but the brink point is distinguishable at 600 masl (GNSS). Deltas at N4 and N3 are extensive (∼0.8 and 1 km², respectively) and have sharp and well-defined brink points at 458 and 345 masl (GNSS), respectively. There are also minor breaks at 584 and 333 masl (GNSS). N2 and N1 deltas are considerably smaller in size and have no distinguishable brink points. In the following sentences we describe outcrops from the N1 (site 4), N2 (site 5), and N3 (site 6) levels. Site 4 features a ∼5 m high and ∼30 m wide outcrop that reveals alternating decimetric layers of sand-pebble and granule-boulder (Figures 6A,B). Sand-pebble layers correspond to matrix-supported pebbles in a matrix of angular medium-grained to coarse-grained metamorphic sand, showing horizontal stratification. Granule-boulder layers are composed of clast-supported subrounded pebbles to boulders (up to ∼30 cm diameter) in a granule matrix. Imbrication is occasionally observed in the deposit. Site 5 consists of a 2.3-m high section composed of clast to matrix-supported pebbles to cobbles of metamorphic origin in a medium-grained sandy matrix (Figure 6C). Clasts are sub-angular to sub-rounded. Well-developed imbrication and parallel lamination are evident in the finer sediments. At site 6 we find a ∼2 m high outcrop of tabular morphology composed of 2 conformable units recognized as topset and foreset beds (Figure 6D). Foreset beds consist of 1.5 m of alternating steeply dipping (∼30°NW) layers of clast-supported pebbles (∼1 cm) and layers of 3–20 cm of fine-grained sand with parallel lamination. We collected two samples for OSL dating from this unit ∼1 m below its upper limit (Table 2; Supplementary Figures S1, S.2), yielding an weighted mean age of 14.8 ± 1.1 ka. The topset is composed of ∼1 m of massive matrix to clast-supported pebbles-to-cobbles in a coarse sand matrix. The topset/foreset contact was also recognizable at N4 (site 7, Supplementary Figure S8C) and N5 (site 8, Supplementary Figure S8D) delta outcrops in the Dos Arroyos sector, but not accessible in situ.

TABLE 2

TABLE 2. OSL dates of delta-front deposits from Dos Arroyos site.

We found two small raised deltas in the Buena Vista sector (Figure 3) falling in the elevation range of the N5 and N4 clusters (Figures 7A,B). These exhibit smooth and inconspicuous brink points at 576 and 443 masl (GNSS) and minor slope breaks at 542 and 431 masl (GNSS). A very distinctive paleo-beach within N3 shows parallel and slightly concave gravel lineations, with its lowest shoreline at 331 masl (GNSS) (Figure 7C). Small deltas within the N3 and N2 elevation ranges are preserved in this sector but their brink points are not well-defined.

FIGURE 7

FIGURE 7. (A–C) Field photographs from the Buena Vista sector. (D) Shorelines between the Dos Arroyos and Buena Vista sectors. Note the remarkable preservation and continuity of N3. See location in Supplementary Figure S7.

We detect several small terraces within the elevations of the N5-N3 clusters between Dos Arroyos and Buena Vista (Figure 7D). The continuity of the N3 shoreline is remarkable, extending almost uninterrupted for ∼5 km with an average elevation of ∼330 masl.

We identify an extensive delta within the N3 range and smaller deltas at N2 elevations in the Río Brown sector (Figures 3, 8A). The brink point in the N3 delta is well defined at 324 masl (GNSS). There is a lateral section in its delta-top that shows alternating layers of clast-supported subrounded pebbles and clast-supported granules (site 9, Supplementary Figure S8E). The N2 deltas are considerably smaller with restricted, underdeveloped delta fronts, and a major slope break at 217 masl (GNSS), and minor breaks at 289 and 278 masl.

FIGURE 8

FIGURE 8. Field photographs from the Río Brown (A) and Paso La Leona sectors (B–E). See location in Supplementary Figure S7.

At Paso La Leona sector (Figure 3) we found geomorphic features that fall in the elevation ranges of clusters N5 to N3 (Figures 8B–E). The N5 landforms are preserved as isolated terraces and paleo-beaches. The lowest break of slope measured in the field has an elevation of ∼530 masl (GNSS). N4 landforms comprise some sharp isolated terraces as well, one of them with an altitude of 414 masl (GNSS), and a small delta with no distinguishable brink point. A large delta with a gentle sloping top surface shows a well-defined brink point at 319 masl (GNSS), and some associated terraces fall in the N3 elevation range. Topset beds were observed at the N3 delta (site 10, Supplementary Figure S8F). They consist of centimetric sub-horizontal layers with alternations between pebble-rich layers and finer sediments. West of Paso La Leona we observe several small terraces belonging to clusters N5 to N2 (Figure 3), whose elevations were extracted from ALOS PALSAR DEMs.

4.2 Isostasy and Paleolake Levels

The elevation dataset used to calculate the glacio-isostasy includes GNSS (sites with in situ field measurements) and ALOS PALSAR DEMs (sites lacking in situ field measurements) data. We first evaluated the accuracy of DEMs by comparing them with our GNSS data. They show a very strong correlation (r = 0.999) with a RMSE of 7.94 m. Although these are good results, we note that ALOS PALSAR systematically underestimates elevations, with a bias of −6.49 m (Supplementary Figure S3; Supplementary Tables S4, S5). Likewise, ASTER G-DEM shows a very strong correlation (r = 0.997) with a RMSE of 10.18 m and a bias of −6.38 m (Supplementary Figure S4; Supplementary Tables S4, S5). The elevation differences between the DEM and our GNSS data range between −18 and 2 m for ALOS PALSAR, and between −22 and 19 m for ASTER G-DEM, resulting in a higher terrain ruggedness (Supplementary Figure S5). Based on these results we favor the ALOS PALSAR dataset.

We plotted the elevation data obtained at the slope breaks of shoreline-related landforms against distance to the Caracoles col to explore trends along the former glacier flowline (Figure 9). Second-order polynomial regressions applied to each cluster (N5, N4, and N3) show very strong correlations (r ≥ 0.89). We observe successively decreasing slopes from the highest to the lowest clusters, which implies that glacio-isostasy has been stronger at higher (and older) than at lower (and younger) levels. Average gradients, between 0 and 100 km (i.e., between the Río Blanco moraines and the western end of the lake) are 1.16, 0.81, and 0.51 m km⁻¹ for N5, N4, and N3, respectively. The y-intercept of each regression represents the lake level during each phase, i.e, the paleolevel (see equations in Figure 9). The orthometric elevations of these are ∼485, ∼385, and ∼295 masl for N5, N4, and N3, respectively. Gradients for N2 and N1 were not calculated owing to the scarcity of morphologies classified within these clusters and the difficulty to identify their brink points. The reconstructions of these phases, therefore, are based on current elevations.

FIGURE 9

FIGURE 9. Glacio-isostatic gradients along the flowline of the former LCP glacier lobe. Polynomial regressions are shown with their equations, coefficients of determination, and the total number of points used (n = DEM/GNSS). Differential GNSS data are shown as white circles. Elevations in this figure are ellipsoidal and in the text are referred in their orthometric conversion (∼20 m lower). Current elevation of LCP shown as dashed line for reference.

Considering the paleolake reconstructions of phases N5, N4, and N3, together with the inferred ice position, which will be discussed in Section 5.1, we estimate water discharges of ∼150–200 km³ for the N5-N4 transition, 100–150 km³ for the N4-N3, and ∼300 km³ during the N3 drainage.

4.3 Lake Sediment Cores

We retrieved overlapping sediment cores from Lago Maldonado (core series PC0904A, B, and C) from a water depth of 680 cm. Core series PC0904A (Supplementary Figure S6) is 13-m long and features a blueish/grey inorganic mud (organic matter <5%) unit between 1,300 and 760 cm depth, with occasional laminations and matrix-supported clasts in the granule and pebble-sized range. This shifts to an organic-rich silty unit at 760 cm through a gradual transition between 760 and 742 cm depth. The organic unit persists with variations in the amount of organic material until the modern water-sediment interface. Seven radiocarbon dates obtained from the lowermost 80 cm of organic sediments yielded a sequence of ages between ∼12.2 and ∼9 ka (Table 3). Four out of five radiocarbon dates obtained in the narrow interval between 759 and 755 cm depth yielded statistically identical ages that produce an error weighted mean age of 9,828 ± 21 ¹⁴C yr BP with a median calibrated age of ∼11.2 ka. One of them (CAMS-144603) yielded an age of ∼12.2 ka that is inconsistent in the context of the stratigraphic sequence, being >3σ older than the enveloping radiocarbon dates UCIAMS-146440 and UCIAMS-146441. Therefore, we consider the age of sample CAMS-144603 as an outlier. A Bayesian age model applied to all the radiocarbon dates in the Lago Maldonado core pinpoints the onset of organic-rich sediments in Lago Maldonado at a median age of ∼11.3 ka (2σ range: 11.226–11.345 ka) (Figure 10).

TABLE 3

TABLE 3. AMS radiocarbon dates from sediment cores and calibrated ages.

FIGURE 10

FIGURE 10. Inorganic density, organic percentage, and age model from the Lago Edita, Lago Augusta, and Lago Maldonado records. The magenta line (lower panel) represents the median probability of the interpolated age, the green lines indicate the 95% confidence interval. The white circles in the lower panels indicate the median probability of the calibrated age for each radiocarbon date, the black horizontal bars represent the 95% confidence interval for each calibrated date. The yellow-filled diamonds represent the calibrated radiocarbon dates from Lago Edita. Notice that the age model from Lago Augusta includes an in-situ obtained radiocarbon date and the correlated basal radiocarbon dates from Lago Edita. The time intervals N5, N4, and N3 are shown as shaded rectangles.

Lake sediment cores from Lago Augusta and Lago Edita (Villa-Martínez et al., 2012; Henríquez et al., 2017) show abrupt transitions from high-density blueish/grey inorganic mud (organic matter <5%) to lower-density organic silts similar to that observed in Lago Maldonado. The exact onset of organic sedimentation in Lago Edita is constrained by two accelerator mass spectrometry (AMS) dates (UCIAMS-133418, CAMS-144454) succeeded by progressively younger dates further up in the stratigraphy. The age model estimate for this irreversible transition is ∼19.3 ka (2σ range: 19.5–19 ka) (Figure 10). In the case of Lago Augusta the shift from glaciolacustrine mud is constrained by one AMS date (CAMS-144600), overlain by one anomalously young date (CAMS-144599) located 8 cm above, and several AMS dates in stratigraphic order further up in the section. A Bayesian age model that excludes CAMS-144599 and incorporates CAMS-144600 along with the basal Lago Edita AMS dates mentioned above, yielded an interpolated age of ∼19.3 ka (2σ range: 19.6–18.7 ka) for the end of the basal inorganic mud. Light brown/yellowish low-density organic silts with plant macrofossils persist in the Lago Augusta record until 343 cm (mean organic content: 10.3%), followed by a shift to higher-density grey organic silts without plant macrofossils until 325 cm (mean organic content: 7.7%). From that point onward the inorganic density of the sediments declines and their organic content increases steadily to a maximum of 60% at ∼11.9 ka. The reversal to grey inorganic mud is bracketed by radiocarbon dates CAMS-144598 and CAMS-146711, and age model-based interpolations of ∼15.8 ka (2σ range: 16.9–15.1 ka) and ∼14.8 ka (2σ range: 15.2–13.9 ka) for its onset and termination, respectively (Figure 10).

5 Discussion

5.1 Evaluation of Paleo-Shorelines and Paleolake Reconstructions

We clustered the landforms in six groups according to their elevation, numbered N5 to N0 from highest to lowest, with N0 corresponding to the modern active deltas (Section 4.1). These clusters attest to at least five ancient phases of stabilization of GLC, which occupied the LCP basin following glacier withdrawal from the Río Blanco moraines.

Landforms grouped in N5 are present on the northern and southern valley walls surrounding the LCP basin (Figures 2, 3), with its westernmost position found at the Dos Arroyos sector. This distribution suggests that the LCP glacier lobe receded ∼95 km from the Río Blanco moraines, stabilized and established an effective dam for the development of GLC west of Dos Arroyos (Figure 11A). The reconstruction of paleolevel N5 shows that LCP and Valle Chacabuco were connected through the Puesto Tejuela col (580 masl), as stated originally by Glasser et al. (2016). This is in agreement with the findings reported by Henriquez et al. (2017), who observed exposed glaciolacustrine beds and lake terrace fragments (590 masl) around Lago Edita, documenting the highest-elevation proglacial lake phase in this area. Drainage of GLC during this phase must have occurred through the easternmost spillway of the basin, the Caracoles col (480 masl), toward the Atlantic Ocean.

FIGURE 11

FIGURE 11. Summary of the main phases proposed for the evolution of GLC. White dashed lines show inferred ice positions. Blue dashed lines show possible drainage routes. Paleolevels are shown in blue text. N5 to N3 phases (A–C) show glacio-isostatically corrected landscapes. N2 and N1 reconstructions (D,E) are based on current elevations.

N4 landforms are widespread from the Río Blanco moraines to the western end of LCP, with its westernmost position at La Ponderosa (Figures 2, 3), which suggests that the LCP Glacier Lobe retreated and stood just west of this sector by this time (Figure 11B). This phase began with a water level drop of ∼100 m from N5 to N4, which may represent an ice dam rupture or, alternatively, thinning and recession of the damming ice front. The N4 paleolevel reconstruction (Figure 11B) shows that GLC ceased to inundate Valle Chacabuco through the Puesto Tejuela col after the N5 phase and, likewise, ceased to drain through Caracoles col. We also note a ∼1.7 km long shoreline on the western side of Cerro Tamango within the elevation range of the N4 cluster (Figure 3), indicating continuity of GLC between Valle Chacabuco and the LCP basin through this short portion of Valle Río Baker.

N3 landforms are the most ubiquitous and extensive around LCP, suggesting that they formed during an extended period of lake stability, and/or the fluvial discharge and sediment load was higher during this phase. Landforms associated with this level are abundant along the adjacent portion of Valle Río Baker, suggesting local ice-free conditions during this interval (Figure 11C). Considering the distribution of N3 landforms, we infer that the damming glacier bodies must have lain downstream from LCP along the modern drainage system, potentially by the expanded Colonia and Calluqueo glaciers (Figure 1).

Geomorphic features clustered in N2 and N1 are very sparse and small, suggesting that these represent minor events during the late phases of GLC. The relatively small difference in their current elevations (N2: 230–220, N1: 180–170 masl) suggests ephemeral stabilizations in lake level during the regressive phase following N3. Considering the slight difference in elevation with the modern lake, it is likely that the current drainage configuration was achieved during these phases (Figures 11D,E).

5.2 Isostatic Adjustments

Glacio-isostatic uplift results are presented in Figure 9, and yield average gradients of 1.16, 0.81, and 0.51 m km⁻¹, for N5, N4, and N3 levels, respectively. The variability of elevations within each level is probably associated with the effects of glacio-isostatic rebound, the erosion over the brink points, the imprecision of ALOS PALSAR DEMs, and the difficulties involved in recognizing brink points in areas not visited during fieldwork. In the Dos Arroyos sector (where morphologies are best preserved), we find that the uplift due to glacio-isostasy reaches ∼115 m for the N5 delta, ∼70 m for the N4 delta, and ∼50 m for the N3 delta (this is considering the paleolevel of each phase vs. their current elevation).

5.3 Timing of Events

In this section we assess the age of the lake marginal features associated to GLC through comparisons with radiocarbon and exposure-age dated lacustrine and glacial features in and around the study area. The innermost Río Blanco moraines (Hein et al., 2010) have a recalculated mean age of 20.7 ± 1.3 ka (n = 3), which constitutes a maximum limiting constraint for the N5 level. A minimum limiting age for this level is afforded by the radiocarbon-dated transition from glaciolacustrine mud to organic lake sediments at ∼19.3 ka in Lago Edita (570 masl, GNSS) (Henríquez et al., 2017), which was flooded by GLC during the N5 phase. Lago Edita became isolated from GLC during the regressive glaciolacustrine phase that led to the N4 phase of GLC, attesting to the end of flooding of the Puesto Tejuela col.

The age for the cessation of glaciolacustrine influence in Lago Augusta at ∼19.3 ka (Villa-Martínez et al., 2012; Henríquez et al., 2017) is statistically identical to the age from Lago Edita. Unlike Lago Edita, however, Lago Augusta is located in Valle Chacabuco at 440 masl (GNSS). These findings imply a synchronous cessation of glaciolacustrine flooding in the highest section of Valle Chacabuco and Puesto Tejuela col by the end of the N5 phase. The stratigraphy from Lago Augusta indicates a subsequent transgressive glaciolacustrine phase in Valle Chacabuco between ∼15.8–14.8 ka. Together these results imply that: 1) the glaciolacustrine regressive phase following N5 must have reached elevations below 440 masl in Valle Chacabuco after ∼19.3 ka, and 2) the subsequent transgressive phase of Valle Chacabuco must have started from elevations below ∼440 masl before ∼15.8 ka.

Considering our N4 reconstruction, it is possible that GLC remained separated from and slightly higher than the Chacabuco glacial lake at the beginning of the N4 phase, probably by an ice dam west to Cerro Tamango (Figure 11B). Glaciolacustrine transgression of Valle Chacabuco during N4 might have occurred through hydraulic continuity with GLC through the damming ice body (englacially) and/or by slight recession and thinning of the ice dam. Rising Chacabuco glacial lake levels reached and surpassed the elevation of closed-basin Lago Augusta, depositing the grey organic silts without plant macrofossils. This persisted until ∼14.8 ka followed by a reversal back to organic-rich lake sediments in the Lago Augusta sediment cores, thus affording a minimum limiting age for the culmination of the N4 phase.

The N3 phase in the evolution of GLC is constrained by two OSL dates that yielded a weighted mean age of 14.8 ± 1.1 ka for the construction of the thickest and most extensive delta in the Dos Arroyos sector. The OSL dates inform us of the timing of construction of the delta front, which does not necessarily imply a close limiting age estimate for the beginning or end of this geomorphic feature. We note, however, that our chronology for the end of the N4 flooding of Lago Augusta is indistinguishable from the weighted mean age of the Dos Arroyos N3 delta. Our radiocarbon dates from Lago Maldonado, located at a current elevation of 325 masl (GNSS) and ∼12 km west of the N3 delta front in Dos Arroyos, afford a minimum limiting age of ∼11.3 ka for a regressive phase of GLC after the deposition of the N3 landforms. This independent minimum age constraint is coherent and compatible with the young confidence age limit of the OSL age. Combining these dates, we interpret that the N3 phase must have taken place between ∼14.8–11.3 ka (Figure 11C). Drainage of GLC after ∼11.3 ka was punctuated by brief stabilizations marked by the N2 and N1 deltas in the Dos Arroyos sector, for which chronologic data is currently unavailable.

5.4 Regional Implications

In this section we compare our reconstruction of the evolution of GLC with previous studies in the LGCBA, Valle Río Baker, and sectors adjacent to LCP.

(1) Our reconstruction of paleolevels in the LCP basin distinguishes three main stages between ∼20.7–11.3 ka, named N5, N4, and N3 phases. These reconstructions consider the elevation of more than 500 landforms and glaciolacustrine deposits along the flowline of the former LCP ice lobe. Although we followed a similar procedure to that reported by Thorndycraft et al. (2019), our results indicate larger-magnitude uplift estimates in the western sector, with average gradients of 1.16, 0.81, and 0.51 m km⁻¹ for N5, N4, and N3, respectively. Thorndycraft et al. (2019), in contrast, reported average gradients of 0.78, 0.65, and 0.51 m km⁻¹ for the equivalent paleolevels. We note an important difference for N5 and N4 which, according to our observations, emerges from the selection of landforms made by the different authors and the distinct elevation datasets used. Our study uses a combination of GNSS data with ALOS PALSAR DEM, which yielded higher precision results than ASTER G-DEM (used by Thorndycraft et al. (2019)), and a more accurate basis for the correct delineation of geomorphic features and extraction of their elevation data.

(2) Glacial withdrawal and rapid recession of the LCP ice lobe from the Río Blanco moraines occurred between ∼20.7–19.3 ka (N5 phase), which differs from the ∼18 ka age estimate for the LGCBA lobe (Bendle et al., 2017). This implies that the LCP glacier lobe abandoned the LGM moraines ∼2,700 years earlier and retreated much faster than the LGCBA lobe, most likely driven by intense calving. The LCP glacier lobe then stabilized in a position ∼50% from its final LGM extent and featured a much slower, step-wise recession through the remainder of T1.

(3) We consider the onset of organic sedimentation in small closed-basin lakes Lago Edita and Lago Augusta at ∼19.3 ka as the decisive marker for the termination of the N5 phase. In addition, this transition constitutes a maximum age for the emplacement of the N4 landforms and deposits. This implies that the 16.2 ± 0.5 ka age from the María Elena moraine (450–600 masl) (Boex et al., 2013), located adjacent to Lago Augusta on the Valle Chacabuco floor, should be considered as a minimum age estimate for local ice-free conditions. Shielding by the proglacial lake, shown in the paleolevel reconstructions, and potential exhumation may account for these anomalously young cosmogenic ages. Thorndycraft et al.’s (2019) model for the lake evolution did not take into account these constraints and interpretations, diverging from our reconstruction.

(4) The elevation, extent, and age of the GLC N3 phase are equivalent to Turner et al.’s (2005) Lower United Lake phase and “Lago Chalenko” in Thorndycraft et al.’s (2019) terminology. Our results suggest that N3 landforms developed between ∼14.8–11.3 ka, an interval within which Soler glacier abandoned the Lago Bertrand moraine LBM3 (14.0 ± 0.5 ka) (Davies et al., 2018), permitting the connection between LCP and LGCBA through the upper Valle Río Baker. We note that the ∼95% confidence interval for LBM3 (15–13 ka) overlaps with the confidence interval for our maximum limiting age estimate for the beginning of N3 (15.2–13.9 ka) and thus, are statistically indistinguishable. This supports that the connection between LGCBA and LCP basins started with the shrinking of the Soler glacier, as proposed by Thorndycraft et al. (2019).

(5) Turner et al. (2005) reported eleven radiocarbon dates aimed at constraining the chronology of the Lower United Lake phase. They deduced a minimum limiting age constraint of ∼12.8 ka, considering radiocarbon dates from shoreline sediments and minimum limiting ages from kettle holes. This estimate is ∼1,500 years older than the transition from glaciolacustrine to organic mud in Lago Maldonado. Unlike the results reported by Turner et al. (2005), the elevation and age estimates from Lago Maldonado are robust and precise, based on geodetic GNSS measurements with a vertical error of ±70 cm, and dating of autochthonous organic matter produced within a small closed-basin lake environment, in the context of a complete stratigraphic record from the onset of ice-free conditions until the present that shows a gradual and comformably lain transition from glaciolacustrine to organic mud. Moreover, our age estimate is constrained by multiple, statistically indistinguishable AMS radiocarbon dates sampled from contiguous and narrow stratigraphic interval immediately above the cessation of glaciolacustrine deposition, as Lago Maldonado detached from a lowering GLC during the abandonment of the N3 level. We attribute the divergence with Turner et al.’s (2005) estimate to imprecisions related to single radiocarbon dating measurements from each site, intersite variability, and differences in the elevation and diagnostic value of kettle holes to constrain lake-level fluctuations which, by definition are located on morainal topography which do not necessarily constitute prime locations for monitoring lake-level changes.

(6) The N3 glaciolacustrine phase provides the timing and context for the subaqueous deposition of the Salto and Esmeralda moraines (Figure 3). Davies et al. (2018) reported a date of 12.5 ± 0.4 ka (recalculated) for the Salto moraine at 359 masl, and a mean age of 13.0 ± 0.4 ka (n = 5, recalculated) (Supplementary Table S1) for the Esmeralda moraines (360–520 masl). These authors also reported paleoshorelines carved on their proximal and distal slopes at ∼350 masl. The depositional setting and ages for these moraines are coherent with our findings from LCP.

(7) According to Thorndycraft et al. (2019) the end of the Lower United Lake/Chalenko phase occurred between ∼12.4–11.8 ka, an inference based on one radiocarbon date reported by Turner et al. (2005) from a kettle hole in El Maitén (geographical location poorly constrained) and six cosmogenic exposure ages, originally reported by Glasser et al. (2012), that they reinterpreted as boulders that were shielded by the water column of the Lower United Lake. This minimum limiting age estimate is significantly older than our interpretation from Lago Maldonado which, as discussed in previous paragraphs, we consider robust and precise.

(8) We estimate water discharges of ∼150–200 km³ for the N5-N4 transition, 100–150 km³ for the N4-N3, and ∼300 km³ during the N3 drainage. The latter estimate assumes that LCP and LGCBA were connected during the Lower United phase. Our estimate is in broad agreement with Thorndycraft et al.’s (2019) Stage 7, with an estimated water volume of 312 km³ during the drainage of their Lago Chalenko stage.

(9) Benito and Thorndycraft (2020) interpreted catastrophic floods from the geomorphological and sedimentological evidence presented in the middle Valle Baker, such as gravel and sand bars, antidune landforms, and erosional evidence. Some of these sites are located downstream from Río Cochrane, at elevations ranging from ∼170–100 masl. We speculate that these floods might be related to transitions between the ephemeral N2 and N1 phases of GLC.

(10) Glasser et al.’s (2016) timing and structure of changes in glacial lake development in this area is incompatible with our findings in terms of chronology. Their stage A, potentially equivalent to our N5 phase, is dated at ∼13.8–11.5 ka on the basis of cosmogenic ages from Lago Tranquilo, Lago Negro, and the valley mouths from Valle Leones and Valle Chacabuco, originally published by Glasser et al. (2006) and Glasser et al. (2012). For reference, our N5 phase is constrained between ∼20.7–19.3 ka. Moreover, Glasser et al. (2016) place the termination of the Lower United Lake at 8.5 ± 0.9 ka, based on OSL samples from sites located around LGCBA. Our estimate for the end of N3, equivalent to Turner et al.’s (2005) Lower United Lake phase, is ∼11.3 ka (2σ range: 11.3–11.2 ka). We observe that their modeled effect of freshwater hosing from the PIS into the SE Pacific Ocean between ∼13.8–8.5 ka must be reassessed in the context of the new evidence reported in this study (timing, duration, and volume).

5.5 Climatic Implications

The pollen record from Lago Edita (Henríquez et al., 2017), along with the L. Unco and L. Mellizas sites (Vilanova et al., 2019; Villa-Martínez and Moreno, 2021) from the Coyhaique area, located ∼200 km north of LCP (Figure 1), afford valuable information for examining the climatic context in which the LCP deltas were built (Figure 12). These pollen sites indicate warming at ∼17.9 ka and an increase in precipitation brought by stronger SWW influence in pulses centered at ∼16.6 and ∼14.8 ka. Cool-temperate and hyper humid conditions persisted between ∼14.8–11.7 ka, spanning the Antarctic Cold Reversal (ACR) and Younger Dryas (YD) chonozones, followed by warming and decline in SWW influence during the early Holocene (∼11.7–9.4 ka). Based on these findings, we interpret that the large volume of the N3 deltas in LCP represents a combination of factors including a ∼3,500-year long period of GLC stability brought by cool-temperate conditions, enhanced sediment supply by increased precipitation, and profuse snow and ice melting driven by air temperatures higher than during the LGM and early T1. This climatic scenario differs from the cold and dry conditions that prevailed during the formation of the N5 deltas in LCP, constrained between ∼20.7–19.3 ka, and the relatively temperate, and intermediate-precipitation N4 phase constrained between ∼19.3–14.8 ka. We note that the end of the N3 phase at ∼11.3 ka occurred shortly after the warm pulse and decline in SWW influence that initiated the Holocene at pan Patagonian scale (∼11.7 ka) (Moreno et al., 2021).

FIGURE 12

FIGURE 12. Palynological record from Lago Mellizas (blue line), Lago Edita (black line), and Lago Unco (red line). Curves represent weighted moving averages. Black and blue dashed lines represent the Conifers and Nothofagus average during the last 10,000 years, respectively. The bottom panel shows the temperature variations recorded in the Dome C ice core, Antarctica. The time intervals N5, N4, and N3 are shown as shaded rectangles.

6 Conclusion

(1) GLC stabilized multiple times above the modern elevation of Lago Cochrane (∼150 masl) during the T1. We recognize five clusters of shoreline features at ∼600–500 (N5), ∼470–400 (N4), ∼360–300 (N3), ∼230–220 (N2), and ∼180–170 masl (N1). We estimate new glacio-isostatic gradients and generated isostatically corrected surfaces of the distinct phases of GLC. From these results we estimate uplifts of ∼150, ∼70, and ∼50 m for the N5, N4, and N3 deltas, respectively, in the Dos Arroyos sector.

(2) During the N5 phase GLC was connected to Valle Chacabuco through the Puesto Tejuela col and drained towards the Atlantic Ocean through the Caracoles col. Cosmogenic dates of 20.7 ± 1.3 ka from the innermost Río Blanco moraines afford a maximum limiting age constraint for the N5 phase. A minimum limiting age constraint comes from basal radiocarbon dates from Lago Edita and Lago Augusta, which afford an age of ∼19.3 ka.

(3) The surface elevation of GLC dropped ∼100 m after ∼19.3 ka, below the elevation of Caracoles and Puesto Tejuela cols. This event could represent an ice dam rupture or a lake level decline driven by thinning and recession of the damming ice front. The distribution of N4 landforms suggests that the ice dam was positioned at the westernmost end of LCP against Cerro Tamango. The Lago Augusta stratigraphic record suggests that the Chacabuco and Cochrane glacial lakes could have remained independent of each other at the beginning of the N4 phase and were eventually reconnected by drainage from the LCP basin, causing a transgressive glaciolacustrine phase in Valle Chacabuco. A subsequent regressive glaciolacustrine phase in Lago Augusta affords a minimum limiting age of ∼14.8 ka for the end of the N4 phase.

(4) The N3 landforms are widely distributed in the LCP basin and Valle Río Baker attesting to a major glacier retreat phase. We constrain its chronology with a maximum limiting age of ∼14.8 ka from Lago Augusta, two OSL dates that yielded a weighted mean age of 14.8 ± 1.1 ka for the construction of the N3 delta at Dos Arroyos, and a minimum limiting age of ∼11.3 ka from Lago Maldonado. Our N3 phase in the LCP basin represents the local expression of the Lower United Lake or Chalenko phase, made possible by recession of the Soler glacier from the Lago Bertrand moraines.

(5) We infer the opening of the lower Valle Río Baker and the irreversible establishment of the current drainage setting shortly before ∼11.3 ka. The subsequent N2 and N1 phases represent ephemeral stabilization events related probably to moraine collapses.

(6) We observe that the N3 landforms are the most ubiquitous, well-preserved, and largest morphologies in the study area, aspects that may represent an extended period of glacial lake stability, enhanced sediment supply by increased precipitation, and profuse snow and ice melting driven by temperate conditions between ∼14.8–11.3 a. This climatic scenario differs from the colder and drier conditions that prevailed during the N4 and N5 phases between ∼20.7–14.8 ka. The end of the N3 phase occurred shortly after a warm pulse and decline in precipitation at ∼11.7 ka that initiated the Holocene at western Patagonian scale.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

AV: Conceptualization, methodology, investigation, visualization, writing—original draft preparation. VF-A: Conceptualization, methodology, writing—review and editing. ES: Conceptualization, writing—review and editing. RH: Methodology. RV-M: writing—review and editing. PM: Conceptualization, methodology, formal analysis, visualization, writing—review and editing. JL: Methodology.

Funding

This work was funded by Iniciativa Científica Milenio NC120066, FONDECYT 1180815 (RV-M), and FONDECYT 1191942 (VF-A).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with one of the authors ES.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank D. Briones, M. Martini, L. Guerra, E. Fercovic, and V. Rivillo for providing field assistance, and A. Arancibia for GNSS data processing. Our appreciation goes to O. Pesce, W. Henríquez, and I. Jara for providing coring assistance. We thank T. Guilderson, R. De Pol, and J. Southon for their contributions to the development of radiocarbon dates from Lago Maldonado.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.817775/full#supplementary-material

References

Anderson, R. F., Ali, S., Bradtmiller, L. I., Nielsen, S. H. H., Fleisher, M. Q., Anderson, B. E., et al. (2009). Wind-Driven Upwelling in the Southern Ocean and the Deglacial Rise in Atmospheric CO 2. Science 323 (5920), 1443–1448. doi:10.1126/science.1167441

PubMed Abstract | CrossRef Full Text | Google Scholar

Balco, G., Stone, J. O., Lifton, N. A., and Dunai, T. J. (2008). A Complete and Easily Accessible Means of Calculating Surface Exposure Ages or Erosion Rates from 10Be and 26Al Measurements. Quat. Geochronol. 3 (3), 174–195. doi:10.1016/j.quageo.2007.12.001

CrossRef Full Text | Google Scholar

Bell, C. M. (2008). Punctuated Drainage of an Ice‐dammed Quaternary lake in Southern South america. Geografiska Annaler: Ser. A, Phys. Geogr. 90 (1), 1–17. doi:10.1111/j.1468-0459.2008.00330.x

CrossRef Full Text | Google Scholar

Bendle, J. M., Palmer, A. P., Thorndycraft, V. R., and Matthews, I. P. (2017). High-resolution Chronology for Deglaciation of the Patagonian Ice Sheet at Lago Buenos Aires (46.5°S) Revealed through Varve Chronology and Bayesian Age Modelling. Quat. Sci. Rev. 177, 314–339. doi:10.1016/j.quascirev.2017.10.013

CrossRef Full Text | Google Scholar

Benito, G., and Thorndycraft, V. R. (2020). Catastrophic Glacial-lake Outburst Flooding of the Patagonian Ice Sheet. Earth-Science Rev. 200, 102996. doi:10.1016/j.earscirev.2019.102996

CrossRef Full Text | Google Scholar

Blaauw, M., and Christen, J. A. (2011). Flexible Paleoclimate Age-Depth Models Using an Autoregressive Gamma Process. Bayesian Anal. 6 (3), 457–474. doi:10.1214/11-BA618

CrossRef Full Text | Google Scholar

Boex, J., Fogwill, C., Harrison, S., Glasser, N. F., Hein, A., Schnabel, C., et al. (2013). Rapid Thinning of the Late Pleistocene Patagonian Ice Sheet Followed Migration of the Southern Westerlies. Sci. Rep. 3 (1), 1–6. doi:10.1038/srep02118

PubMed Abstract | CrossRef Full Text | Google Scholar

Bourgois, J., Cisternas, M. E., Braucher, R., Bourlès, D., and Frutos, J. (2016). Geomorphic Records along the General Carrera (Chile)-Buenos Aires (Argentina) Glacial Lake (46°-48°S), Climate Inferences, and Glacial Rebound for the Past 7-9 Ka. J. Geology. 124 (1), 27–53. doi:10.1086/684252

CrossRef Full Text | Google Scholar

Caldenius, C. C. z. (1932). Las Glaciaciones Cuaternarias en la Patagonia y Tierra del Fuego. Geografiska Annaler 14 (1-2), 1–164. doi:10.2307/519583

CrossRef Full Text | Google Scholar

Clark, P. U., Dyke, A. S., Shakun, J. D., Carlson, A. E., Clark, J., Wohlfarth, B., et al. (2009). The Last Glacial Maximum. Science 325 (59411), 710–714. doi:10.1126/science.1172873

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, B. J., Darvill, C. M., Lovell, H., Bendle, J. M., Dowdeswell, J. A., Fabel, D., et al. (2020). The Evolution of the Patagonian Ice Sheet from 35 Ka to the Present Day (PATICE). Earth-Science Rev. 204, 103152. doi:10.1016/j.earscirev.2020.103152

CrossRef Full Text | Google Scholar

Davies, B. J., Thorndycraft, V. R., Fabel, D., and Martin, J. R. V. (2018). Asynchronous Glacier Dynamics during the Antarctic Cold Reversal in central Patagonia. Quat. Sci. Rev. 200, 287–312. doi:10.1016/j.quascirev.2018.09.025

CrossRef Full Text | Google Scholar

Denton, G. H., Anderson, R. F., Toggweiler, J. R., Edwards, R. L., Schaefer, J. M., and Putnam, A. E. (2010). The Last Glacial Termination. Science 328 (5986), 1652–1656. doi:10.1126/science.1184119

PubMed Abstract | CrossRef Full Text | Google Scholar

Douglass, D., Singer, B., Kaplan, M., Mickelson, D., and Caffee, M. (2006). Cosmogenic Nuclide Surface Exposure Dating of Boulders on Last-Glacial and Late-Glacial Moraines, Lago Buenos Aires, Argentina: Interpretive Strategies and Paleoclimate Implications. Quat. Geochronol. 1 (1), 43–58. doi:10.1016/j.quageo.2006.06.001

CrossRef Full Text | Google Scholar

Durcan, J. A., King, G. E., and Duller, G. A. T. (2015). DRAC: Dose Rate and Age Calculator for Trapped Charge Dating. Quat. Geochronol. 28, 54–61. doi:10.1016/j.quageo.2015.03.012

CrossRef Full Text | Google Scholar

Galbraith, R. F., Roberts, R. G., Laslett, G. M., Yoshida, H., and Olley, J. M. (1999). Optical Dating of Single and Multiple Grains of Quartz from Jinmium Rock Shelter, Northern Australia: Part I, Experimental Design and Statistical Models. Archaeometry 41 (2), 339–364. doi:10.1111/j.1475-4754.1999.tb00987.x

CrossRef Full Text | Google Scholar

Glasser, N. F., Harrison, S., Ivy-Ochs, S., Duller, G. A. T., and Kubik, P. W. (2006). Evidence from the Rio Bayo valley on the Extent of the North Patagonian Icefield during the Late Pleistocene-Holocene Transition. Quat. Res. 65 (1), 70–77. doi:10.1016/j.yqres.2005.09.002

CrossRef Full Text | Google Scholar

Glasser, N. F., Harrison, S., Schnabel, C., Fabel, D., and Jansson, K. N. (2012). Younger Dryas and Early Holocene Age Glacier Advances in Patagonia. Quat. Sci. Rev. 58, 7–17. doi:10.1016/j.quascirev.2012.10.011

CrossRef Full Text | Google Scholar

Glasser, N. F., Jansson, K. N., Duller, G. A. T., Singarayer, J., Holloway, M., and Harrison, S. (2016). Glacial lake Drainage in Patagonia (13-8 Kyr) and Response of the Adjacent Pacific Ocean. Sci. Rep. 6 (1), 1–7. doi:10.1038/srep21064

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, N. F., Jansson, K. N., Harrison, S., and Kleman, J. (2008). The Glacial Geomorphology and Pleistocene History of South America between 38°S and 56°S. Quat. Sci. Rev. 27 (3–4), 365–390. doi:10.1016/j.quascirev.2007.11.011

CrossRef Full Text | Google Scholar

Hein, A. S., Hulton, N. R. J., Dunai, T. J., Sugden, D. E., Kaplan, M. R., and Xu, S. (2010). The Chronology of the Last Glacial Maximum and Deglacial Events in central Argentine Patagonia. Quat. Sci. Rev. 29 (9–10), 1212–1227. doi:10.1016/j.quascirev.2010.01.020

CrossRef Full Text | Google Scholar

Heiri, O., Lotter, A. F., and Lemcke, G. (2001). Loss on Ignition as a Method for Estimating Organic and Carbonate Content in Sediments: Reproducibility and Comparability of Results. J. Paleolimnology 25, 101–110. doi:10.1023/A:1008119611481

CrossRef Full Text | Google Scholar

Henríquez, W. I., Villa-Martínez, R., Vilanova, I., de Pol-Holz, R., and Moreno, P. I. (2017). The Last Glacial Termination on the Eastern Flank of the central Patagonian Andes (47 ° S). Clim. Past 13 (7), 879–895. doi:10.5194/cp-13-879-2017

CrossRef Full Text | Google Scholar

Hogg, A. G., Heaton, T. J., Hua, Q., Palmer, J. G., Turney, C. S., Southon, J., et al. (2020). SHCal20 Southern Hemisphere Calibration, 0-55,000 Years Cal BP. Radiocarbon 62 (4), 759–778. doi:10.1017/RDC.2020.59

CrossRef Full Text | Google Scholar

Hubbard, A., Hein, A. S., Kaplan, M. R., Hulton, N. R. J., and Glasser, N. (2005). A Modelling Reconstruction of the Last Glacial Maximum Ice Sheet and its Deglaciation in the Vicinity of the Northern Patagonian Icefield, South America. Geografiska Annaler: Ser. A, Phys. Geogr. 87 (2), 375–391. doi:10.1111/j.0435-3676.2005.00264.x

CrossRef Full Text | Google Scholar

Kaplan, M. R., Strelin, J. A., Schaefer, J. M., Denton, G. H., Finkel, R. C., Schwartz, R., et al. (2011). In-situ Cosmogenic 10Be Production Rate at Lago Argentino, Patagonia: Implications for Late-Glacial Climate Chronology. Earth Planet. Sci. Lett. 309 (1–2), 21–32. doi:10.1016/j.epsl.2011.06.018

CrossRef Full Text | Google Scholar

Lal, D. (1991). Cosmic ray Labeling of Erosion Surfaces: In Situ Nuclide Production Rates and Erosion Models. Earth Planet. Sci. Lett. 104 (2–4), 424–439. doi:10.1016/0012-821X(91)90220-C

CrossRef Full Text | Google Scholar

Lamy, F., Kaiser, J., Arz, H. W., Hebbeln, D., Ninnemann, U., Timm, O., et al. (2007). Modulation of the Bipolar Seesaw in the Southeast Pacific during Termination 1. Earth Planet. Sci. Lett. 259 (3–4), 400–413. doi:10.1016/j.epsl.2007.04.040

CrossRef Full Text | Google Scholar

Lifton, N., Smart, D. F., and Shea, M. A. (2008). Scaling Time-Integrated In Situ Cosmogenic Nuclide Production Rates Using a Continuous Geomagnetic Model. Earth Planet. Sci. Lett. 268 (1–2), 190–201. doi:10.1016/j.epsl.2008.01.021

CrossRef Full Text | Google Scholar

Moreno, P. I., Henríquez, W. I., Pesce, O. H., Henríquez, C. A., Fletcher, M. S., Garreaud, R. D., et al. (2021). An Early Holocene westerly Minimum in the Southern Mid-latitudes. Quat. Sci. Rev. 251, 106730. doi:10.1016/j.quascirev.2020.106730

CrossRef Full Text | Google Scholar

Moreno, P. I., Videla, J., Valero-Garcés, B., Alloway, B. v., and Heusser, L. E. (2018). A Continuous Record of Vegetation, Fire-Regime and Climatic Changes in Northwestern Patagonia Spanning the Last 25,000 Years. Quat. Sci. Rev. 198, 15–36. doi:10.1016/j.quascirev.2018.08.013

CrossRef Full Text | Google Scholar

Murray, A. S., and Wintle, A. G. (2000). Luminescence Dating of Quartz Using an Improved Single-Aliquot Regenerative-Dose Protocol. Radiat. Measurements 32 (1), 57–73. doi:10.1016/S1350-4487(99)00253-X

CrossRef Full Text | Google Scholar

Murray, A. S., and Wintle, A. G. (2003). The Single Aliquot Regenerative Dose Protocol: Potential for Improvements in Reliability. Radiat. Measurements 37 (4–5), 377–381. doi:10.1016/S1350-4487(03)00053-2

CrossRef Full Text | Google Scholar

R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Available at: https://www.R-project.org/.

Google Scholar

Rabassa, J., and Clapperton, C. M. (1990). Quaternary Glaciations of the Southern Andes. Quat. Sci. Rev. 9 (2-3), 153–174. doi:10.1016/0277-3791(90)90016-4

CrossRef Full Text | Google Scholar

Stone, J. O. (2000). Air Pressure and Cosmogenic Isotope Production. J. Geophys. Res. 105 (B10), 23753–23759. doi:10.1029/2000JB900181

CrossRef Full Text | Google Scholar

Stuiver, M., Reimer, P. J., and Reimer, R. W. (2021). CALIB 8.2. [WWW program]. Available at: http://calib.org (Accessed 7 30, 2021).

Google Scholar

Thorndycraft, V. R., Bendle, J. M., Benito, G., Davies, B. J., Sancho, C., Palmer, A. P., et al. (2019). Glacial lake Evolution and Atlantic-Pacific Drainage Reversals during Deglaciation of the Patagonian Ice Sheet. Quat. Sci. Rev. 203, 102–127. doi:10.1016/j.quascirev.2018.10.036

CrossRef Full Text | Google Scholar

Toggweiler, J. R., Russell, J. L., and Carson, S. R. (2006). Midlatitude Westerlies, Atmospheric CO2, and Climate Change during the Ice Ages. Paleoceanography 21 (2). doi:10.1029/2005PA001154

CrossRef Full Text | Google Scholar

Turner, K. J., Fogwill, C. J., Mcculloch, R. D., and Sugden, D. E. (2005). Deglaciation of the Eastern Flank of the north Patagonian Icefield and Associated continental‐scale lake Diversions. Geografiska Annaler: Ser. A, Phys. Geogr. 87 (2), 363–374. doi:10.1111/j.0435-3676.2005.00263.x

CrossRef Full Text | Google Scholar

Vilanova, I., Moreno, P. I., Miranda, C. G., and Villa-Martínez, R. P. (2019). The Last Glacial Termination in the Coyhaique Sector of central Patagonia. Quat. Sci. Rev. 224, 105976. doi:10.1016/j.quascirev.2019.105976

CrossRef Full Text | Google Scholar

Villa-Martínez, R., and Moreno, P. I. (2021). Development and Resilience of Deciduous Nothofagus Forests since the Last Glacial Termination and Deglaciation of the central Patagonian Andes. Palaeogeogr. Palaeoclimatol. Palaeoecol. 574, 110459. doi:10.1016/j.palaeo.2021.110459

CrossRef Full Text | Google Scholar

Villa-Martínez, R., Moreno, P. I., and Valenzuela, M. A. (2012). Deglacial and Postglacial Vegetation Changes on the Eastern Slopes of the central Patagonian Andes (47°S). Quat. Sci. Rev. 32, 86–99. doi:10.1016/j.quascirev.2011.11.008

CrossRef Full Text | Google Scholar