The SPI is characterized by a marked seasonal temperature variation and spatially variable precipitation patterns. South of 49°S, the precipitation is equally distributed throughout the year, with slight maxima in March and April (Sagredo and Lowell, 2012), yet to the north, there is a marked annual cycle. The longitudinal distribution of precipitation is strongly influenced by the presence of the Andes, which, although presenting relatively low elevations in Patagonia, still generates an extreme precipitation gradient with humid western slopes and arid eastern slopes. Annual and interannual changes in precipitation have been shown to strongly impact the surface mass balance of Patagonian glaciers. In terms of the long-term trend, recent studies suggest possible reductions in the amount of snowfall in this region due to climatic warming. Rasmussen et al. (2007), for example, estimate that there has been a 5% reduction in solid precipitation between 1960 and 1999. This agrees with the findings of García-Lee et al. (2024), who determined an annual upward trend of the freezing level throughout Patagonia between 1951 and 2021.

Wilson et al. (2018) presented the first large-scale inventory of glacial lakes in Chile and Argentina covering the Central Andes, Northern Patagonia and Southern Patagonia, reporting an overall increase in glacial lake area of 27% between 1986 and 2016. This work built upon the findings of Loriaux and Casassa (2013), who reported a 64.9% expansion of glacial lakes in the Northern Patagonia Icefield (NPI) between 1945 and 2011. In this study, we present a multi-temporal inventory of glacial lakes in the vicinity of the SPI (Figure 1), characterizing their physical attributes, water volume, and spatial-temporal distribution using Landsat and Sentinel-2 satellite imagery acquired between 1986 and 2023. This research represents an update of the Wilson et al. (2018) glacial lake inventory for the SPI, reporting changes in glacial lakes on an annual basis between 2015 and 2023 and comparing this to the multi-decadal evolution of these features between 1986, 2000 and 2023. To estimate glacial lake volumes, we applied a newly developed empirical area-volume scaling relationship, representing an advancement over the methods used in previous studies such as O’Connor et al. (2001); Loriaux and Casassa (2013); Cook and Quincey (2015); Shugar et al. (2020). Additionally, this study identifies and discusses past GLOF events originating from the SPI, providing valuable insights for future regional hazard management.

In this study, glacial lakes were mapped only if they were located within or immediately adjacent to the 1945 extent of the SPI outlet glaciers. This search area was chosen as it likely contained lakes that had experienced the most significant changes during the observation period (Loriaux and Casassa, 2013). The boundaries of the outlet glaciers in 1945 were derived from a 1:250,000 scale map created and published in 1954 by the Chilean Geographic Military Institute (Instituto Geográfico Militar, IGM). This map was based on the Trimetrogon aerial photographic survey conducted by the U.S. Army Air Force between December 1944 and March 1945. To enhance the accuracy of the analysis, original Trimetrogon aerial photographs were also utilized, as there were some inaccuracies in the maps.

To complete the multi-temporal glacial lake inventory for the SPI, a total of 93 Landsat and Sentinel-2 satellite images were used, with acquisition dates between 1986 and 2023 (Table 1). The satellite imagery used was selected based on image availability and the presence of snow cover, cloud cover and mountain shadowing. All images were obtained from the United States Geological Survey’s (USGS) Earth Explorer interface (https://earthexplorer.usgs.gov/). Named lakes were identified using 1:50,000 maps available from the IGM and reports from the Argentine Geographic Institute.

Each glacial lake mapped was characterized by several quantitative and qualitative attributes, including an identification ID code, name (if available), lake type, central coordinates, elevation, and surface area. Following the guidelines of Wilson et al. (2018) and Carrivick and Tweed (2013), glacial lakes were categorized into three sub-types based on the nature of their impounding dams (Figure 2): (1) moraine-dammed lakes (encompassing all subtypes of moraine dams); (2) bedrock-dammed lakes (situated within bedrock over-deepenings); and (3) ice-dammed lakes (impounded by ice). This categorization was done through visual analysis supported by geomorphological observations made by Davies et al. (2020). These observations included detailed maps of the moraines surrounding the SPI, which were used to help identify the composition of the dams.

To characterize glacial lake changes at a sub-annual scale between 2015 and 2023, when a significant change was observed, we included additional images to increase the temporal resolution of the inventory. However, due to the temporal resolution of the satellite data used (e.g., 16 days for Landsat eight and 5–10 days for Sentinel-2 data), cloud cover and other data quality issues, often the temporal resolution achieved was not high enough to properly characterize the observed changes. Lakes Argentino (1,368.5 km 2 ), Viedma (1,211.5 km 2 ), and O’Higgins/San Martín (1,042.6 km 2 ), which are located on the eastern margin of the SPI (Figure 1), were excluded from the inventory, as they are several orders of magnitude larger than any other glacial lake considered in this study. Supraglacial lakes were also excluded due to their seasonal nature. A small number of lakes ( < 30 ) were not mapped in some images due to localized cloud cover.

Glacial lakes were mapped in this study using a semi-automated approach, combining the use of the normalized difference water index (NDWI) (McFeeters, 1996) and manual editing. The NDWI is calculated following Equation 1: N D W I = P green − P NIR P green + P NIR (1)

where P green and P NIR denote the green and near-infrared (NIR) bands, respectively.

The Normalized Difference Water Index (NDWI) has been widely used in glacial lake inventories (Huggel et al., 2002; Zhang et al., 2021). The effectiveness of the NDWI for glacial lake mapping can be hindered by image quality issues, including snow, ice, cloud cover, and mountain shadowing. Shadowed areas are common in mid-latitude mountain regions such as Patagonia and exhibit a spectral signature similar to that of glacial lakes with low turbidity, leading to potential misclassification (Gardelle et al., 2011; Loriaux and Casassa, 2013). In this study, shadow areas were discriminated from each other through the use of a slope map derived from the 2000 Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) at a 90 m spatial resolution. This slope map was visually inspected together with the satellite imagery to correct any potential misinterpretations (Figure 3).

