Research ArticlePALEOCLIMATE

Central Europe temperature constrained by speleothem fluid inclusion water isotopes over the past 14,000 years

See allHide authors and affiliations

Science Advances  05 Jun 2019:
Vol. 5, no. 6, eaav3809
DOI: 10.1126/sciadv.aav3809


The reasons for the early Holocene temperature discrepancy between northern hemispheric model simulations and paleoclimate reconstructions—known as the Holocene temperature conundrum—remain unclear. Using hydrogen isotopes of fluid inclusion water extracted from stalagmites from the Milandre Cave in Switzerland, we established a mid-latitude European mean annual temperature reconstruction for the past 14,000 years. Our Milandre Cave fluid inclusion temperature record (MC-FIT) resembles Greenland and Mediterranean sea surface temperature trends but differs from recent reconstructions obtained from biogenic proxies and climate models. The water isotopes are further synchronized with tropical precipitation records, stressing the Northern Hemisphere signature. Our results support the existence of a European Holocene Thermal Maximum and data-model temperature discrepancies. Moreover, data-data comparison reveals a significant latitudinal temperature gradient within Europe. Last, the MC-FIT record suggests that seasonal biases in the proxies are not the primary cause of the Holocene temperature conundrum.


The long-term evolution of hemispheric and global temperatures during the Holocene remains a subject of controversial discussion [for example, (14)], as proxy and model data show diverging temperature trends between the different reconstructions. Global and hemispheric temperature reconstructions inferred from pollen (4) and transient climate simulations (3) challenge previous assumptions of an early to middle Holocene Thermal Maximum [HTM; ~10 to 6 thousand years (ka) BP (before present; BP = 1950)] in the Northern Hemisphere and a subsequent cooling trend toward the present (2), also known as the Holocene temperature conundrum (3). Climate models and new pollen-inferred temperature reconstructions show evidence for lower temperatures during the early Holocene and a continuous warming until approximately 2 ka BP in the Northern Hemisphere (4), a temperature evolution that is supported by qualitative temperature records from Eurasia (5). The ultimate causes for the Holocene temperature conundrum remain ambiguous, and current explanations range from seasonal biases in marine records (3, 4), underrepresentation of pollen records in earlier temperature reconstructions (4), strong influence of records from the North Atlantic, regional cooling effects related to the remnant ice sheets during the early Holocene (3, 5) to incomplete forcing or insufficiently sensitive feedbacks of climate model simulations (6, 7). The first global Holocene temperature reconstruction, for instance, is a stack of 73 temperature records and largely based (~80%) on marine proxy records (alkenones, Mg/Ca ratios, TEX-86) (2). A large proportion of these records are from the North Atlantic realm, leading to the suggestion that the “North Atlantic is the chief source of the temperature conundrum” (8). Note that the vast majority of quantitative Holocene temperature reconstructions are either based on marine or on terrestrial biological proxies (2, 9), which, however, are seasonally biased toward the warm season (9). It is therefore necessary to develop more independent and nonbiogenic quantitative temperature reconstructions without major seasonal biases (1012) to re-examine the global and northern hemispheric temperature evolution and better constrain future warming based on warm-like periods (13).

Quantitative paleotemperature reconstructions can be obtained from speleothems by analyzing their fluid inclusions (12), which are natural repositories of rainwater falling above the cave at the time when the fluid inclusions were formed. Recent analytical advances allow to extract and measure hydrogen (δDfi) and oxygen (δ18Ofi) isotopes in speleothem fluid inclusion water at high precision, enabling us to calculate paleotemperatures (12, 14, 15). Here, we provide a quantitative decadal to multidecadal mean annual temperature reconstruction from two uranium-series dated (230Th hereinafter) stalagmites (samples M6 and M8) from the Milandre Cave [47.49°N/7.02°E, 373 to 508 m above sea level (a.s.l.)] in Northwestern Switzerland (Fig. 1). Since the beginning of temperature monitoring in 1864, Switzerland has experienced a 1.8°C warming, whereas global temperature has increased only by ~0.85°C, underlining the study site sensitivity to global warming. The chronologies of both stalagmites are based on 53 230Th dates covering the past 14,500 years continuously (figs. S1 and S2). A total of ~3340 calcite oxygen isotope (δ18Oc) and 296 δDfi and δ18Ofi analyses were performed on stalagmites M6 and M8 (Fig. 2, B, E, and F).

Fig. 1 Location and spatial correlation maps.

(A) Records discussed in the main text (stars) and the Supplementary Materials (dots) are shown together with the Laurentide and Fennoscandian ice sheet (LIS and FIS, respectively) extents at ~10.5 ka BP (30, 35). Pollen-gridded temperature differences (0.5 ka BP minus 10.5 ka BP; red is warming and blue is cooling) is also given (4). (1) Milandre Cave (red star; this study), (2) Bunker Cave, (3) Mediterranean site MD95-2043, (4) Kinderlinskaya Cave, (5) Greenland Ice Sheet Project 2 (GISP2), (6) North Greenland Ice Core Project (NGRIP), (7) Qunf Cave, (8) Cariaco Basin, (9) Gerzensee Lake, (10) Ammersee Lake, (11) Chauvet Cave, (12) Gemini Lake, and (13) Atlantic site MD99-2275 records. (B) Modern spatial correlation between Basel instrumental temperatures for the interval 1901–2014 and the CRU TS4.01 temperature dataset for Europe (see fig. S3). The red color indicates a high correlation between the temperature variations at the study site and the rest of Europe.

