Response of lacustrine glacier dynamics to atmospheric forcing in the Cordillera Darwin

Abstract

Calving glaciers respond quickly to atmospheric variability through ice dynamic adjustment. Particularly, single weather extremes may cause changes in ice-flow velocity and terminus position. Occasionally, this can lead to substantial event-driven mass loss at the ice front. We examine changes in terminus position, ice-flow velocity, and calving flux at the grounded lacustrine Schiaparelli Glacier in the Cordillera Darwin using geo-referenced time-lapse camera images and remote sensing data (Sentinel-1) from 2015 to 2022. Lake-level records, lake discharge measurements, and a coupled energy and mass balance model provide insight into the subglacial water discharge. We use downscaled reanalysis data (ERA5) to identify climate extremes and track land-falling atmospheric rivers to investigate the ice-dynamic response on possible atmospheric drivers.

Meltwater controls seasonal variations in ice-flow velocity, with an efficient subglacial drainage system developing during the warm season and propagating up-glacier. Calving accounts for 4.2% of the ice loss. Throughout the year, warm spells, wet spells, and landfalling atmospheric rivers promote calving. The progressive thinning of the ice destabilizes the terminus position, highlighting the positive feedback between glacier thinning, near-terminus ice-flow acceleration, and calving flux.

1. Introduction

In recent decades, the temperature increase in South Patagonia has been lower than the global average (NOAA National Centers for Environmental Information, 2023). Remarkably, the ice fields of South Patagonia are one of the fastest shrinking ice bodies (Braun, 2019). Since most glaciers in Patagonia calve either into a fjord system or freshwater lake (Warren and Aniya, 1999), a considerable amount of mass loss is attributed to ice-dynamical adjustment (Sakakibara and Sugiyama, 2014; Mouginot and Rignot, 2015; Braun, 2019; Sauter, 2020; Minowa and others, 2021).

In the southern part of the Andes, three icefields host the majority of South Patagonia's glacier volume. These include the Northern and Southern Patagonia Icefields on the mainland of South America and the Cordillera Darwin Icefield, located in the western part of Isla Grande de Tierra del Fuego, which hosts the southernmost temperate glaciers (Fig. 1). It consists of a main ice body with a total area of 1760 km² centered on Monte Darwin (2207 m) to the east, glaciers that coalesce around Monte Sarmiento (2261 m), and some smaller glaciers separated by fjord systems to the west covering an area of 179 km² (Meier and others, 2018; Rada and Martinez, 2022). As in the whole region of South Patagonia, strong baroclinicity in the westerly wind belt characterizes the prevailing very harsh climatic conditions (Schneider and others, 2003; Garreaud, 2009; Garreaud and others, 2013; Sauter, 2020). Westerly winds transport moist, temperate air masses from the South Pacific (Langhamer and others, 2018) and are responsible for the extremely humid climate, low seasonal temperature ranges, strong winds, and high annual precipitation totals (Miller, 1976; Endlicher, 2000; Garreaud and others, 2013). The Andes act as an effective climatic barrier perpendicular to the main flow, creating a heterogeneous precipitation pattern. Annual precipitation amounts at the Patagonian icefields are estimated to be $6.09 \pm 0.64 \, {\rm m\,a}^{-1}$ (Sauter, 2020). Compared to other regions in South America, the Patagonian Icefields exhibit the highest negative specific mass balance rates (Braun, 2019). These high rates obtained from geodetic mass-balance observations are inconsistent with climate records (Malz and others, 2018; Braun, 2019). Previous estimates of frontal ablation suggest that at least 30% of the mass loss in the Southern and Northern Patagonian Icefields occurs via dynamic processes (Minowa and others, 2021), highlighting the importance of ice-dynamic processes.

Figure 1. Overview of South Patagonia and the icefields. The star highlights the location of Monte Sarmiento in the Cordillera Darwin. The inset map displays the glacier outlines at Monte Sarmiento Massif (RGI Consortium, 2017).

When studying glacier dynamics and their response to atmospheric forcing, it is crucial to consider glacier runoff, defined as the flux of water leaving the glacier system (Cogley, 2010), including meltwater, precipitation, and water input from the margins, which can drain down to the ice-bedrock interface and potentially increase basal lubrication (Weertman, 1972). However, whether runoff increases basal lubrication also depends on the state of the subglacial discharge system. The subglacial discharge system drains runoff water efficiently by channelization (Röthlisberger, 1972) or inefficiently by cavitation (Lliboutry, 1968). In the latter, runoff reaches the ice-bedrock interface, increasing basal water pressure, which reduces basal drag and improves lubrication, accelerating ice motion (Weertman, 1972; Iken, 1981; Iken and others, 1983; Iken and Bindschadler, 1986). Thus, the interplay between runoff and discharge efficiency determines a glacier's characteristic seasonal velocity pattern (Moon and others, 2014). Ultimately, high runoff volumes due to single weather events, such as rapid temperature increases or heavy precipitation events, can affect a glacier's ice dynamics (e.g., Kamb, 1987; Meier and others, 1994; Luckman and others, 2006; den Ouden and others, 2010; Schellenberger and others, 2015; Little and others, 2019; Tuckett and others, 2019) when runoff inputs exceed the ability of the existing englacial channel system to discharge water (Bartholomaus and others, 2008). Generally, thinning at the terminus of a calving glacier reduces the effective pressure at the glacier bed (difference between glacier weight and basal water pressure) and causes instantaneous acceleration and rapid retreat as the glacier approaches flotation (Howat and others, 2008; Stearns and van der Veen, 2018). Because of the often asynchronous behavior of calving glaciers on a regional scale, these glaciers belong to the group of non-representative glaciers with respect to climate variation (Post and others, 2011). The dynamic behavior, particularly the ice-flow velocity, terminus position, and calving, can undergo abrupt changes, and their adjustments are partially decoupled from the climate (Benn and others, 2007). For instance, in South Patagonia, different seasonal velocity patterns can exist side-by-side, and retreat or advance can occur in immediate proximity (Sakakibara and Sugiyama, 2014).

To date, no long-term records are available to study the dynamic response to the atmospheric forcing of Cordillera Darwin's calving glaciers. Glaciological fieldwork in the Cordillera Darwin Icefield is challenging and uncommon due to the remoteness and harsh conditions. Therefore, only remote sensing-based measurements of glacier area, volume, length, and velocity changes and numerical simulations of mass balance are available (Melkonian and others, 2013; Braun, 2019; Weidemann and others, 2020; Temme and others, 2023). However, the latter neglect ice dynamics. The spatial and temporal resolution of repeat-pass synthetic aperture radar (SAR) remote-sensing capabilities limits the ability to detect small isolated calving events and associated frontal changes of lacustrine glaciers.

Since 2013, we have recorded ice-dynamic adjustments at a representative lacustrine glacier, Schiaparelli, located at the Monte Sarmiento Massif (Fig. 2). The surface ice-flow velocity in the glacier's ablation area is estimated using a combined approach of SAR images and images from in-situ monoscopic camera systems. Furthermore, we used the camera systems to record changes in the ice-front position at two-hour intervals from 2015 to 2022, with occasional interruptions. A coupled snowpack and ice surface energy and mass balance model estimates glacier runoff, which we compare to the measured lake discharge to evaluate the efficiency of the subglacial discharge system. We use automatic weather station records to downscale ERA5 temperature and precipitation to the complex mountainous terrain to reproduce the atmospheric forcing over the whole glacier. This study provides a detailed picture of the spatial and temporal variations in ice-front position and surface velocity. We study the relationship between glacier discharge and ice-flow velocity using (time-lagged) correlation analysis to understand the interplay between runoff and ice-flow velocity. In addition, we evaluate recurrent weather extremes and assess their contribution to calving and surface ablation with the aim of quantifying potential atmospheric drivers conducive to controlling the ice dynamics and total mass balance of Schiaparelli Glacier.

Figure 2. Overview of the study area with Sentinel-1 surface ice-flow velocity estimates averaged from 2015 to 2022. White lines indicate terrain height (Pléiades © CNES 2020-02-03, Distribution Airbus D&S; Marti and others (2016); Beyer and others (2018); Deschamps-Berger (2020)) and red lines ice thickness in meters (Farinotti, 2019). Blue areas with contour lines present the lake depth in meters. Black circles (P1 to P7) represent the region of interest for the calculation of surface ice-flow velocity presented in Figures 11, 9, 14, and 16 along the centerline (black dashed line obtained from Maussion, 2019). The black dots illustrate the section of the photo above that captures the glacier terminus, the proglacial Lago Azul, and its discharge into the fjord system. The photo was taken at a ridge $\approx 710\, \rm{m\,msl}$ in March 2020.