Estimating the mapping errors associated with the glacial lake area calculations is challenging without using high-resolution reference data. Several factors influence this estimation, such as the spatial resolution of the imagery (e.g., Landsat = 15–30 m and Sentinel-2 = 10 m), geometric accuracy of the images (e.g., Landsat = 15 m and Sentinel-2 = 10 m), the expertise of the operator performing the classification, and specific image quality issues previously reported (Paul et al., 2013; Wilson et al., 2018).

Following these considerations, we adopted the error estimation approach proposed by Hanshaw and Bookhagen (2014), as adapted by Lesi et al. (2022). This methodology estimates error as a function of the number of edge pixels in a lake polygon, with adjustments to account for duplicated nodes. Key parameters include the total number of nodes in the lake polygon ( N total ) , image spatial resolution ( G ) , lake perimeter ( P ) and a coefficient of 0.6872, which reflects that approximately 69% of pixels along lake edges are subject to classification uncertainty.

The equation used to estimate the error is: Error 1 σ = P G − N int × 0.6872 + G 2 2 (2)

where,

• Error 1 σ is the 1 σ error in square meters,

In the original formulation, the number of edge pixels corresponds to the fraction P / G . However, corner pixels at interior nodes were counted twice, leading to an overestimation of uncertainty. To correct for this, Lesi et al. (2022) introduced the term N int in Equation (2), producing a correct count of edge pixels.

Following Lesi et al. (2022), we use two different expressions to calculate N int depending on whether or not the polygon contains internal islands.

For polygons without islands, we estimate the number of inner nodes following Equation 3 as: N int = N total − 4 − 1 2 (3) where N total is the total number of nodes in the polygon, and the subtractions account for four corner nodes and one duplicated endpoint.

For polygons with islands, we apply Equation 4 as follows: N int = N total − N island + 1 × 5 2 (4) where N island is the number of internal islands. Each island is assumed to contribute five nodes (four corners plus one duplicated node).

The total number of nodes for each lake delineation was computed using the exterior.coords function from the GeoPandas package (Jordahl et al., 2020) in Python 3.9.

Finally, the relative uncertainty D is computed following Equation 5 as: D = Error 1 σ A × 100 % (5)

where A is the lake area in square meters.

Individual glacial lakes were assigned unique identification (ID) numbers. To consistently associate the same ID with corresponding lake outlines identified from different years, an equivalent circular area was first calculated for each lake. This area represents the space a circular lake with the same perimeter as the mapped lake would occupy, allowing for the standardization of irregular lake shapes into a comparable metric. Lakes were then matched by comparing their centroids with those of lakes identified in previous time periods. If the centroid of a lake fell within the equivalent circular area of a previously mapped lake, it was then considered to be the same lake. All lakes that did not meet this criterion were assigned ID letters along with their emergence year (year ID), unless they already had a name. If a lake did not fit within the reference circular area, we evaluated whether the areas of both lakes intersected. If an intersection existed, the compared lake was kept as the reference unless it already had a name. In cases where a name was present, the existing name was preserved, and the lake was renamed only if no reference lake was matched.

