RADIOCARBON CONSTRAINTS ON PERIODS OF POSITIVE CAVE ICE MASS BALANCE DURING THE LAST MILLENNIUM, JULIAN ALPS (NW SLOVENIA)

Abstract

Caves containing perennial ice deposits make up a little-known, but emerging part of the cryosphere under increasing scrutiny from the scientific community. M-17, a sag-type ice cave opening at 1879 m asl in the Tolminski Migovec massif of the Julian Alps (NW Slovenia) contains a perennial underground ice deposit whose paleoclimate sensitivity is poorly understood and whose longevity under current climate change is at risk. The past mass balance of this cave is constrained using wood macro-remains embedded in ice. Accelerator mass spectrometry radiocarbon dating of 18 wood samples embedded in ice provides the largest currently available dataset for a subterranean ice deposit in the southern European Alps. The reconstructed chronostratigraphy reveals three main phases of likely positive ice balance around 900–1100 AD, 1200–1300 AD, and 1700–1800 AD, as well as a period of negative mass balance around 1300–1400 AD. The onset of cave glaciation is deemed to have occurred no later than about 900 AD, with evidence of overall positive ice mass balance during multi-decadal periods characterized by cooler-than-average summers and wetter-than-average springs. Conversely, negative mass balance is recorded during a period warmer-than-average summers and dry springs. The cave has experienced ice mass loss since its discovery in the 1980s.

INTRODUCTION

The observation and systematic investigation of ice cave phenomena in the European Alps largely find their origins in the 18th century (Meyer et al. 2017 and references therein). Variations of seasonal ice level, relationships between cave and outside air temperatures, and ice movement were among the early observations in caves of the Eastern and Western Alps as well as the Jura Mountains (Fugger 1888; Balch 1900). One of the first detailed micro-meteorological studies of an ice cave was performed in Dachstein-Rieseneishöhle (central Austria) in the late 1920s (Saar 1956). With the spread of speleology in the 20th century, the reports of Alpine ice cave phenomena multiplied, bringing ever more potential study sites under scrutiny. To date, many thousands of ice caves are known worldwide and have been surveyed (e.g., Perşoiu and Lauritzen 2018). Monitoring of ice-level fluctuations (Perşoiu et al. 2012), radiocarbon age determination of organic matter included in ice (Luetscher et al. 2007), plant remains (Leunda et al. 2019), trace elements (Munroe et al. 2018) and stable isotopes in ice (Yonge and Macdonald 1999; Bădăluţă et al. 2020) are methods used to address questions related to the long-term dynamics of these poorly understood underground ice deposits and their link to climate.

Ice-bearing caves are commonly classified according to (1) their geometry and attendant ventilation patterns, and (2) the type and mechanism of ice deposition (Luetscher and Jeannin 2004). The present study bears on a specific type of ice caves with a steeply descending morphology. At the bottom of these so-called sag-type caves, firn-derived deposits form by a combination of snow diagenesis, water percolation and freezing (Luetscher and Jeannin 2004). In addition, woody macro remains (twigs, branches, needles, pinecones, etc.) may be wind-borne or avalanche into the cave and thereafter become embedded within underground firn or ice. The process of ice accumulation results in stratified underground ice bodies of limited volume but predictable geometry, from which preserved past climate signals may be recovered (e.g., Perşoiu et al. 2017; Sancho et al. 2018).

The key to unlocking the potential record archived in underground ice is to place it in a robust chronological framework. In ice caves, this task is achieved by various methods, including annual layer counting where practicable (e.g., Rygielski et al. 1995) or radiogenic ²¹⁰Pb age determinations for ice deposits no older than about 150 years (Appleby 2008), as demonstrated by Luetscher et al. (2007). Radiocarbon age determination of wood inclusions is the most frequently applied dating method for underground deposits (Kern 2018, and references therein). The oldest wood macro remains hitherto identified and dated in ice caves of the European Alps place the earliest onset of (preserved) deposition at around cal 3366–3030 BC (2σ, Achleitner 1995; Spötl et al. 2018). Despite the complex stratigraphic relationships often exhibited by underground deposits due to the constrained geometry of the accumulation sites and the various ice accumulation processes such as snow ingress and water freezing (Malyudov 2018), fragments preserved in individually accumulated units date the time at which that unit was deposited. The combination of radiometric dating of inclusions and the evaluation of stratigraphic relationship between ice units yields unique chronostratigraphic histories for each ice cave (Kern 2018).

In the last two decades, several detailed chronostratigraphies of Alpine ice caves have been published. A detailed study was conducted at St-Livres ice cave (Jura Mountains), where wood inclusions in the basal ice were radiocarbon-dated back to about 750 AD (Stoffel et al. 2009). In this study, a combination of radiocarbon dating and dendrochronology was used to reconstruct ice deposition gaps, whose timing and duration was subsequently cross-validated with paleoclimate indices such as the North Atlantic Oscillation (NAO). In the Alps, Hundsalm Eis- und Tropfsteinhöhle (hereafter HETH) a cave situated in Tyrol, western Austria, was the first target of an extensive radiocarbon sampling campaign. The 19-sample dataset from this cave was interpreted as evidence of substantially positive mass balance around 600 AD and to a greater extent from about 1300 AD onwards (Spötl et al. 2014), suggesting a positive mass balance during the so-called “Late Antique” Little Ice Age (LALIA) and the subsequent Little Ice Age (LIA) (Büntgen et al. 2011). At the HETH site, however, the complex stratigraphy and apparent out-of-order ¹⁴C ages precluded both the systematic sampling of a proxy and the production of an age-depth model, which could otherwise be used to derive quantitative mass balance values. Spötl et al. (2014) summed the probability density functions of the calibrated ¹⁴C ages and interpreted maxima (minima) of the resulting curve as alternatively favorable (unfavorable) periods for positive ice mass balance, highlighting for instance the remarkable correlation between glacier advances in the Northern Alps with the summed probability curve. More recently, the detailed sampling of Alpine ice-bearing caves was carried out at sites where Roman Age ice deposits were found to have been preserved (Racine et al. 2020).