Fig. 2 New speleothem water and calcite isotopes and inferred paleotemperatures over the past 14,000 years for Milandre Cave compared with mid- and high-latitude records.

(A) 230Th ages and 2σ errors for M6 and M8 stalagmites. (B) Stalagmite fluid inclusion δDfi records for M6 (black), M8 (red dots), and M8 replicate (dark red circles). (C) New composite (M6 and M8 stalagmites) temperature record (MC-FIT) and uncertainty. (D) Milandre noble gas temperatures (24). (E) Stalagmite fluid inclusion δ18Ofi records for M6 (black), M8 (red dots), and M8 replicate (dark red circles). (F) Stalagmite calcite δ18Oc for M6. (G) NGRIP ice core δ18O (21). (H) Lakes Ammersee (18) and (I) Gerzensee (19) δ18O records. The blue arrow corresponds to monitored cave temperature in 2012–2013 (15). B-A, Bølling-Allerød; YD, Younger Dryas; NGT, noble gas temperature; VPDB, Vienna Pee Dee Belemnite; VSMOW, Vienna Standard Mean Ocean Water.

To reduce uncertainties surrounding the interpretation of calcite and fluid inclusion stable isotope measurements, a comprehensive 2-year monitoring program (2013–2014) was conducted in and above the Milandre Cave (15). Present cave air temperatures at the sampling site of stalagmites M6 and M8 are constant and average 9.8° ± 0.2°C that is almost identical to local mean annual air temperature (MAAT) of ~9.6°C above the cave (15). In addition, precipitation events (n = 332) were sampled between 2012 and 2014 to determine the main climatic controls on δDp and δ18Op in local precipitation and to interpret fluctuations in stalagmite δ18Oc, δDfi, and δ18Ofi values from the Milandre Cave. Temperature changes at the Milandre Cave are representative for central Europe, both with respect to temperature and stable isotopes in precipitation (Fig. 1B and fig. S3) (15).

For calculating paleotemperatures, only Milandre stalagmite δDfi values were used for the following reasons: First, δ18Oc values in speleothems from mid-latitudes are influenced by multiple climatic (e.g., surface air and cave air temperatures and evaporation) and nonclimatic (e.g., kinetic isotope fractionation within the cave) processes and are therefore difficult to interpret (16, 17). Second, postdepositional processes can alter the original δ18Ofi in fluid inclusion water and thus limit the use of δ18Ofi for paleotemperature calculations (17). Third, δDfi is more suited for calculating temperature as it is not affected by isotopic fractionation during calcite precipitation and remains unaltered as there is no hydrogen source once the water is entrapped in the calcite matrix. Furthermore, stalagmite δDfi values are intrinsically linked to drip water δD (δDdw) and precipitation δDp, respectively, where δDp values above the Milandre Cave are strongly governed by local temperature (15).

The Milandre Cave fluid inclusion temperature (MC-FIT) record is based on δDfi values by using a linear δDfi/T transfer function (TF) of 3.84 per mil (‰)/°C (see Materials and Methods for further details) and reflects the MAAT. Before calibration, the Milandre Cave δDfi record was corrected for changes in global ice volume. The errors in the MC-FIT reconstruction are nonparametric and vary between ±0.3° and ± 0.5°C. The MC-FIT record may be slightly biased toward the cold season based on the current water balance, as the epikarst above the Milandre Cave is predominantly recharged during autumn, winter, and spring (15).

New fluid inclusion water isotope records and inferred MC-FIT

The Milandre Cave δ18Oc, δDfi, and δ18Ofi records show strong evidence for the Bølling-Allerød (B-A), Younger Dryas (YD), and onset of the Holocene (Fig. 2 and fig. S4), and both their timing and amplitude are in excellent agreement with other central European δ18O records such as Ammersee (18) and Gerzensee lakes (19) (Fig. 2, H and I, and fig. S5). Furthermore, the 8.2-ka cold event is also marked and characterized by negative shifts in δ18Oc, δDfi, and δ18Ofi (Fig. 2, B, E, and F). After the end of the 8.2-ka cold event, a distinct long-term trend toward more negative isotope values is evident in the Milandre Cave δDfi and δ18Ofi time series but is absent in the δ18Oc record (Fig. 3, E and I). This suggests that δ18Oc is potentially influenced by kinetic fractionation within the cave during calcite precipitation (16, 17). These divergent long-term isotope trends between δ18Oc and reconstructed δDdw records were also observed in Holocene stalagmites from the Bunker Cave, Germany (Fig. 3J) [(20) and references therein]. They found that, despite no δ18Oc trend, the parent drip water feeding the stalagmite should have been enriched by ~0.5‰ in δ18Odw (20), in line with the ~4‰ measured on Milandre δDfi. (Fig. 3, E and J). The hydrogen isotope (δDfi) values are therefore the most suitable temperature proxy as they witness surface air temperature (15), as further discussed in the following. This assumption is supported by the close resemblance between the Milandre Cave δDfi and Greenland ice core δ18O records (21) (Fig. 2, B and G), and the comparison also reveals a direct spatial link between water isotope records in high- and mid-latitude areas in the North Atlantic realm.