2. Study site

The scenery of the Cordillera Darwin Icefield, with its numerous fjords, mountain peaks, and calving glaciers, has fascinated scientists and adventurers since its description by the famous Beagle Expedition in the 19th century. ‘The occurrence of glaciers reaching to the water's edge & in summer, in Lat: 56 $^\circ$ is a most curious phenomenon: the same thing does not occur in Norway under Lat. 70 $^\circ$. From the number of small ice-bergs the channel represented in miniature the Arctic Ocean’ (van Wyhe, 2006, p. 139). These were the words Charles Darwin used to describe the unique climate and environmental conditions at the Cordillera Darwin in his Beagle Diary on January 29, 1833. Lithographs taken during the second Beagle Expedition in 1836 provide the earliest documentation of the glaciers around Monte Sarmiento (Darwin and others, 1839), increasing our awareness of glaciers’ adaptation to climate change in the mid-latitudes of the Southern Hemisphere.

The proximity to the fjord system results in a temperate maritime climate with an annual mean temperature of 1.2 °C and a low seasonal amplitude of 5.9 °C. Located south of the core of the westerlies and upstream of the mountain's main ridge, orographically enhanced uplift leads to high annual precipitation totals of 3914 mm and almost permanent cloud cover. During winter, precipitation sums are slightly lower ($\approx 200 \, {\rm mm\,month}^{-1}$) than in summer ($\approx 330 \, {\rm mm}\,{\rm month}^{-1}$) due to a weakening of the westerly wind belt (Garreaud and others, 2013; Langhamer and others, 2018). According to a non-parametric Mann–Kendall test, there is a significant (p ≤ 0.05) linear trend in both ERA5 temperature and precipitation of 0.09 °C decade⁻¹ and 4.0 mm decade⁻¹, respectively, illustrating sustained changes in climate forcing from 1950 to 2022.

The Monte Sarmiento Massif comprises a glaciated area of around $70 \, {\rm km}^2$ (Barcaza, 2017) with the four largest glaciers named Conway, Lovisato, Emma and Schiaparelli (Fig. 1). Schiaparelli Glacier covers the largest area (24.3 km² in 2016 (Meier and others, 2018)), extending from an elevation of 2013 m above mean sea level (m msl) on the northeast side of Monte Sarmiento (Rada and Martinez, 2022) to the moraine-dammed lake, Lago Azul, at 17 m msl. Schiaparelli Glacier, named by the Salesian missionary, explorer, mountaineer and scientist Alberto De Agostini in 1913 to honor the Italian astronomer Giovanni Schiaparelli (1835-1920) (De Agostini, 1955), is a grounded lacustrine glacier, like most glaciers in South Patagonia (Warren and Aniya, 1999). The glacier calves into Lago Azul, a proglacial lake that formed after the glacier receded from its former terminal moraines in the 1940s (Meier and others, 2019). Through a cut in the moraine ridge, the lake discharges westward down a stream (Rio Azul) into the Magdalena Channel, providing a suitable setup for glaciological studies (Fig. 2).

Since 2015, two automatic weather stations (AWSs), one in the glacier's vicinity and the other close to the centerline, have been recording air temperature, air pressure, precipitation amount, relative humidity, wind direction and velocity, snow height and global radiation (Fig. 2). Fieldwork included measurements of ice thickness (Gacitúa and others, 2021), lake bathymetry and the discharge of Lago Azul. The lake-level and water temperature fluctuations are monitored using a water pressure sensor to determine a functional relationship between lake-level and lake discharge. On-site unoccupied aerial vehicle (UAV) missions provide spatially detailed annual height and length changes of the glacier margin. Ablation stakes have been used to validate surface mass balance models in the lower ablation area (Weidemann and others, 2020; Temme and others, 2023). To capture the most recent changes and their causes, camera systems have recorded the glacier's terminus since 2015.

3. Materials and methods

3.1 Atmospheric data

The ERA5 reanalysis dataset is the highest temporally and spatially resolved global reanalysis dataset to date (Hersbach and others, 2020). Variables at surface and pressure levels are available in hourly time steps at a horizontal resolution of approximately 31 km. Surface variables of temperature, relative humidity, pressure, wind velocity and cloud cover are used to generate the climatic forcing for the coupled snowpack and ice surface energy and mass balance model in python (COSIPY). Therefore, ERA5 data is extracted at the four closest grid cells to the AWS and statistically downscaled to the local conditions at the AWS by quantile mapping (Gudmundsson and others, 2012). Quantile mapping is a technique for downscaling climate-model data via statistical bias correction by adjusting the cumulative distribution function of the model data to the observed data. It has been successfully applied in recent studies in Southern Patagonia (e.g., Weidemann and others, 2018, 2020) and shows a good performance in the Monte Sarmiento Massif (for evaluation results, see supplementary table S1 in Temme and others (2023)). Since wind velocities at the AWS are strongly influenced by the local topography, we decided to take this variable directly from the ERA5 data. Statistically downscaled air temperature and pressure are spatially interpolated from the AWS over the topography using a linear temperature lapse rate and the barometric equation, respectively. By fine-tuning the parameters within both the atmospheric and melt model, we determined a lapse rate of $0.60 \, {\rm K}\,(100\, {\rm m})^{-1}$ (Temme and others, 2023). This estimate is consistent with values previously used in the region (e.g., Strelin and Iturraspe, 2007; Koppes and others, 2009) and the average obtained from the ERA5 dataset. The generated high-resolution dataset of atmospheric conditions at the study site is used to force COSIPY following the approach in Weidemann and others (2020) and for the climatological analysis of extreme events. The precipitation is downscaled by an orographic precipitation model (see the following section). The required input is taken from ERA5 and consists of upwind information on geopotential height, air temperature, wind vectors and relative humidity between 850 hPa and 500 hPa and total precipitation over the study site. We use the horizontal integrated water vapor transport (IVT) to detect atmospheric rivers.

3.2 Orographic precipitation model

Since precipitation events can be highly variable spatially and temporally, downscaling is challenging over complex terrain. Furthermore, observations in southern Patagonia are known to be error-prone due to strong winds (Schneider and others, 2003, 2007). Thus, linear extrapolation of such sparse and uncertain observations with an elevation-dependent lapse rate is critical. Therefore, precipitation is simulated following a physically motivated approach using a linear model of orographic precipitation (Smith and Barstad, 2004; Barstad and Smith, 2005; Sauter, 2020), which has been shown in previous studies to reproduce in-situ measurements and their temporal variability well over mountainous terrain (e.g., Schuler and others, 2008; Jarosch and others, 2012; Weidemann and others, 2013, 2018; Sauter, 2020; Temme and others, 2023). The model is based on the linear steady-state theory of orographic precipitation and calculates the precipitation resulting from forced orographic uplift over a mountain, assuming saturated and stable conditions. The cloud water and hydrometeor density are given from advection, condensation due to terrain-forced uplift, and the conversion of cloud water to hydrometeors, as well as the fallout of those producing precipitation (Smith and Barstad, 2004; Barstad and Smith, 2005; Sauter, 2020; Weidemann and others, 2020). The total precipitation is calculated by adding the orographic precipitation computed by the model to the large-scale precipitation, which is given by removing the orographic component from the ERA5 precipitation (Weidemann and others, 2013, 2018; Sauter, 2020; Weidemann and others, 2020; Temme and others, 2023). The model parameters include a threshold of relative humidity above which precipitation can occur (90%), determined on observations, and the fallout and conversion time of hydrometeors ($1200 \, {\rm s}$), calibrated for the Monte Sarmiento Massif, in combination with surface-mass balance model parameters based on observations of geodetic mass balance (Temme and others, 2023).

3.3 Glacier mass balance and runoff model

