Subglacial carbonate deposits as a potential proxy for a glacier's former presence

The retreat of ice shelves and glaciers over the last century provides unequivocal evidence of recent global warming. Glacierets (miniature glaciers) and ice patches are important components of the cryosphere that highlight the global retreat of glaciers, but knowledge of their behaviour prior to the Little Ice Age is lacking. Here, we report the uranium–thorium age of subglacial carbonate deposits from a recently exposed surface previously occupied by the disappearing Triglav Glacier (southeastern European Alps) that may elucidate the glacier's presence throughout the entire Holocene. The ages suggest the deposits' possible preservation since the Last Glacial Maximum and Younger Dryas. These thin deposits, formed by regelation, are easily eroded if exposed during previous Holocene climatic optima. The age data indicate the glacier's present unprecedented level of retreat since the Last Glacial Maximum and the potential of subglacial carbonates as additional proxies to highlight the extraordinary nature of the current global climatic changes.

1 Introduction

Glaciers respond to climatic changes, making them valuable archives with which to study the effects of past and current climatic changes (Benn and Evans, 2010). Their worldwide retreat over the last century provides evidence in support of current global climate change even though the decrease in summer insolation in the Northern Hemisphere favours climate cooling (Solomina et al., 2015; IPCC, 2018). The uniqueness of this trend can only be understood when compared to past retreats and advances of different types of glaciers throughout the entire Holocene.

Figure 1 Mount Triglav, the Triglav Glacier and the sampling location of subglacial carbonates relative to the years 1932 (a; arrows) and 2017 (b; dots). Close-up of the exposure of subglacial carbonate (c); arrows point to fluted subglacial carbonate, attached to the lee side of the Upper Triassic limestone. Note the colour change between the bedrock surfaces: greys (longer exposure) versus beige and brown on the recently exposed surfaces. This has often been explained by microbial activity (Dias et al., 2018). Photos courtesy of (a) Janko Skerlep (© ZRC SAZU Anton Melik Geographical Institute archive), (b) Miha Pavšek and (c) Matej Lipar.

The data that inform a glacier's dynamics in the past can be provided by subglacial carbonate deposits (Hallet, 1976; Sharp et al., 1990). These are thin carbonate crusts formed between the ice and bedrock by regelation at the glacier base on the lee side of a bedrock protuberance (Hallet, 1976; Lemmens et al., 1982; Souchez and Lemmens, 1985) and can provide information on chemical and physical processes present at the time of their formation (Hallet, 1976).

Unlike other glacial deposits (e.g. moraines), subglacial carbonates may be eroded in a matter of decades (Ford et al., 1970), suggesting their exposure can be recent. They have been reported from deglaciated areas of northern North America (Ford et al., 1970; Hallet, 1976; Refsnider et al., 2012), northern Europe and the European Alps (Lemmens et al., 1982; Souchez and Lemmens, 1985; Sharp et al., 1990; Lacelle, 2007; Gabrovec et al., 2014; Colucci, 2016; Thomazo et al., 2017), Tibet (Risheng et al., 2003), New Guinea (Peterson and Moresby, 1979), and Antarctica (Aharon, 1988; Frisia et al., 2017). The uranium–thorium (U–Th) method remains the main dating technique for these carbonates, whilst the ¹⁴C technique is invalidated by modern carbon contamination (Aharon, 1988) and a possible dead carbon effect.

The deposits recently exposed by the rapid retreat of the Triglav Glacier (Fig. 1) offer a unique opportunity to gain additional knowledge of this glacier's behaviour in the past. The Triglav Glacier extends from the northeastern side of Mount Triglav (2864 m a.s.l.) in the Julian Alps (Figs. 2, S2), Slovenia's highest mountain, consisting of Upper Triassic limestone and dolostone (Ramovš, 2000; Pleničar et al., 2009). This region's montane climate is characterised by precipitation from moisture-bearing air masses from Mediterranean cyclones, typically in autumn and late spring, and by frequent freeze–thaw cycles (Komac et al., 2020).

