Physical Limnology and Sediment Dynamics of Lago Argentino, the World's Largest Ice‐Contact Lake

Proglacial lakes, whose numbers have been growing around the world, may drive accelerated glacier retreat and provide valuable records of past glacier and climatic changes. Despite their importance, few studies have investigated the sedimentary properties and processes acting within large proglacial lakes. Lago Argentino (LArg) is a 1,500 km2 ice‐contact lake on the eastern flank of the Southern Patagonian Icefield. Here, we describe the results from a detailed analysis of 47 sediment cores obtained throughout this lake basin, supplemented with remotely sensed data. We show that: (a) LArg exhibits a seasonal variation in sediment properties (varves); (b) varve formation results from three distinct processes, driven by seasonal changes in glacial sediment input, seasonal changes in fluvial sediment input, and seasonal variations in lake mixing; and (c) distance from glacier calving fronts provides the first‐order control on sediment grain size and accumulation rate. Our findings highlight the exceptional preservation of annual laminations within proglacial lakes, their potential for reconstructing past glacier changes, and their relevance for forecasting future glacier–lake interactions.

Sediment in proglacial lakes is commonly deposited with seasonal variations in grain size and composition. These annually periodic sediment changes, known as varves where preserved in sedimentary sequences, allow for the construction of high-resolution age-depth models and sedimentation-rate time series (e.g., Ojala et al., 2012;Zolitschka et al., 2015). Varves are commonly formed due to strong seasonal changes in temperature, glacial meltwater discharge, and ice cover, while biogenic processes may also be locally important. Glacio-lacustrine sediments from proglacial lakes may be used alongside sub-aerial data to constrain the timing of past ice advances (e.g., Larsen et al., 2011) and to interpret past climatic changes (e.g., Bertrand et al., 2005;Elbert et al., 2015;Leonard, 1997). Carrivick and Tweed (2013) highlighted the lack of detailed limnological and sedimentological assessments of proglacial lakes. Data are particularly lacking for large, ice-contact lakes which exert the greatest controls on lake-terminating glacier behavior. This knowledge gap is critical because of the potential for these lakes to accelerate deglaciation by absorbing solar radiation, transporting heat to the ice margin, and promoting iceberg calving (Matero et al., 2020;Shugar et al., 2020;Sutherland et al., 2020;Warren & Kirkbride, 1998;Watson et al., 2020).
In this study, we provide a detailed assessment of the sediment properties and sedimentation dynamics within Lago Argentino (LArg), the world's largest ice-contact lake. The varied geography and range of sedimentary and lacustrine environments in LArg make it an excellent case study for understanding the processes occurring in ice-contact lakes.
LArg is the largest of several ultra-oligotrophic, ice-marginal lakes located on the eastern flank of the SPI at 50°2′S, 72°4′W ( Figure 1). LArg's maximum depth reaches over 600 m (Magnani et al., 2019;Sugiyama et al., 2016). While precipitation reaches several meters water equivalent per year on the SPI itself, most of the landscape surrounding LArg is semiarid and receives less than 500 mm of precipitation per year (e.g., Garreaud et al., 2012;Lenaerts et al., 2014).
LArg is located within a tectonically and volcanically active region, and three active volcanoes are located within 100 km of the lake. Bedrock geology varies from metasedimentary and volcanic rocks in the Andean fold and thrust belt to the west of LArg, to a forearc basin sedimentary sequence and plateau basalts to the east of the fold and thrust belt (Coutand et al., 1999).
The surface air temperature at LArg varies from a summer (January) maximum of around 10°C to a winter (July) minimum of around − 2°C, as measured at El Calafate airport. The wind speed follows a similar cycle, with a monthly maximum of 14.2 m/s (in January) and minimum of 6.8 m/s (in June). Local winds are dominantly westerlies, with 70% of measured wind directions in the range 230-320° (SW-NW). Glacier melt rates and flow velocities also vary seasonally, with a maximum in the summer (Minowa et al., 2017;Mouginot & Rignot, 2015;Richter et al., 2016). High wind speeds keep the majority of LArg ice-free even during the winter, when temperatures are below freezing.
Six glaciers feed into LArg, of which three (Upsala, Perito Moreno, and Spegazzini) calve directly into the lake and three (Onelli, Mayo, and Ameghino) calve into smaller peripheral lakes (Figure 1). Upsala glacier drains around 400 cubic kilometers, approximately three quarters of the ice volume in LArg's catchment (Carrivick et al., 2016;Millan et al., 2019), and calves into the longest and deepest fjord in the north-western branch of LArg. Two major rivers are located at the east of the main basin of LArg, one of which flows into the lake (Río La Leona), and the other flows out of the lake (Río Santa Cruz). The daily discharge of Río La Leona has been discontinuously monitored since 1956, and its average March discharge (527 m 3 /s) is more than six times its average September-October discharge (83 m 3 /s). Río Santa Cruz' discharge exhibits a similar variability, with In this paper, we describe the changes in sediment grain size, composition, and accumulation rate around LArg through a detailed analyses of 47 sediment cores. In particular, we highlight the presence of alternating light and dark laminations, and demonstrate that they are annual in nature (varves). Finally, we evaluate the controls on  The key processes affecting the formation and evolution of ice-contact lakes, based on Carrivick and Tweed (2013). These effects are two-way, with the lake in turn affecting the glaciers, land systems, and (for very large lakes) the local climate. sediment accumulation in LArg through a conceptual varve formation model, and discuss its implications for the recording of past glacial changes and glacier-lake interactions.