Similar to Weidemann and others (2020), we use an updated version of the COSIPY model (Sauter and others, 2020) to estimate the glacier's total runoff. The runoff, including meltwater and liquid precipitation, represents the liquid water flux leaving the glacier system (Cogley, 2010). The model performance of COSIPY was compared to three different surface mass balance models of varying complexity in the Monte Sarmiento Massif (Temme and others, 2023). Mass balance estimates of COSIPY agree well with the other model results and show good performance with regional and local geodetic mass balance estimates and ablation stake measurements (Temme and others, 2023). COSIPY combines a surface energy balance with a multi-layer subsurface snow and ice model, where the computed surface meltwater serves as an input for the subsurface model (Sauter and others, 2020). The total ablation includes surface melting, which is calculated by solving all energy fluxes at the glacier surface, sublimation and subsurface melting. Accumulation is possible through snowfall, refreezing and deposition. Snowfall is derived from downscaled precipitation, distinguishing between solid/liquid phases with a logistic transfer function scaling around a threshold temperature of 1.0 °C. Furthermore, the model was extended with a basic parameterization of snowdrift (Warscher and others, 2013; Temme and others, 2023), modifying the snowfall distribution based on the topography and the prevalent wind direction. The simulations cover the Monte Sarmiento Massif, including Schiaparelli Glacier, with a 200 m spatial and three-hourly temporal resolution for 2000–2022. The model parameters were partly calibrated (albedo of ice and firn and ice roughness length) and partly fixed based on sensitivity tests and reference values. For a detailed description of the model setup, calibration and evaluation, we refer to Temme and others (2023) and Tables 1 and S2 therein.

3.4 Discharge model

Pressure readings from the water-level sensor are converted to relative water levels of Lago Azul by the hydrostatic equation, accounting for the ambient air pressure. The relative water level of 0 m corresponds to the lowest point of the lake outflow at Rio Azul, when no discharge would be possible. The discharge was measured at the outlet of Rio Azul over a cross-section of $\approx 20 \, {\rm m}$ (Fig. 2). Following the velocity area method, water flow velocity measurements were taken at 1 m intervals (Morgenschweis, 2018). A two-point method was used, measuring the flow velocity at 20% and 80% of the water depth in 2018 and 2019 and a fixed depth approach in 2022. A total of 14 discharge measurements were conducted at different lake levels to establish a relationship between discharge and lake level (Fig. 3). The discharge-lake-level relation follows a power function of the form f(x) = ax ⁿ (Morgenschweis, 2018). We achieve the best fit with a coefficient of determination R ² = 0.87 using the coefficient a = 0.51 and the exponent n = 0.32. We use the 95% prediction interval to provide model uncertainties.

Figure 3. Empirical discharge function (blue line) representing the relationship between lake level and discharge.

3.5 Identifying extreme meteorological events

Environmental responses to intense and frequent climate and weather extremes are complex and globally diverse (Karl and others, 1997). To quantify extremes, the Expert Team on Climate Change Detection and Indices (ETCCDI) suggests a set of 27 temperature and precipitation indices to assess terrestrial climate variability and change (Karl and others, 1999; Peterson and others, 2001, 2002). Six consecutive days of daily maximum temperatures in the upper and lower 10th percentile of maximum temperature centered on a five-day window (1961–1990) are defined as a warm spell duration index (WSDI) or cold spell duration index (CSDI), respectively (Karl and others, 1999; Peterson and others, 2001; Zhang and others, 2005). In southern South Patagonia, hardly any events (eleven warm spells and five cold spells) of six consecutive days could be attributed to WSDI or CSDI from 2015 to 2022. Thus, four instead of six consecutive days, as defined above, are used to identify warm or cold spells in accordance with the 90th percentile of all identified consecutive days. We used a base period from 1991 to 2020 as suggested by the World Meteorological Organization (WMO). Additionally, we count the number of consecutive wet days (CWD) (consecutive dry days (CDD)) of at least (less than) 1 mm daily precipitation sum. Wet and dry spells are identified if the number of consecutive days exceeds the 90th percentile of all identified CWD and CDD (i.e., 17 and 7 consecutive days, respectively).

3.6 Atmospheric river tracking

We detect and track landfalling atmospheric rivers (ARs) with an Image-Processing-based Atmospheric River Tracking method (IPART) (Xu and others, 2020b). While conventional detection methods of ARs typically use either IVT thresholds (e.g., ${\rm IVT} \geq 250 \, \hbox{kg m}^{-1} {\rm s}^{-1}$) or a relative IVT magnitude threshold (e.g., 85th percentile of local climatology) (Ralph and others, 2004; Dettinger and others, 2011; Lavers and others, 2012; Rutz and others, 2014), the IPART algorithm works by the ‘top hat by reconstruction’ technique (Xu and others, 2020b). Additionally, this method can capture even the genesis and decaying stage of single ARs, especially in high-latitudes where less water vapor is present in the atmosphere (Xu and others, 2020a,2020b).

According to Xu and others (2020b), we set a geometrical filter to drop ARs with an area smaller than 5×10⁵ km² or greater than 18×10⁶ km², a length subceeding 800 km or exceeding 11 000 km and having a length-width-ratio of smaller than two. We evaluate all ARs with their centroids inside $120 ^\circ$W and $60 ^\circ$W and $10 ^\circ$S and $80 ^\circ$S. Due to the high spatial resolution of the ERA5 dataset, the kernel size (E) has been adapted to E = [16, 18, 18] (number of grid cells). We assume that any precipitation event in the study area that falls within the contours of a detected AR of at least one reanalysis time step per day is associated with an AR (Viale and others, 2018).

3.7 Lake sounding

In April 2018, bathymetric measurements of Lago Azul were conducted using a Garmin Echomap Plus 42cv echosounder mounted on a touring kayak. At three-second intervals, lake depth measurements were taken along parallel transects spaced 40 m to 60 m apart (Fig. 4). Subsequently, the resulting 14 100 lake depth measurements were interpolated to the entire lake area by the Triangulated Irregular Network method (Peucker and others, 1976).

Figure 4. Bathymetric map of Lago Azul. The red line represents the survey grid of April 2018. The black lines show the glacier retreat, which has been derived from UAV mosaics captured in October 2016, March 2017, March 2018, April 2019, March 2020 and January 2022.

3.8 Photogrammetric processing

Two identical monoscopic camera systems (Canon EOS 1200D, as shown in Fig. 5) recorded changes in the glacier-front position and on the glacier surface (Fig. 2). The lower camera system (example image in Fig. 6) was installed close to the glacier terminus and operated from October 2015 until March 2019 (Weidemann, 2021). The upper camera system has been operating since March 2018, and it is located approximately 200 m above the glacier terminus perpendicular to the main glacier flow (example image in Fig. 5). Based on a time interval of two hours, images with similar illumination and lighting were chosen on a three-to-five-day frequency. To estimate the ice-flow velocity from the selected camera images, we use the Python-based open-source tool PyTrx (How and others, 2020). The tool is adaptable, performing image transformations of 2D images into a real-world coordinate system and feature tracking for ice-flow velocity estimates. PyTrx uses a pinhole camera model that assigns each pixel of the image coordinate system (u,v) into a world coordinate system (X _w, Y _w, Z _w),

(1)$$z_c \left(\matrix{u\cr v\cr 1 }\right) = P \left[\matrix{X_w\cr Y_w\cr Z_w\cr 1 }\right],\; $$

Figure 5. The digital single-lens reflex camera (DSLR) system and time-lapse images. (a) Installed time-lapse camera system. (b) Original image from the upper time-lapse camera on November, 23 2020. The dashed black frame refers to the zoomed-in section captured in the images shown in (c) to (f). The dotted frame refers to the zoomed-in section captured in (g) and (h). This sequence of images presents the largest calving events on record. (c) Ice front image from November 19, 2020 and its fracture line (dashed red line) and (d) from November 20, 2020 immediately after the calving event. (e) Ice front image from November 1, 2021 and its fracture line (dashed red line), and (f) from November 3, 2021 after the calving event. (g) Ice front image from April 27, 2021 and (h) from April 29, 2021 after the calving event.

Figure 6. Relative position of the glacier front and rate of length change were derived from the lower (blue dots) and upper (red dots) time-lapse camera. The corresponding line indicates the centered 30 d rolling mean. The thick black line shows the water level relative to the reference height. The rate of length change is depicted in the panel above as centered 30 d rolling mean in m d⁻¹. The images above show the ice front from the lower camera system in November 2016 (left) and January 2022 (right). The background color indicates the ERA5 daily temperature anomaly with respect to the 2015–2022 mean. The 75% of the longest-lasting identified wet, and dry (top), cold, and warm spells (bottom) are shown separately in the upper panel. The roman letters within the Figure correspond to the following assignment: I – both camera systems were operated simultaneously; II – largest observed calving event in November 2020 (cf. Figs. 5c, d); III – calving event in April 2021 (cf. Figs. 5g, h); IV – calving event in November 2021 (cf. Figs. 5e, f).