Fig. 3 Northern Hemisphere key records and trends for the Holocene.

Temperature records from (A) Greenland GISP2 ice core (10), (B) pollen from high- to mid-latitude areas with Scandinavia (34) (black), 53°N to 62°N/18°E to 27°E area (green) and 43°N to 49°N/5°E to 13°E area (brown) (4), (C) Germany (NGT) (20), and (D) Switzerland speleothems (this study). Water isotope records from (E) Switzerland speleothems (this study) and (F) Greenland NGRIP ice core (21). Precipitation records from (G) Oman speleothem (31) and (H) Venezuela marine sediments (41). Calcite oxygen isotope records from (I) Switzerland (this study) and (J) Germany (27) (yellow) speleothems. On (J), calculated drip water values (δ18Odw) (20) are also shown (gray) and show a similar behavior as Milandre Cave water isotopes. Linear fits over the past 10,000 or 8000 years.

The MC-FIT record reveals that MAATs during the B-A varied around 4.9° ± 0.3°C and then dropped by more than 4.3° ± 0.8°C during the YD when MC-FIT was close or even below 0°C (Fig. 2C). The ~4°C drop in MC-FIT record corresponds with other temperature estimates from Europe [e.g., (9, 22, 23)]. The termination of the YD cold phase and onset of the early Holocene is marked by a rapid increase in MC-FIT by over 5.4° ± 0.7°C (Fig. 2C and fig. S6) within a few decades. The MC-FIT shifts across the B-A, YD, and early Holocene are similar with independent noble gas temperature (NGT) estimates from the Milandre Cave stalagmites (24) and therefore support the validity of the δDfi temperature TF used for the development of the MC-FITrecord (Fig. 2, C and D, fig. S6B, and Materials and Methods). Furthermore, the good comparison between MC-FIT and Milandre Cave NGT estimates indicates that changes in the seasonality of rainfall and storm tracks did not significantly bias the temperature signal in δDp during the B-A, YD, and Holocene transitions. For the same time interval, we also compared the MC-FIT record with another independent technique based on the temperature-dependent oxygen isotope fractionation between calcite and water (25). These temperature estimates are in very good agreement with the MC-FIT record (fig. S6B).

The early Holocene (fig. S7) is characterized by a continuous rise of 1.34° ± 0.16°C ka−1 in MC-FIT until approximately 9.6 ka BP. The timing and rate of warming is very similar to the increasing sea surface temperature (SST) of 1.32° ± 0.17°C ka−1 in the Western Mediterranean Sea (23) (11,620 to 9590 years BP) but significantly lower compared to the temperature increase of 3.73° ± 0.01°C ka−1 in Greenland (11,620 to 9550 years BP) (10). The lowest Holocene MC-FIT of 5.4° ± 0.4°C occurred at 11,370 ± 150 years BP and corresponded within age uncertainties to the so-called Preboreal Oscillation (also termed as the 11.4-ka event) in Greenland ice cores (11,270 ± 30 years BP) based on the new ice core chronology (21). On the MC-FIT time series, the HTM in central Europe lasted from ~9800 to 6400 years BP when MC-FIT was approximately 1.4°C above the preindustrial level (preindustrial level refers to MAAT in Basel between 1755 and 1850). However, the MC-FIT record reveals that the central European HTM was interrupted by a cold episode, lasting from 8800 to 7800 years BP. The 8.2-ka cold event is superimposed on this multicentennial colder period and characterized by an abrupt drop in MC-FIT of 1.2° ± 0.5°C. This temperature pattern is consistent with numerous temperature (and precipitation) records from across the Northern Hemisphere (18, 26, 27). After the 8.2-ka event, the mean rate of continental long-term gradual cooling (~7800 to 1900 years BP) is −0.12° ± 0.01°C ka−1 years with internal centennial- to millennial-scale variations in the order of 0.7°C. The cooling trend for Mediterranean SST is similar (−0.13° ± 0.03°C ka−1, ΔT of ~0.77°C, 7820 to 1900 years BP), whereas it is −0.18° ± 0.01°C ka−1 in Greenland (ΔT of ~1.06°C, 7800 to 1900 years BP). Moreover, the MC-FIT cooling between ~6400 and ~1900 years BP (ΔT of ~0.70°C) is similar to the mean Northern Hemisphere cooling based on geological records (2). The past 2500 years are particularly well documented in the Milandre record with alternating cold and warm phases (fig. S8) and display, namely, a marked long-term 0.4°C “Dark ages” cooling trend between 1400 and 1090 years BP and an abrupt Little Ice Age cooling of 1.5° ± 0.5°C starting at 690 years BP, with timing and amplitude similar to those observed in summer temperatures based on the ice-cap growth from Canada-Arctic and Iceland (28). Last, for the interval of 450 years BP to present, a good agreement is observed between the MC-FIT record and reconstructed absolute mean annual temperature for the grid representing the Milandre Cave site (29), including the instrumental part (fig. S9).

Climate forcing related to temperature variations and Holocene cooling trend