Despite the great number of sag-type ice caves in the Southern Alps and the Dinarides, only relatively few radiocarbon chronostratigraphies have been published (Paar et al. 2013; Staut et al. 2016). Most recently, Kern et al. (2018) sampled woody macro remains in three caves of the Velebit Mountains (Croatia) and reported periods of positive mass balance at around 1500 BC, 1000 AD, and 1500–1800 AD for old deposits in the Vukušić, Kugina, and Ledena caves, respectively. In this study, wood fragments as well as individual trunks were sampled primarily from layers rich in organic matter, interpreted as unfavorable for positive ice mass balance. Substantial ice loss in recent years (Kern and Perşoiu 2013), accelerated by extreme precipitation events (Perşoiu et al. 2021) which are expected to increase under climate scenarios RCP4.5 and RCP8.5 (Púčik et al. 2017) underscores the vulnerability of this unique archive.

From a cryospheric standpoint, the presence of underground ice bodies, recorded by historical cave exploration, represents a way to circumvent the lack of large surface glaciers in the Julian Alps of northern Slovenia (Colucci et al. 2016) and provides an archive sensitive for precipitation, recording past intervals of positive ice mass balance. Results from M-17 archive can be compared to and contrasted with the growing corpus of records of past cave glaciation in Europe. The purpose of the present study is thus to provide a new detailed mass balance record based on cave ice in the Julian Alps, through systematic sampling of woody macro remains in a stratified, subterranean ice deposit. This study underscores the need to sample complex stratigraphies to constrain the onset and duration of periods of past positive mass-balance.

Study site

The Tolminski Migovec Massif

The Julian Alps in northwestern Slovenia are dominated by limestone and dolomite ranges with peaks mostly ranging between 1800–2400 m asl (Figure 1), with a maximum of 2864 m asl at Mount Triglav (Perko et al. 2020). The Dachstein Formation, made up of thick, well-bedded carbonates of Triassic age, is ubiquitous in these mountain ranges and forms the backbone of the Tolminski Migovec massif (Buser 1986; Ogorelec and Buser 1996). The Alpine karst developed therein is characterized by several deep (> 1 km) shaft and canyon cave systems (Racine 2019). The surface expression of this karst comprises many large dolines, Schichttreppenkarst (stepped limestone pavements), and snow- or ice-bearing pits locally referred to as kotliči. Such ice caves are a common feature of Alpine karst in general (Maire 1990), and the Migovec plateau in particular (Figures 1 and 2).

Figure 1 (A) and (B) study area localization. (C) Topographic map of the Tolminski Migovec massif and the cave entrances (orange circles). A multi-directional hillshade of the digital terrain model was extracted from the Slovenia LIDAR dataset (ARSO 2021). (1) M-17, (2) M-15 and (3) M-10 (4) Planika Jama, and (5) Bertijeva Jama. The field of view of photograph (D) is highlighted in yellow. (D) View of the karstic plateau of Tolminski Migovec and the Bohinj Ridge in the background, photo: C. Spötl (2020). (Please see electronic version for color figures.)

Figure 2 Simplified vertical profiles of the studied sag-type ice caves in the Migovec area. Depth is plotted relative to the highest entrance of each site. Dashed blue lines and year of observation indicate past snow levels. Adapted from digital maps provided by ICCC and JSPDT.

The Tolminski Migovec massif rises to 1881 m asl at the western (windward) side of the Bohinj Ridge (Figure 1). The study area mostly receives Mediterranean derived westerly and southwesterly moisture (Vičar and Bertalaničnik 2006; Gómez-Hernández et al. 2013). The annual precipitation is 50–100% higher on Tolminski Migovec than at the closest mountain meteorological station of Vogel (WMO ID: 14006, 1535 m asl, Ogrin and Ortar 2007), which is located on the leeward side and rain shadow of the Bohinj Ridge (Ogrin and Kozamernik 2020). According to observations by the Environmental Agency of the Republic of Slovenia, the mean seasonal snow accumulation in the study area was between 200–420 cm per season between 1971 and 2000 (ARSO 2021); during the same period, mean annual precipitation ranged between 2600 mm and 3200 mm a^–1 (ARSO 2021). Over the period 1970–2019, the Tolminski Migovec area in the Slovenian region that experienced (1) the smallest increase in the maximum temperature of the warmest month, (2) the maximum interannual variability of spring and autumn precipitation, and (3) the highest and least variable snow height (Kozjek et al. 2017). At the Vogel meteorological station, precipitation maxima were recorded in spring (AM) and autumn (SON) between 2010–2019 (suppl. Figure 3, Reliable Prognosis 2021), denoting the strong influence of Mediterranean cyclogenesis in autumn.

The Migovec Plateau hosts several caves containing perennial snow and ice, five of which are of particular interest with regards to their geometry and ice content. Each of the caves presented in Figure 2 is characterized by a simple pit geometry, with one or more entrance dolines, which act as snow funnels in winter. These caves open within a narrow elevation range between 1720 m and 1870 m asl and their entrance diameters span at least 5 m. M-15 and Bertijeva Jama are very similar open pits characterized by a single entrance, and a depth of about 50 m (Figure 2). In the case of M-17, M-10, and Planika Jama, the elevation difference between the multiple openings is very small (Figure 2). The base of the snow/firn/ice deposits coincides with the deepest accessible point in the caves, and each ice deposit rests on a cobble to boulder grade angular clastic fill rather than limestone bedrock (Figure 2).

The local timber line (closed forest limit) extends to about 1600 m asl, and the krummholz zone, containing dwarf specimens, extends up to about 1900 m asl. The predominant vegetation at the cave entrance elevation comprises a patchwork of dwarf pine thickets and grasses of the transitional zone between timber and tree lines. Rare examples of Larix decidua saplings are found in the shelter of dwarf pine thickets. Each of the sag-type caves presented in Figure 2 therefore receives both a substantial quantity of solid winter precipitation (ARSO 2021) and organic debris sourced from the surrounding vegetation.

