1. Introduction
Lake terminating glaciers have been shown to retreat faster (and present higher variability) than land-terminating glaciers (King and others, Reference King, Dehecq, Quincey and Carrivick2018; Tsutaki and others, Reference Tsutaki2019; Sutherland and others, Reference Sutherland, Carrivick, Gandy, Shulmeister, Quincey and Cornford2020). This accelerated retreat reflects the influence of glacial lakes on ice dynamics (Tsutaki and others, Reference Tsutaki, Sugiyama, Nishimura and Funk2013, Reference Tsutaki2019), which can partially decouple glacier response to climate forcing (King and others, Reference King, Dehecq, Quincey and Carrivick2018, Reference King, Bhattacharya, Bhambri and Bolch2019). In Patagonia and other mountainous regions, glacial lakes are currently increasing in number and extent (Loriaux and Casassa, Reference Loriaux and Casassa2013; Wilson and others, Reference Wilson2018; Shugar and others, Reference Shugar2020). Consequently, understanding the sensitivity of glacier−lake systems to climate variability is relevant for the interpretation of past, present and future glacier fluctuations (Cuffey and Paterson, Reference Cuffey and Paterson2010).
Glacial lakes can affect ice dynamics through several processes. Most notably, lake development can reduce effective pressure and therefore increase surface velocity (Carrivick and Tweed, Reference Carrivick and Tweed2013; Tsutaki and others, Reference Tsutaki, Sugiyama, Nishimura and Funk2013). As lakes expand, the formation of a calving front introduces an additional mass-loss mechanism that can be comparable to surface ablation in its contribution to glacier (Boyce and others, Reference Boyce, Motyka and Truffer2007; Minowa and others, Reference Minowa, Schaefer, Sugiyama, Sakakibara and Skvarca2021). Beyond their impact on glacier dynamics, proglacial lakes influence downstream systems. They modulate catchment hydrology and function as sediment traps, altering sediment and nutrient delivery to rivers and fjords (Carrivick and Tweed, Reference Carrivick and Tweed2013; Milner and others, Reference Milner2017). These changes can affect proglacial ecosystems and reshape downstream watercourses (Carrivick and Heckmann, Reference Carrivick and Heckmann2017). Furthermore, glacial lakes pose high risk of glacial lake outburst floods (GLOFs), which can cause substantial socioeconomic impacts (Tweed and Russell, Reference Tweed and Russell1999; Iribarren Anacona and others, Reference Iribarren Anacona, Norton and Mackintosh2014; Carrivick and Tweed, Reference Carrivick and Tweed2016; Harrison and others, Reference Harrison2018; Taylor and others, Reference Taylor, Robinson, Dunning, Rachel Carr and Westoby2023). Although GLOF frequency may increase in the future, their magnitude could decrease (Dussaillant and others, Reference Dussaillant, Benito, Buytaert, Carling, Meier and Espinoza2010; Vandekerkhove and others, Reference Vandekerkhove, Bertrand, Crescenzi Lanna, Reid and Pantoja2020; Zheng and others, Reference Zheng2021; Veh and others, Reference Veh, Lützow, Kharlamova, Petrakov, Hugonnet and Korup2022).
Over the past decades, the Northern and Southern Patagonian Icefields have experienced increasing melt rates (Schaefer and others, Reference Schaefer, Machguth, Falvey and Casassa2013; Dussaillant and others, Reference Dussaillant, Berthier and Brun2018, Reference Dussaillant2019; Bravo and others, Reference Bravo, Bozkurt, Ross and Quincey2021; McDonnell and others, Reference McDonnell, Rupper and Forster2022). Extensive research has documented glacier thinning (Foresta and others, Reference Foresta2018; Dussaillant and others, Reference Dussaillant2019; Hugonnet and others, Reference Hugonnet2021; McDonnell and others, Reference McDonnell, Rupper and Forster2022) and recession (Davies and Glasser, Reference Davies and Glasser2012). These processes often coincide with observed increases in glacial lake area and volume (Wilson and others, Reference Wilson2018; Shugar and others, Reference Shugar2020). This is significant, as frontal ablation can account for up to 20% of the total mass loss for some lake-terminating glaciers in the Southern Patagonian Icefield, such as Viedma and Tyndall, and for most lake-terminating glaciers in the Northern Patagonian Icefield (Minowa and others, Reference Minowa, Schaefer, Sugiyama, Sakakibara and Skvarca2021). In addition, recent studies have shown that lake-terminating glaciers across the Northern and Southern Patagonian Icefields exhibited the most negative geodetic mass balances during 2000–20 (McDonnell and others, Reference McDonnell, Rupper and Forster2022), underscoring the need to improve our understanding of glacier–lake systems.
Climatic conditions across the Patagonian Icefields are shaped by the steep orographic gradient of the southern Andes, resulting in strong west–east contrasts in both precipitation and temperature. The western flanks of the Northern and Southern Patagonian Icefields lie within a hyperhumid zone dominated by the Southern Westerlies (Schaefer and others, Reference Schaefer, Machguth, Falvey and Casassa2013; Bianchi and others, Reference Bianchi, Villalba, Viale, Couvreux and Marticorena2016; Sauter, Reference Sauter2020). In contrast, the eastern flanks experience significantly lower precipitation and higher mean annual temperatures due to a pronounced rain shadow effect. These climatic gradients have a direct impact on glacier mass balance, with western outlet glaciers generally receiving higher snow accumulation (Bravo and others, Reference Bravo2019; Troch and others, Reference Troch, Åkesson, Cuzzone and Bertrand2024). Recent modeling and paleoclimatic reconstructions further show that precipitation has played a dominant role in modulating glacier variability in western Patagonia on multi-centennial timescales, underscoring the importance of precipitation as a key driver of glacier evolution in this region (Bertrand and others, Reference Bertrand, Hughen, Lamy, Stuut, Torrejón and Lange2012; Troch and others, Reference Troch, Åkesson, Cuzzone and Bertrand2024).
Despite these significant climatic and glaciological features, lake-terminating glaciers on the western flanks of the Northern and Southern Patagonian Icefield remain understudied compared to their leeward eastern counterparts. Research has focused mainly on the more accessible eastern outlets glaciers, including Upsala, Viedma and Perito Moreno in Southern Patagonian Icefield, and Leones and Nef in the Northern Patagonian Icefield (Sakakibara and Sugiyama, Reference Sakakibara and Sugiyama2014) and Leones and Nef in the Northern Patagonian Icefield (Warren and others, Reference Warren, Benn, Winchester and Harrison2001; Harrison and others, Reference Harrison2008). This disparity underscores the need for more research on western lake-terminating glaciers and their evolution of ice-to-lake in hyperhumid conditions.
Consequently, this study integrates lake characteristics and ice dynamics to identify the main controls on ice-to-lake evolution. We selected four glaciers located at the northwest flank of the Northern Patagonian Icefield: Gualas, Reicher, Grosse and Exploradores (Fig. 1). These glaciers were entirely or almost entirely covered by ice in 1945, as observed from aerial images. Therefore, we can narrow down the timeline of lake formation using aerial images from 1945 and 1976, along with more recent satellite imagery. Additionally, these glacier−lake systems exhibit distinct geomorphological characteristics, allowing for comparison of different stages and configurations of ice-to-lake evolution over the last 80–50 years. To describe their spatial and temporal evolution, we mapped and characterized the geomorphological and dynamical patterns of lakes and glaciers and compiled available information. In addition, we conducted the first bathymetry surveys of four lakes. By integrating these observations, we aim to identify key features of past glacier−lake responses and assess the main controls on their evolution.
Figure 1. (a) Trimetrogon aerial image (1945) east view: Gualas Glacier (1) and Reicher Glacier (2). (b) Trimetrogon image (1945) east view: Grosse Glacier and an ice-dammed lake (3) and Exploradores Glacier (4). (c) Bing aerial image composite with glacier extent (DGA, 2022). Airplanes icons show the approximate position and orientation of the aerial images. Orange lines correspond to glacier centerlines; orange points to locations where glacier surface velocity was extracted.
2. Study site
The glacier−lake systems of Gualas, Reicher, Grosse and Exploradores Glaciers are located north and northwest of the Northern Patagonian Icefield (Fig. 1). The Northern Patagonian Icefield extends north−south over 125 km (46.5º S to 47.5º S), occupied an area of 3675 km2 in 2016 (Meier and others, Reference Meier, Grießinger, Hochreuther and Braun2018) and is composed of 38 glaciers larger than 0.5 km2. Its elevation ranges from sea level to 4032 m a.s.l. at the peak of Mount San Valentin. The Patagonian Ice fields are significantly influenced by frontal systems associated with mid − latitude cyclones originating over the South Pacific Ocean. The study area is affected by the Southern Annular Mode and El Niño Southern Oscillation. Notably, over the second half of the 20th century the Southern Annular Mode has trended positively, which is associated with higher temperatures (Gillett and others, Reference Gillett, Kell and Jones2006), a strong increase in precipitations in the southern westerlies at higher latitudes and an increase in precipitation in the south of the Patagonian Icefields (Hall and Visbeck, Reference Hall and Visbeck2002). Locally, the Northern Patagonian Icefield is characterized by a west-east negative precipitation gradient.
The four studied glaciers have accumulation area in the flanks of Mount San Valentin and glacier terminus reaching lake elevation that span from 26 m a.s.l. (Gualas Lake) to 207 m a.s.l (Grosse lakes). All glacier termini are situated in a maritime high-precipitation climate, with mean annual precipitations on the Western slope of the Northern Patagonian Icefield averaging 5.83 ± 0.64 m a−1 and precipitation can reach to 15 m a−1 locally in the highest parts of the icefield (Schaefer and others, Reference Schaefer, Machguth, Falvey and Casassa2013; Bianchi and others, Reference Bianchi, Villalba, Viale, Couvreux and Marticorena2016; Sauter, Reference Sauter2020). Measured ablation rates at the terminus of Exploradores Glacier (170 m a.s.l.) can reach up to 20 m w.e. a−1 (Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022), but are typically on the order of 14 m w.e. a−1 (Willis and others, Reference Willis, Melkonian, Pritchard and Ramage2012; Schaefer and others, Reference Schaefer, Machguth, Falvey and Casassa2013). Meteorological data at the glacier termini show mean annual temperatures of 5–7ºC (Aguayo and others, Reference Aguayo2024). A description of each glacier along with the code on the Randolph Glacier Inventory 6.0 RGI60 and Chilean Water Directorate Inventory DGA − CL (DGA, 2022) is provided (Fig. 1).
Gualas Glacier: Gualas Glacier (RGI60-17.15809, CL111440130) has an area of 122.2 km2 and an elevation range of 46–4034 m a.s.l. The glacier is located at the western slopes of Mount San Valentín. The equilibrium line altitude (ELA) is reported at 1087 m a.s.l. (Rivera and others, Reference Rivera, Benham, Casassa, Bamber and Dowdeswell2007) and has a calving terminus in Gualas Lake (∼26 m a.s.l.). Intermediate recession moraines, formed after Little Ice Age (LIA), located at the periphery of Gualas Lake were dated to between 1879 and 1965 by dendrochronology (Harrison and Winchester, Reference Harrison and Winchester1998). This is coherent with trimetrogon images of 1945 that show the glacier occupying the current (2024) lake extent (Fig. 1a,c).
Reicher Glacier: Reicher Glacier (DGA-CL111440111) also spelled Reichert (RGI60-17.15815) has an area of 68.5 km2 and is located on the western slopes of Mount San Valentín. The ELA was not reported in (Rivera and others, Reference Rivera, Benham, Casassa, Bamber and Dowdeswell2007), yet average ELA for Northern Patagonian Icefield is 1190 m a.s.l. The LIA moraine at the south-western lake margin was dated to 1933 by dendrochronology (Harrison and Winchester, Reference Harrison and Winchester1998). In addition, on the 1945 trimetrogon aerial images, the glacier fills up most of the current Reicher Lake extent splitting at the lake center into a northern and southern branch and forming a T-shaped terminus. The lake is described as a fault-controlled lake (Harrison and Winchester, Reference Harrison and Winchester1998) as the lake orientation follows the Liquiñe−Ofqui fault system (Cembrano and others, Reference Cembrano, Hervé and Lavenu1996).
Grosse Glacier: Grosse Glacier (RGI60-17.15825, DGA-CL111421043) is a 63.9 km2 glacier located at the north slopes of Mount San Valentín (4032 m a.s.l.) and ELA at 1096 m a.s.l. (Rivera and others, Reference Rivera, Benham, Casassa, Bamber and Dowdeswell2007). The glacier ablation area is covered by debris and has a calving front terminus in Grosse Lake. The debris-covered area shows semi-circular ice cliffs (Fig. 1b). In addition, the 1945 aerial images show an ice-dammed lake (with an estimated area of 1.8 km2), which caused a glacial lake outburst flood in 1957 (Gorsic and others, Reference Gorsic2025).
Exploradores Glacier: Exploradores Glacier (RGI60-17.15831, DGA-CL111421023) is an 83.4 km2 glacier located at the north slopes of Monte San Valentín with an ELA of 1187 m a.s.l. (Rivera and others, Reference Rivera, Benham, Casassa, Bamber and Dowdeswell2007). Previous research has suggested that Exploradores Glacier is at a tipping point: the formation of the glacial lake will shift the glacier’s current low downwasting rate at the glacier front (<0.7 m a−1) to a fast-retreating phase (Aniya and others, Reference Aniya, Barcaza and Iwasaki2007; Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022). It has a prominent Little Ice Age moraine (Aniya and others, Reference Aniya, Barcaza and Iwasaki2007) and ice-thickness models suggest a maximum bed overdeepening of 388 m below current lake levels (Fürst and others, Reference Fürst2024). Notably, Exploradores Glacier suffered two GLOFs in the last decade. On 2015 a tributary valley El Chileno produced a GLOF that entered the main trunk of Exploradores Glacier (Wilson and others, Reference Wilson2019). The GLOFs entered the east-marginal lake and currently a river discharges its water to the lake. Later in 2018, the tributary valley of Glaciar Bayo produced the sudden drainage of Triángulo Lake with the flood also entering the main trunk (Bañales-Seguel and others, Reference Bañales-Seguel, Salazar and Mao2020).
3. Materials and methods
3.1. Ice front outlines, glacier surface elevation and surface velocity
The ice front positions and lake outlines were manually mapped from satellite images using the Google Earth Digitalization Tool GEEDiT (Lea, Reference Lea2018). Glacier front retreat rates were computed along the glacier’s centerline and statistically significant changes were identified by change–point detection algorithm (Killick and others, Reference Killick, Fearnhead and Eckley2012) similar to applications in other glacier–lake systems (Hill and others, Reference Hill, Carr, Stokes and Gudmundsson2018; Dell and others, Reference Dell, Carr, Phillips and Russell2019). The data set includes Landsat imagery starting in 1985, ASTER and Sentinel imagery until 2023. In addition, the mapping was extended by including aerial images from Hexagon flights in 1976 (McDonnell and others, Reference McDonnell, Rupper and Forster2022) and trimetrogon oblique photographs acquired by the U.S. Air Force in 1945. Glacier surface elevation was obtained from (McDonnell and others, Reference McDonnell, Rupper and Forster2022) for 1976 and from Hugonnet and others (Reference Hugonnet2021) for years 2000, 2005, 2010, 2015 and 2020. The surface elevation models were obtained by adding the observed reported change rates (dh/dt) to a reference DEM of 2000. Glacier bed elevation models were obtained from ice-thickness models from Fürst and others (Reference Fürst2024). Surface velocity datasets were extracted from ITS_LIVE (Gardner and others, Reference Gardner, Fahnestock and Scambos2022) and velocity trends inspected annually by change–point detection algorithms (Killick and others, Reference Killick, Fearnhead and Eckley2012). Lastly, we computed the maximum depth at the glacier front (Hmax) by intersecting the glacier perimeter with the lake bathymetry and extracting the maximum value.
3.2. Lake bathymetry and subaqueous ice volume
Lake bathymetry surveys were performed during the year 2023 using two sonar systems a Hummingbird Helix 7 CHIRP and a Garmin GPSMAP 527xs mounted on kayaks (Exploradores, Grosse and Reicher Lakes) and on a motorboat (Gualas Lake). The sonars were configured to sample at a frequency of 1 s. Horizontal coordinates were acquired by a built-in GPS and are expected in the order of 5–10 m according to the manufacturer.
Bathymetry point measurements for Exploradores, Gualas and Reicher lakes were interpolated in MATLAB using scattered-data linear interpolation, producing a continuous lake-floor surface constrained by the 2023 lake perimeter. The surface was sampled onto a regular 50 m grid for analysis and volume calculations, a resolution chosen to accommodate the non-homogeneous distribution of survey tracks across the lakes and to facilitate integration with other variables, such as glacier extent mapped at finer scales.
The bathymetric survey at Grosse Lake was less dense, covering a single track along the lake’s main axis, a secondary track along its side, and a perpendicular section parallel to the calving front, located 600 m away. Due to data sparsity, traditional interpolation methods were unsuitable, and linear interpolation tends to produce a narrow, V-shaped cross-section, which contrasts with the broader U-shaped geometry typically associated with glacial overdeepening and observed near the calving front. Given the lake’s symmetry along its major axis and the distinct U-shaped morphology observed in bathymetric profiles, we modeled the bathymetry by fitting a polynomial function to depth data and the lake perimeter, ensuring the model represented the U-shaped morphology to approximate lake volume and structure (Supplementary material, Fig. S1). This procedure was chosen only in the absence of denser data and specific lake geometry. Given this, the polynomial interpolation residuals were −0.5 m mean bias and ±12 m standard deviation (Supplementary material, Fig. S1). For Grosse Lake a sensitivity analysis was done with linear, cubic and nearest-neighbor interpolation methods. The U-shaped valley geometry (polynomial fit) yielded systematically larger lake volumes (20−35%) than a V-shaped geometry produced by linear interpolation. In contrast, maximum lake depth at the front is less sensitive to the interpolation method, as the survey passes near the lake centerline and is expected to capture the deepest part of the basin. Overall, given the distribution of survey transects it is not possible to fully assess the interpolation bias.
Lastly, subaqueous ice loss, corresponding to the portion of the glacier submerged within the lake basin, is estimated from changes in lake volume through time. Lake volume is calculated by integrating the difference between lake level and the interpolated bathymetry on a regular 50 m grid, assuming that the glacier front observed at the surface extends vertically downward to the lakebed. The resulting volume change is converted to mass using an ice density ρi = 913 kg m−3 to obtain the subaqueous ice loss. This assumption provides the submerged ice loss in an area where direct observation is nearly impossible (Fig. 2).
Figure 2. Proglacial lake formation. White hatched area shows elevation loss measurable from geodetic studies. Prograde and retrograde slopes are indicated in the upper plot. Light blue-hatched area indicates unaccounted subaqueous ice volume loss in geodetic mass balance studies. Acronyms h refer to elevation variables with the following subscripts: hdem (surface elevation), hw (water level), hb (lakebed elevation) and hab height above flotation.
3.3. Uncrewed aerial vehicle surveys
We deployed an uncrewed aerial vehicle (UAV), specifically a DJI Mavic 3 Enterprise equipped with differential GNSS, to survey the fronts of Gualas (19 March 2024) and Exploradores Glacier (20 April 2024 and 17 November 2024). UAV flight plans were designed to achieve a resolution of 20 cm px−1 for orthomosaics. We established a fixed GNSS base for each glacier to facilitate repeated surveys on different dates. The image processing workflow involved updating image location coordinates by postprocessing kinetics in Emlid Studio. Subsequently, the images were processed photogrammetrically using Structure from Motion-Multiview Stereo workflow via Agisoft Metahshape Professional 1.8 to obtain a digital surface model and orthomosaic. To compute surface velocity for Exploradores Glacier, we co-register the surveys using virtual tie-points on stable terrain from a master survey (Benoit and others, Reference Benoit2019; Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022), then resample the orthomosaic at 1 m px−1 resolution and applied the feature-tracking algorithm ImCORR in SAGA (Scambos and others, Reference Scambos, Dutkiewicz, Wilson and Bindschadler1992). Supraglacial ponds on Exploradores Glacier were identified from the orthomosaic by segmentation using Orfeo ToolBox in QGIS, then validated and classified manually.
3.4. Glacier terminus buoyant conditions
To inspect the glacier tongue buoyant condition or flotation condition (Boyce and others, Reference Boyce, Motyka and Truffer2007) follow a similar approach to (Minowa and others, Reference Minowa, Schaefer and Skvarca2023). We define the flotation height (hf) as the elevation when ice overburden pressure equals buoyancy force,
\begin{equation}{h_f} = {h_b} + \left( {\frac{{{\rho _w}}}{{{\rho _i}}}} \right)\left( {{h_w} - {h_b}} \right)\end{equation}where hb is the bed elevation, hw is the lake water-level elevation, assuming ρw = 1000 kg m−3 and ρi = 913 kg m−3 the density of lake water and glacier ice, respectively.
Then, we calculated the glacier’s buoyant condition (hab) by the difference of the glacier surface elevation (hDEM) to flotation height:
\begin{equation}{h_{ab}} = {h_{DEM}} - {h_f}\end{equation}Consequently, if hab > 0 the glacier is grounded (upper panel in Fig. 2), and hab < 0, the glacier is in buoyant condition (lower panel in Fig. 2).
Buoyancy uncertainty was estimated by propagating errors in hab which estimates to ± 10 m. A vertical uncertainty of ± 5 m was adopted for the DEM, consistent with near-flat terrain (Hugonnet and others, Reference Hugonnet2021). Uncertainties of ± 15 m for bathymetry and ± 5 m for water level were assumed as conservative. Note that for Grosse Glacier, the bathymetry is well constrained in the area used for buoyant conditions analysis, supported by crossing survey tracks, with residuals ± 12 m (Supplementary material, Fig. S1, S2). These uncertainty ranges are comparable to those reported by Minowa and others, (Reference Minowa, Schaefer and Skvarca2023; ∼ 14 m). Buoyancy estimates for Exploradores Glacier are based on bed elevation models. Near the eastern lake front, bed elevation is likely underestimated, as available bathymetry indicates greater depth. Extrapolations across the remaining glacier tongue should therefore be interpreted with caution.
4. Results
4.1. Lake bathymetry
Bathymetric surveys of the four proglacial lakes provide first-order constraints on lake geometry and basin overdeepening, which are essential for assessing buoyant conditions, estimating subaqueous ice volume, and for enabling a comparative evaluation of glacier–lake and glacier dynamic responses (Minowa and others, Reference Minowa, Schaefer and Skvarca2023; Zhang and others, Reference Zhang2023).
Gualas Lake, in November 2023, had an area of 4.4 km2 and a maximum length of 3 km. The bathymetry survey of Gualas lake revealed a maximum depth of 272 m. The bathymetric surveys consistently indicate the shallowing of the lake depth towards the glacier front (Figs. 3a and 4). Glacier bed reconstruction (Fürst and others, Reference Fürst2024) overestimated the maximum lake depth by 250 m, as well as the current glacier bed elevation at the glacier front in 2023.
Figure 4. Gualas glacier−lake evolution. (a) Glacier elevation profiles, glacier bed elevation, lake bathymetry and lake level. Blue vertical line indicates lake depth at the plotted location. Squares mark the glacier front position for each DEM year, with line colors matching the corresponding year. (b) Surface velocity (mean annual in black dots), maximum lake depth at the glacier front (Hmax, from lake bathymetry), lake area and relative front position. (c) Lake bathymetry, glacier extent and glacier buoyant conditions.
Reicher Lake is a T-shaped lake, 8.6 km long, 1 km wide and covers 8.7 km2 as of 2023. The southeastern branch is notably larger and deeper, with the bathymetry survey indicating a maximum depth of 347 m (Figs. 3b and 5). The lake depth decreases as it approaches the current glacier calving front. The glacier bed model (Fürst and others, Reference Fürst2024) correctly replicated the slopes of the glacier bed (lake bottom), but it resulted in an offset of ∼120 m shallower than the measured lake bottom.
Figure 5. Reicher glacier−lake evolution. (a) Glacier elevation profiles, glacier bed elevation, lake bathymetry and lake level. Blue vertical line indicates lake depth at the plotted location. Squares mark the glacier front position for each DEM year, with line colors matching the corresponding year. (b) Surface velocity, maximum lake depth at the glacier front (Hmax, from lake’s bathymetry), and relative front position. (c) Lake bathymetry, glacier extent and glacier buoyant conditions. Year 1976 had no data, and after 2000 no bathymetry is available.
Grosse Lake is situated between steep valley walls, flanked by prominent lateral moraines (∼240 m above the lake surface). In 2023, the lake measured 4 km in length and 1 km in width. The bathymetric survey was conducted along the past glacier center line, with a quasi-perpendicular cross-section ∼600 m from the calving front. Due to adverse weather conditions, the survey was incomplete. However, the bathymetry along the center line indicated a smooth deepening of the lake towards the calving front, reaching a maximum depth of 206 m at the point closest to the calving front (Figs. 3c and 6). The perpendicular-cross section revealed a U-shaped lake morphology (Supplementary material, Fig. S1). The lake bottom slope trend shows that the lake deepens towards the calving front, consistent with bed elevation models that suggest that the bed over deepening extends for more than ∼ 3.5 km from the 2023 glacier front position (Fürst and others, Reference Fürst2024).
Figure 6. Grosse glacier−lake evolution. (a) Glacier elevation profiles, glacier bed elevation, lake bathymetry and lake level. Blue vertical line indicates lake depth at the plotted location. Squares mark the glacier front position for each DEM year, with line colors matching the corresponding year. (b) Surface velocity, maximum lake depth at the glacier front (Hmax, from lake’s bathymetry), and relative front position. (c) Lake bathymetry, glacier extent and glacier buoyant conditions.
Lastly, the bathymetry survey at Exploradores Glacier, conducted in January 2023, encompassed most of the lake’s extent at that time. However, the lake only extends 500 m from the eastern margin towards the glacier’s centerline within the 3 km wide glacier. Despite the limited scope of lake bottom exploration, the results indicate a deepening lake bottom, reaching a maximum depth of 333 m at a distance of 1 km from the centerline. The glacier bed elevation model estimates a maximum lake depth of 338 m at the glacier center line of the glacier (Fürst and others, Reference Fürst2024), suggesting that the glacier bed model may underestimate the future lake’s maximum depth (Figs. 3d and 7).
Figure 7. Exploradores glacier−lake evolution. (a) Glacier elevation profiles, glacier bed elevation, lake bathymetry and lake level. Blue vertical line indicates lake depth at the plotted location. Squares mark the glacier front position for each DEM year, with line colors matching the corresponding year. (b) Surface velocity, maximum lake depth at the glacier front (Hmax, from lake’s bathymetry), and relative front position. (c) Lake bathymetry, glacier extent and glacier buoyant conditions.
Across the four sites, bathymetric surveys reveal pronounced overdeepening, with maximum depths ranging from ∼200 m at Grosse Lake to >300 m at Reicher and Exploradores Lakes. Surveys at Exploradores and Grosse were spatially limited and bed slopes continued to deepen toward the glacier fronts, suggesting that maximum depths may be underestimated. Comparisons with glacier bed elevation models show depth offsets near the glacier fronts despite broadly consistent bed-slope trends. Overall, these results demonstrate discrepancies between modeled and measured lake basins and underscore the importance of bathymetry surveys to constrain proglacial lake geometry. Together, the bathymetries and glacier extents show that the glaciers are at different stages of retreat and lake formation, a distinction that forms the basis for the process-based interpretation developed in the next section.
4.2. Glacier dynamics changes (front retreat, velocity and elevation changes)
4.2.1. Gualas glacier
The first aerial images from 1945 show Gualas Glacier covering the entire extent of the current lake. Hexagon images from 1976 reveal the initial formation of Gualas lake at the outlet, with an area of 0.065 km2 (Fig. 3). Buoyancy analysis based on bathymetry and surface elevation data suggest that the glacier was grounded in the deepest parts of the lake in 1976 (Fig. 4c). However, the 1976 DEM does not cover the full extent of the terminal tongue, and grounding conditions cannot be confirmed along the entire front. Satellite images from 1999 show increased crevassing and fragmentation of the glacier surface. Large icebergs up to 100 × 150 m, are visible floating on the proglacial lake. Between 1976 and 2000, the glacier transitioned from a stable front and homogeneous surface to a crevassed and fragmented partly floating tongue.
The early 2000s has marked a shift in Gualas Glacier dynamics, with much of the tongue entering floating conditions (Fig. 4c). Between 2000 and 2007 Gualas Glacier underwent rapid frontal retreat (420 m a−1, Supplementary Material, Fig. S3), withdrawing from the deepest parts of the basin in 2005 and expanding the lake area by 2.1 km2 (Fig. 4b). Surface elevation model and lake bathymetry indicate that the glacier tongue was floating in 2000, and by 2005 hab values had decreased further, preceding disintegration. At that time, satellite images show the lake densely covered with irregularly sized icebergs (Supplementary material, Fig. S8). After 2005, the glacier front retreated from a reverse to a prograde bed slope. By 2006–07, the lake attained an extent similar to that of 2024. The velocity trends for Gualas indicates a declining trend. Velocity data for 1986 averaged 197 m a−1, while in 2000 velocities averaged 186 m a−1, the highest values observed. Velocity trend then declined, reaching a mean of 65 m a−1 by 2024. Glacier dynamics along the first 12 km of centerline profiles (2000–22) show a stronger velocity decline toward the terminus and reduced surface lowering in the same sector (Supplementary material, Fig. 12).
In November 2023, the UAV survey reveal bedrock outcrops extending ∼700 m beneath the glacier front, with only ∼300 m of glacier still in contact with the lake. This indicates that the glacier is close to disconnecting from the lake and reverting to land-terminating (Supplementary material, Fig. S7).
4.2.2. Reicher glacier
Aerial images from 1945 show that Reicher Glacier had a T-shaped terminus, covering most of the present-day lake area. Lake development began at the southern branch, near the lake outlet. The glacier flowed northwest and split 90° into northeastern and southwestern branches (Fig. 1). By 1976, aerial images showed Reicher Lake divided into two sections, with active outlets to the north and south. The northern lake section covered 1.3 km2 and contained tabular icebergs measuring up to ∼ 350 m in length. The southern branch formed a 0.8 km2 lake between the glacier and the northern moraine margin.
In 1985 and 1987, Reicher Glacier still split the lake in two parts, with the terminal tongue in contact with the opposite shore. A well-defined calving front with small icebergs developed in the north, while the southern lake hosted larger tabular icebergs. Between 1991 and 1994, the southern branch of Reicher Glacier retreated by 3.7 km, expanding the lake to 2.21 km2. From 1998 and 2000, satellite images show that the ice front continues dividing the lake into two sections, a condition that persisted until September 2001, when retreat unified the lake.
Between the 1986 and the early 2000s, Reicher Glacier experienced rapid dynamic changes. Surface elevation data for 2000 indicate that the terminal tongue was in floating conditions (Fig. 5c). After 2001, the ice front retreated 1.2 km, forming a coherent ice front perpendicular to the main flow direction. Elongated tabular icebergs, up to 1100 × 380 m drifted towards shallower parts of the lake. In 2002, a landslide occurred on the southern margin, less than 1 km from glacier front. During 2003–04, retreat continued rapidly. By March 2004, the glacier front had reached a position close to that observed in 2023, resting on a prograde lake floor (Supplementary material, Fig. S9, 2006).
From late 2004, the glacier entered a slow retreat phase, with rates decreasing from ∼332 ma−1 to ∼8 m a−1 (Supplementary material, Fig. S4). Over the next 20 years, the glacier receded 400 m. Surface velocities declined from ∼180 m a−1 in the early 2000s to 120 m a−1 in 2023. A short-lived acceleration occurred in 2015, followed by deceleration, but this did not produce significant frontal advance. After 2023, velocities increased again, and the front advanced ∼50 m.
4.2.3. Grosse glacier
In 1945, Grosse Glacier extended to the Little Ice Age moraine, entirely covering the present-day lake. The glacier surface was covered with debris and showed semi-circular ice cliffs (Fig. 1b). Aerial images also show an ice-dammed lake, which caused a GLOF in the 1950s (Gorsic and others, Reference Gorsic2025). By 1976, Hexagon images reveal no evidence of the lake. The 1976 DEM does not cover the first 2.3 km of glacier front, and only a 1 km section intersects with lake bathymetry. This indicates that the glacier was not afloat in that section, but conditions at the first kilometer remain unknown.
The first evidence of lake formation appears in 1986 Landsat images. The glacier tongue was still connected to the main glacier body, but small supraglacial ponds disrupted surface (Supplementary material, Fig. S10). From this image, funnel-shaped and crescent-shaped ice cliffs typical of debris-covered glaciers can be observed (Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021; Kneib and others, Reference Kneib2023). Between 1987 and 1999, these ponds expanded and coalesced. By 1999, a small proglacial lake had formed, while supraglacial ponds persisted on the debris-covered tongue. A 2000 Aster image reveals that supraglacial ponds are connecting with the lake and the glacier front is becoming increasingly fragmented. Buoyancy analysis indicates that the first hundred meters of the glacier front was in floating conditions. The next 1.5 km consisted of patchy debris-covered ice that disintegrated into isolated blocks, likely stagnant ice.
In 2003 and 2004, the proglacial lake expanded further, and the glacier terminus became increasingly disrupted with numerous ponds. By 2008, the tongue showed severe fragmentation, and the front appeared disconnected from stagnant dead ice. This is the first evidence of a calving front perpendicular to the centerline (Supplementary material, Fig. S10). Buoyancy maps from 2000, 2005, 2010 and 2015 show consistent conditions: a narrow section near the centerline afloat but pinned laterally against the valley walls.
By 2016, the proglacial lakes had expanded and most of the dead ice had melted. In 2020, satellite images depict a larger proglacial lake, and a coherent calving front perpendicular to the main flow. The most recent images (2024) show that Grosse Glacier’s front has continued to retreat, with a well-defined calving front. Icebergs are present but relatively small. Throughout the entire period of lake development, the glacier front has remained on a reverse bed slope (retrograde slope).
4.2.4. Exploradores glacier
Oblique aerial images from 1945 provide the first evidence of proglacial lake development at the northwest margin. Additionally, a few supraglacial ponds were observed near the glacier terminus. At that time, the glacier had two active outlets in the northwest and northeast, with the latter still active today. In 1976, Hexagon images depicted the northwest lake with an area of 0.37 km2. Ice−cliffs were observed within the first 300 m of the glacier, though the presence of supraglacial ponds is unclear from the image. Since then, the west marginal lake has continued to expand southward between the lateral moraine and the glacier to the present day. At the east margin, lake development began later, where the first evidence appeared in 2005 (Aniya and others, Reference Aniya, Barcaza and Iwasaki2007).
Buoyant conditions were calculated using glacier bed models (Fürst and others, Reference Fürst2024) due to limited bathymetric data. Uncertainties in these models are large. Fürst and others (Reference Fürst2024) reports errors up to 50− 150 m per pixel in the lower section of Exploradores Glacier. Bathymetry from the east marginal lake measured a maximum depth of 330 m, ∼1 km from the centerline, where bed model is shallower. At the centerline, bed models indicate a maximum depth of 338 m. With low confidence, this suggests that bed models may underestimate the depth of the overdeepening near the east marginal lake. Taking this into consideration, results indicate the glacier was grounded in 1976. By 2000, parts of the front were buoyant, and this zone expanded as thinning continued. These interpretations remain limited by uncertainties in the bed models.
For Exploradores Glacier, the front position is not shown because the front has not retreated along the centerline. However, glacier retreat can be inspected by the increase in area of the east marginal lake. By 2015, it covered 1 km2, with a calving front parallel to the glacier’s main flow direction. In October 2023, a large calving event produced tabular icebergs up to 500 × 250 m. By December 2023, all lakes and supraglacial ponds together covered 4.14 km2. In May 2024, the east marginal lake expanded to 2.5 km2, primarily through calving.
From 2000 to 2017 the velocity remained between 100 to 150 ma−1. After 2014, velocities increased, peaking in 2015 before decreasing in 2017. Change point algorithm detected an increase in velocity trend at the end of 2016 and another increase at the end of 2022 reaching maximum values at the end of the time series.
Two GLOFs affected the Exploradores system in recent years. In December 2015, a GLOF occurred in El Chileno Valley, a tributary of Exploradores Glacier (Wilson and others, Reference Wilson2019). Velocities had already increased since 2014, peaked in 2015, and declined afterward. Wilson and others (Reference Wilson2018) suggested that part of the drained volume may have originated from the east marginal lake of Exploradores, although no direct evidence of lake level change was reported and/or lake expansion associated with the event. In 2018, another GLOF occurred from Laguna Triángulo at Bayo Glacier, a tributary from Exploradores Glacier (Bañales-Seguel and others, Reference Bañales-Seguel, Salazar and Mao2020). At Exploradores, no immediate expansion of the marginal lake was observed following the event; however, later that year calving activity increased and contributed to lake enlargement (Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022). Although GLOFs are known to influence glacier dynamics (Huss and others, Reference Huss, Bauder, Werder, Funk and Hock2007), in this case their specific impact on dynamics or lake development cannot be confirmed.
Two UAV surveys were conducted on 20 April 2024 and 17 November 2024 covering an area of ∼13 km2 (Fig. 8). The first 2 km of the glacier lies <25 m above lake level. Marginal lakes extend >4 km along the lateral moraines, while supraglacial ponds are concentrated in the first 1.5 km. Elevation profile A−A’ shows a crevassed area near the centerline (black arrow in Fig. 8), with water levels matching the lake. The digital surface model (DSM) colorbar was stretched to 175–200 m a.s.l. to highlight the elevation difference with the lake level (∼175 m a.s.l.). Elevation inside ponds and lakes may contain artifacts due to poor texture matching, but the DSM was left uncropped to display iceberg presence. Image pair velocities resulted in maximum velocities of 180 m a−1. Errors from stable areas showed biases and standard deviations of 0.2 ± 0.5 m on the east axis and − 0.4 ± 0.4 m on the north axis. Histograms and spatial distribution of the error are provided in the Supplementary material, Fig. S11. Notably, the velocity east and west of the crevassed area (back arrow) differs in both magnitude and direction.
Figure 8. (a) DSM (17 November 2024), marginal lakes and supraglacial ponds (black line polygons) and glacier contact with moraines (thick black line). (b) Orthomosaic (17 November 2024), surficial velocity from 30 April 2024 to 17 November 2024 (blue to yellow colors) and velocity direction (black lines). (c) Elevation profile (A−A’) on (a), Blue band indicates lake level.
Inspection of the centerline velocities and elevation changes from 2000 to 2020 reveals an inflection point at ∼ 4 km from the terminus. Downstream of this point, velocities for 2020−22 exceed the 2000–22 mean, whereas upstream (4–7 km) remain below average. Elevation changes display the opposite pattern, with rates decreasing toward the terminus and increasing between 4 and 7 km (Supplementary Fig. S12, S13). Bed elevation maps indicate a topographic high near the 4 km mark, which may influence glacier dynamics (Minowa and others, Reference Minowa, Schaefer and Skvarca2023). Previous work similarly suggests dynamic thinning at Exploradores Glacier, with pronounced surface lowering between 4 and 6 km from the terminus (Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022).
4.2.5. Comparative glacier-lake evolution
A comparative assessment of the four glacier–lake systems reveals partial associations between front retreat rates, changes in maximum lake depth at the glacier front (Hmax), surface velocity trends, and lake area growth. Despite differences in glacier geometry, similar changes in the trends of individual variables can be observed across sites in accordance to lake development (Supplementary material Fig. S3-S6 and Table S1). Change–point detection identified intervals with coincident trend changes across multiple variables.
Periods of limited lake influence are characterized by relatively low retreat rates (∼50–85 m a−1), small lake area growth (<0.11 km2 a−1), and modest negative Hmax changes (−1.8 to −12 m a−1), indicating gradual deepening of glacier front. Velocity changes during these intervals are small or weakly negative where data are available.
In contrast, intervals of rapid glacier–lake interaction show a marked increase in system activity across all sites. During these periods, retreat rates increase to > 300 m a−1 at Gualas and Reicher and exceed 100 m a−1 at Grosse, accompanied by accelerated lake expansion (0.04–0.42 km2 a−1). These changes coincide with large-magnitude variations in Hmax, including rapid deepening rates of up to −65 m a−1 at Gualas and −71 m a−1 at Exploradores, reflecting glacier retreat into deeper lake basins. Velocity responses during these intervals vary among glaciers, with recent acceleration observed at Exploradores (up to 12 m a−2) and Grosse (0.6–0.9 m a−2), while Gualas and Reicher exhibit declining velocity trends following peak retreat.
At Gualas and Reicher, subsequent intervals are characterized by a sharp reduction in retreat rates (<10 m a−1), minimal lake area growth (∼0.01 to 0.02 km2 a−1) together with consistently negative velocity changes (−2 to −4 m a−2). Such behavior is not observed at Grosse and Exploradores, where retreat rates remain elevated, lake area continues to increase, and velocity trends remain neutral or positive.
5. Discussion
5.1. Glacier dynamics during terminus transition from land to lake (to land)
The transition from land-terminating to lake-terminating glacier involves marked shifts in ice dynamics (Tsutaki and others, Reference Tsutaki, Sugiyama, Nishimura and Funk2013; Tsutaki and others, Reference Tsutaki2019). However, direct observations of subaqueous glacier–lake processes are inherently limited (Truffer and Motyka, Reference Truffer and Motyka2016). Here we relate changes in glacier dynamics based on observed lake and glacier characteristics. We first use empirical evidence from bathymetry, terminus position and buoyancy conditions to interpret controlling mechanisms on the glacier–lake systems. We then examine how these observations relate to changes in glacier state and rates, including retreat, surface velocity, and lake depth at the glacier front. Processes not directly constrained by the available data, such as lake thermal structure and upstream, glacier-wide dynamics, are not explicitly considered (see Section 5.4, Limitations and outlook).
Based on the four study sites, we identify three characteristic phases: (i) initial lake formation, (ii) glacier retreat along retrograde slopes and into deeper lake area, and (iii) retreat over prograde bedrock and disconnection (Fig. 2). While broadly consistent with conceptual models proposed in earlier studies (Kirkbride, Reference Kirkbride1993; Tsutaki and others, Reference Tsutaki2019; Minowa and others, Reference Minowa, Schaefer and Skvarca2023), our analysis integrates observations from western Patagonia and highlights the distinct patterns of transitional behaviors across neighboring glaciers. In the following table, we summarize the main processes and observations associated with each phase Table 1. Note that proposed dates indicate when observable changes occurred based on available data and do not necessarily represent abrupt phase transitions.
Table 1. Main observations during land to lake transition.
5.1.1. Lake initiation
The initial lake formation is characterized by slow lake area increase rates. We identify two patterns attributed to the presence or absence of a debris cover on the glacier tongue. Lake initiation in debris-covered tongue of Grosse and sections of Exploradores occurred through the coalescence of supraglacial ponds associated with ice cliffs backwasting. The development of supraglacial ponds which can enhance melt rates (Reid and others, Reference Reid, Carenzo, Pellicciotti and Brock2012; Miles and others, Reference Miles2017; Steiner and others, Reference Steiner, Buri, Miles, Ragettli and Pellicciotti2019; Buri and others, Reference Buri, Miles, Steiner, Ragettli and Pellicciotti2021) is attributed to low-gradient (less than 2º) tongues Alps (Quincey and Glasser, Reference Quincey and Glasser2009). These conditions are present at Grosse and Exploradores, closely resembling descriptions by Kirkbride (Reference Kirkbride1993) in the New Zealand Alps. In contrast, the initial stages of lake development at Gualas and Reicher glaciers were primarily associated with outlet expansion, tributary river confluences, and lateral widening between the ice and the moraines (Supplementary material, Fig. S8 and S9). Exploradores Glacier exhibits features of both patterns. Inspection of buoyancy conditions during initial lake development is limited; however, flotation may have begun at Gualas and Grosse between 1976 and 2000. We interpret that floating conditions did not play a major role in lake initiation at these glaciers; however, their onset (if reached) may represent a threshold for phase 2.
5.1.2. Glacier recession along retrograde bed topography and into deepest lake section
This stage is characterized by rapid frontal retreat as the terminal tongue, or part of it, attains buoyant conditions and calving of large tabular icebergs, retreating through retrograde slopes and the deepest lake sections.
Gualas and Reicher glaciers showed a sharp transition into stage 2, whereas Grosse Glacier exhibited only moderate features of this stage. We interpret that Exploradores Glacier has now entered the early phase of stage 2. During this stage, Reicher and Gualas Glaciers experienced the highest retreat rates and surface velocities of the record, coinciding with their advance into retrograde and deeper sections of their lakes. Retreat was principally driven by the detachment of large tabular icebergs that drift in the lake after detachments. This can be attributed to form when a portion of the tongue is in floating conditions. These conditions persisted until the glacier fronts reached prograde slopes steeper than 10°, marking the onset of stage 3. Interpreting the evolution of Grosse Glacier is more challenging. Satellite imagery suggests that, instead of producing large tabular icebergs, the debris-covered tongue disintegrated through the expansion and coalescence of supraglacial ponds, which surrounded isolated ice − cliff peaks. After that, buoyant conditions for Grosse for 2000‒2010 indicate that a portion of the glacier near the centerline is on floating conditions whilst it is pinned on the lateral borders (Fig. 6c). Grosse Glacier has developed a coherent calving front but has not yet advanced to prograde slopes. The glacier continues to show an increasing velocity trend, and further deepening of the lake may enhance frontal ablation in the coming years, accelerating mass loss (Supplementary Material, Fig. S5).
Glacier modelling studies have pointed out that glacier retreats on retrograde slopes are faster. For example, (Sutherland and others, Reference Sutherland, Carrivick, Gandy, Shulmeister, Quincey and Cornford2020) modelled glacier retreat under lake’s influence and identified a period of ‘rapid collapse’ from the lip of the over deepening basin into the commencement of reverse slope and deeper water (>200 m). In the case of Grosse its geometry and lateral pinning appear to be critical to stabilize the calving front on a low angle (<3º) retrograde lake bottom slope. Overall, this evidence underscores that lake geometry exerts control on retreat dynamics. Once buoyant conditions are reached, rapid disintegration often follows (Tsutaki and others, Reference Tsutaki, Sugiyama, Nishimura and Funk2013, Reference Tsutaki2019; Minowa and others, Reference Minowa, Schaefer and Skvarca2023). Similar patterns have been documented at Upsala and O’Higgins Glaciers, where increases in flow speed preceded rapid retreat, attributed to flotation of the tongue and buoyancy − driven calving (Minowa and others, Reference Minowa, Schaefer and Skvarca2023).
Other factors that could alter the dynamics of the glacier are landslides that descend onto the glacier (Van Wyk de Vries and others, Reference Van Wyk de Vries, Wickert, MacGregor, Rada and Willis2022). Such events can lead to temporary speedups, changes in flow patterns, or even long-term modifications in glacier behavior. In 2002, a rockslide of ∼300 m extent fell onto Reicher Glacier, covering a 30–50 m wide sector of ice relative to the ∼1 km wide tongue. This event coincided with a period of rapid frontal disintegration (Supplementary Fig. S9, 2002). However, due to limited data density, the specific dynamics and causal relationships surrounding this event remain unclear.
5.1.3. Formation of a coherent calving front, prograde glacier bed and lake disconnection
This stage is characterized by stabilization of glacier front after rapid retreat, associated with the transition onto prograde bed slopes and culminating with lake disconnection.
At Glaciers Gualas and Reicher, stabilization occurred once the fronts reached prograde glacier bed slopes, accompanied by a decline in surface velocities. Reichter, however, displays a more complex pattern. Its front has alternately advanced and retreated, with surface velocities rising before 2015 and then declining, followed by renewed acceleration since 2023. Note that surface velocity estimates for Reicher and Gualas were derived from measurements located upstream of a knickpoint and may therefore not fully represent dynamic conditions at the glacier tongue; interpretations of velocity changes during this stage should therefore be treated with caution.
5.2. Reclassification of Exploradores glacier as lake-terminating
Previous studies and glacier inventories have classified Exploradores Glacier as land-terminating (RGI, Reference RGI2017; Minowa and others, Reference Minowa, Schaefer, Sugiyama, Sakakibara and Skvarca2021; DGA, 2022). However, it has been suggested that the glacier is undergoing a transition towards a lake-terminating state (Aniya and others, Reference Aniya, Barcaza and Iwasaki2007; Irarrazaval and others, Reference Irarrazaval, Dussaillant, Vivero, Iribarren-Anacona and Mariethoz2022). In line with these observations, our results provide new evidence showing that marginal lake development and glacier dynamics now support reclassifying Exploradores Glacier as lake-terminating.
First, marginal lakes now surround most of the terminal tongue. On the western margin, a narrow lake (50–200 m wide) extends for over 4.5 km, covering an area of ∼0.7 km2. On the eastern margin, a larger lake extends for more than ∼4.3 km and covers ∼2.4 km2, with a prominent calving front. At its downstream end, the glacier remains in contact with ∼1.5 km of the terminal moraine (Fig. 8a).
Second, bathymetric surveys of the eastern marginal lake confirm that the terminal tongue of Exploradores Glacier lies within an overdeepened basin, reaching a maximum measured depth of 330 m. This finding is consistent with ice-thickness model predictions (Fürst and others, Reference Fürst2024). Although the bathymetric coverage is limited and the bed elevation models involve considerable uncertainties, the presence of an overdeepening can be confirmed, even if its full shape and extent cannot yet be resolved (Fig. 3).
Third, surface dynamics provide additional evidence. Crevasse openings and calving events on the east margin involve large tabular icebergs detachment (Fig. 8a). A discontinuity in the velocity field further shows that the eastern section of the glacier is drifting northeast, while the main trunk flows northward (Fig. 8b). This pattern is interpreted as evidence of a large tabular iceberg detaching, resembling the dynamics observed in other glaciers during proglacial lake development (Dell and others, Reference Dell, Carr, Phillips and Russell2019). Since 2020, and up to the most recent data available in September 2024, surface velocities at the glacier tongue have increased (Fig. 7). The recent acceleration is consistent with observed lake development and reduction of effective pressure and basal drag (Tsutaki and others, Reference Tsutaki, Sugiyama, Nishimura and Funk2013, Reference Tsutaki2019; Minowa and others, Reference Minowa, Schaefer and Skvarca2023). Collectively, these findings indicate that the behavior and morphology of Exploradores Glacier are characteristic of a lake-terminating glacier.
5.3. Ice mass loss underestimation from unaccounted subaqueous ice loss
Geodetic mass balance methods use surface elevation models to derive ice volume loss between two dates. However, these methods do not account for ice volumes below lake water levels. This discrepancy has led to an underestimation of mass loss of up to 10% in the Central Himalaya between 2000 and 2020 (Zhang and others, Reference Zhang2023), for example. In this study, subaqueous ice loss is estimated assuming a vertical downward extension of the glacier front to the lakebed. This simplification can have biases as studies have shown that the glacier can develop an underwater foot or have floating tongues (Truffer and Motyka, Reference Truffer and Motyka2016). In those cases, subaqueous volume loss might be overestimated or underestimated, respectively.
Here, ice loss volume estimates from previous geodetic mass balance studies are compared with the unaccounted subaqueous ice computed in this study. Data from McDonnell and others (Reference McDonnell, Rupper and Forster2022) for the periods 1976–2000, Dussaillant and others (Reference Dussaillant, Berthier and Brun2018) for 2000–18 and Hugonnet and others (Reference Hugonnet2021) for 2000–20 are used to integrate or accumulate the total ice volume loss over a given period. Subaqueous ice loss is obtained by the increase in lake volume corrected by ice density over the specific period. This value is subsequently divided by the total ice mass loss (sum of geodetic and subaqueous components) over the same period and expressed as a percentage to quantify the relative contribution of subaqueous ice loss. In Table 2 the ratio of subaqueous ice loss to the combined subaqueous and geodetic volume loss is presented.
Table 2. Contribution of subaqueous ice loss relative to total ice mass loss derived from combined geodetic and subaqueous estimates.
a McDonnell and others (Reference McDonnell, Rupper and Forster2022).
b Dussaillant and others (Reference Dussaillant2019).
c Hugonnet and others (Reference Hugonnet2021).
Volume loss from all glaciers is underestimated as all glaciers formed lakes during the study period. Gualas glacier experienced most of the lake expansion during the period 2000–18 and therefore its maximum unaccounted subaqueous ice volume loss is up to 10%. Between 1976 and 2000, 40% of the ice loss from Reicher Glacier occurred subaqueously. This is attributed to the significant expansion of the lake during this period, with water depths reaching more than 300 meters. Between 2000–18 and 2000–20, Reicher Glacier continued to experience a higher volume of unaccounted subaqueous ice loss. Note that the difference between 14% and 25% is due to the difference in estimated annual mass balance of −0.59 m a−1 (McDonnell and others, Reference McDonnell, Rupper and Forster2022) and −0.94 m a−1 (Dussaillant and others, Reference Dussaillant2019) calculated from different data and processing techniques. Grosse Glacier rapidly disintegrated between the 1980s and 2000s, leading to an unaccounted subaqueous ice volume loss of up to 19%. Retreat rates have remained steady, resulting in up to 12% of unaccounted ice volume loss between 2000–18 and 2000–20. The front of Grosse Glacier lies on a retrograde slope, indicating that underestimated subaqueous ice loss is likely to increase for the coming years. Lastly, subaqueous ice loss at Exploradores Glacier accounted for 9% to 10% of total ice loss between 2000–18 and 2000–20. However, this underestimation could increase as the lake continues to expand.
5.4. Limitations and outlook
The interpretation and analysis presented here are constrained by data uncertainties, particularly in bed elevation models, bathymetric surveys, velocity records which are sparse in earlier periods and improve only in recent years. Bathymetric coverage is especially limited in recently deglaciated areas not included in the surveys, as approaching the glacier front is hazardous. Bed elevation models also remain coarse and require continuous improvement to better characterize glacier dynamics and forecast future evolution (Carrivick and others, Reference Carrivick, Tweed, Sutherland and Mallalieu2020; Sutherland and others, Reference Sutherland, Carrivick, Gandy, Shulmeister, Quincey and Cornford2020).
Another limitation is the restricted assessment of upstream sectors and glacier-wide dynamics. Previous studies have found that dynamic changes near the calving front can propagate upstream. High frontal ice loss can trigger dynamic thinning that extends into upper glacier sections, amplifying glacier-wide mass loss. Future work should therefore investigate the upstream impacts of near-terminus acceleration and their role in shaping long-term glacier mass balance (Liu and others, Reference Liu2020). Subglacial topography indicates additional overdeepenings upstream at Gualas and Reicher Glaciers, including a retrograde slope beginning just 1.2 km from (Gualas, Reference Gardner, Fahnestock and Scambos2022) front (Fig. 4). Combined with its low surface slope (2–4°) over the first 3 km, this suggests potential for new lake formation under continued thinning. Future work should refine bed data targeting these zones to anticipate further glacier–lake evolution.
This study is also limited by the lack of analysis of the lake’s thermal. This regime is critical for heat transfer at the ice–water interface, where it can accelerate melt. However, cold subglacial inputs to the deeper sections of a lake may suppress melting in those areas (Sugiyama and others, Reference Sugiyama2016, Reference Sugiyama2021). In some systems, underwater sills act as barriers that constrain the distribution of subglacial meltwater, isolating parts of the lake with cold water, and thereby with the potential to reduce underwater melting and calving (Sugiyama and others, Reference Sugiyama2021). At the three sites with available bathymetric data (Gualas, Reicher, and Grosse), no sills were identified. However, the T-shaped geometry of Reicher Lake may justify a closer inspection of its thermal regime. In this study, we assume that the thermal regimes of these lakes were not strongly compartmentalized by basin morphology. Instead, the thermal regime evolved from early stages dominated by ice contact with limited atmospheric influence on later stages increasingly controlled by atmospheric forcing. The extent and influence of this process remain uncertain. Further studies should collect thermal profile data from the onset of lake formation through subsequent development, as the complex interplay between atmospheric conditions, lake topography and thermal processes remain poorly understood but are critical for advancing knowledge of glacier–lake evolution (Truffer and Motyka, Reference Truffer and Motyka2016).
6. Conclusions
This study analyzes land-to-lake transitions at four western Patagonian glaciers using imagery, bathymetry, velocity, and geodetic data. The record spans the sequence from initial lake formation to near disconnection, providing comprehensive evidence of glacier–lake evolution. Although broadly consistent with previous conceptual models the land-to-lake transition can be described in three phases: (i) thinning with slow frontal retreat, (ii) fast disintegration associated with buoyant-driven calving, and (iii) the formation of a stable calving front on prograde slopes. The four glacier–lake systems examined illustrate different points along this trajectory. While Gualas and Reicher have advanced through all three stages within the past five decades. Grosse, remain in an intermediate configuration, with retreat moderation by bed geometry and lateral pinning. These contrasts underscore local topography control on glacier–lake evolution.
Exploradores Glacier is undergoing dynamics changes linked to the commencement of fast retreat phase. Combined with the ongoing expansion of marginal lakes, this evidence supports reclassifying Exploradores Glacier as lake-terminating glacier in future glacier inventories. The case of Exploradores also raises concerns for tourism, as nearly 9000 visitors undertake ice-walking tours each year. Lake-terminating glaciers can disintegrate abruptly, underscoring the need for updated safety guidelines.
Rapid transitions from land-to lake-terminating have direct implications for GLOF assessment. Current inventories lag behind the rapid pace of lake development, risking the omission of emerging hazardous sites. Anticipating new lakes and updating records in near real time is therefore essential. More broadly, these findings highlight the need for more accurate bed-elevation models which are critical for predicting future lake formation and, in turn, for improving GLOF assessment at newly developing sites.
Subaqueous losses accounted for up to 40% and 19% of total ice loss at Reicher and Gualas Glaciers, respectively, during fast retreat (1976–2000). These estimates remain uncertain due to simplifications in representing the unobserved underwater ice geometry. This highlights the need for new bathymetric data to refine mass−loss calculations and assess additional processes such as thermal influence, particularly in hyperhumid glacier systems. Discrepancies between measured bathymetry and modeled beds highlight the need for refinement. Future work should focus on quantifying subaqueous ice loss during rapid−retreat phases and enhancing bedrock models through targeted data acquisition, particularly for smaller glaciers (<100 km2).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10144.
Acknowledgements
NatGeo Grant EC−95830R−2: Changing Landscape: From Glacier to Lakes, FONDECYT Postdoc ANID 3230337, Ecosystem, Climate Change and Socio-environmental linkages along the continent−ocean continuum: Long-Term Socio-Ecological Research in Patagonia PATSER, ANID R20F0002, FONDECYT ANID 1230433, FONDEF IDeA I+D 2021 ANID ID21I10094 and travel grant from NTNU. We are thankful to Camilo Rada G. from SAGAZ Project (www.SAGAZ.org), Francisco Croxatto and Yerko Ortega who operated the packrat-mounted echosounder in Exploradores Lake, and Benjamín Sotomayor for piloting the UAV during the April 2024 survey. Lastly, we thank the Editor and the three anonymous reviewers for their valuable contributions, which considerably improved the manuscript.
Author contribution statement
The study was conceptualized by I.I., M.S. and I.D. Data compilation, glaciological mapping and GIS analysis was done by I.I., B.M. and E.L. Bathymetry data were collected by I.I., M.S., P.E., I.D. and B.R. UAV data were collected by I.I. and B.S. The article was written by I.I. with input and edits from all coauthors.
Open access