At present, the Triglav Glacier is one of only two remaining ice masses in Slovenia since the last extensive Pleistocene glaciation (Bavec and Verbič, 2011; Ferk et al., 2017; Triglav-Čekada et al., 2020). The retreat of the glacier has exposed a glaciokarst environment (Figs. 2, S2) comprising a range of erosional (shafts (i.e. vertical caves), karren, polished surfaces and roches moutonnées) and depositional (moraines, boulders, till and carbonate deposits) features (Colucci, 2016; Tičar et al., 2018; Tóth and Veress, 2019). The known extent and behaviour of the Triglav Glacier spans from the Little Ice Age (LIA), the cool-climate anomaly between the Late Middle Ages and the mid-19th century (Grove, 2004; Nussbaumer et al., 2011), to the present (Colucci, 2016; Colucci and Žebre, 2016) and is based on geomorphological remnants, historical records and systematic monitoring since 1946 (Gabrovec et al., 2014). Over the last 100 years, the glacier retreated from ca. 46 ha (extending between 2280 and 2600 m a.s.l.) to ca. 0.5 ha (between 2439 and 2501 m a.s.l.) (Fig. S1), with a downwasting rate of around 0.6 m/yr (1952–2016) (Triglav-Čekada and Zorn, 2020). The glacier has thus evolved from the plateau type through glacieret to the present ice-patch type.

Figure 2 Locality map of the Triglav Glacier, the sampling sites discussed in the text and general geomorphology.

Glacierets are defined as a type of miniature (typically less than 0.25 km²) glacier or ice mass of any shape persisting for at least 2 consecutive years (Cogley et al., 2011; Kumar, 2011). Serrano et al. (2011) discriminate them from ice patches, arguing that glacierets are “the product of larger ancient glaciers, still showing motion or ice deformation, although both very low. They have a glacial origin, glacial ice and are never generated by new snow accumulation”; ice patches on the other hand are “ice bodies without movement by flow or internal motion”. Despite their small size, glacierets occupy a significant ice volume fraction at regional scales (Bahr and Radić, 2012) and can therefore be considered an important target for palaeoclimate studies. Accordingly, many present-day glacierets are closely monitored and studied (Gądek and Kotyrba, 2003; Grunewald and Scheithauer, 2010; Gabrovec et al., 2014; Colucci and Žebre, 2016), but the peculiarity of current global climate change requires more evidence from different proxies and from the past when current ice patches and glacierets were still glaciers.

Here, we present preliminary geochemical and petrological data of subglacial carbonate deposits recently exposed by the retreat of the Triglav Glacier. The aim is to highlight the occurrence of deposits in terms of their possible preservation since the Last Glacial Maximum (LGM), discuss the complexity of the deposit, and validate the results in the context of the present climate regime of rising temperatures and the global retreat of glaciers.

2 Methods

The extent of the Triglav Glacier has been measured annually since 1946 and systematically photographed since 1976 (Meze, 1955; Verbič and Gabrovec, 2002; Triglav Čekada and Gabrovec, 2008; Triglav-Čekada et al., 2011; Triglav-Čekada and Gabrovec, 2013; Gabrovec et al., 2014; Del Gobbo et al., 2016), using a panoramic non-metric Horizont camera. The photos were transformed from a panoramic to a central projection in order to allow the calculation of the area and estimation of the volume (Triglav-Čekada et al., 2011; Triglav-Čekada and Gabrovec, 2013; Triglav-Čekada and Zorn, 2020). The early measurement technique was by measuring tape and compass, which enabled measurement of the glacier's retreat from coloured marks on the rocks around the glacier (Meze, 1955). Accurate and continuous geodetic measurements began in the 1990s: standard geodesic tachymetric measurements, GPS measurements and lidar (Triglav Čekada and Gabrovec, 2008; Gabrovec et al., 2014). In addition, several extensive field campaigns were conducted throughout 2018 with the focus on the central part of the Triglav Plateau at the side of the present and former glacier (Figs. 1, 2). During this time, five subglacial carbonate samples (Fig. 3) were collected at multiple localities (Figs. 1, 2, S2) for laboratory analyses. These were collected 50 to 100 m from the current ice patch. A hand drill was used to obtain powdered samples for geochemical analysis; sampling was restricted to areas consisting of the densest and thickest sections. Prior to and after each drilling, the surface was cleaned with deionised water and dried under clean air. Drilling was performed with a 0.6 mm drill bit in a clean environment. A stainless-steel spatula and weighing paper were used to harvest the drilled powder and transfer it into clean 0.5 mL sterile vials.