Of these five caves, M-17 contains the largest volume of underground ice deposit in the Tolminski Migovec area (Figure 2), as well as the most promising layered outcrop containing wood fragments, described in greater detail below.

Ice Occurrences in M-17

M-17 cave (coordinates: 46.253580°N,13.760417°E, cadastral number: 5878), opens at an elevation of 1879 m asl. It is an 88-m-deep cave (Figure 3) developed in limestones of the Dachstein Formation (Buser et al. 1982; Ogorelec and Buser 1996), which was first explored in 1983 and later surveyed and partly photographed in 1994 (ICCC and JSPDT 2007). A more recent survey was conducted with a laser distance-meter in 2018–2019, the result of which is presented in Figure 3 and discussed below.

Figure 3 Detailed vertical profile of M-17 cave, adapted from digital maps provided by Imperial College Caving Club (ICCC) and Jamarska Sekcija PD Tolmin (JSPDT). Depth is plotted relative to each of the cave’s main entrance (a/b). The position of the 1994 snow line is indicated. Markers depict the position of data logger. The inset refers to the well-exposed stratigraphic section (see Figure 5). Survey data courtesy of ICCC and JSPDT.

Three shafts lead to the main ice chamber of M-17 (Figure 3). Entrances (a) and (b), separated by a slender rock bridge, form a pit with an average diameter of 5 m. This constitutes the historical access route to the main ice chamber (Figure 3). The connection of this ice chamber to the open pit (c) was first noted in 2018, by crawling through a short tunnel developed in the snow cone. Pits (c) and (d) connect above the main ice chamber and given their average diameters, and the total catchment area of their respective dolines should contribute more snow to the ice body than the historical entrance pit a/b. Entrance pit a/b remains open to air circulation all year long, whereas c/d is only seasonally open; this was first observed between August 2018 and November 2018, as well as on return visits in August 2020. Given sufficient seasonal snow accumulation, the feeder doline c/d may close again (as observed by the author in June 2019), and effectively turn M-17 into a single entrance cave.

In M-17, as in the other Migovec sag-type sites, the input of solid precipitation during winter leads to a substantial height increase at the snow cone entrance, in the order of 2–5 m, depending on the year. However, by late summer, a combination of slightly positive cave air temperature (ca. 0.5–1.5°C) near the cave entrance, relatively warmer rain, and to a limited extent solar radiation promote loss of snow near the entrance. Seasonal snow-height variations reached an amplitude of at least 5 m at entrance a/b, with a smaller amplitude in doline c/d. The a/b snow cone may disappear entirely and reveal a steep talus ofcryoclasts and underlying ice.

Within the main ice chamber, a gently sloping ice floor covering about 10 × 20 m reaches a thickness of ca. 2 m locally. Frost shattering on the limestone bedrock provides an ample of source of cryoclasts, many of which form a talus cone at the base of the historical entrance (Figure 4b) or end up embedded in the ice floor. This layered ice body is attributed to the dripping and regelation of water that has traveled a short distance through the ca. 20 m rock overburden. Indeed, at the location of major fractures in the bedrock, drip water forms a variety of seasonal ice speleothems including stalactites, stalagmites and pillars. One such 3-m ice stalagmite depicted in the 1994 survey (Figure 4), probably represented several years of seasonal growth under an active drip. It has altogether disappeared since, replaced by a hole in the ice floor several meters across. Recent observations (summers 2017 to present) suggest that most seasonal speleothems, which grow during the spring season, degrade entirely by mid-summer. Water percolating through the fissures of the epikarst drips into the main chamber and enhances floor-ice melting. The formation of small 10-cm-deep bédières (surface water channels) advects infiltration and melt water to the ice-free outskirts of the main chamber. This process is significant with respect to the annual ice mass-balance, since it removes ponded water otherwise available for refreezing during the cold season, as noted in the case of Scărişoara ice cave by Perşoiu et al. (2021).

Figure 4 (A) View of the M-17 main chamber and one of the snow cones, looking towards entrance a/b. Some of ice accumulation processes discussed in the text are highlighted. Photo: Richard Anderson (1994) (B) Similar view of the main ice chamber and snow cone with reduced area, revealing a scree slope and layered ice underneath. Note the absence of a (previously) perennial 4-m-high ice stalagmite. Photo: C. Spötl (2020).

The 2018 survey (Figure 3) indicates that the ice body occupies most of the rock-cave volume at depths ranging between 45 and 85 m. The ice body contains mostly small organic inclusions such as leaves, needles, pinecone, and twig fragments. Such inclusions highlight the stratigraphy in well exposed and very steeply dipping (∼80°) strata (suppl. Figure 1). The deep ice pit outcrop (Figure 5) is the most promising section of the cave for reconstructing past mass balance. The deepest point in the cave contains the base of the ice body, which is predominantly made up of coarse (up to several cm) crystalline, anisotropic ice, overlying a scree slope deposit. This deep part is characterized by the absence of stratal architecture, compounded by the lack of any woody macro remains. Using the topographic survey, the 2019 ice volume and available exchange surface were estimated, assuming a cylinder with a height of 38±2 m and a radius of 4±0.5 m. This calculation yields a total volume of 1.9 ± 0.5 × 10³ m³ and a total exchange surface (ice-air and ice-rock interfaces) of 1.1±0.2 × 10³ m².

Figure 5 Ice stratigraphy in the deep ice pit of M-17, (A) detailed sketch of the exposed outcrop (cf. Figure 3), (B) location of the outcrop within the cave, and (C) summary of the stratigraphy recorded along the ice and sampling transect for woody macro remains. Filled circles denote sampled macro-remains. (D) Rescaled photographs of selected wood fragments that were used for radiocarbon dating.

Methods

Cave Geometry