To assess area changes, we compared the calculated area of each lake across different time periods. If the most recent area measurement fell within the error range of the earlier measurement, we considered the lake’s area unchanged. If the most recent area exceeded the previous measurement and fell outside the error margin, we concluded that the lake had increased in size. Otherwise, we determined that the lake’s area had decreased.

To estimate glacial lake volume, we used the empirical area–volume relationship proposed by Shugar et al. (2020) as a starting point. This approach employs a mixed model with different equations for small and large lakes, using a threshold area of 0.5 km² to divide the two groups. As this method is widely applied at the global scale, it provides a useful benchmark for evaluating its applicability to glacial lakes in the Southern Patagonian Icefield (SPI). The equations used by Shugar et al. (2020) are described in Equations 6–8. For small lakes A ≤ 0.5 km 2 : ln V = β 0 + β 1 ⋅ ln A + ε (6) Back - transformed for volume prediction: V = e β 0 ⋅ A β 1 ⋅ e ε (7) For large lakes A > 0.5 km 2 : V = k 1 ⋅ A k 2 + ε (8)

where V is lake volume in km³;, A is lake area in km²;, and ε is the residual error.

Given the high proportion of small lakes ( < 3   km 2 ) in our inventory, we found it necessary to adjust the model to the characteristics of Patagonian glacial lakes. Therefore, we constructed a new mixed model (V_pow–poly3), which combines a power-law function for small lakes and a third-order polynomial for larger lakes. This model is an empirical fit to the 120 lakes of known volume used by Shugar et al. (2020) and aims to improve volume estimation accuracy, particularly for smaller lakes which the model used by Shugar et al. (2020) tends to overestimate.

The best model was found by minimizing the misfit using the coefficient of determination R 2 , the sum of squares due to error (SSE) and the root mean squared error (RMSE), these quantities were defined in Equation 9 as: R 2 = 1 − S S E S S T (9)

where SST is the total sum of squares around the mean. The SSE is given by Equation 10 as: S S E = ∑ i = 1 N V i − V ̂ i 2 (10) and the SST by Equation 11 as: S S T = ∑ i = 1 N V i − V ̄ 2 (11)

where V i is the observed volume of the i -th lake, V i ̂ is the estimated volume, and V ̄ is the mean volume of the N data points evaluated.

The root mean squared error (RMSE) was calculated by Equation 12 as: R M S E = S S E N (12)

The fitting process was carried out using the robust least absolute residuals (LAR) algorithm in Matlab Center (2021a), and applied to the entire dataset for each polynomial.

Lakes have a natural water level variation resulting from changes in water input. However, on ice-dammed lakes, these variations are much larger due to changes in the ice dam and GLOFs that can fully drain a lake due to the formation of a subglacial or englacial channel. When a significant decrease in area was observed in an ice-dammed lake, we assumed it was due to a GLOF, and we chose to timestamp the event using the date of the image in which the area reduction was detected. However, it might have happened anytime between the date of that image and the previous one in the inventory. For each GLOF, we used our area-volume empirical relationship to calculate the water volume evacuated. O’Connor et al. (2001) developed an empirical model correlating this evacuated volume with the peak discharge of the GLOF. We used this model to estimate the peak discharge of all drainage events identified. Both evacuated volume and peak discharge must be considered lower bounds, as it is likely that our images captured the lake before it was fully drained or after it had partially refilled. The peak discharge was estimated using Equation 13 as follows: Q max = 0.054 ⋅ V max 0.66 (13)

where V max is the total volume drained in cubic meters m 3 . The processes that trigger the sudden drainage of ice-dammed lakes are complex and have been related to changes in ice thickness, lake bathymetry, calving activity, subglacial water pressure, and thermal conditions, among other factors (Roberts, 2005; Clague and Mathews, 1973; Clague and O’Connor, 2021). With this in mind, the GLOF peak discharge estimations presented here are likely to include large uncertainties and, therefore, should only be considered first-order approximations.

Our database shows that 85% of the inventoried lakes have areas between 0 and 3 km 2 . Various methods exist to estimate lake volume, including the mixed model proposed by Shugar et al. (2020). Shugar et al. (2020)’s model is innovative because it applies different area-volume scaling relationships for lakes of different sizes. However, 84% of the lakes they studied had areas below 0.5 km 2 , a point where the mixed model applied exhibited a discontinuity, resulting in similar volume estimates for lakes with significantly different areas. This issue was especially serious for lakes with areas of around 2.82 km 2 .