Figure 3 Subglacial carbonate samples used for this study. Upper and lower sides refer to the surface (facing open air before being collected) and underside of the carbonate samples, respectively. The depths were measured relative to the drilling starting point of each sample as they are shown in this figure (lower-side perspective for sample T.03 and upper-side perspective for all the other samples); only the T.03 sample had precipitates thick enough to allow us to obtain several dates at three different depths.

Five thin sections (30–50 µm) were examined using an Olympus BX51 polarising microscope equipped with an Olympus SC50 digital camera. Due to the fragility of the samples, they were embedded in epoxy resin under vacuum before being cut and polished.

The mineralogy of 12 (sub)samples was determined by X-ray diffraction (XRD) using a Bruker D2 PHASER diffractometer equipped with an energy dispersive LYNXEYE XE-T detector, located at the ZRC SAZU Karst Research Institute, Slovenia. Powdered samples were scanned from 5 to 70^∘ at 2θ at a 0.02^∘2 θ/0.57 s scan speed. The diffractograms were interpreted using EVA software by Bruker (DIFFRACplus 2006 version).

Five (sub)samples were dated by the U–Th method at the University of Queensland, Australia. To ensure that samples were suitable for U–Th analysis, they were first measured by inductively coupled plasma mass spectrometry (ICP-MS) for their trace element concentrations. U–Th age dating was carried out using a Nu plasma multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS) in the Radiogenic Isotope Facility (RIF) at the School of Earth and Environmental Sciences, The University of Queensland. Ages were corrected for non-radiogenic ²³⁰Th incorporated at the time of deposition assuming an initial ²³⁰Th / ²³²Th activity ratio of 0.825±50 % (the bulk-Earth value, which is the most commonly used for initial and detrital ²³⁰Th corrections). Age errors are reported at the 2σ level. Full details of the method are provided in the Supplement – Table S1.

The stable-isotope composition (δ¹³C and δ¹⁸O) of five (sub)samples was measured at the School of Geography, The University of Melbourne, Australia. Analyses were performed on CO₂ produced by the reaction of the sample with 100 % H₃PO₄ at 70 ^∘C using continuous-flow isotope-ratio mass spectrometry, following the method previously described in Drysdale et al. (2009) and employing an Analytical Precision AP2003 instrument. Results are reported using the standard δ notation (per mille ‰) relative to the Vienna Pee Dee Belemnite (V-PDB) scale. The uncertainty based on a working standard of Carrara Marble (NEW1) is 0.05‰ for δ¹³C and 0.07‰ for δ¹⁸O.

3 Results

The subglacial carbonate deposits occur on the lee sides of small protuberances on a bare polished and striated limestone bedrock surface near the Triglav Glacier. They are most abundant on the bedrock recently uncovered by the retreating ice, and their occurrence rapidly decreases with the distance from the edge of the present glacieret margin. The fluted and furrowed crust-like deposits are brownish, greyish or yellowish in colour. The deposits do not exceed 0.5 cm in thickness and are occasionally internally laminated.

3.1 Mineralogical and petrographic data

X-ray diffraction (XRD) analysis shows that the carbonate deposits mostly consist of calcite and of mixtures of calcite with small amounts of aragonite. XRD also confirms the calcite composition of the host rock (Fig. S5). The petrographic study has allowed the identification of different fabrics (Fig. 4). Due to the similarities of subglacial carbonate textures to those of speleothem deposits, we have used, where possible, the formal terminology of Frisia and Borsato (2010):

Figure 4 Petrography of the carbonate crusts. (a) Primary columnar calcite crystals growing perpendicularly to the substrate on the steeper side of the bedrock irregularity; (b) alternation of aragonite (fibrous) and calcite crystals forming layered textures. The base of the aragonite fans of crystals nucleates in dark micritic aggregates. Both images were taken under plane-polarised light.