The caves were surveyed using a calibrated laser distoX2 with an accuracy better than 0.05 m for lengths less than 30 m, and less than 0.5° for azimuth and inclination. A comparison between the original cave survey of 1994 as well as photographic evidence from 1994–1997 and 2017 (ICCC and JSPDT 2007) were used to assess the decadal variation in snow cone height in the main ice chamber (Figure 4).

The cave geometry plays a critical and heavily localized role in promoting or preventing snow accumulation in specific parts of the cave. The bottlenecks of a typical shaft-and-canyon series may easily be blocked by snow or falling debris, greatly diminishing the exchange of air mass between the outside and the cave interior. As a sag-type cave, the descending morphology of M-17, with its two separate vertical entrances, also has a profound influence on its meteorology and airflow dynamics and consequently on its annual energy balance. Upon each visit, the state of the connection between entrance a/b and c/d was noted and photographed.

Ice and Water Sampling

Ice, snow, and drip water were sampled in M-17 over the course of three field campaigns. Ice samples were extracted with the aid of an ice screw and transferred to polyethylene bags, where ice was allowed to melt with minimal head space, before a final transfer to standard laboratory 2-mL vials. The δ ²H and δ ¹⁸O values of these samples were measured using a PICARRO L2140–I cavity ring down spectroscopy analyzer (2019) at the University of Innsbruck, following the procedures outlined by Geldern and Barth (2012). All stable isotopic compositions are reported in permil (‰) against the Vienna Standard Mean Ocean Water (VSMOW). Calibration of the instrument was accomplished using VSMOW2, GISP2, and SLAP standards. The long-term 1-σ analytical precision of the δ ²H and δ ¹⁸O values is 0.05 and 0.1‰ (1σ), respectively. Monthly precipitation collected at the Ljubljana Global Network of Isotopes in Precipitation (GNIP; Vreča et al. 2014), 60 km from the study site, was used to compute a local amount-weighted meteoric waterline from 336 monthly observations (Figure 6), using the formula of Hughes and Crawford (2012).

Figure 6 Stable isotope data of three types of firn, congelation ice and drip water sampled at M-17. For each plot, the solid line is the Local Meteoric Water Line at the Ljubljana GNIP station (see text for details) and the dashed line is a least squares linear regression, whose equation is given in each of the firn and congelation ice diagrams. No regression was performed on the drip water samples due to the paucity of samples and sporadic nature of their collection.

Wood Remains Sampling Strategy

Wood samples were selected and taken at 17 locations from the firn deposit of M-17 cave (Figure 5). A total of 13 samples from the bottom of the ice shaft in M-17 (Figure 3 and 5), comprising a mixture of small needles, branch, twig, or pinecone fragments were radiocarbon-dated. The sampling locations shown on Figure 5 are based on a photo-mosaic (suppl. Figure 1) where the ice strata dip very steeply (∼80°) to the east. Woody macro remains are mostly derived from mountain pine (Pinus mugo), with many living specimens currently surrounding the cave entrance. This species dominates an altitude belt spanning 1650–1900 m asl around Tolminski Migovec. Sampling of woody material was performed during three visits to the cave, with preceding dating results guiding the subsequent sampling strategy. The first six samples suggested a younging direction from right (west) to left (east), whereby samples MIG-1, MIG-3 and MIG-5 were selected from prominent debris and organic rich layers. The other samples were taken from otherwise organic-poor section of the ice body. The second sampling campaign focused on the western (a priori oldest) part of the stratigraphic section, in order to test the above hypothesis. Samples MIG-7 to MIG-9 are a mixture of very small twigs or needle fragments embedded in rather organic-poor ice, while sample MIG-10 was sampled in the same organic rich layer as MIG-1. Snow retreat between 2019 and 2020 revealed a previously obscured 3-m-wide section of the stratigraphy between the MIG-3/MIG-4 layer and MIG-5 sample location. The third sampling focused on this part of the stratigraphic section, as an attempt to construct a monotonic growth model for this exposure.

The stratigraphy presented in Figure 5 could be sub-divided into individually accumulated units of cleaner ice where a total of 12 organic fragments were sampled. Samples MIG-3 and MIG-4 were found in the prominent organic-rich layer in the center of the exposure. Samples MIG-5 and MIG-10 was taken from separate irregular organic-rich layers on either side.

Sample MIG-11 was taken between the main ice chamber and the main outcrop in the ice shaft (see Figure 3). This sample was found in organic-poor, sub-horizontal layers exposed just underneath the main ice chamber level. Located 20 m above the main exposure, MIG-11 could not be correlated to any of the units of the underlying outcrop.

M-15 cave hosts a 3-m-high and rather organic-poor ice sequence (Figure 2) at the bottom of a 50-m-deep vertical pit. One wood fragment (MIG-C15, Table 2) was sampled at the base of the ice body.

Radiocarbon Dating

The ¹⁴C/¹²C and ¹³C/¹²C ratios were measured by accelerator mass spectrometry (AMS) at the ¹⁴CHRONO Centre, Queen’s University Belfast, following acid-alkali pre-treatment and graphitization procedures described in Reimer et al. (2015). Ages were calculated according to Stuiver and Polach (1977) and calibrated using the IntCal20 curve (Reimer et al. 2020), further modelling was carried out using OxCal 4.4 (Bronk Ramsey 2009). Conventional radiocarbon ages are reported with a 1σ error (Table 2). Calibrated and modeled ages are reported with two standard deviations (2σ).