After the initial early Holocene continental warming, a long-term cooling trend starts, which we attribute to decreasing boreal summer insolation forcing (Fig. 4A) (3). The influence of the high-latitude ice sheets is perceptible in mid-latitudes of Europe, with a seesaw warming until 9.6 ka BP, coinciding not only with the complete disappearance of the Fennoscandia ice sheet (FIS) in Northern Europe at 9.7 ka BP (30) but also with the end of increasing monsoon precipitation (31). Between 9.8 and 6.4 ka BP (HTM), central Europe is relatively warm and stable and possibly less affected by the Laurentide ice sheet (LIS) that acts mostly on northern Eurasia (32), with abrupt cooling events resulting from the LIS discharge in North America well imprinted in the events at 11.4, 9.3, and especially 8.2 ka BP. On the MC-FIT records, there is a marked decrease in between 8.8 and 7.8 ka BP and during the 8.2-ka cold event caused by a weakening of the Atlantic meridional ocean circulation and reduction in northward heat transport caused by increased flux of meltwater from the LIS into the North Atlantic (26). A longer cold interval from ~9 to ~8 ka BP resulting from a rapid increase of the sea level rise due to LIS meltwater input (33) induces continental cooling (and reverse) until the influence of the LIS disappeared at around 7.8 ka BP, the latter coinciding again with more stable continental Europe temperatures between ~8 and ~6 ka BP (Fig. 4). It shows how continental European temperatures were influenced by ice sheet discharge and meltwater during the early and middle Holocene on a multidecadal to multicentennial timescale. Once the FIS disappeared, the LIS continues to act on high-latitude records [Scandinavia (34), Eurasia (5), and North America (4)], making the region colder than usual, while mid-latitudes were less affected and highlight the ice sheet’s influence on the occurrence of a Northern Hemisphere HTM (Figs. 1 and 4). At around 5.5 ka BP and once the LIS disappeared (35), all temperatures—except of the Russian speleothem—follow the same cooling trend until 3 to 2 ka BP, regardless of whether they are of biological or physical nature, whereas climate models show inconsistency with the data (3). Furthermore, for the past 2.5 ka, our record suggests that, for continental mid-latitudes, major volcanic eruptions influenced decadal temperature with MC-FIT minima, coinciding most of the time with major eruptions in Central Europe (fig. S8). This observation is related to what was argued on centennial- to millennial-scale Holocene temperature variability in high latitudes (10). The synchronicity with the new volcanic eruption chronology (36) documents that for the past 600 years, major eruptions coincide with MC-FIT minima with a coherence of ±7 years. Moreover, the mean MC-FIT cooling anomaly of around 0.4°C is in agreement with a modeled surface temperature decrease of 1.5°C over less than 5 years to volcanic eruptions (37), considering our sample resolution of 10 to 20 years.

Fig. 4 Holocene temperature records and forcing.

Similar temperature reconstructions between the Greenland, Mediterranean Sea, and central Europe areas (A) that show a different early to middle Holocene trend compared to high-latitude proxies and model simulations (B). (A) Switzerland MC-FIT (black; this study), Mediterranean SST (23) (pink) (both smoothed with a seven-value running mean) and Greenland GISP2 temperature (10) (blue; 401 years running mean) are shown together with July insolation at 30°N (51)(gray) and LIS sea level rise rate (33) (red). (B) European pollen record (4) (green) and Eurasia δ18Oc (as winter temperature indicator) (5) (yellow). Also shown are Switzerland MC-FIT (black), regional CCSM3 (3) (orange), and mid-latitude band (30° to 60°N) LOVECLIM (10) (blue) transient simulations of mean annual temperatures. Anomalies are shifted for comparison.

The MC-FIT record is characterized by an HTM between 9.8 and 6.4 ka BP, although temperatures begin to decline shortly after 7.8 ka BP. The long-term cooling trend toward the present is a distinct feature of the MC-FIT record (Fig. 3D) and supported by several other temperature records from the North Atlantic realm. The Bunker Cave NGT record from Germany (20)—albeit much lower resolved—shows strong evidence for a long-term decline in MAATs during the middle and late Holocene (Fig. 3C). Moreover, the MC-FIT evolution is regionally coherent with changes in the size of alpine glaciers that showed not only a smaller size than in the late 20th century but also no significant glacier expansion between 10.5 and 3.3 ka BP, except during rare brief intervals (38). The MC-FIT cooling trend conforms with Northern Hemisphere and European temperature records such as (i) the Greenland Summit (GISP2) temperatures based on combined nitrogen and argon isotopes (10), (ii) the Mediterranean SST of the south Iberian Margin based on alkenones (23), and (iii) reconstructed long-term trends of multiproxy mean temperature [e.g., (2)] (Figs. 2 to 4 and fig. S7), implying a coherent temperature response in mid-latitude Europe.

On a more global scale, in the oceanic domain, an integrated physically based mean ocean temperature using noble gas isotopes shows a slight Holocene cooling trend (11). In the Pacific Ocean, intermediate-depth water temperatures inferred from Mg/Ca ratios in benthic foraminifera are linked to the high-latitude source waters (likely related to winter temperatures) and display high values during the HTM and a continuous cooling afterward (39). The North Atlantic shows a similar pattern (40). Furthermore, the similarity with precipitation proxy records from South America (41) and the Arabian Peninsula (31)—itself linked to the Asian monsoon (42)—highlights together with Pacific and northern Atlantic temperature trends the interdependencies of temperature and precipitation at the Northern Hemisphere scale (Fig. 3 and fig. S7) and stresses the hemispheric signal recorded in the Milandre Cave.