Primary calcite fabrics are composed of transparent crystals with uniform extinction. The first crystals to form directly over the bedrock are short columnar (length-to-width (L∕W) ratios $<\mathrm{6}:\mathrm{1}$) crystals from 50 to 200 µm long (Figs. S6 and S7). Columnar (L∕W ratios ∼ 6:1) and elongated columnar (L∕W ratios $>\mathrm{6}:\mathrm{1}$) crystals up to 2 mm long and 0.5 mm wide constitute the most abundant fabric. Calcite crystals grow perpendicularly to the substrate on the steeper, nearly vertical areas of the bedrock (Fig. 4a), while in the less steep, nearly horizontal areas, crystals grow inclined, oriented downslope and presumably in the sliding direction of the forming glacier (Fig. S6). In some areas, the younger crystals cross-cut the main direction of the crystal growth of the previous layer, resulting in crystal boundaries resembling unconformities (Fig. S6).

Primary aragonite fabrics consist of acicular crystals (L∕W ratio $\gg \mathrm{6}:\mathrm{1}$) generally growing outwards from a common point in the shape of fans. In some areas of the crusts, these fans are aligned in bands interlayered with very dark, dense micrite and transparent equidimensional calcite crystals, forming layered textures (Figs. 4b, S6d).

Aragonite-to-calcite replacement fabrics are characterised by calcite crystals of variable size and patchy extinction patterns. They contain abundant fibrous inclusions interpreted as aragonite relicts that are either aligned in layers or unevenly distributed throughout, which results in very irregular textures (Fig. S6d, e and f). Micrite and microsparite are often associated with aragonite relicts.

3.2 Geochemical data and ages

Stable carbon and oxygen isotope ratios (δ¹³C and δ¹⁸O) of the subglacial carbonate yielded average values of 1.35‰±0.05‰ for δ¹³C and $-\mathrm{4.32}\mathrm{‰}\pm \mathrm{0.07}\mathrm{‰}$ for δ¹⁸O (relative to the V-PDB; Fig. 3).

Whilst the ages of subglacial carbonates were expected to belong to the LIA, the U–Th ages yielded considerably older ages: 23.62 kyr ± 0.78 kyr, 18.45 kyr ± 0.70 kyr and 12.72 kyr ± 0.28 kyr, respectively (Table S1). The results indicate that these subglacial carbonate dates fall within the LGM and the Younger Dryas (YD). In addition, two U–Th ages corresponding to samples drilled in more surficial calcite layers of the thickest obtained sample yielded ages of 3.85 kyr ± 0.09 kyr and 1.96 kyr ± 0.04 kyr, indicating the presence of a glacier of sufficient thickness to also cause regelation during these periods. For details of where the samples were obtained from, the reader is referred to Fig. 3.

4 Subglacial carbonate deposits

The fluted and furrowed morphology of carbonate deposition parallel to the apparent former ice-flow direction, some of the crystal textures, the location of the samples, and the age data imply a subglacial origin of the carbonate crusts. The presence of ice-flow-oriented calcite crystals suggests that precipitation was strongly influenced by the mechanical force of the ice movement. Similar fabrics occur in different types of tectonic veins, for instance, slickensides, where crystals grow obliquely to the sheets in shear veins, indicating the general shear direction (Bons et al., 2012). Previous studies on subglacial carbonates have considered the orientation of calcite parallel to the ice flow to be evidence of its growth in a thin water film confined by sliding regelation ice (Hallet, 1976) and deformation structures such as folds and fractures to be the result of glacially imposed stress (Sharp et al., 1990).

The variability in fabrics of the studied crusts and, especially, the presence of aragonite coexisting with calcite indicate spatial and/or temporal variability in local subglacial water chemistry. The most important parameters controlling the precipitation of aragonite versus calcite appear to be the CaCO₃ saturation state, the Mg∕Ca ratio in the waters (De Choudens-Sánchez and González, 2009) and/or rapid CO₂ degassing rates (Fernández-Díaz et al., 1996; Jones, 2017). In freshwater systems like spring deposits (Jones, 2017), and especially in speleothems, Mg∕Ca ratios seem to be the main factor controlling aragonite vs. calcite precipitation (Frisia et al., 2002; Wassenburg et al., 2012; Rossi and Lozano, 2016). The studied deposits grow over carbonate bedrock with prevailing limestone and some dolostone (Jurkovšek, 1987; Ramovš, 2000), so high Mg∕Ca ratios in the water may be partially responsible for the precipitation of aragonite. Similarly, precipitation of aragonite in subglacial deposits at Vestfold Hills, Antarctica, seems to have been favoured by the presence of Mg²⁺ in the waters, mobilised from the pyroxenes of the gneiss bedrock (Aharon, 1988).