Core Acquisition
Sediment cores enable recovery and preservation of unlithified stratigraphy across diverse depositional settings in lacustrine and marine systems. LArg is a challenging coring environment: the lake is large (>1,450 km 2 and over 100 km long), deep (over 600 m deep in the fjords: Skvarca et al., 2002;Sugiyama et al., 2016;Magnani et al., 2019), and subject to strong winds throughout the year. In addition, ice surface velocities of over a kilometer per year on Upsala and Perito Moreno glaciers (e.g., Mouginot & Rignot, 2015;Moragues et al., 2017;Sakakibara et al., 2013) result in a steady flux of icebergs in the ice-proximal fjords. To overcome these challenges, we conducted the coring in the winter, when wind speeds are at their lowest, and adapted the Continental Scientific Drilling (CSD) Facility Kullenberg piston-coring system and a gravity surface coring device for deployment onboard a local vehicle ferry. We identified coring sites ( Figure 3) through a preceding seismic reflection survey to cover all major LArg depocenters and depositional environments (Magnani et al., 2019).

Core Analyses
After collection, cores were transported to the CSD Facility for analysis. We performed initial whole-core scans of magnetic susceptibility and gamma-beam density, and split and imaged the cores with a high-resolution optical line-scanner (50 pixels per millimeter). We then conducted near-UV to near-IR spectrophotometric and magnetic susceptibility scans on each split core face. We described the key sediment properties, lithological units, and sediment structures of the 47 cores in logging software PSICAT (Reed, 2007).
We designed the detailed sedimentological analyses described below based on this initial core logging, noting the laminated nature of the sediment (Figure 3). We investigated selected core intervals using scanning X-ray fluorescence (XRF) and laser-diffraction particle size analysis (LD-PSA). Grain-size analyses were conducted on a Horiba LA-920 LD-PSA, using optical model 116a010I. LD-PSA analyses were run directly on core sediment with no sample preprocessing due to the very low organic matter content. XRF analyses were conducted using a Cr source ITRAX XRF scanner.
We selected 24 intervals from 18 characteristic sites for sub-lamina scale (0.2-0.5 mm) XRF analyses, representing the major lake depositional settings. High-resolution XRF scans of 15-30 laminae are nested within whole core-length scans at lower spatial resolution (5-10 mm), providing broader geochemical context. We analyzed XRF data as counts-per-second and count ratios, without attempting external calibration to absolute concentrations. Where possible, we coupled lamina-scale XRF with LD-PSA grain size measurements. We conducted 360 LD-PSA measurements, with the majority of results based on three replicate sample runs. We collected several hundred petrographic smear slides, and used petrographic descriptions of smear slides to supplement grain size and elemental analyses.
To account for multi-modal LD-PSA results, we used a Gaussian mixture model (e.g., McLachlan & Peel, 2004) to separate the signal into multiple normally distributed components. We first fit the distribution with a single Gaussian peak, and measure the coefficient of determination (R 2 ) between this fit and the original data. If the R 2 value is less than 0.95, we discard this fit. We then fit the data with two Gaussian peaks, and measure the coefficient of determination against the data. This process is repeated with additional Gaussian peaks until the R 2 value is greater than 0.95, or the preset maximum threshold of 4 Gaussian contributions is reached. The resulting fit to the LD-PSA results is given by Equation 1 below, with n equal to the number of Gaussian contributions. This allows a single grain size analysis to be separated into a series of constitutive endmembers E g , each described by a mean (b i ), a standard deviation (c i ), and an amplitude (a i , describing the relative contribution of each endmember to the total): (1) We counted the number and thickness of alternating light and dark laminations in our cores using sliding-window autocorrelation of digital core scans (CountMYvarves, see a full description in Van Wyk de Vries et al., 2021). Manual counting of close to 60,000 mm-scale laminations (across the 47 cores) is time consuming, subjective, and challenging to reproduce (e.g., Sprowl, 1993;Tian et al., 2005). CountMYvarves conducts multiple simultaneous lamina counts, which allows for uncertainties to be propagated into resulting lamina-depth-count and lamina-thickness curves. Despite the software-package name, the lamina-counting model requires no assumption about lamination periodicity, which must be verified externally prior to chronological interpretation. Mean paired-lamina thickness (adjacent light and dark layers) is obtained at each core locality by taking the depth-average of all available lamina thicknesses in each core. This method is more robust than dividing total lamina count by core length, as it is not sensitive to disrupted or missing sediment, such as at the break between core sections.

Cesium Chronology
We dated two cores from the main lake basin and one from the brazos ("arms": lacustrine glacial fjords) using the 1963-1965 137 Cs peak (Caffau et al., 2021;Piret et al., 2022). We measured 137 Cs counts at Macalester College using a Canberra high-resolution, hyperpure germanium coaxial gamma detector (GC1518). We estimated a relative detector efficiency of ∼40% and a resolution of 0.27 keV using a Spectrum Technologies 10 Ci 137 Cs standard. Blank runs of u-channels and plastic boxes used to hold the sediment showed no measurable 137 Cs activity. We placed samples in the detector for 12-48 hr, depending on sediment volume. We estimated an analytical precision of ±5% based on duplicate runs.
We converted counts to 137 Cs activity (a, in Bq/kg): in which c T are total counts, c N is background noise, t is time in s, and M is sample dry mass in kg. We calculated c T as the sum of full peak-width at half maximum counts (660.9-662.8 KeV) and c N as the average of off-peak spectrum counts (630.1-631.6 KeV and 690.1-691.4 KeV). We measured dry sediment mass (M) by oven-drying samples at 110 • C for 24 hr.