The dated samples were first modeled together as part of an OxCal P_Sequence (Bronk Ramsey 2009) but building a monotonic age-depth model from this set of samples requires ignoring more than half of the radiocarbon dates. In this study, we use instead the Kernel Density Estimate model (KDE_Model) provided by OxCal 4.4 (Bronk Ramsey 2017), which returns an estimate of the underlying distribution function from which radiocarbon dates were sampled (Figure 7). Therefore, the radiocarbon age distribution of fragments issued from otherwise clean ice may be thought of as a marker of ice cave net ice-mass growth, in the very same way that the time-aggregated probability density functions (PDF) of radiocarbon dates is a common method for summarizing a radiocarbon dataset, often interpreted as a proxy for activity (Bronk Ramsey 2017). To this end, Spötl et al. (2014) and Kern et al. (2018) displayed the sum of probability densities of the radiocarbon datasets collected at HETH and Velebit ice caves, respectively. The KDE method has the added benefit of smoothing some of the noise introduced by summation of marginal posterior PDFs. We interpreted the KDE of radiocarbon dates from fragments found at organic-rich unit boundaries as evidence for past negative mass balance. Conversely, the KDE of radiocarbon-dated wood fragments found in clean ice as evidence for past positive mass balance.

Figure 7 Conceptual model for the use of radiocarbon-dated remains to constrain past periods of positive mass balance in ice caves. (a) Example history of ice dynamics, with alternating periods of positive (blue) and negative (orange) mass balance; organic macro remains (indicated by leaf symbols) are assumed to have been deposited at regular time intervals; those entering the cave during times of negative mass balance are shown in red. (b) Corresponding changes in ice body height with time (grey line) and (c) the expected stratigraphy consisting of comparably clean ice layers (A, B, C) and intervening horizons rich in macro remains marking periods of negative mass balance. The solid red line denotes the best possible age-depth model if all macro remains in the sequence were radiocarbon-dated. Flat sections of this line denote hiatuses, which overestimate the duration of mass loss in the cave as indicated by horizontal red arrows. The vertical gray arrows symbolise the vertical displacement of macro remains during intervals of negative mass balance. (d) Kernel Density Estimate of the age distribution based on all macro remains found in the comparably clean ice layers A, B, and C, where the peaks appropriately highlight past periods of positive mass balance with high confidence. Low values of the KDE denote the lack of information for certain time periods and therefore cannot be interpreted as proof of mass loss.

Results

Stable Isotopes of Ice

Statistics on the stable isotopic composition of different types of ice occurrences are presented in Table 1, where firn (densified snow with an infiltration component sensu Luetscher and Jeannin 2004) is distinguished from congelation ice of seasonal speleothems and the floor of the main ice chamber. In total, 26 samples of firn ice, 9 samples of congelation ice and 8 samples of drip water collected in M-17 are presented in Table 1 and Figure 6.

Table 1 Basic statistics on the stable isotope composition of ice and water sampled in M-17 ice cave.

Firn and ice δ¹⁸O values in M-17 range from –10.89±0.06‰ to –8.14±0.03‰ (mean:–9.6±0.8‰, n = 26) along the section sampled. A Wilkes-Shapiro test on the firn ice samples failed to reject the null hypothesis of normality at a significance level of 0.055 (W-statistic: 0.949, p-value: 0.219). Firn ice samples define a Local Meteoric Water Line (LMWL) of δ²H = 8.42±0.24×δ¹⁸O + 17.17±2.30. Congelation ice δ¹⁸O values range from –10.57±0.06‰ to –7.39±0.03‰ (mean: –9.25±1.2‰, n = 10). The Wilkes-Shapiro test on the congelation ice samples fails to reject the normality assumption at a significance level of 0.05 (W-statistic: 0.897, p-value: 0.234). Congelation ice samples define a LMWL of δ²H = 7.75±0.28×δ¹⁸O + 9.67±2.65.

Drip water δ ¹⁸O values range from –10.50±0.08‰ to –8.89±0.04‰ (mean: –9.3±0.6‰). A recent global isotopic analysis by Baker et al. (2019) concluded that European cave sites show a strong positive correlation between the recharge-weighted δ ¹⁸O of precipitation and that of drip water. High total annual precipitation (> 1750 mm yr^-1), mean annual Air temperature < 10°C, and low annual potential evapotranspiration (PET < 800 mm yr^-1) were found to maximize this correlation (Baker et al. 2019). M-17 cave meets the first two criteria, suggesting that the recorded cave drip waters provide a reliable estimate of the local amount-weighted precipitation. The calculated drip water d-excess is greater than 10, varying between 11.8 and 13.6. No regression was performed on the drip water samples due to the paucity of samples and sporadic nature of their collection. Samples were taken during the summer visits to the cave only and could not possibly capture the full seasonal range of stable isotope values in precipitation at the site.

The d-excess values of firn-derived ice, congelation ice, and drip-water are consistent with a western Mediterranean moisture source (Gómez-Hernández et al. 2013). The isotopic range of 2.8‰ for δ ¹⁸O in the ancient firn-derived ice is inferior to the seasonal range recorded at the nearest GNIP station of Ljubljana (suppl. Figure 11). This suggests at least some mixing between winter snow with lower δ ¹⁸O values and summer precipitation with high δ ¹⁸O. The δ ¹⁸O–δ ²H relationship in firn suggests that minimal disequilibrium fractionation in involved in the diagenesis of snow in M-17.

Dating of Wood Macro Remains

Stratigraphic Context

A summary of the calibrated dates is given in Table 2. The samples span the 10th to 15th century AD within a 95.4% confidence interval. Even in this preliminary analysis, the independently calibrated (but unmodelled) radiocarbon dates are not correlated to their stratigraphic position. This outcrop had been previously hypothesized as younging from right (west) to left (east) (Figure 5). Indeed, samples found within the same organic-rich layers exhibit ages with no overlap, e.g., MIG-3 spanning cal 1267–1388 AD (2σ) and MIG-4 spanning cal 1032–1203 AD (2σ). One possible interpretation of this complex age progression is that organic material remnants apparently too old for their given stratigraphic position (e.g., MIG-2, MIG-C2, and MIG-C4) have been re-mobilized from an earlier position in the ice body and reincorporated in relatively younger ice.

Table 2 Radiocarbon data of wood samples from M-17 and M-15 ice caves. Values in parentheses denote the area under the probability density distribution curve.

^a Calibrated using the CALIBomb software (Reimer et al. 2004).

