3.1 Drilling and stable isotope analyses

The FV ice core (4.87 m long and 10 cm diameter) was drilled in May 2016 from the Great Hall of FV (Fig. 1d) using a modified PICO electric drill (Koci and Kuivine, 1984) manufactured by Heavy Duties S.R.L, Cluj-Napoca, Romania. A rock embedded in the ice at 4.87 m below the surface stopped the drilling effort, but previous work in the cave has shown that the thickness of the ice block exceeds 15 m (Orghidan et al., 1984; Kern et al., 2004; Perșoiu and Onac, 2019). The ice core was cut into 1 cm long pieces (taking the annual layering into account), and each piece was subsequently sealed in plastic bags, allowed to melt at room temperature, transferred to 20 mL HDPE (high-density poly ethylene) scintillation vials and stored at 4 ^∘C prior to analysis.

Precipitation samples were collected monthly between March 2012 and December 2018 at Gheţar (GT, 46^∘29^′28.45^′′ N, 22^∘49^′26.02^′′ E; 1100 m a.s.l.; ∼ 13 km southeast of the FV) using collectors built according to International Atomic Energy Agency (IAEA) specifications.

Water samples were analyzed for stable isotope composition at the Stable Isotope Laboratory, Ştefan cel Mare University (Suceava, Romania), using a Picarro L2130-i CRDS (cavity ring-down spectroscopy) analyzer connected to a high-precision vaporizing module. All samples were filtered through 0.45 µm nylon membranes before analysis and were manually injected into the vaporization module multiple times until the standard deviation of the last four injections was less than 0.03 for δ¹⁸O and less than 0.3 for δ²H. The average of these last four injections was normalized on the SMOW–SLAP scale using two internal standards calibrated against VSMOW2 and SLAP2 standards provided by the IAEA and further used in our interpretation. A third standard was used to check the long-term stability of the analyzer. The stable isotope values are reported using standard δ notation, and the precision is estimated to be better than 0.1 ‰ for δ¹⁸O and better than 0.5 ‰ for δ²H based on repeated measurements of an internal standard.

3.2 Radiocarbon dating and age–depth model construction

The wide opening to the surface in the ceiling of the Great Hall (Fig. 1c) allows for a large volume of organic matter to fall into the cave and subsequently become trapped in the ice, including large pieces of wood that tend to cut through layers of ice, with the latter possibly encompassing decades of ice accumulation. Furthermore, ice melting and water freezing processes usually result in inclined ice surfaces and slow tipping of any heavy materials sitting on top of the ice (see the position of tree trunks in Fig. 1). As a consequence, wood with an age older than that of the newly forming ice can be incorporated in the ice block, resulting in sample ages much older than the ice layers (“old wood effect”). However, the challenge of identifying such organic material contrasts with the desire for precise chronologies, which rely on a high number of data points. Thus, all possible organic samples were recovered and dated, and unreliable ages that were identified at the stage of age–depth modeling were removed. Out of 14 samples recovered from the ice core and potentially suitable for radiocarbon dating, two were not datable due to their extremely small carbon yield. Accelerator mass spectrometry (AMS) radiocarbon analyses were performed at the Institute of Physics, Silesian University of Technology, Poland (Piotrowska, 2013). All samples were precleaned with standard acid–alkali–acid treatment, dried and subjected to graphite preparation using an AGE-3 system (Ionplus, Switzerland) equipped with an Elementar vario MICRO cube elemental analyzer and automated graphitization unit (Wacker et al., 2010; Nemec et al., 2010). The ¹⁴C concentrations in graphite produced from unknown samples, Oxalic Acid II standards and coal blanks of comparable carbon masses were measured by the DirectAMS laboratory, Bothell, USA (Zoppi et al., 2007). The results are reported in Table 1. The radiocarbon dates were calibrated using OxCal v4.3 (Bronk Ramsey, 2009) and the IntCal13 calibration curve (Reimer et al., 2013). The NH1 curve (Hua et al., 2013) was used for one post-bomb date.

Table 1Radiocarbon data from the Focul Viu Ice Cave. Agreement indices for individual samples based on the P_Sequence algorithm (Bronk Ramsey, 2008) are provided for accepted dates. Modeled ages for all dated depths are given as mean and σ values, rounded to the nearest five.