The MC-FIT record and the Holocene temperature conundrum

The Milandre Cave temperature reconstruction from the Atlantic realm adds new physical and nonbiogenic evidence for the occurrence of a HTM in Europe and consequently for the existence of a Holocene temperature conundrum (Fig. 4). The temperature conundrum was raised as a data-model discrepancy between proxy records (2) and model simulations (3) that was recently extended to a data-data mismatch with the North America and Europe pollen-based (4) and Eurasia speleothem-based (5) temperature records. The correspondence between the pollen record and model simulations challenges the existence of a temperature conundrum and leaves open the question of seasonal biases (4). Our MC-FIT record is in good agreement with the multiproxy record, which is essentially based on marine records (2), but contrasts with global annual mean pollen-based temperature reconstructions (4) and model simulations (3, 6). Using these (dis-)similarities, we address hereafter both data-data and data-model discrepancies essentially in the light of our MC-FIT record.

Data-data comparison

Both the MC-FIT and pollen-based trends—used as MAAT estimates—are in agreement with regional gridded pollen time series (43°N to 49°N and 5°E to 13°E) (4), showing that the HTM is followed by a cooling trend from 7.8 ka BP onward but with a two to three times smaller amplitude of cooling for MC-FIT (fig. S7A). The MC-FIT record, however, is in contradiction with the general trend inferred for European pollen data (4) (Fig. 4B). The records indicate strong contrasting MAAT trends with warming in northeastern Europe and Eurasia compared to cooling in central-western Europe until 5.5 ka BP for reasons that still need to be investigated (Fig. 4). This statement reveals the spatial complexity of temperature patterns in continental Europe (high versus mid-latitude) and therefore would question the significance of overall Europe (and North America) temperature reconstructions, especially regarding the overrepresentation of pollen records from northeastern Europe that would likely be influenced by the ice sheets (Fig. 1). A process that may be involved in this high-latitude warming trend is the lowering of summer temperatures and a lagged timing of the HTM up to 6 ka BP due to ice sheet forcing (1). The decrease in summer temperatures may then affect the vegetation evolution in boreal region close to the ice sheets.

Furthermore, the fact that the MC-FIT long-term cooling trend is in good agreement with the regional speleothem-based temperature from Germany (20), the Northern Hemisphere multiproxy reconstruction (2), the reconstruction for Greenland (10), the benthic foraminifera from intermediate depths of the North Pacific and Antarctic (39, 40) and, to a lesser extent, with the slight cooling trend in the mean ocean temperature recorded by noble gases in ice cores from Antarctic (11) helps to shed light on the influence of seasonal biases. Since the mean ocean temperature obtained by noble gas measurements from ice cores corresponds to a pure physical (solubility), seasonality-unbiased temperature, it acts as an ideal test record. Although only a few measurements have been made so far within the Holocene, this record shows a small temperature decrease despite the influence of the Southern Hemisphere oceans. A similar trend is observed for the Greenland temperature reconstruction (Figs. 3 and 4) (10), which is also physically constrained and believed to be seasonally unbiased. It shows a two to three times stronger cooling compared to MC-FIT, which is reasonable due to polar amplification. The SST reconstructions and marine data–dominated reconstructions (2) could potentially exhibit strong seasonal biases, a possibility supported by the fact that mean alkenones temperatures are higher than Mg/Ca temperatures [see database in (2)].

Data-model comparison

Recent comparisions of Holocene temperature simulations for the extratropical Northern Hemisphere performed by four different models (6) have shown that the trend behavior for summer and the MAATs for Greenland, northern Canada, and northeastern and northwestern Europe is consistent between models under investigation. As discussed above, the MC-FIT, the multiproxy (2), and other—namely, physical—proxy reconstructions are in contradiction with all MAAT simulations for the early to middle Holocene but do agree with the summer model output. When looking at regional simulation for northwestern Europe (6), the winter simulations exhibit variabilities from strong warming [~+4°C from 9 to 6 ka BP, Community Climate System Model version 3 (CCSM3) simulation] to cooling trends [~−1 to −2°C, FAst Met Office/UK Universities Simulator (FAMOUS) simulation], the latter being similar to our MC-FIT reconstruction.

As stated earlier, the MC-FIT record may only be biased toward the cold season due to a lack of precipitation infiltrated in the epikarst in summer, while the biological proxies used in the multiproxy global temperature reconstruction (2) are thought to be biased toward the summer season (3, 4). Confirmation of this statement would mean that both winter-biased and summer-biased proxies express the same pattern, suggesting in turn that seasonal bias in the biological records (which include North Atlantic SST records) is not the source of the Holocene temperature conundrum. Moreover, a summer-season bias for the multiproxy reconstruction (2) would also mean that, given the similarities between the different reconstructions (including regional pollen data, as shown in fig. S7A), the other supraregional to global proxies (global ocean temperature composite, ocean temperatures at intermediate depths, mean ocean temperatures, and Greenland temperatures) were also summer biased, which seems unlikely, especially for the physically constrained proxy. The contradiction between this array of global cooling evidence with warming in the high-latitude Eurasian, European, and American pollen and speleothem records (4, 5) sheds light on the influence of the LIS and FIS that would act on high-latitude continental proxies, as modeled for summer temperatures (1), while mid to low latitudes are less affected. This observation associated to the similarity between MC-FIT and pollen trend for the study site (fig. S7A) suggests a strong difference in temperature between northern and central-western Europe, which may hamper a global mean reconstruction valid for every grid cell of the Northern Hemisphere, which may also be an issue for models and composites of hemispheric reconstructions.