where z _c is an arbitrary scaling factor and P the camera perspective projection matrix (Xu and Zhang, 1996). P comprises internal camera properties denoted as the intrinsic camera matrix K and the external camera matrix such as orientation and location,

(2)$$P = {\bf K} [ {\bf R}\vec{T}] ,\; $$

with the rotation matrix R and the translation vector $\vec {T}$ (Xu and Zhang, 1996). The pinhole camera model does not consider distortion effects due to the camera lens (Xu and Zhang, 1996). PyTrx does evaluate the radial and tangential distortion coefficients and the intrinsic camera matrix using a technique based on the OpenCV toolbox by using black-white chessboard camera images (How and others, 2020). The extrinsic camera matrix is determined by the exact location and camera pose (How and others, 2020). Therefore, the camera position must be measured by the differential global positioning system (DGPS), and the pose must be estimated using a stereo reference camera technique. For the latter, we used the photogrammetric processing software AgiSoft PhotoScan Professional (Version 1.7.0) (2020). We estimated the camera pose based on camera images from 17 different DGPS-measured locations and common features in the image plane with the corresponding DGPS-measured position markers, so-called ground control points (GCP). Appropriate GCP are stable features in the camera field view of known world coordinates, such as boulders, quartz veins in rocks, or mountain peaks.

Because very small variations (${\cal O} \approx 10^{-4\, \circ }$) of the estimated pose lead to large deviations (${\cal O} \approx 10 {\rm m}$) in the transformation (Equation (1)), PyTrx has included a camera pose optimization routine based on GCP to refine the estimated camera model (How and others, 2020). Combined with a digital elevation model (DEM) retrieved from multiple UAV surveys during subsequent field campaigns in 2016, 2017, 2018, 2019, 2020 and 2022, a world coordinate is assigned to each pixel of the image coordinate system.

3.9 Ice-flow velocity estimates

PyTrx uses the sparse feature-tracking approach that produces a spatially more detailed velocity map than the alternative dense feature-tracking method. Sparse feature tracking identifies patterns in the image plane with unique pixel-intensity distributions between image pairs (How and others, 2020). Static features visible in the fore- and background of the image plane correct the camera motion. When the feature-tracking or motion correction fails, it produces unrealistic movement patterns that are too high or against the main direction of the ice flow. Thus, we filtered for derived velocities within the 95 % confidence interval and for the direction of movement with orientation along the main flow of the glacier. As the final estimated ice-flow velocity derived from one image pair, we use the mean of all estimated ice-flow velocities that fall inside a circle ($r_{\rm P1} = 400 \, {\rm m}$) close to the terminus (Fig. 2).

However, the camera images only produce a detailed temporal estimate of the ice flow near the terminus. Ice-flow velocity estimates from SAR obtained from the Sentinel-1 mission provide a spatial view of the ice flow over the entire glacier. SAR data is independent of daylight, cloud cover, weather and the sun's illumination of Earth's surface and thus appropriate to study the ice flow in heavily cloudy regions such as southern South Patagonia (Jawak and other, 2015). Worldwide post-processed Sentinel-1 velocity maps of 12 glacierized regions outside the large polar ice sheets are available on the online application RETREAT (2021) (Friedl and others, 2021). RETREAT presents an open-access interactive web interface providing spatial velocity maps with a spatial resolution of 200 m and a time interval of 6–48 days since 2014 (Friedl and others, 2021). We calculated the ice-flow variation along the centerline in 50 m intervals. Here, we consider the mean ice flow within a circle of 400 m (Fig. 2).

3.10 Glacier terminus changes and calving flux

The terminus line T is defined as the dividing line between the glacier front and lake water. We determine the terminus line in the image plane and project it into the world coordinate system using PyTrx. Projecting a terminus retreat on a DEM causes errors in assigning the respective world coordinate whenever the terminus line is projected on the glacier front or surface. Thus, we use a modified DEM instead, where the glacier height is equal to the lake surface. The resulting terminus line on the lake surface is separated into 1 m intervals. For each interval i, we calculate the distance δT _i to the newly estimated terminus position parallel to the centerline (Fig. 7). By considering their means ($\overline {\delta T}$), we retrieve a time-dependent terminus length change relative to the initial terminus position.

Figure 7. Schematic representation of the calculation of glacier length changes. $T_{t_0}$ is the initial terminus position and $T_{t_0 + \Delta t}$ the position at any time step Δt. Length changes are calculated parallel to the centerline (dashed line) in 1 m increments.

To evaluate the reliability of the estimates of glacier-length change of both camera systems, we calculated changes along the ice front using a sample image pair from both the lower and the upper camera system. This process was repeated 20 times, with the position of the ice front determined manually each time, resulting in a standard error of 0.06 m d⁻¹ and 0.08 m d⁻¹ for the lower and upper camera estimates within the camera coordinate system, respectively. In addition, we georeferenced the camera model projections onto the DEM to align with the terminus line determined by the UAV missions (Fig. 4) to reduce the offset to the real world coordinate system.

We estimate the calving flux using a combined approach of a fixed flux gate located upstream of the calving front position and the temporal glacier terminus position (Evans and others, 2022). Let $T_{t_0}$ be the ice-front position at the time t ₀ and $T_{t_0 + \Delta t}$ the new observed ice-front position after any time step Δt. Without any calving, the ice-front position would be pushed forward by the glacier motion and reach a new simulated position,

(3)$$T_{{\rm sim}, t_0 + \Delta t} = T_{t_0} + v_{\Delta t} \cdot \Delta t ,\; $$

where we assume the ice-flow velocity (v _Δt) is the observed mean velocity close to the terminus occurring during Δt (Fig. 8). The calving flux (q _c) is calculated by the difference between the simulated and observed ice-front position,

(4)$$q_{c, t_0 + \Delta t} = ( T_{{\rm sim}, t_0 + \Delta t} - T_{t_0 + \Delta t}) \cdot w \cdot h \cdot \rho_{\rm ice},\; $$

with the width w of the ice-front. The ice thickness h is estimated by the sum of the bathymetry measurements of the lake depth close to the terminus (Fig. 2) and yearly DEMs based on the UAV missions. Intra-annual changes are assumed to be linear. By considering the density of ice $\rho _{ice} = 900 \, {\rm kg.m}^{-3}$, we calculate the water equivalent calving flux (Equation (4)). This approach neglects subaqueous melting and assumes that any change along the ice front is spread across the entire ice depth.

Figure 8. Schematic illustration of calving flux calculation between two terminus positions T at two time steps, t = t ₀ (black line) and t = t ₀ + Δt (grey line). The simulated terminus position is shown in green. Similar to the estimation of glacier length changes (Fig. 7), changes parallel to the centerline (dashed line) are calculated in 1 m increments.

3.11 Time-series analysis

Separation of the spatial and temporal signatures of glacier velocity allows us to identify the primary driving force acting on ice-dynamical processes in response to atmospheric forcing. In many cases, several temporal signals overlap, with seasonal fluctuations often superimposed on low-frequency signals. To separate the signal in its short- and long-term contribution, we emulate the original signal by linear combinations of elementary functions. Riel and others (2021) suggest splitting the signal into a short-term seasonal contribution by third-order B-spline functions, which may vary in amplitude from year to year. Long-term changes, including non-steady and non-periodic signals, are represented by time-integrated B-spline functions (Riel and others, 2021). This methodology is especially suited for discontinuous and heterogeneous time series. Its performance was tested on Greenland's surface ice-flow velocity estimates based on Sentinel-1 data (Riel and others, 2021).

To investigate calving's response (a non-linear dataset) to an extreme event (non-linear phenomena), we first apply a peak detection method (Du and others, 2006). The peak detection identifies local maxima of the original data set, in our case the calving flux. If there is a peak in calving flux within five days after an extreme weather event, the event is considered a potential cause. We test the resulting frequency's significance by comparing this result to surrogate datasets. Surrogates are used to generate pseudo-observations that approximate the calving flux by having the same power spectrum and amplitude distribution. We generated 10 000 surrogates by an iterative amplitude adjusted Fourier transform (IAAFT) algorithm (Kantz and Schreiber, 2004; Venema and others, 2006). When the original dataset's peak contribution of events such as landfalling AR or the onset of warm, cold, wet and dry spells is greater than the 95th percentile of the surrogate samples, the results are considered statistically significant.