The U–Th ages of LGM and YD are in good accordance with the glacier's history and when it was expected to be thick enough to cause regelation. However, aragonite-to-calcite replacement fabrics raise the difficult question of whether all subglacial carbonate crystals were primarily aragonite and were subjected to complete recrystallisation to calcite without leaving traces of replacement fabrics, which may provide a false identification of calcite fabrics as primary. This, in turn, may lead to inaccurate interpretation of the U–Th ages (Bajo et al., 2016). Previous studies on corals and speleothems have shown that the diagenetic transformation of aragonite into calcite can affect the accuracy of U–Th dating (Ortega et al., 2005; Scholz and Hoffmann, 2008; Lachniet et al., 2012; Bajo et al., 2016; Martín-García et al., 2019) due to uranium loss, which usually leads to older-than-true U–Th ages. However, the extent of uncertainty can be highly variable depending on several factors such as the initial amount of aragonite and the timing of its transformation to calcite (Bajo et al., 2016), the possible additional redistribution of Th (Ortega et al., 2005; Lachniet et al., 2012), or the degree of opening of the system during diagenesis (Ortega et al., 2005; Martín-García et al., 2019) produced by the chemical exchange with younger subglacial meltwater.

Notably, the U concentration (in ppm; Table S1) in the youngest sample (2 kyr; T.03_b1) within this study is 1.77 ppm, whereas it is around 0.41 and 0.46 ppm in two of the old samples, T.01_a1 and T.03_a1, which date to the LGM and YD, respectively. This could indicate either U loss during aragonite-to-calcite recrystallisation and consequently provide older-than-true ages or, to the contrary, that these two samples represent primary calcite with non-preferential incorporation of U in calcite with respect to aragonite and thus proving to be the most reliable. Assuming the first possibility, the oldest sample (24 kyr; T.05) has a relatively high U concentration (1.33 ppm), which would indicate less loss of U but still provide data indicating that subglacial carbonates may date back to the last glacial period. Assuming the second possibility, the ages of T.03_a1 and T.03_a2 in correct stratigraphic order would strengthen the reliability of the results. Aharon (1988) used the U–Th method for dating two samples with a combination of aragonite and calcite fabrics with high but also variable U content (23.2 and 41.4 ppm); however, the indistinguishable dates compared to those from pure aragonite led to the conclusion that samples provide reliable age estimates. Nevertheless, due to the external Th incorporated into the samples, Aharon regarded the dates as “maximum ages”, which is also the appropriate approach with the dates discussed in this study. Similar U–Th ages (19–21 kyr) were reported also from northern Canada (Refsnider et al., 2012); however, the carbonate fabrics were not described. The U–Th dating from Antarctica by Frisia et al. (2017) and radiocarbon dating from the French Alps by Thomazo et al. (2017) were performed on calcites. More studies should therefore be carried out in this direction to positively confirm the LGM and YD ages of carbonates and especially to gain the high-resolution dates needed to construct the whole timeline of subglacial carbonate precipitation. Nevertheless, since three of our ages fall in the period of 12–24 ka, we will proceed with the discussion of their susceptibility to weathering and the implication of this for the glacier's existence.