One sample located nearest to the ice chamber, at a depth of 50 m (MIG-C11), indicates positive annual mass balance between 1600–1800 AD, i.e., during the LIA. It is the only sample of this period found inM-17. The stratigraphic context for this sample is an ice section located under the doline c/d feeder (Figure 3), lacking macro remains and exhibiting near horizontal stratification.

The sample from M-15 cave was dated and a “modern” age (post 1950 AD) was obtained. This age of wood remain was also calibrated using the CALIBomb software (Reimer et al. 2004) and the dataset for post-bomb tropospheric radiocarbon spanning 1959–2012 (Levin and Kromer 2004; Levin et al. 2013). This yielded a calibrated age spanning cal 1963–1967 AD (2σ), further highlighting the relative youth of the M-15 ice depositcompared to that hosted in M-17.

Summarizing the Radiocarbon Dataset

Dating organic remains preserved in the ice provides constraints both on the mean annual mass balance of sections showing a regular age progression and on the duration of hiatuses (periods of negative ice mass balance) in the case of sequences characterized by a step-like age progression. We dated a set of 18 samples in M-17, including 13 samples in a continuous exposure located in the ice shaft (Figures 3 and 5). However, the challenging geometry of this stratified deposit, with its steeply dipping strata and apparent unconformable contacts between strata of variable organic content precludes a simple monotonic age-depth model.

Instead, the Kernel Density Estimates (KDE_Model, Bronk Ramsey 2017) of the dated samples from this exposure were used to compute a probability of past positive mass balance over time for M-17. The KDE indicates probable positive mass balance between 900 and 1300 AD (Figure 8). This suggests at least one phase of positive mass balance during the Medieval period. Despite detailed sampling, no wood inclusions significantly older than 900 AD were found. This by no means precludes a positive ice mass balance prior to this date. The complete lack of woody macro remains at the very base of the ice body, however, documents that the bottom of this shaft was ice-free for at least a few decades prior to the onset of the snow and ice deposition commencing at around 900 AD, allowing for the full decay of any organic material that had fallen into the cave beforehand. Samples derived from organic rich layers indicate a peak of probably mass loss at around 1350 AD (Figure 8). A final mode of the probability of positive mass is centered around 1700 AD, based on one sample only (MIG-C11, Table 2). This single sample suggests strongly, however, that M-17 experienced strongly positive mass balance during the LIA, as it is embedded in a thick section of relatively clean ice overlying the Medieval period ice section.

Figure 8 Comparison of the radiocarbon-dated records from selected caves of the Dinarides and Eastern Alps. Kernel Density Estimates of the age distribution of ice sections are plotted in blue (orange) when they denote periods of positive (negative) mass balance. Low KDE values denote the absence of information for that period. Sections marked with an asterisk denote calibrated posterior probability distributions for individual samples. Red markers denote the individual posterior median ages underlying the KDE estimates.

Discussion

The superposition of clean stratified ice overlain by “dirty” layers with a slight unconformity is a common feature in sag-type ice caves (e.g., Spötl et al. 2014; Kern et al. 2018). One key to understanding the mass balance history of underground firn and ice deposits is the interpretation of prominent wood-rich layers such as the layer containing MIG-5 or MIG-3 and MIG-4 (Figure 5), with respect to past ice mass balance. In the exceptional case of seasonally freezing pools, the seasonality of organic infall and ice loss leads to the formation of an alternating sequence of summer-autumn organic-rich and clean winter-spring ice layers (Perşoiu et al. 2017; Kern 2018) resulting in an annually resolved stratigraphy. This is not the case in M-17, where clean, rather organic-poor ice of multi-annual to multi-decadal origin are present.

The time gap between wood death and its subsequent transport and eventual inclusion into the ice is a recognized bias in dating wood inclusions rather than the ice directly (Luetscher et al. 2007; Perşoiu 2018). Indeed, this lag between wood death and its incorporation into the ice deposit may be substantially longer than the error in radiometric dating of that very sample (Gradziński et al. 2016). M-17 is largely devoid of such large wood samples, with MIG-3 being the only sub-sample from a Pinus mugo branch about 3 cm in diameter. The time gap between wood death and its transport into the ice cave, however, comprises no more than a few decades for fragile fragments such as needles or small twigs, between one order of magnitude lower and the same order as the uncertainty in calibrated radiocarbon ages.

Once transported into the cave, fragile woody macro remains such as needles and small twig fragments are stored within the preserving matrix of the ice body. Wherever such organic remains occur as sparse inclusions in otherwise clean ice, their presence suggests that the organic sampled were quickly embedded in the ice deposit due to e.g., rapid ice buildup in the cave (e.g., Gradziński et al. 2016; Spötl et al. 2014). Many of the samples taken from M-17 cave fall in this category, except for MIG-1, MIG-3, MIG-5 and MIG-10, which were sampled in prominent organic- and cryoclast-rich layers. These striking layers, often in unconformable stratigraphic contact with the underlying ice, have been interpreted as markers of prolonged negative multiannual or multidecadal mass balance (Kern et al. 2018; Sancho et al. 2018). Organic material present in such layers can be interpreted as mixture of materials dating opposite events: (1) material previously deposited and embedded in the ice, but released by the ice retreat and remobilized, thus dating a period of positive ice mass balance or (2) material transported into the cave during the period of ice retreat, linked therefore to a period of negative ice mass balance (Figure 7).