We use a Spearman rank-order test to assess the statistical dependence between the rankings of two variables. Throughout the study, we only mention significant correlations (i.e., p ≤ 0.05).

4. Results

4.1 Terminus position

Over the seven years, we observed a terminus retreat of about 200 m with an identifiable seasonal cycle and with ice advances during the cold season. Due to a malfunction in the camera systems (an empty timer or camera battery and no memory access due to defective SD cards), no data was collected from October 2016 to March 2017 and April 2019 to April 2020. We used two different camera locations to determine the position of the ice front. Both camera systems operated simultaneously from March 2018 to August 2018 (highlighted with I in Fig. 6). We use this period to verify accuracy and to evaluate any system-related differences. The estimations of the terminus position of both camera systems are consistent with a maximum deviation of 1% in the distance estimates to the reference line of 2015 (Fig. 6). On average, the calculated distance of the ice front position from the reference position deviates by 3% from the world coordinate system comparing the estimates of the camera system and the UAV measurements.

4.2 Ice-flow velocity

The functional representation of Sentinel-1 ice-flow estimates shows a seasonal cycle, with the smallest velocities at the end of the melt season and maxima in the middle of the cold season (Fig. 9). Figure 10 shows the functional representation of the ice-flow velocity estimates along the centerline. A phase shift in the seasonal ice-flow velocity pattern is apparent, indicating a time delay that increases with increasing distance from the terminus. Variations in velocity are inversely related to variations in melt and runoff water. Daily temperature anomalies affect the total glacier runoff but not the ice-flow velocity (Fig. 9). These observations are confirmed using the results of ice-flow velocity estimates based on the time-lapse camera images (Fig. 11). Considering the monthly rolling mean, the ice-flow velocity reaches its peak in late cold winter ($\approx 0.45 \, {\rm m\,d}^{-1}$), while a minimum is reached during the melt season ($\approx 0.33 \, {\rm m\,d}^{-1}$), when melt and runoff water input is relatively high. Figure 11 considers only ice-flow velocity estimates from the upper camera system, which continuously worked from May 2020 to January 2022. Difficulties in obtaining appropriate ice-flow velocity estimates arise when the fraction of the ice surface visible in the camera image is too small and when the oblique angle of the camera to the glacier flow reduces the signal-to-noise ratio.

Figure 9. Ice-flow velocity estimates at circles 1, 4 and 7 (cf. Fig. 2) based on Sentinel-1. Solid lines represent the superimposed functional representation of the ice-flow velocity after Riel and others (2021), shaded areas the uncertainty, and the colored markers the individual estimates (P1: crosses, P2: circles, P3: triangles). The black solid and dashed lines indicate the modeled glacier runoff and meltwater (monthly rolling mean) across the whole glacier. The background color shows the daily mean temperature anomalies averaged over the glacier. Velocity estimates are missing from April to November 2018.

Figure 10. (a) Sentinel-1 ice-flow velocity and (b) modeled glacier runoff along the centerline of the glacier. The gray line presents the air temperature. Velocity estimates are missing from April to November 2018.

Figure 11. Ice-flow velocity estimated by the upper time-lapse camera (orange dots). The black line represents their monthly centered rolling mean, and the orange line represents their superimposed functional representation after Riel and others (2021) and its uncertainty (gray shaded area). The blue solid and dashed lines indicate the modeled daily glacier-wide runoff and meltwater contribution (monthly rolling mean), respectively, as daily sums in meter water equivalent (m w.e. d⁻¹). The background color shows the daily mean temperature anomalies averaged over the entire glacier.

Concerning the estimates obtained from the time-lapse camera, we observe a greater variance in the ice-flow velocity during the cold season. In particular, days with low illumination and visibility limit the ability to track features along the glacier surface and fail when there is snow cover. In general, velocity estimates from the time-lapse camera are 15% higher and show a greater variance than Sentinel-1 estimates. We cross-checked the results by DGPS measurements close to the glacier terminus (the cyan cross in Fig. 2) from January 16–21, 2022. We measured an average ice-flow velocity of 0.36 m d⁻¹. During the same period, time-lapse camera-based estimates of flow velocity were 16% higher (0.42 m d⁻¹). It was not possible to validate the velocity estimates of Sentinel-1 due to data unavailability. Nevertheless, similar velocity patterns and relative changes are found in both approaches. However, absolute values should be handled with care.

An early season peak is visible in the monthly rolling mean ice-flow velocity estimates obtained from the time-lapse camera in December 2020 (Fig. 11). Gaps in the Sentinel-1 observations may explain the missing early season peaks in December 2020 and other years (Fig. 9). For instance, at the beginning of the warmest summer spanning from December 2019 to January 2020, a time when the glacier runoff exceeded that of other years by over 20%, the smaller reduction in velocity observed near the glacier's terminus (P1) might suggest the occurrence of an early season peak in ice velocity.

4.3 Calving flux, runoff and discharge

To estimate the calving flux from 2015 to 2022, we used a combination of both Sentinel-1 mean velocity estimates close to the terminus in P1 (Fig. 9) and time-lapse camera terminus changes (Fig. 6). The average calving flux is estimated to 0.137 m³ s⁻¹ (Table 1)), which corresponds to an annual mass loss of 4.33 Mt a⁻¹. The calving peaks during the melt season, when the lake temperature 2 m below the surface ($\overline {T}_ {\rm melt-season} = 1. 9 \, ^\circ {\rm C}$) is on average 0.9 °C warmer than in the cold season. During the winter months, the calving flux is reduced by more than half (0.083 m³ s⁻¹) as the glacier advances. The presented uncertainties of the calving flux include the standard deviation of the ice-front position estimates and the errors in the surface ice-flow velocity. In addition, we consider an accuracy of $\pm 10 \, {\rm m}$ for the ice thickness estimates to account for neglecting subaqueous melting.

Table 1. Mean m³ s⁻¹ values for glacier runoff, surface ablation, calving flux and mass balance (MB = accumulation - ablation) over the entire study period whenever estimates of the calving flux are available (cf. Fig. 12). The numbers in bold indicate their corresponding seasonally detrended means. The data is also grouped by specific atmospheric events. The number of days n defined as warm season or cold season is given in parentheses (n _warm,n _cold)

Since 2020, the mass loss from calving has increased by 8%. Greater ice-front retreat in summer, as presented in Figure 6, increases the calving flux by 60%. Conversely, since 2020, the calving flux in winter has decreased by 30% compared to the previous winter months. This was especially evident during 2020, the coldest winter of the study period ($\approx 1 \, ^\circ {\rm C}$ colder than the other winter periods), when the lake was partially frozen, combined with the terminus acceleration, resulting in a glacier advance of 25 m (Fig. 6).

There is a significantly high correlation between the modeled glacier runoff and lake discharge (R = 0.89), which reaches its maximum with a time lag of three days (R = 0.92) (Figs. 12 and 13). The temporal evolution indicates similar extremes, and the modeled runoff and lake discharge indicate high positive correlations to air temperature. Note that the lake discharge represents all water reservoirs (e.g., non-glacial snowmelt, groundwater, precipitation on non-glacierized areas of the catchment) contributing to the catchment of Lago Azul. Thus, it is reasonable that lake discharge exceeds modeled runoff. During the warm season the total glacier runoff correlates significantly (R > 0.77, p ≤ 0.05) with the lake discharge, with a leading time lag of two days (Fig. 13). In comparison, there are significant correlations of R > 0.71 up to a four-day lag during the cold season.

Figure 12. Solid lines indicate the five-day centered rolling mean of glacier runoff, lake discharge, calving flux, lake temperature and downscaled ERA5 air temperature. The solid black line in the upper panel indicates the mean height of the terminus derived from UAV missions. The vertical lines indicate the onset of landfalling ARs and warm spells. For clarity, we only present ARs that lead to large temperature increases according to the 75th-percentile (i.e., $\Delta T \geq 2. 7 \, ^\circ {\rm C}$).

Figure 13. Spearman's rank correlation coefficient (R) and leading time-lag. ‘Runoff’ denotes the combined modeled rain and melt runoff from the glacier. Note that the lake level is, by definition, related to the lake discharge due to the discharge-lake-level relation (Fig. 3).