The cold Alpine environment during the glacial period with low biological respiration rates could be indicated in the relatively high δ¹³C signal, also reported by others (Lemmens et al., 1982; Fairchild and Spiro, 1990; Lyons et al., 2020), but the (re)freezing of the subglacial water causes supersaturation with respect to carbonate, and the non-equilibrium conditions produced by this process can affect the stable isotopic composition of the subglacial carbonate, usually leading to isotopic enrichment in the carbonate minerals (Clark and Lauriol, 1992; Courty et al., 1994; Lacelle, 2007). In any case, the climate during the precipitation of subglacial carbonate had to be “warm” enough to produce subglacial water at the glacier bed. Kuhlemann et al. (2008) suggested a lowering of summer temperature of about 9–11 ^∘C at the LGM peak, meaning that the summer temperatures at Triglav Glacier should have been around −4.7 to −6.7 ^∘C (based on the present summer temperatures at around +6 ^∘C (Slovenian Environment Agency, 2020a) and considering the recent warming of up to +2.03 ^∘C in the area; Colucci and Guglielmin, 2015; Hrvatin and Zorn, 2020), which would have provided conditions for the cold-based glacier and the absence of subglacial water flow. Whilst small-scale detailed palaeoclimate conditions in the southeastern Alps are still uncertain, the current U–Th ages of subglacial carbonates within this research fit with the temporary Garda Glacier (Italy) withdrawal phases reported by Monegato et al. (2017): the 23.6 kyr age of subglacial carbonate relates to the withdrawal phase between 23.9 and 23.0 ka, and 18.5 kyr could relate to the glacier retreat from around 19.7 to 18.6 ka or even to the final collapse of the glacier around 17.7 to 17.3 ka. Both time periods could also relate to commencement of subglacial water flow at the Triglav Glacier and consequently the precipitation of subglacial carbonates. In addition, the 12.7 kyr age marks the early phase of YD cooling (12.9–11.5 kyr ago) (Renssen and Isarin, 1997; Alley, 2000; Broecker et al., 2010).

The δ¹⁸O_PDB of subglacial carbonate can be transformed into δ¹⁸O_SMOW using the equation ${\mathit{\delta }}_{\mathrm{SMOW}}=\mathrm{1.03037}{\mathit{\delta }}_{\mathrm{PDB}}+\mathrm{30.37}$ (Faure, 1977). The mean values of subglacial carbonate δ¹⁸O_SMOW would therefore be 25.92‰, ranging from 27.44‰ to 24.75‰. Using the fractionation factor of 1.0347 for calcite and water at 0 ^∘C (Clayton and Jones, 1968) and the δ_SMOW values of subglacial carbonate, we obtain δ¹⁸O_SMOW values ranging from −7.02‰ to −9.62‰. On the other hand, if assuming calcite crystals could all have originated primarily as aragonite crystals, we can use the fractionation factor of 1.0349 for aragonite and water at 0 ^∘C (Kim et al., 2007) and obtain relatively similar δ¹⁸O_SMOW values ranging from −7.21‰ to −9.81‰. This relates to the average Triglav Glacier ice meltwater values (δ¹⁸O${}_{\mathrm{SMOW}}=-\mathrm{9.3}$‰) measured in summer 2018 (Carey et al., 2020). However, the isotopes are slightly heavier when compared to glacier ice samples, which range between −10.0‰ and −12.7‰. The reason could be that water in equilibrium with the growing ice is progressively impoverished in heavy isotopes (Jouzel and Souchez, 1982) or simply that the remnants of residual ice have a different isotopic ratio than the basal ice at the time of carbonate precipitation, which itself could represent a large range in values (Lemmens et al., 1982). In addition, δ¹⁸O differences of a few per mille in the carbonate precipitate can also arise due to variations in subglacial hydrology shifting from closed to open geochemical systems (Hanshaw and Hallet, 1978). However, a geochemical study of the present Triglav ice area (Carey et al., 2020) shows that ice samples within the cave of Triglavsko Brezno (see Fig. 1 for location) are much lighter in deuterium than those from the Triglav Glacier, which suggests that they were deposited during colder times, indicating they are remnants of older glacier ice.

Table 1 Calculation of subglacial carbonate existence under a meteoric environment based on the carbonate's thickness, years of exposure and denudation rate.

5 Implications for a continuously existing glacier during the Holocene

The LGM and YD ages of the Triglav subglacial carbonates provide the first physical evidence that the Triglav Glacier persisted through the Holocene to the present day. Glacierets of southern Europe, including the Triglav Glacier, have generally been viewed as relicts of the LIA, suggesting a discontinuous presence following the Holocene climatic optimum (HCO) (Grunewald and Scheithauer, 2010), a period of high insolation and generally warmer climate between 11 000 and 5000 years ago (Renssen et al., 2009; Solomina et al., 2015).