Download XLSX

Because organic material can fall into the cave decades to centuries before being trapped in the ice (see Fig. 1b), we have carefully screened the radiocarbon results prior to age–depth modeling with the aim of selecting the most reliable dates forming a chronological sequence. In total, four dates were selected for age–depth modeling. For the top of the ice core, a uniform age distribution from 1991 to 2016 CE was assigned, allowing for the possibility of surface ice melting. The model was constructed using the OxCal P_Sequence algorithm (Bronk Ramsey, 2008) with a variable prior k parameter (k=1, $U(-\mathrm{2},\mathrm{2})$; Bronk Ramsey and Lee, 2013) and extrapolated to a depth of 4.86 m. The agreement index of the model was 85 %, confirming a good statistical performance when the threshold of 60 % is surpassed. The sections between the dated depths were assumed to have a constant deposition rate. The complete age–depth model is shown in Fig. 2. The mean age derived from the model was used for further analysis and is also reported in Table 1. The constructed age–depth model was compared with that of Maggi et al. (2008), plotted in green in Fig. 2, and a broad agreement was found between both chronologies. All of our rejected ages are older than those of Maggi et al. (2008), and we suspect that these were based on dating old wood that had already been in the cave for decades before being incorporated in the ice (see the “old wood effect” discussion above and Fig. 1c). Further, the same authors identified several potential markers of volcanic eruptions, and their ages agree within ±40 years with those of our model. Another evaluation was obtained by comparing our stable isotope record with those of Kern et al. (2004) and Forizs et al. (2004); although the latter lack a precise chronology, a simple visual correlation between the records indicates a satisfactory match. In order to avoid any circular reasoning, we decided against the use of other records, e.g., the regional temperature record, to further anchor our chronology.

Figure 2Age–depth model of the Focul Viu ice core. The calibrated age range of samples used in the model is indicated in blue; rejected samples are indicated in red. Samples in green are from the ice core drilled in 2004 (Maggi et al., 2008). Dark and light blue shading indicates the 95 % and 68 % confidence ranges of the model, respectively.

Download

3.3 Climate data

The sea surface temperature (SST) is extracted from version 5 of the Extended Reconstructed Sea Surface Temperature data (ERSSTv5) of Huang et al. (2018). This dataset covers the period from 1854 to present and has a spatial resolution of 2^∘ × 2^∘. The AMO index used in this study was obtained from https://climexp.knmi.nl/data/iamo_ersst_ts.dat (last access: 19 October 2020) and is also based on the ERSSTv5 dataset. Station-based meteorological data were provided by the Romanian National Meteorological Administration for three stations (Baia Mare, Sibiu and Timișoara) that have some of the longest instrumental records in Romania and bracket the location of the study site. To remove the short-term variability and retain only the multidecadal signal in our data, both the temperature time series and the SST data were smoothed with a 21-year running mean filter prior to correlation analysis.

4.1 Ice accumulation in Focul Viu Ice Cave

The results of the radiocarbon analyses performed on organic matter recovered from the ice are shown in Table 1, and the age–depth model is displayed in Fig. 2. The maximum age of the ice is 1000 ± 20 cal BP at 4.45 m below surface, based on direct dating of organic remains (Table 1) and extrapolated to 1100 cal BP at 4.86 m below surface

High accumulation rates were recorded between 850 and 950 CE (0.39–0.41 cm yr⁻¹) and between 1220 and 1970 CE (0.36–0.44 cm yr⁻¹). Between 950 and 1220 CE, the net accumulation rate dropped to between 0.29 and 0.34 cm yr⁻¹. The highest net accumulation rates recorded after 1970 CE (0.56 cm yr⁻¹) contradict recent findings from other ice caves in the Carpathian Mountains (Kern and Perșoiu, 2013), which all register record melting. However, this value might be an artifact of the age–depth modeling (see above) as well as of the very short time span considered; thus, it is unreliable for further interpretation. The low accumulation rates spanning the MWP are similar to those recorded in the Scărișoara Ice Cave (Perșoiu et al., 2017, Bădăluţă, 2019), the Hundsalm ice cave in Austria (Spötl et al., 2014) and the ice caves in the Velebit Mountains in Croatia (Kern et al., 2018), potentially suggesting a regional signal of climatic conditions unfavorable for ice accumulation. Ice can melt as a result of either warm summers with enhanced conductive heat transfer to the cave or wet summers with rapid ablation resulting from water flowing across the top of the ice block. Subsequently, ice growth is influenced by the amount of water present at the onset of freezing, the timing of this onset and its duration. The low accumulation rates during the MWP were likely the result of enhanced melting during warm (see Sect. 4.2 below) and wet (Feurdean et al., 2015) conditions. After 1450 CE, the climate in the region was dominated by dry summers with frequent storms and cold winters (Perșoiu, 2017). These conditions led to reduced summer melting and enhanced winter growth, which are conditions favorable for net ice accumulation.