The repeated sampling of organic rich layers may thus yield remarkably different ages, tied both to periods of positive and negative mass balance. The maximum age difference between radiocarbon ages from a single layer may help bracket the length of ice ablation period, as proposed by Kern et al. (2018) for the Vukušić deposit. Such a method may however overestimate the length of hiatus given that the older date may reflect remobilized material tied to an earlier positive mass balance (Figure 7). While the targeted sampling of very prominent organic-rich layers only might constrain periods of negative mass balance (e.g., Kern et al. 2018), the sampling of sparse organic fragments embedded in ice over the entire stratigraphy as demonstrated by the present study is expected to constrain periods of past positive mass balance (Figure 7). The compilation of radiocarbon dates from the Velebit Mountains (Croatia) shows a remarkable coherence with the data from HETH (Kern et al. 2018, Figure 4), with maxima and minima in the summed probability distributions of the dated samples broadly agreeing before the 11th century AD and during the 15th century AD. However, while Spötl et al. (2014) interpret maxima in the summed probability distributions of the ice deposit at HETH as periods where positive mass balance was highly likely, Kern et al. (2018) adopt an apparently contradictory approach, whereby those same maxima are interpreted as indicated of time periods unfavorable for positive ice mass balance. Considering that the HETH samples are found in relatively clean ice and that those from Velebit mountains originate from organic rich layers, the two sampling approaches yield complementary records of past positive and negative mass balance, respectively. At M-17, we also combine the two sampling methodologies and present the results in Figure 8.

To better understand the key parameters influencing past mass balance in ice caves, the M-17 record is contrasted with historical temperature and hydroclimate variability in the NE Adriatic, the Alps and Central Europe. Temperature–forced records include that of Grotta Savi, where the calcite δ ¹⁸O values of stalagmite SV1 (Figure 9d), were calibrated against a gridded temperature dataset provided by Luterbacher et al. (2004), itself based on mixture of homogenized instrumental records and proxy temperature reconstructions derived from ice cores and tree-rings. The LIA sensu Frisia et al. (2005) is recorded from 1450 to 1800 AD with a declining trend in δ ¹⁸ O, reaching a minimum at about 1800 AD (Figure 9d), whereas the onset of the LIA in the Northern Alps is placed earlier at 1260–1270 AD (Nicolussi and Patzelt 2000; Büntgen et al. 2011).

Figure 9 Comparison of the M-17 radiocarbon dataset with published Alpine paleoclimate proxy records, (a) tree-ring based JJA temperature, (b) annual AMJ precipitation totals, including a 60–year low-pass filter (Büntgenet al. 2011), (c) subset of the annually resolved record of past hydroclimate from the Old World drought dataset covering 13–14°E and 46–47°N, averaged over 2 × 2 cells, centered over the NE Adriatic, including a 60-year low-pass filter (Cook et al. 2015), (d) δ ¹⁸O record of Grotta Savi (Frisia et al. 2005), (e) KDE of M-17 negative ice-mass balance derived from organic rich layers. (f) KDE of M-17 positive mass balance. (g) KDE of Unit B (h) KDE of Unit C. MC1(2): Medieval Cold period 1(2), MW1(2): Medieval Warm period 1(2), LIA: Little Ice Age (Frisia et al. 2005).

In the last 2000 years, the earliest phase of positive mass balance occurred at HETH, from 600–700 AD (Spötl et al. 2014, Figure 8), which corresponds to MC1. Identifying ice sections lacking wood at Kugina ice cave, Kern et al. (2018) interpret the period between 600–700 AD as favorable for positive ice mass balance. No counterpart is observed in the M-17 record. Positive mass balance may indeed have taken place at the M-17 site but if correct, the ice would have since been lost. At M-17, the first mode of positive mass balance is found at around 1000 AD, a peak which not shared by HETH, but which also corresponds to a period of positive mass balance according to the Kugina record. During MC1, the 60-yr filtered Palmer Drought Severity Index (PDSI) in Northwestern Slovenia get progressively wetter reaching a maximum during the middle of MW1 at 750 AD. Several Alpine glaciers also expanded in the early 9th century AD (Nicolussi and Patzelt 2000), broadly during the transition between MW1 and MC2 sensu Frisia et al. (2005). There is a parallel trend toward greater spring precipitation totals and early summer temperatures in the Northern Alps (Büntgen et al. 2011). Between 900 and 1000 AD, the KDEs of radiocarbon samples from M-17 suggest positive mass balance, at a time where temperatures were dropping in the southern Alps, the northern Alps experienced wetter than average spring climate. Thus, the strongest evidence for positive mass balance at HETH falls during MC1, while at M-17, this occurs at the transition between MW1 and MC2.

A transient shift towards wetter climate is recorded in northwestern Slovenia between 1100 and 1150 AD (Cook et al. 2015), corresponding to the transition between MC2 and MW2. The tree-ring based temperature reconstruction for the Alps also suggests that the peak temperature of the Medieval Warm Period is reached shortly after 1000 AD (Büntgen et al. 2011). At that time, the Velebit dataset suggests negative mass balance at Kugina. This also coincides with very low values of the KDE of past positive mass balance at HETH and at M-17. In other words, while there is strong evidence of mass loss in the Velebit dataset, the M-17 and HETH datasets lack strong evidence of positive mass balance, and vice-versa. The second minor peak in annual mass balance at M-17 occurred in the middle of MW2 and suggests that snow and ice were probably still accumulating despite higher summer temperatures, again lagging the Alpine glacier re-advance by about a century (Nicolussi and Patzelt 2000). In addition, a late Holocene record from a Croatian stalagmite indicates drier conditions in the Central Adriatic during this period (Rudzka et al. 2012). The shift towards drier summers during MW2 may have counter-balanced the effect of elevated summer temperatures altogether by mitigating the advection of sensible heat from warm rain events, as reported by Perşoiu et al. (2021).

The Velebit and M-17 datasets of negative mass balance both indicate with high probability that caves experienced mass loss at around 1300 to 1400 AD. This period sees the KDE values at HETH increase sharply to their second peak at 1300 AD, before decreasing to low values again during the period 1400–1600 AD (Figure 8). The peak at HETH corresponds to the start of the LIA in the northern Alps, which is picked out by the start of a multi-century decreasing trend in temperatures (Büntgen et al. 2011). At M-17, the radiocarbon-dated remains found in organic rich layers suggest probable mass loss during the peak of MW2. This retreat coincides with elevated air temperatures from SV1 and a wetter than average hydroclimate.