The MC-FIT record from Switzerland adds new nonbiogenic evidence for the occurrence of an HTM in central Europe and therefore supports the existence of the so-called Holocene temperature conundrum. Our stalagmite data suggest that seasonal biases are unlikely to be the primary cause of the temperature conundrum. The proxy-data records document a significant temperature gradient between central and northeastern Europe, which may affect global temperature reconstruction through data aggregation. Therefore, further modeling and experimental investigations are required to address the observed Holocene temperature gradients across Europe and data-model disagreement, in particular, the combined influence of biases (spatial, time lag, and seasonal) induced by ice sheets in high-latitude temperature reconstructions and climate model simulations.


234U/230Th dating and samples’ age

For uranium-series dating (230Th hereinafter), approximately 170 to 200 mg of calcite powder was drilled along discrete growth horizons (fig. S1). The 230Th ages were determined on a multicollector inductively coupled plasma mass spectrometer (Thermo-Finnigan Neptune) at the Department of Earth Sciences, University of Minnesota. Details of the methods, including standards used for mass fractionation and yield correction, are found elsewhere (43). The chronology of the stalagmites is based on 42 (M6) and 11 (M8) 230Th ages (figs. S1 and S2 and table S1). The samples have low 232Th concentration (between 139 and 5028 parts per trillion) and low 238U concentration (between 63 and 430 parts per billion). All ages are quoted as years before 1950 (BP = 1950).

Age model

In both samples, the measured ages are in chronological order within their 2σ uncertainty except one age reversal at 56.75 cm in sample M6 (M6_32), which could not be resolved by the statistical program used, and therefore was excluded. Age models were built using the COPRA algorithm (44). The age model of our record continuously covers the time interval of −62 to 14,705 ± 288 years BP (considering 50 CE to be 0 years BP). More precisely, the calcite of the sample M6 was deposited between −62 and 14,546 ± 155 years BP, and the 1.5- to 13.25-cm depth in sample M8 was deposited between 8352 ± 437 and 14,705 ± 288 years BP. On the interval of 10,985 ± 223 to 14,705 ± 288 years BP, the chronology of M8 sample is more reliable, in terms of dated subsamples and dating precision, compared with the one of the sample M6. To improve the chronology of the M6 calcite on its 135.05- to 145.18-cm depth, we applied an eye-fit “event stratigraphy” using the M8 as “master chronology.” On the basis of the assumption that climatic changes above the cave should be recorded simultaneously by each sample, distinct shifts in δ18O were tuned within the error intervals of each record using 25 fix points (table S2). To each fix point, we attributed its corresponding 2σ error given by the M8 COPRA chronology. This approach accounted both for the error within M8 chronology and for the one introduced by the temporal resolution of our calcite records (7 years for M6 and 25 years for M8). No significant growth hiatuses have been identified—macroscopic or by dating—on the measured interval. Growth rates range from 0.007 to 0.192 mm/year (overall mean, 0.053 mm/year) in sample M6 and from 0.004 to 0.073 mm/year (overall mean, 0.028 mm/year) in sample M8. Age ranges from the B-A (14,705 ± 288 to 12,715 ± 140 years BP), the YD (12,715 ± 140 to 11,615 ± 40 years BP), to the Holocene (past 11,615 years).

Stable isotope analysis

Speleothem fluid inclusion water isotopes were analyzed at the Physics Institute, University of Bern, Switzerland, following a recently developed laser cavity ring-down spectroscopy (CRDS)–based method (14), with further improvements (15). Interlaboratory comparison with isotope ratio mass spectrometry–based speleothem fluid inclusion isotope measurements showed a very good agreement for hydrogen (δDfi) values (12). Oxygen (δ18Ofi) and hydrogen (δDfi) isotope ratios are given in per mil (‰) using the standard delta notation and are reported against Vienna Standard Mean Ocean Water (VSMOW). To sum up the measurement procedure, a calcite block was placed in a copper tube connected to the line that was continuously heated to 140°C and crushed using a hydraulic press. The released sample water was instantaneously vaporized and directly flushed with nitrogen gas to a Picarro L2140-i (280 samples) or L1102-i (16 samples) and measured with the CRDS technology without any prior water treatment. We measured a total of 296 calcite blocks that were extracted at the central growth axis of stalagmites M6 and M8 (255 samples, 15 duplicates, and 4 modern for M6; 15 samples and 7 duplicates for M8) with rough dimensions: a length of ~30 mm (perpendicular to the growth axis), a thickness of ~4 to 5 mm (parallel to the growth axis), and a weight between 0.44 and 1.88 g (mean, 1.35 g). Analytical precision of measurements is typically 0.2 and 1.0‰ for δ18Ofi and δDfi, respectively, when using an L2140-i instrument and for water amounts larger than 1 μl, whereas it is 0.4‰ for δ18Ofi and 1.5‰ for δDfi when using the L1102-i instrument and for released water amounts below 1 μl. Including the natural variations and potential tortuous laminae in all replicate samples, the overall deviation for the duplication of 19 Holocene samples was 0.2‰ for δ18Ofi (0.0 to 0.7‰) and 0.5‰ for δDfi (0.0 to 1.6‰). Crushing of calcite samples released up to 7 μl of water (mean, 2.9 μl). The mean temporal resolutions achieved is ~10 years for the last 800 years BP, ~25 years between 800 and 2400 BP, ~50 years until 10,000 BP, and ~70 years after that for M6, whereas M8, which covers mostly the YD, has a lower resolution (77 to 480 years).