Given the area distribution in our inventory and the issues caused by this discontinuity, a new relationship was developed in this study to improve model accuracy, particularly for lakes with areas smaller than 3 km 2 . To achieve this, two polynomial models are proposed based on the same dataset used by Shugar et al. (2020): V poly 2 = f 2 A = p 21 A 2 + p 22 (14) V poly 3 = f 3 A = p 31 A 3 + p 32 (15)

where V poly 2 and V poly 3 represent lake volume estimates in km³; for a given lake area A in km 2 , and p 21 , p 22 , p 31 , p 32 are the regression parameters. Both models display similar behaviour to that shown for other glacial lake power area-volume relationship found in the literature (e.g., Iturrizaga, 2014; Cook and Quincey, 2015; Wilson et al., 2018; Buckel et al., 2018), f 2 ( 0 ) = 0 and f 3 ( 0 ) = 0 , with strictly increasing slopes for A > 0 , where f 2 ′ ( A ) > 0 and f 3 ′ ( A ) > 0 when p 21 , p 31 > 0 and p 22 , p 32 ≥ 0 .

The second-order polynomial model, f 2 , has also been used by O’Connor et al. (2001) for area-volume relationships in moraine-dammed lakes. With the third-order polynomial defined in Equation 14 and the conditions on p 31 and p 32 , we ensured f ″ 3 ( A ) > 0 to fix the inflection point at A = 0 , guaranteeing the convexity of f 3 for all A ≥ 0 .

The resulting models are shown in Figure 4, and the fitting parameters in Table 2. Both models show a reasonable fit to the lake volume data, although V poly 3 provides better R 2 and RMSE values than V poly 2 , along with a lower S S E for lakes with areas smaller than 3 km 2 . However, from the residual analysis, a systematic overestimation was observed for lakes larger than 0.57 km 2 in V poly 3 , while V poly 2 showed an overestimation for lakes larger than 0.35 km 2 , which represent 94% of the lakes in our database.

Even though both models introduce some bias in estimating smaller lakes, the statistics for the third-order polynomial suggest that V poly 3 may be more useful for prediction. However, the overestimation of the small lakes by V poly 3 is addressed by using a mixed model, as proposed by Shugar et al. (2020). This model avoids discontinuities in the glacial lake volume estimation and better fits lakes with areas smaller than 3 km 2 .

The mixed model developed consists of two components that require an intercept area value of at least 0.57 km 2 . In this mixed model, V poly 3 represents the second part, while the first part is based on the classic power area-volume relationship approach, which we established as follows: V pow = f 1 A = k 1 A k 2 (16)

where V poly 3 is the lake volume estimation in km 3 , A is the lake area in km 2 , and k 1 and k 2 are the regression parameters.

Based on Equations 15, 16, the mixed model that was evaluated was set out as follows: V pow−poly3 = k 1 A k 2 A ≤ A 13   km 2 p 31 A 3 + p 32 A A > A 13   km 2 (17)

where ( A 13 ) represents the mixed model transition point, which corresponds to the intercept between ( V pow ) and ( V poly 3 ) . This can be achieved for lake areas greater than zero by selecting a specific subset of the data.

In total, 313 glacial lakes were detected in the 2023 inventory (Figure 6; Table 4), covering an area of 639.09 km 2 and storing an estimated 34.84 km 3 of water. Of the lakes identified, moraine-dammed lakes accounted for the largest portion (165; 52.72% of total), followed by bedrock-dammed lakes (76; 24.28% of total) and ice-dammed lakes (72; 23% of total). In terms of area, moraine-dammed lakes accounted for 56.43% of the total, followed by ice-dammed lakes (34.87%) and bedrock-dammed lakes (8.70%). Despite their smaller area, ice-dammed lakes were found to have a similar total volume (15.65 km 3 ) to moraine-dammed lakes (16.77 km 3 ). In terms of their morphology, most lakes were found to be elongated and aligned with their associated glacier tongue due to the influence of glacial topography and the direction of ice flow. This was particularly the case for moraine-dammed lakes. However, ice-dammed lakes, which typically occur at higher elevations, were found to exhibit more irregular shapes, both spatially and temporally. This is due to the inherent instability and variable geometry of the ice dams that impound them, which, combined with the large fluctuations in water level observed for this type of lake, cause significant changes in lake shape between observations. Ice-dammed lakes primarily develop in areas that have recently undergone deglaciation, causing them to change their orientation and shape along with fluctuations in the surrounding glaciers. The most significant ice-dammed lakes, such as Bernardo and Pío XI, are located in the northwest part of the SPI and parallel to the main mountain range along the SPI.