Remote Sensing
We use satellite observations to measure relative suspended sediment concentration (rSSC) trends and to assess the presence of internal waves at LArg. We assessed rSSC trends with a modified version of the Ulysses Water Quality Viewer algorithm (Zlinszky & Padányi-Gulyás, 2020). This algorithm measures rSSC directly from band 5 (wavelength 705 nm) of optical satellite Sentinel 2. Where the ratio of suspended sediment to chlorophyll is high (as is the case in ultra-oligotrophic LArg), the radiance of the ∼700 nm band scales close to linearly with water suspended sediment concentration (SSC; Nechad et al., 2010). Due to the limited in-situ SSC measurements to calibrate remotely sensed measurements, we interpret satellite-derived rSSC as a relative measure of near-surface lake water SSC and do not attempt a calibration to absolute SSC values. We obtained remotely sensed suspended sediment values from all 77 cloud-free Sentinel-2 images from 2015 to 2020, and threshold each rSSC image to exclude clouds and icebergs.
Internal waves are oscillatory gravity waves that form in lakes with high vertical density gradients (e.g., at the thermocline: Antenucci & Imberger, 2001;Filatov et al., 2012). We assess the presence of internal waves in LArg using synthetic aperture radar (SAR) satellite Sentinel-1 following the approach described by Klemas (2012). The propagation of internal waves modulates the lake's surface wavefield, causing long-wavelength zones of alternating rough and smooth surface water. Rough water surfaces produce a brighter radar backscatter than smooth surfaces, which we measure using Sentinel-1's C-band (wavelength of 5.5 cm) SAR in vertical-vertical (VV) polarization. We manually confirm the presence or absence of internal waves in all 284 Sentinel-1 radargrams of LArg (as of Feb 2020).

Core Acquisition
We recovered a total of 108 m of core from 47 sites in August-September 2019 ( Figure 3). Core-recovery length ranged up to 6.4 m at shallow water (<275 m depth) sites and to 5.5 m at deep water (>275 m depth) sites. The most ice-proximal core was collected 4.3 km from the front of Glaciar Upsala in a region of the lake only deglaciated in 1995 (Moragues et al., 2017).

Core Analyses
Lago Argentino's sediment is dominated by alternating light and dark laminae, with lamina pairs ∼1 mm thick in the main basin of the lake and scaling up to 5-10 mm in the more ice-proximal brazos. The laminated sediment is punctuated by several tephra layers with thicknesses ranging up to 10 cm. Lamina counting using countMYvarves reveals several thousand lamina in the main basin cores and several hundred lamina in the brazos (Van Wyk de Vries et al., 2021). Sediment from the main basin of the lake shows continuous and uninterrupted settling of micron-scale particles, with only minor variation with depth. Sediment from the brazos is coarser grained, with laminated sediment interspersed with cm-scale, non-laminated sand layers. Cores from the easternmost zone of the main lake basin near the source of the Río Santa Cruz and mouth of the Río La Leona also exhibit thicker laminae and coarser grain sizes. There is a strong color contrast from light colored sediment in the main lake basin to dark sediment in the glacial fjords and eastern margin of LArg ( Figure 3). Full resolution core scans and associated stratigraphic information are available at 10.5281/zenodo.5815107.
Analysis of petrographic smear slides shows two dominant grain size fractions throughout the lake, one in the sub-μm range and a second in the 1-10 μm range. The coarsest grains in each slide have diameters of 20-200 μm. The coarser-grained fraction is primarily composed of plagioclase, quartz, one to two pyroxenes, biotite, opaque oxides, and volcanic glass (fresh or devitrified). The concentration of volcanic glass varies from ∼100% in macroscale tephra layers to an average of 10%-20% in background material. Rare pennate and centric diatoms are present, with a higher abundance in the main lake basin than the brazos. The fine-grained fraction is too fine to identify mineral composition using a standard petrographic microscope.
LD-PSA-derived grain-size histograms are consistent with petrographic observations, and show that LArg's sediment is dominantly bimodal. Gaussian mixture model separation shows that 80.5% of samples are composed of two endmembers. The fine endmember ( Figure 4a) exhibits a sharp peak in the grain size histogram, with more than 60% of measurements within the range 0.50-0.56 μm. The coarse endmember ( Figure 4b) exhibits a broader peak, with 60% of values in the range 1.8-3.4 μm. Grain sizes fine with distance from the front of Upsala glacier ( Figure 4c). The fine fraction remains constant at approximately 0.5 μm throughout the entire lake basin, while the coarse fraction declines from around 3-5 μm in the western fjords to 1-2 μm in the main lake basin. LD-PSA results were also grouped based on the lamina coloration. Both laminae sets are composed of similar overall grain sizes, with "light" laminae composed of slightly finer grain sizes, and "dark" laminae of coarser grain sizes (Figures 5a and 5b).
Lamina thicknesses show a similar trend to the grain size data: values are highest close to the glacier front, and low in the main basin of the lake (Figures 4e and 4f). The median paired (light + dark couplet) lamina thickness is 4.8 mm (interquartile range 2.5-6.3 mm) and 1.5 mm (interquartile range 1.2-2.0 mm) in the brazos (glacial fjords) and main lake basin, respectively. Sedimentation rate away from glacier front averaged over the past 50 years shows a moderate fit to an exponential decay (Boldt, 2014;Koppes & Hallet, 2002), with an R 2 of 0.37 and 0.40 for Upsala and Perito-Moreno respectively.
XRF results show that sediment deposited in the eastern-most margin of LArg is compositionally distinct from the majority of the lake basin ( Figure 6). Sediment from near the E margin of the lake exhibits higher Ca/K and Ti/K than the main basin and western portion of the lake. Coupled sub-mm resolution XRF results reveal that the sediment composition also varies between light and dark lamina. Dark laminae exhibit high Ca/K ratios (generally in the range 0.5-0.7) and light laminae low Ca/K ratios (in the range 0.35-0.45; Figure 5c).