Low discharge estimates should be treated cautiously because there is only one discharge measurement $< 7 \, {\rm m}^3{\rm s}^{-1}$ (Fig. 3). A higher uncertainty for the discharge estimates $> 15 \, {\rm m}^3{\rm s}^{-1}$ results from applying a fixed-depth approach instead of a two-point method in 2020 (Morgenschweis, 2018).

Meltwater controls the lake level and shows a seasonal cycle with the maximum water level during the melt season. The air temperature correlates (R = 0.92) with the lake level with a leading time lag of three days (Figs. 6 and 13). Individual peaks in the lake level can be directly related to daily temperature anomalies. Precipitation amounts are less important in controlling the lake level (there is no significant correlation; thus, it is not shown in Fig. 13). These findings are similar to a previous study, which considered a period from October 2015 to April 2016 (Weidemann and others, 2020).

Calving correlates positively with lake temperature (R = 0.57), air temperature (R = 0.39), and lake level (R = 0.44) (Figs. 12 and 13). Warm spells, wet spells and ARs are associated with increases in the average daily lake-level changes at Lago Azul from 3.6 cm d⁻¹ to 8.8 cm d⁻¹. Since 2020, the impact of daily air and lake temperature and lake level on the calving flux increased to R = 0.52, R = 0.62 and R = 0.49, respectively.

4.4 The link between extreme meteorological events and calving

Considering all the days in the study period, 25% of the days are associated with wet spells, 5% with dry spells, 8% with warm spells and 5% with cold spells. While the 75% of the longest-lasting dry spells occur mainly in the cold season, wet spells occur almost only in the melt season (Fig. 6). The longest wet spell lasted 90 days, from November 2020 until February 2021, with a total precipitation of 1153 mm, which corresponds to 34% of the annual precipitation total. The periods covered by wet spells are identical with AR periods and identified warm spells in 36% and 32% of all cases, respectively. From 2015 to 2022, landfalling ARs contributed, on average, 28% to the total precipitation, with a seasonal minimum contribution during winter months and increased the daily mean temperature by $1.0 \, ^\circ {\rm C}$, on average, after onset.

The inherent noisiness of the data, the superposition of multiple events, and their non-linearity and complex coupling pose a challenge in isolating leading processes that increase the calving activity. Thus, we use surrogate datasets to assess the significance of recurring extremes, such as warm, cold, dry and wet spells and landfalling ARs on calving. We consider a leading process to be significantly controlling the calving flux if the contribution is greater than the 95% confidence interval provided by the surrogates. This ensures that any dependencies are not arbitrary. While there is no evidence of a response of ice motion to specific atmospheric events, we identified in 50% of the cases a peak in calving flux after the onset of a warm spell, in 38% of the cases when a precipitation event is associated with an AR, and in 46% after the onset of a wet spell. The seasonally detrended signal is used in further analyses to ensure that the retrieved averages of runoff or MB-related terms are not solely the result of their intra-annual occurrence. Detrended runoff and surface ablation (surface melt plus sublimation) time series show that warm spells, wet spells and ARs are characterized by a 12% to 43% increase in runoff and an 8% to 45% increase in surface ablation (Table 1). Calving flux increases during AR (14%) but hardly increase during wet spells. Warm spells indicate an increase in calving flux by an average of 29%. In contrast, dry and cold spells are characterized by low runoff, low surface ablation and a low calving flux, while only cold spells contribute to a positive mass balance. Notably, calving flux decreases by 43% during cold spells and by 21% during dry spells.

5. Discussion

5.1 Seasonal ice-flow velocity