In terms of their spatial distribution, smaller glacial lakes ( < 3 km 2 ) are shown to be uniformly distributed around the periphery of the SPI, occupying previously deglaciated glacial valleys (Figure 6). In contrast, larger lakes ( > 3 km 2 ) are more concentrated along the southeast and northwest periphery of the SPI. With respect to area, 78.27% of the glacial lakes detected (n = 245) are smaller than 1.3 km 2 . Lake Greve stands out as the largest lake in the 2023 inventory, representing 28% of total area and covering 180.52 km 2 (13.84 km 3 ). Lake Greve is an ice-dammed lake located in the northwest section of the SPI and is dammed by Pío XI Glacier. The elevation analysis shows that 91% of the glacial lakes mapped are located below 800 m. a.s.l (Figure 7).

The assessment of the temporal evolution of glacial lakes, based on a comparison of the 1986, 2000, and 2015–2023 inventories (Table 5), showed that the number of glacial lakes increased from 233 in 1986 to 313 in 2023, reflecting a 34% rise. The total surface area expanded from 496.56 km²; in 1986 to 639.09 km²; in 2023, representing a 29% increase over the same period. In comparison, the total lake volume increased by 31%, from 26.68 km 3 in 1986 to 34.84 km 3 in 2023, with slight interannual variations. Throughout the observation period, 2019 was identified as having the largest total glacial lake area of 646.07 km 2 (Table 5). This total area shrank by 1%–639.09 km 2 by 2023.

Overall, the mapped glacial lakes exhibited various changes, including expansion, coalescence, shrinkage, disappearance, and detachment from their parent glaciers. As expected, bedrock-dammed lakes demonstrated greater stability, as they are typically smaller and situated in geomorphologically stable basins. In contrast, ice-dammed lakes showed high areal variability, while the most significant growth was observed in moraine-dammed lakes (Figure 8).

Throughout the 37-year observation period (Table 5), a total of 352 lakes were mapped (233 lakes in 1986 and a further 199 lakes mapped between 2000 and 2023). Between 2000 and 2023, the newly emerged lakes accounted for a total area and volume of 15.82 km 2 and 0.44 km 3 , respectively. Interestingly, the rate of glacial lake emergence was observed to triple between 2015 and 2023 (an average of 6.7 lakes emerging per year) compared to the 1986–2015 period (an average of 2.2 lakes emerging per year). Within this most recent period (2015–2023), the contrast between 2021 and 2022 is notable, with 17 new lakes emerging, which is far greater than the average of 5.8 lakes emerging per year over the entire observation period.

The analysis of the evolution of lakes that existed in 1986 shows that, in 2023 44% maintained their area (n = 155), 29% increased in area (n = 102), 15% reduced their area (n = 52) and 12% disappeared (n = 44). All the disappeared lakes were dammed by ice. Overall, moraine-dammed lake tended to increase their area while bedrock-dammed lakes tended to maintain their area (Figure 9).

Between 2015 and 2023, several large GLOF events were observed in the SPI (Figure 13; Table 6). These events were primarily associated with ice-dammed lakes, which exhibit dynamic changes due to the formation of drainage channels through the ice that impounds them (Figure 14). One of the most significant GLOF events during the observation period occurred at the lake dammed by Bernardo Glacier (48.59°S, 73.80°W) on 30 May 2018, when a total of ∼ 0.59 km 3 of water drained, reaching a peak discharge of ∼ 33.3 m 3 /s. A similar rapid drainage event was recorded on 15 October 2019, when ∼ 0.66 km 3 drained, accompanied by a peak discharge of ∼ 35.9 m 3 /s (Figures 13, 14). This lake experienced periodic emptying each year and underwent a 53% reduction in area between 2015 and 2023. These changes are likely associated with the recent thinning of Bernardo Glacier, as geodetic measurements from 2000 to 2015/16 indicate a mean elevation change of −2.08 m y − 1 . Another ice-dammed lake that has recently generated a large GLOF is Lake Greve. Dammed by Pío XI Glacier, Lake Greve was relatively stable until 27 June 2020, when ∼ 4.0 km 3 of water drained from the lake with a peak discharge of ∼ 116 m 3 /s, making it the largest GLOF event to have occurred during the observation period. This lake also experienced new drainage events on 3 March 2021 and a smaller one on 4 June 2022.