4.2 Stable isotopes in Focul Viu cave ice – a proxy for summer air temperatures and AMO variability

The variability of δ¹⁸O and δ²H in precipitation at Gheţar (∼ 13 km south of the cave's location and at the same altitude), assessed for the 2012–2017 CE period, follows that of temperature (Fig. 3a), with the maximum values (−3.6 ‰ and −26 ‰ for δ¹⁸O and δ²H, respectively) in July and August and the minimum values (−19.8 ‰ and −140 ‰ for δ¹⁸O and δ²H, respectively) in January. Similar results were found by Bojar et al. (2009) and Ersek et al. (2018) for the same region, suggesting that the ¹⁸O ∕ ¹⁶O and ²H ∕ ¹H ratios in precipitation register temperature changes on a regional scale. The Local Meteoric Water Line, defined by the equation δ¹⁸O = 7.4 × δ²H + 6.1 (Fig. 3b), has a slope and intercept very similar to those found by Ersek et al. (2018). Precipitation in the region is mainly delivered by weather systems carrying moisture from the Atlantic Ocean, with the Mediterranean Sea contributing moisture during autumn and winter (Nagavciuc et al., 2019b). The deuterium excess in precipitation (d-excess or d), defined as d=δ²H − 8 × δ¹⁸O (Dansgaard, 1964) allows for a clear separation of the air masses: the Atlantic Ocean (d-excess close to the global average of 10; Craig, 1961) and the Mediterranean Sea (d-excess between 12 and 17, resulting from the high evaporative conditions in the eastern Mediterranean Sea). Similarly high values of d-excess have been measured in precipitation in southwestern Romania (Drăgușin et al., 2017) and in precipitation in northeastern Romania (Bădăluţă et al., 2019) and have been linked to air masses originating in the strongly evaporated Mediterranean and Black seas, respectively.

Figure 3(a) Temporal variability of δ¹⁸O and δ²H in precipitation and air temperature at Gheţar (10 km south of the Focul Viu Ice Cave and at the same altitude). (b) The Local Meteoric Water Line (LMWL) of precipitation at the same station, plotted against the Global Meteoric Water Line (GMWL).

Download

Observations on the dynamics of cave ice during the past 18 years have shown that it starts to grow in early autumn due to the freezing of water accumulated during summer. As the ceiling of the cave is opened to the outside, precipitation directly reaches the site of ice formation; thus, the stable isotope composition of precipitation is not modified in the epikarst above the cave, and the original δ¹⁸O and δ²H values of summer (June–July–August, JJA) precipitation are preserved in the cave water. However, while freezing processes in caves could alter the original δ¹⁸O (and δ²H) values in cave ice, several studies have shown (e.g., in the nearby Scărișoara Ice Cave; Perșoiu et al., 2011b) that the original climatic signal embedded in the stable isotope composition of cave ice is preserved and can be used as a proxy for external climate variability.

Figure 4Temporal variability of the Atlantic Multidecadal Oscillation instrumental index; summer (JJA) air temperature (anomalies with respect to the 1961–1990 CE period) recorded at the Baia Mare (BM), Timișoara (TM) and Sibiu (SB) weather stations; and FV δ¹⁸O and δ²H (‰) during the instrumental period. In the lower panels, the positive (red) and negative (blue) anomalies are shown against the 1850–2000 CE averages for the FV δ values.

Download

Overall, our observations of cave ice genesis and dynamics and the stable isotope monitoring data clearly indicate that summer air temperatures are registered and preserved in the ice block in the FV. In order to test the long-term stability of these connections, we have analyzed the links between the FV δ¹⁸O_ice record and instrumental data from three nearby meteorological stations over the 1851–2016 CE period. On multidecadal timescales, summer air temperature changes in the region are mainly controlled by the dynamics of the Atlantic Multidecadal Oscillation (Ionita et al., 2012). Figure 4 shows the JJA air temperature at the Baia Mare (BM), Sibiu (SB) and Timișoara (TM) stations, the AMO index and the FV δ¹⁸O_ice. The instrumental temperature data indicate large multidecadal variability, with a cold period between 1890 and 1920 CE, followed by a warm period between 1921 and 1960 CE, a slightly colder period between 1960 and 1980 CE and enhanced warming after 1980 CE, which all follow the AMO variability. The δ¹⁸O_ice values show a similar temporal evolution, with the slight offsets between the observational data and δ¹⁸O likely being due to the dating uncertainty (20–35 years). Further, we have computed the correlation map between the summer mean air temperature at SB station (with the longest instrumental record) and the summer SST as indicator of AMO variability (Sutton and Dong, 2012). Figure 5 clearly shows that positive (negative) temperature anomalies over the analyzed region are associated with positive (negative) SST anomalies over the North Atlantic Ocean, resembling the SST anomalies associated with the positive (negative) phase of the AMO (Mesta-Nuñez and Enfield 1999; Latif et al., 2004; Knight et al., 2005). The strongest correlations (Fig. 5) are found between local air temperature and SSTs in the North Atlantic and the eastern Mediterranean Sea, which are the main sources of moisture feeding local precipitation (Bădăluţă et al., 2019). These results are also in agreement with the results of Della-Marta et al. (2007), who showed that summer positive temperature anomalies and heat waves over Europe are triggered, at least partially, by the phase of the AMO. A recent study of δ¹⁸O variability in oak tree rings in northwestern Romania (∼ 75 km north of our site) also indicated the influence of the AMO on summer temperatures and drought conditions (Nagavciuc et al., 2019a). Thus, we suggest that, during summer, strongly meandering Rossby waves (Ionita et al., 2015, 2017) result in blocking conditions over central Europe that lead to the persistence of high-pressure systems and the occurrence of regional heat waves. These, in turn, favor regional recycling of moisture, resulting in positive δ¹⁸O anomalies in precipitation that are further recorded by climate proxies: stable isotopes in cave ice (this study) and tree ring parameters (width and cellulose stable isotopes) (Popa and Kern, 2009; Nagavciuc et al., 2019a). The increased contribution of recycled moisture to precipitation has been reported from several locations in central Europe, further supporting our inference (Goìmez-Hernaìndez et al., 2013; Kern et al., 2020).

Figure 5Spatial correlation map between sea surface temperature (SST) and average summer (JJA, June–July–August) air temperature at Sibiu (60 km south of Focul Viu Ice Cave, indicated by the green arrow) over the 1850–2011 CE period.

Combining all of the abovementioned data, we find that, on timescales ranging from years to decades, prolonged periods of positive temperature anomalies throughout the summer months, linked to prolonged warm SSTs in the North Atlantic Ocean (and thus a positive AMO index), could be preserved by the δ¹⁸O_ice in FV.

Figure 6Summer climatic conditions recorded by δ¹⁸O and δ²H from the FV ice core (panels e and f, bottom) and comparison with proxy indicators from the Northern Hemisphere: (a) central European summer temperature anomalies (against the 1901–2000 CE mean; Buntgen et al., 2011); (b) Northern Hemisphere air temperature anomalies (against the 1961–1990 CE mean; D'Arrigo et al., 2006), (c) tree ring width index from Albania, southeastern Europe (Seim et al., 2012); (d) summer temperature anomalies in Romania (against the 1961–1990 CE mean; Popa and Kern, 2009). Blue and yellow shaded areas indicate cold and warm periods, respectively.

Download

The FV δ¹⁸O_ice and δ²H_ice records span the 850–2016 CE period (Fig. 6). Decadal- to multidecadal-scale oscillations occur over the entire record, but no discernable long-term trend has been identified. Several periods of excursions towards low δ¹⁸O values (defined as δ¹⁸O_ice below the long-term average), indicating low summer temperatures have been observed at 875–930 CE, 1050–1080 CE, 1260–1330 CE, 1430–1480 CE, 1520–1550 CE, 1710–1750 CE, 1820–1870 CE and 1880–1930 CE. Significant maxima in the FV δ¹⁸O_ice record occurred during 850–870 CE, 1000–1050 CE, 1080–1260 CE, 1350–1390 CE, 1480–1520 CE, 1625–1710 CE and 1950–1970 CE (Fig. 6). Both the δ¹⁸O_ice and δ²H_ice records display a remarkable similarity throughout the entire period; thus, we have relied solely on the δ¹⁸O record in our discussion. We have compared the FV δ¹⁸O_ice with a tree ring width-based reconstruction of summer (JJA) temperature anomalies from the Eastern Carpathian Mountains (Popa and Kern, 2009). The highest similarities between the FV ice core and summer temperature records were found for the cold periods during 1260–1330 CE, 1430–1480 CE, 1520–1550 CE and 1820–1870 CE and the warm periods during 1080–1260 CE, 1625–1710 CE and 1950–1970 CE (Fig. 6). Given the very different nature of the two archives (trees vs. cave ice), of the proxies (tree ring width and δ¹⁸O) and of the chronologies (annual tree ring counting vs. ¹⁴C dating with a ±30-year error), the two records agree remarkably well, further supporting the hypothesis that δ¹⁸O and δ²H values in the FV ice core register both summer air temperature variability during the past ca. 1000 years in east-central Europe and, on a broader spatial scale, the variability of the AMO.

Figure 7Temporal variability of the FV δ¹⁸O (blue), the reconstructed AMO index (Wang et al., 2017), and the ¹⁴C measurement uncertainty between 850 and 2000 CE. Shading indicates the offset (in years) between the FV δ¹⁸O and AMO index values: green denotes less than 20 years, yellow denotes between 20 and 50 years and orange denotes above 50 years.

Download

Similar to the 1850–2016 CE interval described above, the stable isotope record closely mirrors the AMO variability over the entire studied interval (Fig. 7). The relationship between the FV δ¹⁸O_ice and AMO records (Wang et al., 2017) is strongest between 1125 and 1525 CE and 1750 and 2016 CE, with decadal-scale variability in the two records being synchronous. Peak-to-peak matching of the FV δ¹⁸O_ice and AMO records shows that the two records correlate well within the dating uncertainty (±30 years) during the past 1000 years. However, peak-to-peak matching indicates a difference between the two records of up to ∼ 50 years between ∼ 1600 and 1750 CE (Fig. 7), which is likely the result of high (> 50 years) uncertainty in the ice core chronology between 1525 and 1750 CE (Fig. 7). The combined correlations between (1) FV δ¹⁸O_ice and instrumental summer temperature reconstruction over the past 150 years, (2) the instrumental (this study) and proxy-based (Nagavciuc et al., 2019a) summer air temperatures and AMO variability, and (3)  the FV δ¹⁸O_ice and AMO records over the past 1000 years suggest that changes in the North Atlantic are transferred, likely via atmospheric processes, towards the wider Northern Hemisphere, resulting in hemisphere-wide climatic responses to these changes.

The FV δ¹⁸O_ice record is in agreement with other summer temperature reconstructions (e.g., Buntgen et al., 2011) at both regional and hemispheric scales (Fig. 6). Further, regional summer temperature (e.g., Popa and Kern, 2009) and summer temperature-sensitive drought (Seim et al., 2012) reconstructions show warm peaks around 1320 CE, 1420 CE, 1560 CE and 1780 CE and cold ones at around 1260 CE, 1450 CE and 1820 CE, similar to reconstructions and models at the global level (Neukom et al., 2019) and the FV temperature reconstruction (this study). Contrary to the summer season temperature reconstructions, a late-autumn through early-winter season temperature reconstruction from the nearby Scărișoara Ice Cave (Perșoiu et al., 2017) shows that the MWP was rather warm and also wet (Feurdean et al., 2011), whereas the LIA was cold and likely dry, with erratically distributed precipitation. Together, these data suggest a complex picture of climate variability in the wider Carpathian region, with much of the yearly temperature variability during the past 1000 years being attributed to the influence of winter conditions.