The ice-flow velocity shows a seasonal pattern, with higher ice-flow velocities during the cold season (MJJAS) than during the melt season (NDJFM) (Fig. 14). Thereby, the air temperature regulates the availability of runoff (Fig. 13) which could potentially increase basal lubrication. Counter-intuitively, there is an inverse seasonal relationship between runoff and ice-flow velocity, suggesting a transition from an efficient to inefficient subglacial drainage system during the melt season (Andrews, 2014; Moon and others, 2014; Vijay and Braun, 2017). During the melt season, fast lake level response and slower ice motion support that the subglacial drainage system evolves seasonally to efficiently channel summer runoff (Fig. 13). As runoff decreases toward the end of the melt season, velocity increases (Figs. 9, 11, and 14), and the response time of lake level to runoff doubles (Fig. 13), indicating a change to an inefficient englacial drainage system during the cold season. Despite the relatively low availability of runoff, basal meltwater and ongoing meltwater percolation from ice and firn may further increase basal water volume and pressure, thereby enhancing basal lubrication and causing higher ice-flow velocities during the cold season (Brinkerhoff and O'Neel, 2017). In the transition month of December, the ice-flow velocity peaks in the early melt season, considering the monthly rolling mean in Figure 11. At the same time, the greatest annual increase in runoff occurs when the drainage system is still inefficient.

Figure 14. (a) Seasonal ice-flow velocity and (b) modeled glacier runoff anomalies along the centerline of the glacier (from P1 to P7 as indicated in Fig. 2). The gray line presents the seasonal temperature pattern.

The mechanism controlling the development of the subglacial drainage system can be illustrated by the so-called ‘Röthlisberger channel (R-channel) theory’ (Weertman, 1972; Röthlisberger, 1972; Mathews, 1973). An efficient englacial drainage system develops during the melt season when the runoff water widens the channels by frictional heat. Glacier runoff drains through pressurized tubular englacial R-channels, which tend to form main veins as meltwater increases (Röthlisberger, 1972; Fudge and others, 2008; Pohle and others, 2022). This reduces the basal water pressure during the melt season and enhances friction at the glacier bed, which in turn decreases ice-flow velocity during the summer (e.g. Bartholomaus and others, 2008; Bartholomew and others, 2010; Hoffman and others, 2011; Sundal and others, 2011; Andrews, 2014). Once the ambient ice pressure exceeds the water pressure in the channels, they collapse within a few days under viscous deformation (Röthlisberger, 1972; Vieli and others, 2004; Bartholomaus and others, 2008; Pohle and others, 2022). The breakdown of the R-channel system reduces subglacial water drainage capability, causing an increase in the englacial residence time of the runoff.

Over the course of the year, the capability of englacial discharge alternates between efficient and inefficient. This seasonal behavior has been reported primarily in regions where the melt season is long, and runoff rates are high, such as in southeast Greenland (Moon and others, 2014; Solgaard and others, 2022). When the seasonal velocity follows a similar temporal pattern as the runoff, the discharge is inefficient throughout the period (Moon and others, 2014; Solgaard and others, 2022). These glacier types are typically found in colder regions with a shorter melt season and primarily determined by limited meltwater availability. There are also glaciers controlled by meltwater in the Southern Patagonia Icefield, such as Glacier Jorge Montt and the northern front of Glacier Pío (Sakakibara and Sugiyama, 2014). However, the velocity patterns of glaciers in South Patagonia and the controlling mechanism behind them have not yet been thoroughly evaluated.

Seasonal velocity changes propagate up-glacier by approximately 3 km month⁻¹ (Fig. 14a). Such upstream propagating waves have been studied on tidewater glaciers in Greenland at Jakobshavn Isbræ where the discharge is mostly inefficient over the course of the year (Moon and others, 2014), and described as kinematic waves (Riel and others, 2021). Kinematic waves are primarily driven by ice mass redistribution caused by calving or thinning (Riel and others, 2021), a typical feature of glaciers that follow a seasonal velocity pattern in phase with runoff availability (Moon and others, 2014). In the case of Schiaparelli Glacier, we do not observe any event-related velocity increase that we could explain by calving activity. Indeed, glacial runoff follows a similar seasonal pattern as ice-flow velocity, but with a 5-month phase shift and seasonal runoff variations propagating up-glacier at least twice as fast as velocity changes (Fig. 14b). As the melt season progresses, the changing phase of discharge efficiency may shift upstream when channels develop further up-glacier. Thus, we assume an upward propagating wave of basal water pressure changes to explain the upward propagation of the seasonal speed variation. Owing to the underlying bedrock topography, the magnitude of seasonal velocity variations peaks between 1.5 km to 2.2 km from the terminus until an inverse bedrock slope is reached at 2.5 km (Fig. 15). A more pronounced englacial channel system attenuates seasonal speed-up near the terminus, especially where the glacier reaches below lake level.

Figure 15. Cross-section along the centerline with the mean surface ice-flow velocity (2015 to 2022). Dashed white lines show the center point of the black circles (P1 to P7 presented in Fig. 2), and the dotted black line is the lake level.

5.2 Multi-annual ice-flow velocity changes

Short datasets ($< 25 \, {\rm years}$) do not capture long-term trends or variability and make it impossible to assess their response to the underlying climate forcing. The 2015 to 2022 analysis offers a snapshot to evaluate recent changes and discuss their potential atmospheric drivers. Similar to outlet glaciers in Alaska and on the Antarctic Peninsula (Burgess and others, 2013; Tuckett and others, 2019), preceding cold-season temperature or runoff in the years 2015 and 2016 correlate negatively with the subsequent melt-season velocity (Fig. 16). However, multi-annual speed-ups seem to be spatially decoupled. Low temperatures until the onset of the melt season result in reduced runoff water availability, especially at higher elevations where temperatures remain below freezing. Reduced runoff water availability degrades the formation of the englacial drainage system in the upper area and causes the onset of a multi-year acceleration. This can be seen in 2015 at distances $> 1500 \, {\rm m}$ from the terminus, where the bedrock elevation of the glacier is higher than the lake level (cf. Fig. 15).

Figure 16. (a) Multi-annual ice-flow velocity and (b) modeled glacier runoff anomalies along the centerline of the glacier. The black line presents the multi-annual air temperature anomalies. Velocity estimates are missing from April to November 2018.

Since 2019, a speed-up close to the terminus can be observed (Fig. 16a). The acceleration is inherent to thinning (Fig. 12). Once the loss of bed traction outweighs the loss of driving stress, a gradual change in dynamics occurs (Howat and others, 2007; Pfeffer, 2007). Basal drag is reduced, and parts near the front accelerate (Pfeffer, 2007), a common feature of sudden glacier acceleration (Thomas, 2004). This process has positive feedback similar to the marine instability hypothesis, which explains a sudden increase in the dynamic instability of tidewater glaciers (Truffer and Motyka, 2016). Progressive thinning causes a gradual decrease in effective pressure at the glacier bed and an acceleration of the parts of the glacier close to the front (Amundson and Truffer, 2010; Post and others, 2011; Stearns and van der Veen, 2018). The acceleration causes further stretching-induced thinning close to the terminus (Venteris, 1999; Howat and others, 2008). Moreover, the combination of thinning and acceleration often results in an increase in mass loss due to calving (e.g. Meier and Post, 1987; Howat and others, 2007; Post and others, 2011). Between 2019 and 2020, there was an enhanced decrease in the terminus height (Fig. 12), but there was no enhanced mass loss due to melting (Fig. 9). We suspect an additional dynamic thinning near the terminus with acceleration-induced ice stretching as the driving force.

5.3 Changes of ice-front position

Apart from the long-term trend of the terminus retreat, there are seasonal variations in the ice-front position, with a retreat during the melt season when ice motion is slower (Fig. 6). Retreat and advance can be linked to seasonal changes in air temperature. The transition between seasonal advance and retreat coincides with temperature anomalies. Thus, thinning may intensify seasonal variations (Luckman and others, 2006; Howat and others, 2007; Pfeffer, 2007) and become highly unstable when the terminus thins to floating (Howat and others, 2008). In addition, larger seasonal temperature anomalies, warmer (longer) summers and colder (shorter) winters supported the stronger terminus advance during winter and enhanced the retreat in summer from 2020 to 2022.

The underlying bed topography is a major controlling factor in the duration and extent of glacier retreat once frontal instability begins (Meier and Post, 1987; Howat and others, 2007). Schiaparelli Glacier tends to retreat rapidly due to a shallow and slightly reversed slope forming within 1 km of the ice-front resulting in a over-deepening of the basin (Fig. 15). This intensification may progress or remain at a point determined by the bedrock topography as basal drag increases and the height-above-buoyancy criterion is reestablished (Benn and others, 2007). Since 2020, the mean terminus height has dropped below 10 m above lake level (Fig. 12), with a section near the centerline reaching the flotation thickness ($< 5 \, {\rm m}$; cf. Figs. 5 and 6), considering a lake depth in this section of 50 m (Fig. 1). Simultaneously, the stability of the ice-front position has decreased. Daily ice-front variations have more than doubled from 0.26 m d⁻¹ (September 2015 to April 2019) to 0.61 m d⁻¹ (May 2020 to January 2022) (Fig. 17). This period is characterized by a substantial ice-mass loss related to individual single calving events with average terminus retreats of up to 20 m, as observed in November 2020, April 2021 and November 2021 (highlighted with II, III and IV in Fig. 6). If thinning at the terminus section continues at the rate of recent years from 2015 to 2022 of −10.6 m a⁻¹, larger parts of the terminus could become ungrounded, leading to increased instability of the ice front position during the next decade.

Figure 17. Mean location and daily changes of the ice-front position separated in 1 m intervals along the ice front for the two periods, September 2015 to April 2019 and May 2020 to January 2022. The dashed blue and red lines represent the extreme positions of 2015–2019 and 2020–2022, respectively. The background image is the mosaic from the UAV mission in 2017 and matches the mean ice-front position from 2015 to 2019.

As the glacier advances, areas near the centerline get thinner until they become buoyant, where the most pronounced variations along the ice front are observed. Note that the maximum advance in 2020 and 2021 was reached one month before the largest retreats in November 2020 and 2021, highlighting the increased vulnerability when the ice-front position becomes ungrounded (Fig. 5). Changes in lake level increase the effective principal stress in the terminus due to buoyancy torque and can cause ablation at the ice front (Benn and others, 2007). These observations are similar to those reported from lacustrine glaciers in southeast Alaska (Boyce and others, 2007; Trüssel and others, 2013).

5.4 Calving

In total, calving contributes 4.2% to the total mass loss of Schiaparelli Glacier and has a pronounced seasonal calving pattern (Fig. 12). Despite the minor contribution to the mass loss, calving nevertheless makes an important contribution of 14% to the total mass balance of −0.99 m³ s⁻¹ (Table 1). A comparison of the surface mass balance (MB) with the geodetic MB at Schiaparelli Glacier resulted in a calving flux estimate of 4.26 Mt a⁻¹ (Temme and others, 2023), which is in very good agreement with our estimate (4.33 Mt a⁻¹).

When the glacier reached its maximum advance in November 2020 and 2021, an almost-flat glacier surface with a backward tilted block near the ice front can be temporally observed (Figs. 5c, e). These records suggest that parts of the once grounded terminus become buoyant as the ice thins, creating torque and tensile stress in the ice front area. This is a highly unstable and transient situation, making floating termini on temperate glaciers a rare observation (Walter and others, 2010; Trüssel and others, 2013). Surface mass balance, dynamic thinning and lake-level changes are essential components controlling the buoyancy close to the ice front of lacustrine glaciers (Benn and others, 2007; Truffer and Motyka, 2016). Estimates from the MB model indicate ablation-induced thinning as the main controlling factor, causing the ice front to become partially ungrounded. Ongoing thinning and a rising lake level force the terminus out of hydrostatic equilibrium, creating substantial bending forces in the vicinity of the non-buoyant ice junction (Benn and others, 2007). During the largest calving events of 60 m x 160 m in November 2020 and of 40 m x 150 m in November 2021, the fracture developed along a transverse crevasse separating the buoyant and non-buoyant area close to the glacier's centerline (Figs. 5c, e). The collapse of ice rotated outward, indicating an imbalance of outward forces (How and others, 2019). In April 2021, a 50 m x 90 m undercutted thermo-erosional notch collapsed on the non-buoyant left edge as the lake-level dropped, marking the third-largest calving event (Figs. 5g, h). Immediately prior to separation, the non-buoyant area bent forward toward the lake surface, increasing shear as the transverse crevasse widened along the rupture line.

By analyzing the time-lapse camera images, we noticed areas along the glacier front where certain calving styles predominate: Marginal sections of the subaerial ice front overhang due to flaking larger lamellae (Figs. 5 and 6). The ice seracs rotate outwards and collapse into the lake above well-developed waterline notches, mostly observed on the right margin. Stack topples (one event depicted in Figs. 5g, h) are limited to the margins, where the lake bed is steepest (cf. Fig. 4), and the ice is farthest above the water level (Fig. 6). Subaqueous melting near the water surface creates an unstable ice front above the waterline notch along the entire ice front that promotes subaerial calving (How and others, 2019). Since 2019, this observation has been limited to the higher margins. Occasionally, subaqueous calving is found when ice blocks have appeared at the ice front while the ice-front above the waterline has not changed (How and others, 2019).

Such observations show similarities to a marine-terminating glacier in Svalbard (How and others, 2019). However, for lacustrine glaciers, the lack of density-driven circulation combined with the cold water results in very low subaqueous melt rates and inhibited melt undercutting (Truffer and Motyka, 2016). Vertical lake-water profiles of temperature and salinity in summer 2022 show very small differences in temperature and density between lake and runoff water. The lake water was thermally stratified and there was no evidence of meltwater pulses or glacier runoff perturbing the water column near the ice front. The seasonal near-surface water temperature of the lake varies between $0.5 \, ^\circ {\rm C} \hbox{ to } 3.0 \, ^\circ {\rm C}$ (Fig. 12). These observations are consistent with records from proglacial lakes of other lacustrine glaciers (Boyce and others, 2007; Trüssel and others, 2013; Truffer and Motyka, 2016; Sugiyama and others, 2019, 2021). Thus, at lacustrine glaciers, subaqueous melt plays a subordinate role. In addition, in deeper regions along the ice front, where the water is near freezing, a subaqueous ice terrace or ‘ice foot’ may form due to the warmer water in the upper layer (Benn and others, 2007; Sugiyama and others, 2019).

5.5 Atmospheric extremes can trigger calving

In general, ARs advect warm, moist air masses from the subtropics and can be responsible for both a temperature increase and enormous precipitation sums. Almost a third of the total precipitation is associated with a land-falling AR at Schiaparelli Glacier. On average, the occurrence of an AR leads to an increase in mean daily temperature, which reduces the ratio of solid to liquid precipitation from 3/4 to 2/3. In addition, we observe an increase in sensible and latent heat flux at the glacier surface and an increase in incoming longwave radiation (not shown), which may enhance melting on Schiaparelli during AR events, similar to findings from a case study on Brewster Glacier, New Zealand (Kropač and others, 2021). Within five days after the onset of a landfalling AR, calving flux reaches a local maximum in 38% of such cases. We observe an increase in calving flux (14%) and an increase in surface ablation (13%) for a water equivalent height change of 16 mm d⁻¹, resulting in increased mass loss during the days associated with ARs (Table 1). Combined with lower mass accumulation due to greater contribution from liquid precipitation, the MB is 9% more negative. In general, the effect of landfalling ARs on the surface mass balance is controversial. While it can enhance melting by increasing the energy fluxes from liquid precipitation, turbulent fluxes, latent heat release from the advected air mass, and enhanced downward longwave radiation, the glacier may also gain mass by snowfall where temperatures are below freezing (Little and others, 2019; Saavedra and others, 2020; Francis and others, 2021; Wille and others, 2021). At Antarctic ice shelves, strong winds associated with atmospheric extremes and ARs caused exceptionally large calving events triggered by storm surges and a wind-driven ocean slope (Francis and others, 2021, 2022). In contrast, for lacustrine glaciers, the wind-driven circulation controls the temperature stratification of large lakes and determines the calving mechanism (Sugiyama and others, 2021). At Schiaparelli Glacier, the generally low lake temperature and the comparatively small Lago Azul suggest that the temperature stratification in the lake, possibly influenced by atmospheric extremes, plays a minor role in controlling calving.

During wet spells, the daily mean air temperature and incoming longwave radiation increase by 2.5 °C and 12%, respectively, at Schiaparelli Glacier. After the onset of a wet spell, runoff into Lago Azul increases and the calving flux peaks in 46% of such cases. However, results highlight that a response to wet spells does not necessarily imply an increase in calving flux (Table 1). Warm spells increase the runoff by 43%, promoting rapid lake-level rise. In 5% of all cases, a peak in the calving flux can be related to the onset of a warm spell. Here we register an average 29% increase in calving flux. These events also have the largest impact on surface ablation, increasing the total mass loss by 45%, resulting in a tripling of the negative MB term. Dry spells occur only during the cold season, but the mass balance terms are strongly influenced by the lack of precipitation. The resulting lack of accumulation outweighs the inhibited ablation, which doubles the negative MB (Table 1). The lowest surface ablation and calving flux can be attributed to cold spells. Here we record a positive MB.

The largest calving events occurred when the grounded frontal center became buoyant. These calving events are accompanied by rapid lake-level changes (Fig. 6). Prior to the largest calving event on November 20, 2020 (Fig. 5), one of the most rapid changes in lake level occurred. A landfalling AR with a total precipitation sum of 89 mm and a daily mean temperature increase of 4.7 °C caused an intensive lake-level rise of 0.73 m within 7 days. Almost a year later, the second-largest calving event was subject to a similar mechanical process. The sudden ice loss on November 3, 2021, was preceded by a lake-level rise of 0.25 m. This event was accompanied by an AR, which caused a temperature increase of 3.5 °C and 27.5 mm accumulated precipitation within three days. Additionally, this event can be classified as a five-day warm spell. The ice front collapsed on April 29, 2021, while the average daily air temperature dropped by 2.8 °C over the previous three days, reducing glacier melt and thus lake level by 0.29 m.

There are no incidences of short-term acceleration triggered by those calving events, as commonly observed at tidewater glaciers (e.g. Nick and others, 2009; Moon and others, 2014; Schellenberger and others, 2015; Kehrl and others, 2017). More generally, no evidence exists that extreme events such as warm, cold, dry or wet spells affect short-term variations in ice velocities. The ice flow's internal variability seems greater than the event-related potential increase in basal water pressure and its effect on basal lubrication.

6. Conclusion

We analyzed records from in-situ camera systems and remote sensing data based on Sentinel-1 from 2015 to 2022 to estimate ice-front position, surface velocity and the resulting calving flux of Schiaparelli Glacier, a grounded lacustrine glacier in the Cordillera Darwin Icefield. In this study, changes were quantified, underlying concepts were illustrated, and the mechanisms and possible connections with atmospheric factors were shown.

At Schiaparelli Glacier, calving accounts for 4.2% of cumulative ice loss, which is in the range of the lower percentile when comparing estimates for lacustrine glaciers from the Patagonian icefields (Minowa and others, 2021). The complex interactions at the lake-glacier interface make an accurate assessment of purely atmospheric-driven effects difficult. However, we recognize three recurring atmospheric extremes responsible for the increased calving activity: warm spells, wet spells and landfalling ARs. Rapid changes in lake level trigger the largest calving events due to temperature variations and associated glacial meltwater input to the lake. The calving flux increase during warm spells and AR, together with the higher surface ablation, contributes to a considerably greater mass loss during the occurrence of such events. Cold spells favor conditions that reduce both surface ablation and calving, resulting in a positive short-term MB during cold spells on average. Although dry spells tend to occur during the cold season, they contribute on average to a doubling of negative MB as a result of reduced precipitation, while calving flux and surface ablation are below average.

Seasonal ice-front advance in winter can be linked to sustained negative daily air-temperature anomalies characteristic of the cold season. The advance during the cold season is more pronounced when higher ice-flow velocities occur. Changes in the ice front are in sync with the calving flux, with higher calving rates during the retreat phases of the melt season.

Observed changes in Schiaparelli Glacier provide insights in line with established dynamic concepts of marine-terminating and lacustrine glaciers. These include:

1. – Progressive melt-induced thinning is inherent to an acceleration of the terminus section of the glacier. This process destabilizes the ice-front position, which includes an intensification of retreat and seasonality of the ice-front position, while the largest calving events are observed simultaneously.

2. – The seasonal velocity pattern is inversely related to the seasonal runoff. We hypothesize that the englacial discharge system alternates between efficient (melt season) and inefficient (cold season).

3. – Seasonal changes in velocity propagate up-glacier, as do changes in glacier runoff.

The mass loss due to calving response to atmospheric extremes must be considered to evaluate and predict the rapid changes in the lacustrine glaciers of southern Patagonia. Neglecting calving in the total MB would result in a 14% difference in the case of Schiaparelli Glacier. Our results show that individual weather conditions can trigger calving, which in turn can increase the calving flux. Individual glaciers may respond very differently to similar atmospheric forcing. However, the potential atmospheric drivers outlined in this study could also explain recent ice dynamic changes at other sites in Patagonia or elsewhere.