The onset of this LIA phase is well bracketed by radiocarbon dated remains at HETH. It is dated to 1300 AD (Figure 8, Spötl et al. 2014). Despite the presence of prominent layers of the Velebit dataset, positive mass balance during the LIA is postulated at Vukušić ice cave and Ledena pit (Kern et al. 2018). One sample in the middle of the layered, but otherwise organic-poor, upper ice section of M-17 (MIG-11) also indicates that ice also accumulated at the time. Based on radiocarbon-dated remains and stratigraphic constraints, we find that ice cave datasets in the southern Alps point to overall positive ice mass balance during the coldest phase of the LIA, at a period when the tree-ring evidence suggests average hydroclimate conditions in the Alps (Cook et al. 2015). In addition, Croatian stalagmite record indicates that the LIA was on average wetter there (Rudzka et al. 2012). This corroborates an earlier reconstruction of Mediterranean winter precipitation by Luterbacher and Xoplaki (2003), suggesting a period of elevated precipitation and colder temperatures during the Maunder minimum (1645–1715 AD), relative to the average during the period 1500–1995 AD. The apparent lack of further LIA age wood samples in M-17 and the Velebit dataset can be explained by the dilution of organic fragments during periods of strongly positive ice mass balance. The northern Alps experienced two transient shifts towards a wetter climate at 1650 AD and 1810 AD, corresponding to the Maunder and Dalton solar minima, respectively (Büntgen et al. 2011), when periods of likely positive ice mass balance are also recorded in HETH (Spötl et al. 2014). Ice caves of both sides of the Alps appear to be highly responsive to multi-decadal periods of colder and wetter winters, which led to substantial ice gain.

Regarding the past availability of wood samples, while the elevation of the closed forest timberline may have migrated in response to colder temperatures during the LIA, this would not have necessarily resulted in a depression of the dwarf pine belt to the same extent. The krummholz specimens of dwarf pine survive winter frost mainly thanks to a protective snowpack layer (Troll 1973). Assuming a negligible change in the input of wood into the cave, the overall wood-free appearance of the upper ice body can be explained by a higher ratio of snow to organic input into the cave at that time.

In addition, the geometry of the cave may put a fundamental limit on the snow accumulation rates, as has been invoked in the case of cave A294 in the Central Pyrenees (Sancho et al. 2018). At several points during the LIA for instance, decades of excessive snow input within M-17, driven by wet and cold winters (Luterbacher and Xoplaki 2003) may have fully blocked both entrances with an ice plug preventing the throughput of organic material as well as the deposition of additional snow due to the resulting lack of accommodation space.

The paleoclimatic context of mass balance changes in the southern Alps, recorded by the M-17 and Velebit datasets suggest that changes of air temperature and hydroclimate together are required to explain ice loss or gain in alpine caves. The phases MC1 and MC2 coincide with well bracketed periods of positive mass balance in HETH and M-17, respectively. There is also evidence of ice gain during the first half of MW2 at M-17, despite elevated air temperatures. The wetter, second half of MW2, however, sees stronger evidence for mass loss, which can be explained by the shift to greater amounts of warmer, liquid precipitation entering the cave. The colder air temperatures of the LIA coincide with generally positive mass balance in caves north and south of the Alps, despite their differing hydroclimates. The relative scarcity of wood dating back to the LIA in M-17 precludes a more detailed discussion of this ice cave’s response to the several distinct phases of the LIA, which incidentally saw repeated glacier advances. Additional age-depth chronologies from simpler ice exposures are required to tackle these questions.

Conclusion AND Outlook

We presented the first detailed chronostratigraphic investigation of a site from the Julian Alps (northwestern Slovenia). Strategically sampled wood remains dated by radiocarbon record peaks in mass balance at 900–1100 AD, 1200–1300 AD and to a lesser extent at 1700–1800 AD. These intervals coincide with periods when the 60-yr filtered local PDSI (JJA) was near zero or slightly negative, corresponding to average or dry summers. The 1200–1300 AD peak in likelihood of positive mass balance is also recorded in underground ice deposits in both the Northern Calcareous Alps (Spötl et al. 2014) and the Dinarides (Kern et al. 2018). Positive mass balance is also recorded during periods when the regional air temperature anomaly was negative or positive (Frisia et al. 2005). Stable isotopic values of cave ice plot along a local meteoric water line consistent with values recorded at the Ljubljana GNIP station. The double funnel geometry of the M-17 cave likely plays an important role in the distribution of woody material and the development of the ice stratigraphy. The opening or closing of one or both entrances fundamentally alter the geometry of the cave, affecting snow transport and redistribution and thereby controls the influx of organic remains into the cave.

Photographic documentation suggests that overall, M-17 experienced a negative mass balance during the three last decades. The loss is mostly confined to the upper part of the cave. Within the entrance shaft, an ice-plug located at about –20 m depth is now gone and replaced by a seasonal snow cone. Multi-annual ice speleothems, like the 3-m-tall ice-stalagmite from the main chamber, also melted to a significant extent and have since been replaced by depressions and holes on the ice surface where ceiling drips have contributed to the ice decay. Recent ice retreat is not limited to M-17 nor nearby ice-bearing caves. Such loss has been recognized in many ice caves worldwide with a lost ice volume between 1908 and 2010 AD estimated at about 23,000 m³ (Kern and Perşoiu 2013). Substantial ice loss was recently reported from ice caves of southeastern Europe, with the extreme rainfall events of summer 2019 likely playing a key role in the accelerated demise of these ice bodies (Perşoiu et al. 2021). In the specific case of M-17, we propose that exceptional rainfall at the beginning of November 2019 played a key role in opening the connection from the ice chamber to c/d doline and allowing forced convection between the two entrances. This result highlights a dual vulnerability of M-17 to the expected frequency increase of similarly extreme rainfall events over the 21st century (Púčik et al. 2017). Ice melt due to sensible heat advected by liquid water into the cave notwithstanding, the permanent opening of the doline c/d connection caused by a drop in the ice level would allow forced convection between the two entrances, whose immediate consequences are not obvious, but may for instance enhance ice sublimation.