Being prone to fast weathering (Ford et al., 1970), subglacial carbonate deposits are generally found only on recently deglaciated areas (Ford et al., 1970; Sharpe and Shaw, 1989). The Triglav Glacier area is no exception. The fact that the deposits had not been reported in the literature (Meze, 1955; Šifrer, 1963, 1976, 1987) until the year 2005 (Hrvatin et al., 2005; Gabrovec et al., 2014) reinforces the likelihood of their recent exposure. Whilst Refsnider et al. (2012) discussed the preservation of the subglacial carbonates by a very cold and dry arctic climate, this is unlikely to be the case for the plateau of the Triglav Glacier, either in the last century or during any possible period of absence of an ice body, because of the large mean annual precipitation at the site (2038 mm; recorded at Kredarica hut, 2515 m a.s.l.) (Slovenian Environment Agency, 2020a) (see Fig. 2 for location).

This is supported by theoretical numerical data for the chemical denudation of the subglacial carbonate, which are based on denudation rates of limestones (Table 1). Chemical denudation rates on carbonate rocks can vary from ca. 0.009 to 0.14 mm/year (Gabrovšek, 2009). Taking the low and high extreme values for, e.g., 6 kyr during the HCO, surface lowering would be between 54 and 840 mm, so the exposed 5 mm thick subglacial carbonate would have been denuded by this time.

In addition, carbonate surfaces in periglacial areas are exposed not only to chemical weathering but also to intensive frost weathering, promoting the physical disintegration of minerals (Matsuoka and Murton, 2008). Therefore, had the subglacial carbonate been exposed in the past, it would be expected that it would have been eroded by dissolution or frost weathering. This indicates that subglacial carbonate was constantly covered with glacier ice and, as Steinemann et al. (2020) have demonstrated, that glacial abrasion on carbonate bedrock is minimal compared to the abrasion on crystalline bedrock. This would promote the preservation of subglacial carbonates on limestone plateaus in the Alps. Moreover, since the subglacial carbonates form in lee positions of bedrock protuberances, it is likely that abrasion would vanish as abrading rock fragments diverge from the bed at sites of subglacial precipitation due to regelation ice growth.

Organic matter (charcoal or wood) from a non-vegetated, scree-covered moraine ca. 300 m below the main ice patch of the Triglav Glacier (as it stood in the year 2006) was analysed in the radiocarbon (¹⁴C) laboratory in Erlangen, Germany. The ¹⁴C result yielded a 5604–5446 kyr BP age, which provides additional evidence of pre-LIA and post-LGM and post-YD ice cover (unpublished analysis by Karsten Grunewald and his team). Similarly, the two younger U–Th ages obtained within this study (3.85 and 1.96 kyr) also provide evidence of a pre-LIA ice cover.

6 Glacier variations and palaeoclimatic implications

The complex global pattern of Holocene glacier fluctuations indicates the influence of multiple climatic mechanisms and that individual glaciers may not respond uniformly to a particular set of climate forcings (Kirkbride and Winkler, 2012; Solomina et al., 2015). In addition, glaciers are also influenced by topographic conditions (DeBeer and Sharp, 2017), as well as their size and flow dynamics (Sugden and John, 1976; Nussbaumer et al., 2011). The Alps has experienced several glacial advances and retreats during the Holocene (Nussbaumer et al., 2011), with reports that some glaciers were even smaller than today or absent (Leemann and Niessen, 1994; Hormes et al., 2001; Ivy-Ochs et al., 2009; Solomina et al., 2015). On the other hand, certain regions show evidence of an unprecedented modern retreat of glaciers beyond their previous Holocene minima (Koerner and Fisher, 2002; Antoniades et al., 2011; Miller et al., 2013).

The current climate near the Triglav Glacier is characterised by rising temperatures in the ablation season from May to October (Fig. S3) and a descending trend in the highest seasonal snow elevation (Fig. S4). If subglacial deposits indicate its ongoing existence throughout the Holocene, the recent retreat of the Triglav Glacier suggests the regional climate during previous Holocene optima was cool enough to sustain the glacier, unless an additional (presently unknown) forcing component prevailed which was not important in the past. Exceptionally rapid 21st-century melting, for example, has also been reported from Barnes Ice Cap in northern Canada (Gilbert et al., 2017) and in the Alps, where the exceptional retreat has been attributed to deposition of industrially sourced black carbon (Painter et al., 2013), and is also evident from a large number of archaeological and palaeontological finds due to recent melt-out (Mol et al., 2001; Basilyan et al., 2011), including the famous “Iceman” Ötzi in the Tyrolean Alps (Baroni and Orombelli, 1996; Solomina et al., 2015).