Other important lakes that have produced GLOFs during the observation period are the lakes dammed by Chico, Viedma and Occidental glaciers. At the lake dammed by Chico Glacier (48.98°S, 73.13°W), cycles of lake filling and emptying are observed, where the lake during emptying periods becomes two smaller lakes. The lake dammed by Viedma Glacier (49.54°S, 73.05°W) has also undergone a number of filling and emptying cycles. This lake recorded its smallest area on 7 July 2020 (0.37 km 2 ) after a short period of expansion and filling, but has generally undergone a period of shrinkage, having reduced its area by 91% between 2015 and 2023. Lastly, the lake dammed by the Occidental Glacier (48.82° S, 74.04° W) was also observed to have experienced a number of drainage events (6 in total). The most notable of these events occurred on 8 January 2016, when ∼ 0.22 km 3 of water was drained at a peak discharge of ∼ 17.3 m 3 /s, leaving only a small body of water behind, surrounded by a large number of grounded icebergs demarcating the former lake extent. Subsequent to this event, the lake continued with its filling and emptying cycles, with similar near-complete drainage events occurring in 2018 and 2019 (Figure 13).

It is worth mentioning that, as Figure 14; Table 6 show, all observed drainage events associated with ice-dammed lakes are indistinctly referred to as GLOF events hereafter. However, we have not identified the exact timing of each emptying episode.

Our methodological approach opened the possibility to produce a large and detailed dataset for studying glacial lakes and their evolution through time at decadal and annual scales. However, some of the processes we have studied, such as GLOFs, have annual periodicity and can exhibit significant changes at a daily scale. Consequently, our estimates of GLOF frequency may be influenced by the temporal availability of imagery. For example, years with more images have a higher likelihood of capturing a GLOF event. Similarly, the timing of image acquisition within a year may introduce bias, as ice‐dammed lakes are more or less likely to be in a filled or drained state depending on the month. In the case of our lake volume calculations, for simplicity and due to the scarcity of validation data, we have treated all glacial lakes as equal, even though different processes have formed, which could result in differing area-volume relationships (Cook and Quincey, 2015). Nevertheless, as Shugar et al. (2020), this study neglects such differences and, therefore, the volume estimations reported must be considered first-order approximations.

Between 1986 and 2023, the number, area, and volume of glacial lakes surrounding the SPI increased by 34%, 29%, and 31%, respectively. By 2023, this study identifies the existence of 313 lakes with a total volume of 34.84 km 3 . The growth of glacial lakes in this area was primarily driven by the expansion of moraine-dammed lakes, resulting from the retreat of many of the SPI’s outlet glaciers. This retreat has led to the formation of new proglacial lakes and the enlargement of existing ones. These findings align with observations from other glaciated regions (Carrivick and Tweed, 2013). For the SPI, the results observed have three key implications. Firstly, it demonstrates that the meltwater storage capacity of the periglacial environment continues to increase. This storage not only has implications for global sea level assessments (Shugar et al., 2020), but may also affect the hydrological regime, sediment transport, water quality, and aquatic ecology of the many proglacial rivers originating from the SPI (Miserendino et al., 2023; Milner et al., 2017; Tweed and Carrivick, 2015; Tiberti et al., 2019). Secondly, the findings underscore the significant role that glacial lakes play in the glacier dynamics of the SPI. When in contact with their source glacier, glacial lakes can negatively impact glacier mass balance through subaqueous melting and the promotion of mechanical calving (King et al., 2019). The latter process is particularly relevant for the calving outlet glaciers of the SPI, some of which, such as Jorge Montt Glacier, have experienced periods of stability and rapid retreat in response to individual calving cycles that are decoupled from the local climate signal (Bown et al., 2019). These influences on glacier change in the SPI emphasize the need for continued monitoring of glacial lakes in this region into the future. Thirdly, southern Patagonia has been identified as a hotspot for GLOFs in the Andes, many of which have been generated by ice-dammed lakes and moraine-dammed lakes surrounding the SPI (Wilson et al., 2018). The continued emergence and expansion of lakes in this region has implications in terms of GLOF hazard and risk management, as discussed in the following sections, with up-to-date glacial lake inventories representing an important baseline dataset for subsequent GLOF hazard assessments (e.g., Iribarren Anacona et al., 2014).