Speleothem calcite oxygen stable isotope (δ18Oc) samples were measured using a Finnigan Delta V Advantage mass spectrometer equipped with an automated carbonate preparation system (Gas Bench II) at the Geological Sciences Institute, University of Bern, Switzerland. Stalagmite M6 was continuously micromilled at increments of 0.2 mm (from 136.6 to 145.18 cm from the top) and 0.5 mm (0 to 136.6 cm from the top) and stalagmite M8 at increments of 0.5 mm (from 2 to 13.25 cm from the top). Calcite powder δ18Oc was measured in ~3340 powder samples in both stalagmites (average resolution, 4 years). About 150 to 200 μg of powdered sample was placed in 12-ml vacutainers, flushed with He, and reacted with five drops of 99% P5O5 acid at 70°C for 80 min. Results are reported relative to the international Vienna Peedee Belemnite (VPDB) standard. The analytical reproducibility based on repeated measurement of the internal standard for δ18O is lower than 0.07‰ VPDB (1σ) (45).

Isotope temperature conversion

The composite paleotemperature record (MC-FIT) from M6 and M8 stalagmites was reconstructed for 267 samples based on the robust regional water isotope/temperature relationship (15, 46, 47). Before temperature reconstruction, δD was corrected for ice volume effect with a gradient obtained for δ18O (48), converted by a factor of eight to δD of −0.064‰ per meter of sea level rise (49). Furthermore, we used δ18O in modern precipitation and corresponding temperature time series for the GNIP (Global Network of Isotopes in Precipitation) station in Basel (1986–2017) that is most representative for the cave site and the GNIP station in Bern, which offers the longest interval (1971–2017) to estimate short (monthly), seasonal, and long-term δ18O/ΔT relationships. Seasonal correlations for an individual station such as Bern show the highest degree of covariations with temperature (95%), whereas monthly correlations suffer from strong and long-term correlations from low variability (17%). When pooling many stations, as done for Europe, the long-term correlation increases significantly as documented by (46).

We used a linear δD/T TF to calculate paleotemperatures (table S3). Because of the abovementioned consideration in the selection of TF, we reconstructed temperature based on monthly values with a slope of 0.36‰/°C (Basel, 0.35 ‰/°C, R2 = 0.57; Bern, 0.37 ‰/°C, R2 = 0.60) and 0.60‰/°C (R2 = ~0.5) for a spatial relation of long-term trends (46), as well as on the mean TF of 0.48‰/°C (fig. S6A), the latter being also the mean value for Basel based on MAATs and isotope values (0.47‰/°C, R2 = 0.17), and the time (number of years)–weighted mean seasonal TF of intervals 1971–1990 (46) and 1994–2002 (47) at the Bern station with a value of 0.48‰/°C. The same value was obtained on the basis of a forward modeling method from multiple sites of borehole temperatures in Greenland (50). These δ18O values were further transposed in δD values based on the equilibrium fractionation factor of eight. We discuss temperature based on the mean relationship of 0.48‰/°C (i.e., for δD 3.84‰/°C) anchored at the most recent sample (mean depth, 0.9 cm; i.e., 1928–2011 CE) set to 8.3°C, corresponding to the 1928–2011 MAAT for Basel MeteoSwiss station corrected for the elevation of Fahy, that is, 596 m a.s.l. Temperature reconstructions were not possible for the interval of ~127 to 62 years BP (six samples) due to the low amount of water extracted from those samples with correspondingly low confidence in their obtained isotope values. For the YD interval, due to the tortuous shape of laminae in samples, we based the temperature reconstruction only on the three samples from M8 for which we managed to follow the laminae contour most precisely. Other measurements of this time interval likely suffered from mixing of inclusion water from side laminae that smooth the isotope values, as seen on Fig. 2B. We intentionally anchored the MC-FIT at the MAAT because of (i) the lack of information regarding temporal changes in the seasonal distribution of water infiltration to the karst and through it and (ii) seasonal water balance values with rather high uncertainties combined with only 3-year records of parallel precipitation and drip water isotope measurements at our sampling site (15). Yet, the temperature could also be anchored at a seasonally biased temperature excluding summer based on the seasonally biased water balance values (15). This would just result in constant offsets of the presented mean annual temperatures (see also the “Error estimation and potential biases” section).

Error estimation and potential biases