Local physiographic influences can insulate small glaciers from the warming effects of the regional and global climate (Grunewald and Scheithauer, 2010); it has been shown that some small glaciers are less sensitive to climate fluctuations than previously thought (Colucci and Žebre, 2016). The natural resilience of the Triglav Glacier is due to its relatively high elevation, its bright limestone substrate with higher albedo, the vertical water drainage through karstified rocks and the consequent lack of subglacial lakes as heat collectors, and avalanche feeding (Grunewald and Scheithauer, 2010; Gabrovec et al., 2014). Assuming these factors were relatively constant throughout the Holocene, the possible unprecedented retreat may highlight the consequences of direct anthropogenic forcing (Solomina et al., 2015). On the other hand, further comparative research on other small glaciers is needed to test palaeoclimatic data due to the geomorphological peculiarity of glacier regions. For example, despite the apparent resilience of the Triglav Glacier until recently, its retreat has been more evident than those of the glacierets and ice patches of Kanin (Italy), Montasio West (Italy) and Skuta (Slovenia) (all in the southeastern Alps), which are all lower in elevation than the Triglav Glacier but with large differences in mean annual precipitation (e.g. in the Triglav area water equivalent precipitation (2038 mm) is 62 % of that in Kanin area (3335 mm)); potential annual solar radiation also plays a major role in controlling for differences in glaciers' dynamics (Colucci, 2016).

7 Conclusions

Three U–Th ages of subglacial carbonate exposed by the retreating Triglav Glacier fall within the LGM and YD, providing the first direct evidence of the Triglav Glacier at that time. The high erodibility of these deposits, once exposed, strongly suggests continuous glacier cover since their deposition throughout all but the most recent part of the Holocene, including the HCO. This defines the subglacial carbonates as a complex, but an important, palaeoenvironmental proxy and as a research subject for further analyses with great potential that may further highlight the extraordinary nature of the current global warming in the context of the Holocene.

Several separate lines of additional evidence could strengthen this assertion:

• A

  The identification and geochemical analyses of the primary aragonite would contribute data on primary U concentration in samples (where there is evidence of primary aragonite), which can be used to detect U loss. In situ U-series dating by laser ablation would represent the best approach as the samples (and laminae) are small and the aim is to target different crystal fabrics, but this is challenging due to low U concentrations.

• B

  Use the ³⁶Cl nuclide dating method on the limestone hosting the subglacial carbonate deposits to extract the exposure time of the limestone surfaces hosting the carbonate crusts (i.e. no glacier cover) and provide additional data as to whether or not those areas have experienced multiple episodes of exposure during the HCO. The production rate at 2500 m a.s.l. (the mean altitude of the upper part of the present Triglav Glacier) is relatively high, which is an advantage for dating young exposures. In addition, calcite and limestone comprise one of the best mineral and rock systems for ³⁶Cl because they often have low ³⁵Cl abundances, so there is less uncertainty regarding contributions from factors that affect the thermal neutron production of ³⁶Cl, such as water or snow shielding (Fabel and Harbor, 1999; Marrero et al., 2016).

• C

  The analyses of material on the small remaining ice masses and, if possible, ice cores can highlight the existence of a present-day input to accelerate melting. For example, black carbon and other mineral dust particles in glaciated regions, which accumulate on the ice, can accelerate the melt of glaciers by reducing albedo (Ramanathan and Carmichael, 2008; Ming et al., 2013; Painter et al., 2013; Gabbi et al., 2015).

• D

  Numerical data on frost weathering can be used. There is a relatively extensive literature concerning the chemical denudation and glacial erosion of the limestone but scarce numerical data for frost weathering and no direct measurements concerning subglacial carbonates. This can be experimentally studied by direct freezing–thawing experiments with controlled humidity and additional control for the type of porosity and pressure gradient (Ducman et al., 2011).