Cesium Chronology
All three cores analyzed exhibit above-background 137 Cs activity (see profiles in Figure S1 in Supporting Information S1). Core 11A has a maximum 137 Cs peak of 0.39 ± 0.02 Bq/kg at a depth of 127 ± 5 mm, and negligible 137 Cs concentrations below 147 ± 5 mm. The top of core 11A is located at a depth of 90 mm, with the core barrel filled with Zorbitrol (Tomkins et al., 2008) above this level. We merge the samples from cores 20A and 20B into a single record using marker horizons in both cores. Core 20A has a maximum 137 Cs peak of 0.31 ± 0.02 Bq/kg at a depth of 157 ± 6 mm, and negligible 137 Cs concentrations below 183 ± 5 mm. The top of core 20A is located at a depth of 105 mm, with the core barrel filled with Zorbitrol (Tomkins et al., 2008) above this level. Core 20A exhibits a second 137 Cs peak near the top of the core at a depth of 108.5 ± 2.5 mm mm, with a 137 Cs peak of 0.37 ± 0.02 Bq/kg. Core 27A has a maximum 137 Cs peak of 0.57 ± 0.03 Bq/kg at a depth of 29.85 ± 2.5 mm, and negligible 137 Cs concentrations below 34.85 ± 2.5 mm. The top of core 27A is located at a depth of 105 mm, with the core barrel filled with Zorbitrol (Tomkins et al., 2008) above this level. . Summary of laser-diffraction particle size analysis (LD-PSA) and lamina thickness results. Panels (a and b) show the grain size histograms for the finest and coarsest endmembers, respectively. Panel c shows the change in grain size with distance from the largest glacier (Upsala), and (d) shows the change in coarse endmember contribution with distance from glacier front + symbols represent the coarse endmember, while o symbols represent the fine endmember. Panels (e and f) show the change in lamina thickness (averaged over the past 50 years) with distance from the front of Upsala glacier and Perito Moreno respectively. The dashed lines represent an exponential best fit to the data.

Remote Sensing
The highest rSSC is found in the northwest fjord of the lake (Brazo Upsala-Norte), close to the calving front of Upsala glacier (Figure 7a).
Step-wise decreases in rSSC with distance from glacier front are observed at bathymetric highs and fjord constrictions, particularly at the boundary between the northern fjord arm and main lake basin ("Boca del Diablo") where the lake is less than a kilometer wide. A high rSSC plume is observed on the eastern margin of the lake, as well as locally high rSSC on the south margin of the lake. We divide the Sentinel-2 for three cores, and panel (c) shows the variation in Ca/K ratio between light and dark laminations in a segment of core 24a. data set into a "summer" (November-April) and a "winter" average (May-October). Lake surface rSSCs are generally the highest during winter months, with the most pronounced seasonal changes proximal to the glacier front ( Figure 7c). This trend is inverted for the eastern-most 10 km of the lake, in which rSSC is higher in summer than in winter months.
We identify internal waves in multiple Sentinel-1 images. Long wavelength (∼10 km), alternating bands of rough and smooth lake water surface are present in 34.5% of radargrams (Figure 7d). Bands are most common in the austral summer and autumn, and are rare in austral winter and spring. Internal waves are most common in January and February (52% and 50% of radargrams) and least common in May and June (14% of radargrams). A full list of positive internal wave identifications is given in Table S1 in Supporting Information S1.