Using geodetic methods, several studies have reported large ice mass losses for the SPI since the 1980s (Rignot et al., 2003; Jacob et al., 2012; Willis et al., 2012; Dussaillant et al., 2019) with Malz et al. (2018) reporting an overall specific mass balance of −0.941 ± 0.19 m w. e. a − 1 between 2000 and 2015/16. It is therefore unsurprising that the number and size of glacial lakes have increased over this period (Figure 11). However, the results presented reveal that the rate of glacial lake area expansion for the SPI has reduced considerably since 2015, reducing from an increase of 0.9% per year between 2000 and 2015 to 0.1% per year between 2015 and 2023. In contrast, the rate of glacial lake emergence was shown to more than triple between 2015 and 2023. However, this increase was driven by the appearance of new ice-dammed lakes, which tended to be smaller and undergo repeated phases of emptying and refilling, with no new moraine-dammed or rock-dammed lakes being observed since 2016. Given this, the reduction in the overall rate of glacial lake area expansion between 2015 and 2023 may suggest that (1) the availability of low gradient ice areas that facilitate lake growth during glacier retreat is reducing in the SPI; (2) glaciers are beginning to detach from their proglacial lakes (in doing so limiting their growth) and retreat to higher gradient portions of their basins and (3) most glaciers have retreated from their main (maybe change for LIA) frontal moraines. Therefore, no new moraine-dammed lakes can be formed. A similar but more pronounced reduction in the rate of glacial lake area expansion was also observed by Wilson et al. (2018) for Northern Patagonia Andes (36°–45°S) but for the earlier period of 2000–2016, with this region having been identified as having undergone more extensive rates of ice loss compared to the Southern Patagonia Andes (45°–55°S) (Wilson et al., 2018; Paul and Mölg, 2014).

From 1986 to 2023, 119 newly emerged lakes were recorded. Taking 1986 as a reference, we looked into the factors promoting the formation of new glacial lakes and the history of their development in relation to their dam type. Among these, bedrock-dammed lakes are usually smaller and offer relatively less variability with respect to area and volume after the peak achieved in 2017 (Figure 12).

On the other hand, ice-dammed lakes are much more variable. Although year after year, new ice-dammed lakes develop, their number is smaller, and they do not last long as opposed to the moraine and bedrock-dammed lakes (Figure 12). This characteristic reflects the observed variability in ice-dammed lakes driven by the drainage and refilling cycles associated with their ice-dam dynamics.

The non‐linear nature of the area–volume relationship used means that a relatively small number of large lakes can contain a disproportionately high proportion of the total lake volume. On the other hand, small lakes are more numerous and, therefore, can have a large water volume. However, we found that small lakes are not numerous enough to compensate for their much smaller individual volumes, and a small number of large lakes monopolizes the water storage reported in our inventory. This was also true for ice-dammed lakes in the periods of 1986–2015 and 2019–2023, when large changes in lake size occurred in the two biggest ice-dammed lakes, dammed by the Pio XI and Bernardo glaciers. This is especially relevant as ice-dammed lakes pose the highest GLOF hazard. This highlights the importance of establishing continuous monitoring of these lakes, especially because of their sustained increase in number and size. New lakes formed between 1986 and 2023 account for around 15.82 km 2 out of an overall increase in lake area of 142.55 km 2 . Moreover, this trend will likely persist in the coming decades, buoyed by the continued ice loss.

As discussed in section 4.2, lakes that have a larger surface area tend to have a disproportionately large impact on water storage. This relationship can also be found when using our mixed approach and has been pointed out in previous studies (e.g., Shugar et al., 2020; Loriaux and Casassa, 2013; Huggel et al., 2002; O’Connor et al., 2001; Cook and Quincey, 2015). In Figures 4, 5, we can see how the model calibrated by Shugar et al. (2020) and by us show how both the model calibrated by Shugar et al. (2020) and our own calibration produce increasingly steep curves for larger lakes, highlighting that lake volume grows disproportionately faster than area as lake size increases. This explains why two ice-dammed lakes, dammed by Pío XI and Bernardo glaciers, make the most significant contribution to the volume increases shown for 1986 to 2015 and from 2019 to 2023.

The Southern Patagonian Icefield (SPI) experienced significant variability regarding the formation and evolution of new glacial lakes, particularly in the northern section (Figure 11). The retreat of glaciers and climatic variations contributed to the dynamics of these areas, alongside the periodic creation and rapid disappearance of ice-dammed lakes due to drainage events likely associated with GLOF phenomena. This disappearance of ice-dammed lakes can be attributed to drainage that occurred by the thinning of the glacial dam because the increasing water depth can cause ice margin flotation, flexure or fracture, and jökulhlaups (Carrivick and Tweed, 2013). Globally, these lakes represent the most common source of glacier outburst floods (Carrivick and Tweed, 2016). However, some of these disappearances may be due to the fact that the image used for the mapping turned out to be from a moment when the lake was empty, but then could have been filled up again. Because of this uncertainty, we cannot be sure that all accounted drainage events could be associated with GLOF phenomena.