We followed a similar error analysis approach as used recently for temperature reconstruction inferred from ice core measurements (50). Uncertainties for the MC-FIT record come from the isotope measurement error and the error resulting from the selected TF. Uncertainties that encompass hydrogen isotope measurement uncertainty of 1‰ are ±0.21°, ±0.26°, and ± 0.35°C (225 samples) for TFs of 0.60, 0.48, and 0.36‰/°C, respectively (linear regression). It is ±0.31°, ±0.39°, and ±0.52°C (42 samples) for a 1.5‰ measurement uncertainty. The domain error arising from the TF choice is nonlinear and asymmetric and displayed in Fig. 2 and fig. S6 as upper and lower boundaries (blue-shaded interval) of each single TF temperature reconstruction. Whereas the change in TFs will induce an offset in MC-FIT and result in a change of MC-FIT amplitudes, the variability remains robust and is mainly dependent on measurement uncertainty. Thus, the record moves in the uncertainty band, but the variability will hardly be affected. Errors given in the text for the interval means were calculated with the following formula: σinterval=σmeas2+σmean2n. The uncertainty of a change is given as the sum of the σinterval. The good δD interlaboratory reproducibility (12) documents that the loss of water by diffusion through the calcite or during measurements is negligible, leaving only evaporation before enclosure to alter δDfi that would tend to produce enriched isotope values and thus higher temperature estimates. Drip water δDdw is enriched compared to the correspondent value of direct precipitation by ~10‰. Yet, evidence that δDfi provides a direct record of past precipitations is given through its modern calibration and high correlation between δDfi and δ18Ofi (slope, 7.92 ± 0.07; R2 = 0.98) (15). We would further like to mention that the recent water balance (precipitation − evaporation) (15) potentially indicates slightly positive values in summer, implying reduced water infiltration and thus bias toward negative temperature values. If this would be true, then the reconstructed temperature would exhibit an offset of ~−3°C, resulting from autumn, winter, and spring (AWS) water infiltration based on instrumental temperature data. The comparison between the regional CCSM3 model simulations for the MAAT and the AWS for the past 14,000 years also results in a mean offset of ~−3°C. Nevertheless, the assumption is made that the water balance is constant over time. We also have no control on seasonality changes of infiltration and precipitation in the past, and we assume that annual infiltration is fairly constant over time. Potential changes in the large-scale atmospheric pattern, i.e., North Atlantic Oscillation or changes in the moisture source, were not taken into account.

Model results used for comparison

TraCE paleoclimate simulation data were provided by the National Center for Atmospheric Research (NCAR) that uses a coupled atmosphere-ocean general circulation model, i.e., the CCSM3 (3). The forcing and boundary conditions include orbital parameters, greenhouse gases in the atmosphere, ice sheets, paleogeography, and meltwater fluxes. Mean annual surface temperatures obtained from CCSM3 TraCE paleoclimate run were extracted for the area: longitudes 2°E to 12°E and latitudes 45°N to 52°N. Further, Northern Hemisphere mid-latitude (30°N to 60°N) average temperature was obtained using the LOVECLIM model version 1.3. Information about the forcing descriptions can be found in (10).


Supplementary material for this article is available at

Supplementary Text

Fig. S1. Stalagmites M6 and M8 with U-Th age locations.

Fig. S2. Depth-age models for M6 and M8 stalagmites.

Fig. S3. Spatial correlation maps for temperatures and isotopes.

Fig. S4. Original δ18O time series for M6 and M8 stalagmites.

Fig. S5. Timing of Milandre Cave calcite δ18Oc shifts compared with other oxygen records.

Fig. S6. Temperature reconstruction (MC-FIT) and uncertainties from speleothem fluid inclusions.

Fig. S7. Comparison of MC-FIT (black lines) with independent temperature records for the Holocene.

Fig. S8. Influence of volcanism on decadal temperature average.

Fig. S9. Comparison between reconstructed (MC-FIT), instrumental, and multiproxy temperature records.

Table S1. 230Th ages of M6 and M8 stalagmites from the Milandre Cave.

Table S2. List of fix points.

Table S3. TF impact for different time interval.

References (5258)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank P. Nyfeler, S. Brechbühl, and T. Aebi for technical issues; S. Breitenbach for inputs with age model; Y. Krüger, W. Tinner, and C. Raible for comments; and Ph. Häuselmann for fieldwork support. We thank M. Messmer for providing NCAR and T. Kobashi for LOVECLIM simulation outputs. From the NCAR community, we thank B. L. Otto-Bliesner, Z. Liu, F. He, E. C. Brady, R. Tomas, P. U. Clark, and A. Carlson. Funding: This study was part of the “STALCLIM” and “STALCLIM 2” Sinergia projects funded by the Swiss National Science Foundation (grant nos. CRSI22-132646 and CRSII2_147674 to D.F. and M.L.). H.C. acknowledges the support by the National Natural Science Foundation of China (NSFC 41888101). Author contributions: S.A., M.L., and D.F. co-designed the study. For the experimental work, S.A. performed fluid inclusion water isotope measurement, A.H. and D.F. performed oxygen calcite measurement, and H.C. and R.L.E. performed the uranium-series dating. S.A., A.H., and D.F. performed field work. S.A., M.L., A.H., and D.F. contributed to the data interpretation. S.A., M.L., and D.F contributed to writing the manuscript, with input from A.H. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Data will be made available at the PANGAEA data repository (

Stay Connected to Science Advances

Navigate This Article