Discussion
To explore past glacial changes through the interpretation of proglacial lake sedimentary records, we must understand the first-order controls on sediment accumulation and redistribution throughout the lake basin. We explore this at LArg through the interpretation of our lake core sequence and satellite imagery. We first test whether the alternating light and dark lamina are varves. A variety of processes may give rise to cyclically laminated sediment, including the seasonal cycle (e.g., Zolitschka et al., 2015), tides (e.g., Kvale et al., 1995), orbital cycles (e.g., Cramp & O'Sullivan, 1999), and regular influx of turbidity currents (e.g., Carrivick & Tweed, 2013). Other sub-annual and multiannual cycles occur at LArg, such as an 85-min period fundamental seiche (Richter et al., 2016) and multiannual El-Niño Southern Oscillation and Southern Annular Mode (Ariztegui et al., 2007;Kaplan et al., 2020;Li et al., 2013).
We review the lamina structure, core composition, remotely-sensed data, and independent geochronological evidence collected in this study to investigate the causes of cyclical sediment lamination: 1. LArg's sediment is primarily composed of alternating light and dark laminations (Figure 3). Periodic laminations are present throughout the entire lake basin, and persist for thousands of cycles in the main basin and at minimum several hundred cycles in the brazos. The thickness and grain size of these laminations decreases with distance from glacier fronts (e.g., Figures 4d-4f). The thickness of the lamina couplets (light + dark layer) varies within the range of 1-10 mm across the majority of the lake 2. LArg's sediment is bimodal (Figure 4), composed of one very fine endmember (∼0.5 μm) and one coarse endmember (∼2-5 μm). This bimodality is present in both light and dark laminae (Figures 5a and 5b) 3. Light laminae have a finer grain size and lower Ca/K ratio, whereas dark laminae have a coarser grain size and higher Ca/K ratio 4. LArg's sediment composition is similar across the main lake basin, except for the easternmost core site, which is higher in Ca and Ti ( Figure 6) 5. LArg's rSSC is highest at the glacier front and lowest in the main lake basin. In addition, rSSC is higher in the winter than in the summer, except in the easternmost region of the lake 6. Internal waves are common at LArg, and are more common in the summer than in the winter. Climatic data from the shores of LArg shows strong seasonal cycles of temperature and wind speed To further test the hypothesis that the light and dark sediment laminae are seasonal, we compared the down-core position of the 1963-1965 137 Cs peak to the number of light-dark lamina couplets. The tops of cores 11A, 20A, and 27A are well preserved -we estimate that at most 10 laminae are disrupted and uncountable at the sediment-water interface for 11A and 20A, and at most 5 laminae are disrupted for 27A (where laminae are thicker). Cores were collected in 2019, meaning that 54-56 annual layers should have been deposited since 1963-1965. Accounting for up to 10 disrupted lamina, 44-56 annual layers should be present (49-56 annual layers for 27A). We take the position of the 1963-1965 137 Cs horizon as between the peak in 137 Cs activity and the point below this peak with background activity ( Figure S1 in Supporting Information S1), with depths of 47 ± 10 mm, 65 ± 13 mm, and 218.5 ± 25 mm for cores 11A, 20A, and 27A respectively. The depth of the 1963-1965 137 Cs horizon is within uncertainty of the lamina-count derived depth for the 1963-1965 horizon in all three cores dated ( Table 1).
The LArg laminations are too numerous and too regular to result from longer period (e.g., multi-annual, multi-decadal) climatic changes and too thick to result from 1.5 hr period seiches. Climatic and remotely sensed data supports seasonal variation in sediment dynamics. In conjunction with independent 137 Cs dating, we conclude that the laminations in LArg are annual-they are varves. Lamina thicknesses therefore provide an annual resolution sedimentation rate record.
Unlike the original light summer layer and dark winter layer varve norm (DeGeer, 1912), summer layers in LArg are darker than winter layers. At LArg, coarser-grained dark layers represent summer settling and finer-grained light layers represent winter settling. The color of clastic glacial varves depends on the bedrock of the lake catchment area, and mineral composition of the summer and winter laminae (O'Sullivan, 1983;Ridge & Larsen, 1990;Ridge et al., 2012;Zolitschka et al., 2015). At LArg, the color of dark layers is contributed by a seasonally higher fraction of coarser grained, dark colored mineral such as pyroxenes rather than organic matter. High K/Ca in light-colored winter layers is consistent with settling of light-colored, K-rich clays (Figures 5a-5c).
We next build a conceptual varve-formation model for LArg and explore the mechanisms that give rise to these seasonal sediment changes. Compositional ( Figure 6) and sedimentological data suggest that at least two distinct processes generate varves in LArg: one driven by fluvial changes affecting only the eastern portion of the lake, and another driven by glacial changes in the remainder of the main basin and brazos of the lake. Hydrodynamic data in LArg is sparse, limited to the summer surveys conducted by Sugiyama et al. (2016) close to the glacier calving fronts and isolated water measurements in the east of the lake for the Río Santa Cruz hydroelectric dam project. To overcome this data paucity, we incorporate remote sensing results into the conceptual varve-formation model.

Varve Formation: Fluvial Influx
XRF compositional data, LD-PSA grain size data, and sedimentation rate data from semi-automated lamina counting all show a strong difference between the eastern-most ∼10 km of LArg and the remainder of the lake. The results show that the majority of sediment in this eastern margin of LArg is derived from Río La Leona. Discharge near the outlet of Río La Leona from Lago Viedma varies seasonally by almost an order of magnitude. The greater late-summer discharge translates to greater fluvial erosion and sediment transport potential. Remotely sensed rSSCs from the easternmost LArg basin show a higher rSSC plume in the summer relative to the winter (Figure 7c). Río La Leona's sediment load is primarily derived from the non-glaciated catchment between Lago Viedma and LArg, while Río La Leona's discharge is primarily controlled by ice melt in its catchment area (e.g., Figures 1 and 3). The higher summer fluvial sediment influx into LArg may represent both increased summer sediment erosion from this non-glaciated region, and increased summer sediment transport potential along Río La Leona due to the higher discharge. The varve-formation mechanisms in this most ice-distal portion of LArg are more comparable to those in non-ice-contact proglacial lakes, in which the lake is connected to the glaciers by a Core name Depth ( 137 Cs counts) Depth (lamina counting) 11A 47 ± 10 mm 52.5 ± 9.5 mm 20A 65 ± 13 mm 60.5 ± 7.5 mm 27A 218.5 ± 25 mm 202.5 ± 14.5 mm Table 1 Comparison Between Depth of 1963-1965Horizon From 137 Cs Counts and Depth of 44-56 Laminae for 11A and 20A, and 49-56 Lamina for 27A (1963-1965 fluvial system -such as that described in nearby Lago Plomo, on the flank of the Northern Patagonian Icefield (e.g., Elbert et al., 2013Elbert et al., , 2015.

Varve Formation: Glacial Influx
In the remainder of the lake, sediment is primarily sourced from glaciers. Figure 5 shows the seasonal variation in sediment properties, with dark-colored lamina enriched in coarser grains and light-colored lamina enriched in finer grains. Two conditions must be fulfilled for varves to form: first the quantity, composition, or other properties of sediment accumulated must vary throughout the year; and second, any disturbances at the sediment-water interface must remain small enough to preserve annual-resolution sediment changes (O'Sullivan, 1983;Zolitschka et al., 2015). We use these sub-annual resolution grain-size measurements to investigate two processes that may control varve formation: a seasonal cycle in glacial sediment delivery and a seasonal cycle in the degree of lake mixing.

Seasonal Cycle in Sediment Delivery
Wind speed, surface temperature, and glacial melt rates all vary seasonally at LArg (e.g., Minowa et al., 2017;Richter et al., 2016). We do not have in-situ measurements of seasonal glacial sediment inflow from LArg, but data from other glaciers around the globe (e.g., Haritashya et al., 2006;Riihimaki et al., 2005) show a seasonal variation in sediment export from the glacier front, with sediment flux highest over the summer months. We examine the effect of a seasonal cycle on sediment accumulation at LArg.
The majority of grain size diameters in LArg range from 0.5 to 5 μm (Figure 4), and are sufficiently fine to not induce turbulence while settling. Stokes' Law can therefore be used to compute settling velocities, w s : where ρ s is the particle (sediment) density, ρ is the fluid (water) density, μ is the dynamic viscosity, g is acceleration due to gravity, and D is the particle diameter. We use ρ s = 2,650 kg/m 3 , ρ = 1,000 kg/m 3 , and μ = 0.001 N⋅s/m 2 , respectively, which are standard values for silica particles in freshwater. Particles with a grain size of 5 μm will have a settling velocity of 708 m/yr, while 0.5 μm particles have a settling velocity of only 7 m/yr.
We bracket possible settling timescales using the endmembers of fully non-mixed water ("plug flow") and fully mixed water ("stirred reactor"; see Appendix A). In the 200 m deep main basin of LArg, 90% of 5 μm particles will settle in 0.25 years for non-mixed water and 0.65 years for fully mixed water. Holding all other factors constant, the settling time for 0.5 μm particles will be two orders of magnitude longer than for 5 μm particles: 25 and 65 years, respectively.
The presence of >5 μm scale particles in the center of the lake basin requires that horizontal mixing of the entire lake occurs on a timescale shorter than the settling timescale of these grains (3-8 months). The long settling timescale of sub-micron particles results in them being almost insensitive to any seasonal cycle in sediment influx, while the sub-annual settling timescale of 5 μm particles allows their settling rate to vary in response to any seasonal cycle in sediment influx (Figure 8). A seasonal variation in sediment influx may therefore lead to the settling of sediment with a seasonally varying grain size distribution, even where the distribution of grain sizes within the sediment source is constant.
A seasonal cycle of sediment influx, with higher summer sediment input and lower winter sediment input, may explain the formation and characteristics of varves in LArg based on sedimentary data alone. However, available ice-velocity data shows only a minor seasonal variation in ice flow speed (Minowa et al., 2017;Mouginot & Rignot, 2015), suggesting that the sediment influx also may not experience a strong seasonality. In addition, remotely sensed observations show that rSSCs are higher in the winter than in the summer. This is inconsistent with a seasonal cycle in sediment influx alone, for which we would expect higher rSSCs in the summer during times of greater sediment delivery. We therefore consider a second mechanism for varve formation, relating to a seasonal cycle in lake mixing.

Seasonal Cycle in Lake Mixing
The timescale of vertical mixing in a lake affects the settling of sediment, and provides another potential varve forming mechanism. The degree of vertical mixing in a lake is primarily determined by a balance between density stratification (due to temperature, suspended-sediment, or salinity contrasts) and wind mixing (e.g., Wüest & Lorke, 2003). Density stratification lowers the center of mass of the lake, reducing its potential energy. Wind adds kinetic energy to the lake system, acting to disrupt this stratification and mix the lake (e.g., Figure 9). A variety of parametric indicators have been used to describe the stability of stratification and degree of mixing in lakes, including the Schmidt stability (St), Lake Number (LN), Wedderburn Number (W), and Richardson Number (Ri). The Richardson number is the ratio of buoyancy (Brunt-Väisälä frequency) to shear flow, which may be modified to describe the stratification regime (e.g., Fischer et al., 1979;Robertson & Imberger, 1994), with , ρ h , ρ e , and ρ a are the densities of the lake average, hypolimnion, epilimnion and air respectively; C D is the drag coefficient of the water surface; and U 10 is the wind speed at 10 m altitude. We calculate Ri for likely density-contrast ρ h − ρ e (Sugiyama et al., 2016); and wind-speed (U 10 ) values for LArg (Figure 10). ; and Regime C, unstable stratification with possible full-lake vertical mixing . Regimes A, B and C are referred to as "Stratified", "Partially stratified", and "Mixed" respectively for the remainder of this paper. LArg is likely in the "Partially stratified" regime the majority of the time, with possible excursions into the "Stratified" regime during calmer than average summer periods (Figure 10b) and 'Mixed" regime during windier than average winter periods (Figure 10c). Higher median summer Ri suggests that the stronger summer thermal gradient promotes lake stratification despite Figure 9. State of the water column for different degrees of mixing. In a well-stratified lake (a), a strong vertical gradient in water properties exists and significant vertical mixing only occurs at the upper and lower boundaries (with surface mixing due to wind action the strongest). Particle concentration in the water column will decrease close to linearly. In an intermediate-mixed lake (b), a diffuse metalimnion will separate surface and bottom waters, and internal waves may form. For a wellmixed lake (c), full-lake-depth eddies will homogenize vertical temperature and density gradients. Particle concentration in the water column will decrease through time as an exponential decay. the higher wind speeds. This is consistent with more frequent observations of internal waves in the summer than in the winter, as these require an internal density gradient to form (Fischer et al., 1979).
Particles will settle more rapidly if their settling timescale is shorter than the local mixing timescale of the lake, and more slowly if their settling timescale exceeds that of the local mixing timescale. It follows that a temporary decrease in the lake mixing timescale (more vigorous mixing, lower Ri) will temporarily reduce sedimentation. Particles with high settling velocities are more sensitive to this effect as the temporary decrease in the lake mixing timescale will account for a larger proportion of their total settling time. A seasonal cycle in lake mixing can thus  Fischer et al. (1979), Robertson and Imberger (1994), and representative values for the main basin of LArg (red shading). The black grids divide the parameter space into the mixed (regime A), partially stratified (regime B), and stratified (regime C). (b and c) show histograms of Ri for typical summer (high wind speed, high density contrast) and winter (low wind speed, low density contrast) values. The arrows show the median value for each season. generate a seasonal cycle in sediment deposition and seasonal fractionation of grain sizes, even in the absence of a seasonal cycle in sediment influx or its particle-size distribution (Figure 9).
We calculate the approximate degree of lake mixing, or mixing timescale, necessary to generate the observed seasonal grain size fractionation. The settling timescale of a sediment particle (in non-mixed water) is given by the ratio of lake depth to particle settling velocity: and the mixing timescale is given by: where K is the eddy diffusivity, itself a measure of the local degree of mixing due to turbulence within the water. Layers of sediment deposited in the summer are enriched in grains larger than around 2-4 μm, and conversely winter layers are depleted in grains larger than around 2-4 μm (Figures 5a and 5b). Using the above equations for particle settling velocity, particle-settling timescale and lake-mixing timescale, we calculate the eddy diffusivity for which the mixing timescale is equal to the settling timescale. For 2-4 μm grains, this requires an eddy diffusivity of K = 1.8-7.2 × 10 −4 m 2 /s. This is many orders of magnitude larger than the particle molecular diffusivity of ∼10 −12 − 10 −13 m 2 /s for micron-scale grains, estimated from the Stokes-Einstein equation.
A seasonal cycle of lake mixing, with lower summer lake mixing and greater winter lake mixing, may also contribute to the formation and characteristics of varves in LArg. Greater mixing of the water column in the winter would bring deeper, sediment-rich water to the surface, explaining the observation of higher winter rSSCs. We do not, however, have any direct data on LArg's seasonal changes in lake mixing timescale, and must rely on indirect observations of internal waves ( Figure 7) and lake Richardson number ( Figure 10).
Overall, varves at LArg may form via either seasonal cycle in sediment delivery, a seasonal cycle in lake mixing, or a combination of both of these processes. Both mechanisms would enhance summer settling of coarser grained sediment when compared to deposition during winter. Seasonal changes in glacier flow (e.g., Moragues et al., 2017;Mouginot & Rignot, 2015;Sakakibara et al., 2013) suggest that some seasonal cycle in sediment delivery is likely. Remote sensing observations and local climatology support a seasonal cycle in lake mixing. We therefore consider that varves likely result from a combination of the two processes across most of LArg.

Synthesis of Processes in Large Proglacial Lakes
Previous studies have investigated the processes occurring in proglacial lakes, such as Rubensdotter and Rosqvist (2009) in Sweden; Schiefer and Gilbert (2008); Richards et al. (2012) in Canada; Larsen et al. (2011) in Iceland;and Boes et al. (2018) in the USA. Carrivick and Tweed (2013) also provide a literature review describing the formation of and processes occurring in proglacial lakes, noting that "It is absolutely crucial to understand sedimentation within proglacial lakes for deciphering the Quaternary record of glaciation and glacier dynamics". Prior studies are, however, strongly biased toward small alpine glacier lakes two to three orders of magnitude smaller than LArg. Our results show that many processes occurring in LArg are comparable to those in smaller lakes, but that the relative importance of these processes may differ ( Figure 11): 1. Excluding regions proximal to fluvial deltas, sedimentation in LArg is almost entirely sourced from glaciers (e.g., Figures 6 and 7a) 2. Settling of micrometer to sub-micrometer scale particles dominates the sediment budget, even within much of the ice-proximal freshwater fjords 3. Sedimentation rate away from the glacier fronts averaged over the past 50 years shows a moderate fit to an exponential decay (R 2 of ∼0.4), consistent with that observed in ice-proximal marine settings (Boldt, 2014;Koppes & Hallet, 2002). The sediment grain size and rSSCs also show a similar decrease away from the glacier fronts, with the decrease from the front of Upsala glacier being the most marked (Figure 4) 4. LArg lies close to the boundary between stable density stratification and full lake mixing, with a higher degree of winter lake water mixing 5. Seasonal variations in glacial and fluvial sediment influx and lake mixing generate cyclical changes in sediment properties (varves). Particle size distribution is bimodal, with a coarser fraction settling on monthly timescales and a sub-micron scale fraction which may remain suspended in the water column for decades A positive correlation has been identified between lake temperature and glacier calving rate (e.g., Warren & Kirkbride, 2003), and persistence of warm waters at the glacier calving front can drive rapid retreat. Consequently, the degree of lake mixing is important for assessing glacier-lake interactions, as it determines the distribution of warm water at the glacier front. Unlike in many marine systems, subglacial meltwater is denser than ambient lake-water (due to its low temperature and high suspended sediment load). If the lake is well-mixed, glacial meltwater will be rapidly diluted into the water column and the entire ice front may experience seasonally high water temperatures. If the lake is density-stratified, this cold meltwater will remain trapped in the hypolimnion, and warming will be concentrated in the epilimnion. In this latter case, only the portion of the glacier front above the thermocline will experience seasonally high temperatures. Sugiyama et al. (2019) identified an ice terrace at the front of Gray glacier, Patagonia which they interpret as the result of this vertically localized glacier frontal melt in a stratified lake water column.
In the absence of basin-wide hydrodynamic and water-column temperature data, our sedimentary and remotely sensed datasets provide preliminary constraints on the degree of lake mixing in LArg. The presence of >5 μm scale particles in the center of LArg, over 50 km from the glacial sediment source, suggests that horizontal mixing of the entire lake occurs on a sub-annual timescale (shorter than the settling timescale of these grains). Decreased lake-surface rSSCs and more abundant internal waves in the summer are consistent with a stratified water column, Figure 11. Summary diagram of the processes occurring in Lago Argentino.
with vertical mixing of suspended sediment inhibited by the vertical density contrast. High influxes of suspended sediment at depth from glacial meltwater likely reinforces this stratification by increasing the sediment concentration and bulk density of lake-bottom waters.
Our results show that mixing and stratification processes in LArg are complex, and likely modulate the seasonal distribution of warm water at the glacier calving fronts. LArg's sediment provides a salient, high-resolution record of glacial sedimentation over the past hundreds to thousands of years. Despite all the complexity of the lake with multiple branching freshwater fjords, variable bathymetry, and rapidly evolving glaciers, sedimentation rates can be approximately described as an exponential decay away from the largest glacier in the basin, Upsala. Varves form by three different processes, the first linked to an annual cycle in lake mixing, the second to an annual cycle in glacial sediment delivery, and a third to changes in fluvial sediment transport (in the far East of the lake). Our data set is limited in the most ice-proximal regions of the lake, where future work is required to evaluate the coupling between the lake, sediment, and glacier front.

Conclusions
We collected 47 deep-water lake cores from LArg, the world's largest ice-contact lake. We then conducted detailed XRF spectroscopy, particle-size analysis, and semi-automated lamina counting to evaluate the processes occurring at LArg. We identify and validate the presence of varves in LArg, with dark summer layers exhibiting coarser grains and higher Ca/K ratios, alternating with light-colored winter layers exhibiting finer grains and lower Ca/K ratios. The grain size of LArg's sediment is bimodal, made up of a 0.5 μm grain size fraction alongside a 1.8-3.4 μm fraction that fines with distance from the glacier front. Varve-count-derived sedimentation rates decrease exponentially with distance from the glacier front. Combining our sedimentological evidence with remotely sensed rSSCs, we describe a model for varve formation at LArg. Three processes give rise to the annual laminations: seasonally varying vertical water-column mixing and seasonally varying glacial sediment input dominate in the majority of the lake, and seasonal changes in fluvial sediment input dominates in the eastern 10% of the main lake basin. These results highlight the importance and the complexity of large proglacial lakes and their relationship to the glaciers that calve into them.

Appendix A: Sediment Settling With and Without Turbulent Mixing
The time required for particles to settle depends on the depth of the lake, the settling velocity, and the degree of vertical mixing, which can counteract settling by redistributing sediment through the water column. Because we lack direct measurements of water column turbulence and sediment-concentration profiles at Lago Argentino, we investigate particle settling and deposition under two endmembers cases: fully unmixed "plug flow" and fully mixed "stirred reactor" conditions. These approximate the limits of fully laminar flow with negligible molecular diffusion (plug flow) and fully turbulent conditions in which vertical turbulent mixing is much more vigorous than particle settling (stirred reactor).
In general, the sediment concentration, C, and depth over which sediment is distributed, h s , co-evolve as: where w s is settling velocity. In plain terms, this equation means that the amount of sediment in a column of water will decrease as sediment settles out of the column and deposits on the bed, and that this can be accomplished through either a change in concentration, C, or a change in the thickness of the water column occupied by the sediment, h s . We define a coordinate system in which w s is positive downwards, and for simplicity consider the concentration profile to be uniform across h s .
In the laminar "plug flow" endmember case, we hold concentration constant. As a result, the entire sediment column settles at a rate w S : We assume that at t = 0, the sediment is uniformly distributed across the height of the water column, h. For times larger than ℎ , sediment concentrations will be zero. We can algebraically rearrange this equation to solve for the time t required to deposit a certain fraction of this initial sediment.
For a fully turbulently mixed water column (i.e., that with a vertical eddy diffusivity much higher than the particle settling velocity), the sediment concentration will remain uniform but will decrease with time as particles settle and mix. We represent this by holding h s = h constant in Equation A1, which leads to the solution: where C 0 is the initial sediment concentration at t = 0. As before, we can rearrange this equation to solve for the time required for a given fraction of C 0 to settle onto the bed.

Data Availability Statement
Data on the LArg cores are available from the CSD Facility at https://cse.umn.edu/csd/projects. All Sentinel-1 and 2 images are freely available online from several sources, including the Sentinel EO-Browser (apps.sentinel-hub.com/eo-browser