As our results have shown, ice-dammed lake drainage events are particularly common in the northern part of the SPI, where lakes have formed in marginal positions around several of the large outlet glaciers, where the greatest changes in glaciers were observed. In particular, in the glacial systems of Pío XI, Bernardo-Témpanos, O’Higgins-Chico, Viedma and Amalia glaciers (Minowa et al., 2021; Mouginot and Rignot, 2015). Periodic or episodic GLOFs were registered for five lakes, dammed by the Bernardo, Chico, Pío XI, Occidental and Viedma glaciers. These events were marked by changes in lake size and/or the presence of newly exposed lake basins, which, in some cases, were scattered with grounded icebergs (Wilson et al., 2018), like the Occidental and Bernardo glacier events (see Figure 14). In that regard, the analysis of ice-dammed lake change between 2000 and 2023 in this study (Figure 12) highlights the dynamic nature of this type of lake, with observations revealing multiple periods of filling and emptying for individual lakes over relatively short time periods. Several drainage events were observed from five ice-dammed lakes surrounding the SPI, some of which reached peak discharges of up to ∼ 116,647 m 3 / s (Table 6), demonstrating the frequency and magnitude of such events in this region.

In particular, the years 2017 and 2020 recorded the highest frequency of GLOF occurrences, as shown in Figure 13. This peak was associated with a rapid glacier retreat. The frequency of GLOFs was particularly high in areas with widely distributed negative ice elevation changes (Dussaillant et al., 2019). The largest of these drainage events occurred at Lake Greve on 27 June 2020, where an 11% reduction in area in 2 months was observed. We estimate that 3.96 km 3 of water was drained during this evident which coincides with that reported by Hata et al. (2022). In terms of discharge volume, the magintude of this GLOF is amongst the largest reported for any glacial lake world (Hata et al., 2022). The relationship between Lake Greve and Pio XI Glacier has previously been investigated (Carrión et al., 2010a; 2010b; Wilson et al., 2018) due to the hazard posed by a potential drainage event into Eyre Fjord if Pio XI Glacier experiences a prolonged period of retreat in the future. Interestingly, Figure 13 shows that post-GLOF, the area of Lake Greve has undergone variations larger than those observed before the event. This suggests that its drainage events might not only be associated with the permanent erosion of a rocky dam, as reported by Hata et al. (2022). Instead, it might be experiencing small drainage events through sublacial, englacial or marginal channels. If that is the case, we consider that a marginal channel is the most likely cause, as a subglacial channel would experience a runaway enlargement and a near-complete drainage.

The second most important event is the drainage of the lake dammed by Bernardo Glacier, in which the area reduction wasover 12 km 2 ( ∼ 60%). However, superimposed in its large areal variations, Figure 13 shows a long-term area reduction trend, most likely associated with the significant retreat and thinning observed at Bernardo glacier (Minowa et al., 2021). Such large variations in an ice-dammed lake can influence glacier dynamics and mass loss, due to: (i) promoting calving, (ii) raising englacial water temperature, (iii) reorganizing the subglacial drainage system due to changes in water supply and pressure, (iv) varying marginal melting rates, (v) promoting ice margin flotation, (vi) raising the englacial water table, (vii) generating ice flexure and fracture during the draining and filling processes, (viii) flushing of sediments from the glacier bed, and (ix) increasing aggradation of sediment at a glacier terminus during GLOFs (Carrivick and Tweed, 2013).

Overall, our results highlight the need to monitor ice-dammed lakes in SPI in particular. Globally, ice-dam failure is responsible for the majority of GLOF events (Carrivick and Tweed, 2016). Although the socio-economic vulnerability to drainage events and GLOFs sourced from the SPI is relatively low, due to the low population density of the surrounding region, ice-dammed lakes, in particular, have the potential to threaten tourism in the region. This was demonstrated in October 2023 when the popular ice trekking routes on the Exploradores Glacier (outlet of the Northern Patagonia Icefield) were temporarily closed to tourists. Located in the San Rafael Lagoon National Park in the Aysén region of Chile, the decision to close Exploradores Glacier by the National Forestry Corporation of Chile (CONAF) was partly due to the rapid expansion of an ice-dammed lake located along the eastern flank of the main glacier trunk. As GLOF frequency is likely to increase because of natural and anthropogenic climate change (Emmer et al., 2022), it is crucial to implement monitoring systems on glacial lakes in order to adapt and take a preventive approach. Currently, one such system is being developed by the SAGAZ Project in the Aysén and Magallanes regions (Rada et al., 2024). Further work should also consider GAPHAZ guidelines (Glacier and Permafrost Hazards in Mountains Group), which have already been implemented in the Peruvian Andes (Allen et al., 2022).