Stalagmite carbon isotopes suggest deglacial increase in soil respiration in western Europe driven by temperature change

. The temperate region of western Europe under-went signiﬁcant climatic and environmental change during the last deglaciation. Much of what is known about the terrestrial ecosystem response to deglacial warming stems from pollen preserved in sediment sequences, providing information on vegetation composition. Other ecosystem processes, such as soil respiration, remain poorly constrained over past climatic transitions but are critical for understanding the global carbon cycle and its response to ongoing an-thropogenic warming. Here we show that speleothem carbon isotope ( δ 13 C spel ) records may retain information on soil respiration and allow its reconstruction over time. While this notion has been proposed in the past, our study is the ﬁrst to rigorously test it, using a combination of multi-proxy geochemical analysis ( δ 13 C, Ca isotopes, and radiocarbon) on three speleothems from the NW Iberian Peninsula and quantitative forward modelling of processes in soil, karst, and cave. Our study is the ﬁrst to quantify and remove the effects of prior calcite precipitation (PCP, using Ca isotopes) and bedrock dissolution (using the radiocarbon reservoir effect) from the δ 13 C spel signal to derive changes in respired δ 13 C. The coupling of soil gas p CO 2 and δ 13 C via a mixing line describing diffusive gas transport between an atmospheric and a respired end-member allows the modelling of changes in soil respiration in response to temperature. Using this coupling and a range of other parameters describing carbonate dissolution and cave atmospheric conditions, we generate large simulation ensembles from which the results most closely matching the measured speleothem data are selected. Our results robustly show that an increase in soil gas p CO 2 (and thus respiration) is needed to explain the observed deglacial trend in δ 13 C spel . However, the Q 10 (tem-perature sensitivity) derived from the model results is higher than current measurements, suggesting that part of the signal may be related to a change in the composition of the soil respired δ 13 C, likely from changing substrate through increasing contribution from vegetation biomass with the onset of the Holocene.

Abstract. The temperate region of western Europe underwent significant climatic and environmental change during the last deglaciation. Much of what is known about the terrestrial ecosystem response to deglacial warming stems from pollen preserved in sediment sequences, providing information on vegetation composition. Other ecosystem processes, such as soil respiration, remain poorly constrained over past climatic transitions but are critical for understanding the global carbon cycle and its response to ongoing anthropogenic warming. Here we show that speleothem carbon isotope (δ 13 C spel ) records may retain information on soil respiration and allow its reconstruction over time. While this notion has been proposed in the past, our study is the first to rigorously test it, using a combination of multi-proxy geochemical analysis (δ 13 C, Ca isotopes, and radiocarbon) on three speleothems from the NW Iberian Peninsula and quantitative forward modelling of processes in soil, karst, and cave. Our study is the first to quantify and remove the effects of prior calcite precipitation (PCP, using Ca isotopes) and bedrock dissolution (using the radiocarbon reservoir effect) from the δ 13 C spel signal to derive changes in respired δ 13 C. The coupling of soil gas pCO 2 and δ 13 C via a mixing line describing diffusive gas transport between an atmospheric and a respired end-member allows the modelling of changes in soil respiration in response to temperature. Us-ing this coupling and a range of other parameters describing carbonate dissolution and cave atmospheric conditions, we generate large simulation ensembles from which the results most closely matching the measured speleothem data are selected. Our results robustly show that an increase in soil gas pCO 2 (and thus respiration) is needed to explain the observed deglacial trend in δ 13 C spel . However, the Q 10 (temperature sensitivity) derived from the model results is higher than current measurements, suggesting that part of the signal may be related to a change in the composition of the soil respired δ 13 C, likely from changing substrate through increasing contribution from vegetation biomass with the onset of the Holocene. 1904 F. A. Lechleitner et al.: Stalagmites suggest deglacial increase in soil respiration climate change, with repercussions for the terrestrial carbon cycle and the release of greenhouse gases to the atmosphere (Clark et al., 2012). The temperate region of western Europe was particularly affected by large and latitudinally diverse environmental changes during the last deglaciation, driven by its proximity to the Scandinavian Ice Sheet and the North Atlantic (Moreno et al., 2014). Over the entire region, terrestrial paleo-climate records indicate a transition from colder to warmer climatic conditions, punctuated by millennial-scale events which closely match the Greenland ice core record (Genty et al., 2006;Moreno et al., 2014). Pollen records from western Europe reveal a general deglacial trend from grassland steppe and tundra ecosystems towards landscapes dominated by temperate forest and provide evidence for a remarkably rapid ecosystem response to temperature changes on millennial scales over the last glacial (Fletcher et al., 2010).
Speleothem carbon isotope (δ 13 C spel ) records from the temperate region of western Europe are often clearly correlated to regional temperature reconstructions during the last glacial (Genty et al., 2003) and the deglaciation (Baldini et al., 2015;Denniston et al., 2018;Genty et al., 2006;Moreno et al., 2010;Rossi et al., 2018;Verheyden et al., 2014) (Fig. 1), pointing towards a regionally coherent mechanism driving the response to the temperature increase. Early on, Genty et al. (2006Genty et al. ( , 2003 suggested that the temperature sensitivity of δ 13 C spel in western Europe was likely related to the response of vegetation and soil respiration to climate warming. Higher concentrations of respired CO 2 in the soil gas lower its δ 13 C signature, due to the increase in strongly fractionated organic carbon in the system. Speleothems can capture this change as they are fed by drip water, which equilibrates with soil gas pCO 2 before proceeding to the dissolution of carbonate bedrock. This mechanism could lead to the observed transitions from higher δ 13 C spel during colder periods to lower δ 13 C spel during warmer periods and may provide a means to quantify past changes in soil respiration, an elusive parameter in the global carbon cycle (Bond-Lamberty and Thomson, 2010). However, formal testing of this mechanism has so far not been attempted, mainly because of the numerous and complex processes that influence δ 13 C spel (Fohlmeister et al., 2020).
Speleothem carbon can originate from three sources: atmospheric CO 2 , biogenic CO 2 from autotrophic (root and rhizosphere) and heterotrophic (soil microbial) soil respiration (from here onwards jointly referred to as "soil respiration"), and the carbonate bedrock itself (Fig. 2). Recent research has additionally suggested that deep underground reservoirs of carbon ("ground air"; Mattey et al., 2016) or deeply rooted vegetation (Breecker et al., 2012) may play a significant role in the karst carbon cycle. The relative importance of these different sources on δ 13 C spel is modulated by hydroclimate and temperature. This can occur as a propagation of a biosphere response to climate change, e.g. changes in vegetation composition (Braun et al., 2019), changes in soil respiration (Genty et al., 2003), and changes in soil turnover rates (Rudzka et al., 2011). Secondly, δ 13 C spel can be modulated by changes in karst hydrology, i.e. the carbonate bedrock dissolution regime (Hendy, 1971). Thirdly, compounded changes in hydrology and cave atmospheric pCO 2 can lead to prior calcite precipitation (PCP) during carbonate precipitation (Fohlmeister et al., 2020). Altitudinal transects in caves in the European Alps have shown that changes in soil respiration, vegetation, and temperature have a traceable effect on speleothem fabrics, stable oxygen isotope ratios, and δ 13 C spel (Borsato et al., 2015). So far, it has not been possible to disentangle these effects and quantify their relative importance on δ 13 C spel records, but this quantification is a crucial step towards evaluating the potential of δ 13 C spel as a paleo-soil respiration proxy. Here, we generate a multiproxy dataset from three stalagmites from the NW Iberian Peninsula and use quantitative forward modelling to show that changes in soil respiration can explain much of the observed deglacial trend in western European δ 13 C spel . Our approach is the first to leverage differing proxy sensitivities to quantitatively model key environmental parameters, in particular soil gas pCO 2 , allowing us to estimate the total temperature sensitivity of soil respiration (Q 10 ), including the effect of changing vegetation communities.

Study site and samples
El Pindal and La Vallina caves are located ∼ 30 km apart on the coastal plain in Asturias, NW Iberian Peninsula, at 23 and 70 m a.s.l., respectively (43 • 12 N, 4 • 30 W; Fig. 1). Both caves developed in the non-dolomitic, Carboniferous limestones of the Barcaliente formation, with an overburden of 10-35 m of bedrock for El Pindal Cave and 10-20 m for the gallery in which samples were collected in La Vallina.
Current climate in the NW Iberian Peninsula is characterised by temperate maritime conditions, with clear precipitation seasonality but no summer drought (Peinado Lorca and Martínez-Parras, 1987). The region is strongly affected by North Atlantic climate conditions, in contrast to the rest of the Iberian Peninsula, where North Atlantic and Mediterranean influences persist (Moreno et al., 2010). Both caves are affected by similar climatic conditions, with ∼ 1250 mm/yr annual precipitation (Stoll et al., 2013) and maximum precipitation occurring in November (140 mm/month) (AEMET meteorological stations at Santander andOviedo, period 1973-2010;AEMET, 2020). Due to the proximity to the coast, temperature exhibits a clear but modest seasonality, with averages of 9 • C for winter months (December-March) and 20 • C for summer months (June-September) (AEMET meteorological station at Santander, period 1987Santander, period -2000AEMET, 2020). For the last deglaciation, quantitative estimates of temperature can be derived from marine records from the western and southern Iberian margins. These likely give a reasonable estimate of the deglacial temperature change in caves on the coastal plain, as the re-Figure 1. Speleothem δ 13 C records covering the last deglaciation in temperate western Europe. (a) δ 13 C spel vs. age, colour-coded by cave. Villars Cave -stalagmites Vil-stm11 (Genty et al., 2006) and Vil-car-1 (Wainer et al., 2011); Chauvet Cave -stalagmite Chau-stm6 (Genty et al., 2006); El Pindal Cave (studied here)-previously published record from stalagmite Candela (Moreno et al., 2010); La Garma Cave -stalagmite GAR-01 (Baldini et al., 2015); El Soplao Cave -stalagmite SIR-1 (Rossi et al., 2018); Père Noël Cave -stalagmite PN-95-5 (Verheyden et al., 2014); Buraca Gloriosa -stalagmite BG6LR (Denniston et al., 2018). All stalagmite data were extracted from the SISAL database, version 2 (Comas-Bru et al., 2020b, a). Shown here is the millennial-scale trend in the records, calculated using a Gaussian kernel smoother (nest package in R, Rehfeld and Kurths, 2014). Please note that the record from Villars Cave is very low resolution but no hiatus is reported between 18-14 ka. Time slices at the top of the figure are as defined for the modelling in this study. GS: Greenland Stadial, as defined in Rasmussen et al. gion's modern seasonal cycle displays a similar amplitude to sea surface temperatures (Stoll et al., 2015). Minimum average temperatures are reconstructed for Greenland Stadial 2.1a (GS-2.1a, Heinrich event 1, 18-15 ka; Rasmussen et al., 2014) and are ∼ 8 • C cooler than those of the Early Holocene (∼ 8 ka; Darfeuil et al., 2016).
Previous monitoring data from the two caves reveals seasonal variations in cave air pCO 2 driven by external temperature variations (Moreno et al., 2010;Stoll et al., 2012). Both caves are well ventilated in the cold season with close to atmospheric pCO 2 values but feature elevated CO 2 concentrations during the warm summer season (Stoll et al., 2012). The caves are covered by thin (< 1 m deep) and rocky soils, and modern vegetation is strongly impacted by Late Holocene land use change, including deforestation of native Quercus ilex (evergreen oak) for lime kilns above El Pindal Cave, and discontinuous pasture maintained by cycles of burning above both caves. At present, the vegetation above the two caves includes pasture and gorse shrub (Ulex), but in some areas above El Pindal Cave, the recent abandonment of pastures has permitted the return of patches of native Quercus ilex forest. Above La Vallina Cave, pastures are interspersed with native oak (Quercus) and planted groves of Eucalyptus, the roots of which penetrate the cave in points directly beneath the tree groves.
Candela is a calcitic stalagmite that grew ∼ 500 m inside El Pindal Cave and was not active at the time of collection (Moreno et al., 2010). Previous investigations revealed that the stalagmite grew between ∼ 25-7 ka and provide highresolution stable isotope and trace element records (Moreno et al., 2010), as well as radiocarbon ( 14 C) measurements between 15.4-8.8 ka (Rudzka et al., 2011). Growth of Candela is strongly condensed between 18-15.5 and 11-9 ka (Stoll et al., 2013). Stalagmite Laura is from El Pindal Cave, while Galia grew in La Vallina Cave. Both Laura and Galia are also composed of calcite. Previous U-Th dating on Galia revealed intermittent growth between 60 and 4 ka (Stoll et al., 2013), including a short growth phase at 26 ka, which together with the Holocene growth is sampled here. Laura grew between 16.1-14.2 ka, covering the GS-2.1a-Greenland Interstadial 1 (GI-1) interval.

Geochemical measurements
To minimise sampling bias, samples from all three stalagmites were drilled from the same locations for all geochemical analyses using either a handheld drill or a semi-automated high-precision drill. An aliquot each of the collected powder was used for U-Th dating, δ 13 C, 14 C, and Ca isotopes. In the case of Candela, where a few U-Th dates were available from previous investigations (Moreno et al., 2010), powders for the remaining proxies were drilled from the same sampling holes. Additional paired multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS) U-Th dates from all three stalagmites are detailed elsewhere (Stoll et al., 2021).
For stable carbon isotopes, an aliquot of powder was analysed on a ThermoFinnigan GasBench II carbonate preparation device at the Geological Institute, ETH Zurich, following the procedure by Breitenbach and Bernasconi (2011). Measurement runs were evaluated using an in-house standard (MS2) that has been linked to NBS19 and the external standard deviation (1σ ) for δ 13 C is smaller than 0.08 per mil (‰). Isotope values are expressed in per mil and referenced to the Vienna Pee Dee belemnite standard (VPDB).
Radiocarbon measurements were performed at the Laboratory for Ion Beam Physics, ETH Zurich, using a MI-CADAS (mini radiocarbon dating system) accelerator mass spectrometer (AMS; Synal et al., 2007) coupled to a gas ion source (GIS; Fahrni et al., 2013). Carbonate powders (∼ 1 mg) were dissolved in 85 % H 3 PO 4 and the resulting CO 2 gas was directly injected into the GIS. Quality control of the AMS measurements was ensured by measuring oxalic  . Schematic representation of the main processes modulating δ 13 C, δ 44/40 Ca, and 14 C in the soil-karst-cave system. The small inserts show the evolution of δ 13 C spel in response to changes in soil pCO 2 , DCF, and PCP. The dissolution process can be constrained using DCF, while PCP is constrained using δ 44 Ca. Soil pCO 2 affects all three proxies but can be constrained further using the coupled relationship with δ 13 C from the mixing lines. Adapted from Day et al. (2021). acid II (NIST SRM 4990C), IAEA C-2 as a carbonate standard, and IAEA C-1 as carbonate blank, and measurement precision was better than 10 ‰. We use the 14 C reservoir effect ("dead carbon fraction", DCF), which quantifies the amount of fossil carbon incorporated in the speleothems and serves as a tracer for changes in karst hydrology or mean soil carbon age (Genty et al., 2001). The DCF is calculated as the normalised difference between the atmospheric 14 C activity (F 14 C;Reimer, 2013) at the time of speleothem deposition (defined through the independent U-Th chronology) and the speleothem 14 C activity corrected for decay. Using paired U-Th and 14 C ages has the advantage of minimising uncertainty from age modelling interpolation techniques. To account for the uncertainty in matching the speleothem chronology with the atmospheric 14 C record (IntCal13 calibration curve; Reimer et al., 2013) the atmospheric record was interpolated to a yearly resolution and matched to 10 000 simulated speleothem ages for each U-Th dating point. The average and standard deviations from these ensembles were then used for the final DCF calculation and uncertainty propagation.
Samples for Ca-isotope analysis were taken from the stalagmites and from three pieces of bedrock overlying both caves. Combined bedrock and stalagmite Ca-isotope analyses allow the reconstruction of the Ca-isotopic composition of the initial growth solution and therefore of the fraction of Ca remaining in solution (f Ca ) at the point of stalagmite growth, a quantitative measure for PCP (Owen et al., 2016). Aliquots of CaCO 3 (200-650 µg) were dissolved in distilled 2 M HNO 3 . The Ca was purified using an automated Ca-Sr separation method (PrepFAST MC, Elemental Scientific, Omaha, NE, USA). This process separates Ca from Sr, Mg, and other matrix elements to avoid isobaric interferences during MC-ICP-MS. Ca-isotope ratios were analysed at the University of Oxford using a Nu Instruments MC-ICP-MS, following the method of Reynard et al. (2011). All solutions were at 10 ± 1 ppm concentration, and the samples were measured with standard-sample bracketing. Each sample was analysed a minimum of 5 times. δ 44/40 Ca was calculated using δ 44/40 Ca = δ 44/42 Ca · ((43.956 − 39.963)/(43.956 − 41.959)) (Hippler et al., 2003) and is reported normalised to NIST SRM 915a. Secondary standards HPSnew (in-house standard) and NIST-SRM-915b (purified alongside the samples) were used to determine accuracy and external precision. Measured values for our purified SRM 915b were δ 44/40 Ca = 0.71 ± 0.06 ‰ (2se, n = 12), which match values obtained by TIMS (thermal ionization mass spectrometry), δ 44/40 Ca = 0.72 ± 0.04 ‰ (2se; Heuser and Eisenhauer, 2008). Uncertainty on Ca-isotope data is quoted as the t-distributionderived 95 % confidence interval on the mean of repeat measurements calculated using either the standard deviation on all repeat measurements on each sample or the standard deviation on all secondary standard analyses, whichever is greater.

Process modelling and sensitivity analysis
Forward modelling of processes occurring in the soil, karst, and cave allow us to investigate the combination of parameters which would simultaneously simulate δ 13 C spel , δ 44/40 Ca, and DCF for each time period sampled. Using δ 44/40 Ca and DCF to quantify changes in PCP and bedrock dissolution conditions (open vs. closed system), respectively, we can remove these effects from δ 13 C spel and derive soil gas pCO 2 and δ 13 C (Fig. 2). We employ the PHREEQCbased, numerical model CaveCalc (Owen et al., 2018), a tool that enables us to evaluate and combine the effects of PCP and bedrock dissolution quantitatively and systematically. We generate large ensembles of simulations from which we then choose the solutions best fitting the measured proxy data. CaveCalc simulates the equilibration between meteoric water and soil CO 2 , the subsequent dissolution of the host carbonate rock by this solution, and the degassing of CO 2 from the solution in the cave environment that leads to the formation of speleothem carbonate. Key model inputs (Table 1) are the concentration and isotopic composition of soil CO 2 and the degree to which isotopic exchange during carbonate dissolution occurs under open/closed or intermediate conditions (gas volume relative to solution volume), which set the initial saturation state and isotopic composition of the drip water. Together with the soil gas pCO 2 , the pCO 2 of the cave environment is modelled to set the degree of oversaturation the solution will have in the cave, and it determines the amount of carbonate which can precipitate before the solution reaches equilibrium. Constraints on the model parameters are given by δ 13 C spel , δ 44/40 Ca, and DCF.
Our primary interest is evaluating constraints on soil respiration, soil gas pCO 2 , and its isotopic composition. Soil CO 2 is a mixture of carbon from respired, atmospheric, and bedrock sources, with its concentration depending mainly on temperature, water content, porosity, and soil depth (Amundson et al., 1998;Cerling et al., 1991). Global regressions find growing season soil gas pCO 2 strongly positively correlated with temperature and actual evapotranspiration (Borsato et al., 2015;Brook et al., 1983), and water balance is responsible for a steep gradient in soil pH globally (Slessarev et al., 2016). As soil gas pCO 2 is typically much higher than atmospheric pCO 2 , CO 2 diffuses from the soil along concentration gradients, and its concentration and δ 13 C value can be approximated using a mixing line between an atmospheric and a soil respired end-member using the Keeling plot approach (Amundson et al., 1998;Cerling et al., 1991;Pataki et al., 2003). Here we use this relationship to test whether changes in soil respiration can realistically explain the observed deglacial δ 13 C spel trend. The likely range of values for the soil respired end-member was constrained through monitoring of cave air pCO 2 and δ 13 C cave-air at La Vallina Cave, supplemented by measurements of local atmospheric pCO 2 and δ 13 C over 1 year. This is possible because the cave, like the soil, is defined by a two end-member mix-ing system (soil gas and atmospheric air), driven by seasonal ventilation. Monthly CO 2 measurements reveal a strong correlation between cave air pCO 2 and δ 13 C cave-air , in particular during the summer, when soil respiration is highest (Fig. 3a). Our estimation of the modern respired end-member, defined along a mixing line which includes the modern global atmospheric end-member, is −26.9 ± 0.8 ‰. This end-member may be more negative than the preindustrial end-member (which characterised the Early and mid-Holocene growth periods of the stalagmites in this study) because atmospheric δ 13 C has decreased by 2 ‰ over the last century due to anthropogenic activities (Suess effect). Assuming modern and decadal age soil carbon pools contributing dominantly to the respired end-member, the preindustrial respired end-member may have been as much as 2 ‰ heavier (∼ −25 ‰). If a significant fraction of the respired pool is older, then the preindustrial respired end-member may fall between −27 ‰ and −25 ‰, but this is unlikely since actively growing stalagmites from the cave show rapid post-bomb spike decrease in 14 C.
Using this respired end-member, we define mixing lines for relevant periods of the late glacial (LG, before 16.5 ka), deglaciation (DEG, 16.5-11.7 ka), and Early Holocene (EH, after 11.7 ka), using atmospheric CO 2 compositions as published from ice core measurements (Table 1, Fig. 3b). The isotopic composition of the atmospheric CO 2 remains within a few tenths of a per mil of the Holocene value (Schmitt et al., 2012). Therefore, assuming a constant respired endmember, the slope of the mixing line is reduced during periods of lower atmospheric pCO 2 (Fig. 3b). The mixing line may also vary if the respired end-member changes. Although there is variation in the respired end-member both within and among biomes, the mean respired end-member for the potential biomes which may have characterised this site over the last 25 ka -temperate broadleaf, temperate conifer, and boreal -feature mean δ 13 C of respired end-members which differ by only 1 ‰ (Fig. 3b;Pataki et al., 2003). This suggests that we cannot predict a systematic change in the δ 13 C of the respired end-member with changes in the biome. Moreover, the fact that deglacial trends in δ 13 C spel across western Europe are very similar also indicates that highly localised factors that may lead to a strong change in respired δ 13 C without a biome change are unlikely. Consequently, we address the potential for variation in the respired end-member by completing a sensitivity analysis of mixing lines which encompass 3 ‰ heavier and lighter respired end-members (Supplement Table S1). For the modelling, we use a maximum soil CO 2 concentration of 8000 ppmv, consistent with predictions based on modern climatology and global regressions of pCO 2 from climatic factors (e.g. Borsato et al., 2015;Brook et al., 1983). Current vegetation density and soil gas pCO 2 may underestimate Holocene conditions that preceded significant land use alteration, but they provide the best available constraints on the end-member. While cave conservation efforts did not permit extensive monitoring of El Pin-dal Cave, the proximity and similar conditions to La Vallina Cave allow us to use this end-member for both sites.
The sensitivities of the measured speleothem proxies (DCF, δ 44/40 Ca, and δ 13 C) to different processes in the soilkarst-cave system allow us to use them to assess the most realistic coupling between measured δ 13 C spel and soil gas pCO 2 . For each combination of soil gas pCO 2 and δ 13 C calculated from the mixing lines, changes in mean soil 14 C concentration, dissolution conditions (termed "gas volume" and indicating the amount of gas that 1 L of groundwater solution interacts with; Owen et al., 2018), and cave air pCO 2 were allowed to vary within realistic bounds (Table 1). These boundary conditions were set based on the available monitoring data, e.g. cave air pCO 2 was left to vary between atmospheric and the maximum soil gas pCO 2 , modelling the effect of cave ventilation dynamics on the proxies. To test whether the system can also be described without invoking changes in soil gas δ 13 C, we performed a second set of experiments ("sensitivity analysis") where all parameters (soil gas pCO 2 , soil 14 C, gas volume, cave air pCO 2 ) were allowed to vary as before but soil gas δ 13 C was kept constant at −18 ‰ ( Table 1).
The model solutions were compared to the measured data from Candela (the stalagmite with the most complete deglacial record) and all solutions matching the measured DCF, δ 44/40 Ca, and δ 13 C spel within a defined interval were extracted. For DCF, the confidence interval of the proxy was chosen, while for δ 13 C spel and δ 44/40 Ca, where measurement uncertainties are much smaller, we defined the threshold at ±1.5 ‰ VPDB and ±0.2 ‰, respectively. Model solutions were filtered sequentially for all three proxies, and each possible permutation of the sequences (e.g. DCF → δ 44/40 Ca → δ 13 C spel ) was calculated. The median and 25/75 % quantiles of all filtered solution ensembles are used as final model result. To avoid too many solutions without matches to the data, we selected the 5 % simulations closest to the measured proxy value for the sensitivity analysis.

Geochemistry
Both Candela and Galia record a substantial decrease in δ 13 C spel between the LG and the EH (Fig. 4). For Candela, δ 13 C spel is highest (−2.48 ‰ and −4.43 ‰ VPDB) at 24.9-15.4 ka and then decreases by about 2 ‰ with the onset of GI-1 (14.4-12.9 ka). After a short-lived increase back to values of ∼ −3 ‰ VPDB at 12.3 ka (corresponding to Greenland Stadial 1, GS-1, Younger Dryas), δ 13 C spel decreases further to −10 ‰ to −7.7 ‰ VPDB in the EH (8.5-7.9 ka). In Galia, δ 13 C spel is −3.88 ‰ VPDB during the LG (26.8 ka) and between −9.78 ‰ and −8.79 ‰ VPDB in the EH (8.7-4.2 ka). Laura covers the time period between 14.3-16.1 ka, where the δ 13 C spel decreases from ∼ −1.7 ‰ to −7.8 ‰ VPDB. Importantly, the absolute values and the magnitude Figure 3. (a) Keeling plot of cave and local atmospheric CO 2 used to define the respired end-member of soil gas. The respired endmember is defined through linear regression of the entire dataset. (b) Mixing lines defined for the model simulations of past soil gas pCO 2 and δ 13 C. We define three mixing lines based on the changes in the atmospheric composition (EH, DEG, LG). All three mixing lines use the same respired end-member, with a variability of ±3 ‰ to account for changes in respired substrate. Mean respired δ 13 C of soil CO 2 across relevant biomes, adapted from Pataki et al. (2003), is shown on the left. of changes in δ 13 C spel in all three stalagmites are comparable over the study period.
The DCF is relatively low in the younger part of the record (∼ 16-4 ka) of all three stalagmites (averages of 6.7 %, 7.2 %, and 13 % for Candela, Galia, and Laura, respectively, Fig. 4). DCF in Candela is slightly higher in the LG portion of the record (∼ 11 %-15 %, 18-20 ka), while values at 24 ka are again comparable with the EH. We disregard the one negative (and physically impossible) DCF value at 24.9 ka, as  Table S2). The DCF in the LG sample from Galia is much higher (23 %) than any in the three stalagmites, but there is no indication of alteration or other reasons why this sample should not be trusted. While the absolute δ 44/40 Ca values in the individual stalagmites are very different, probably reflecting variations in drip path length and drip interval, leading to different amounts of PCP, their temporal variation is remarkably small. In Candela, a slight tendency towards less negative δ 44/40 Ca values can be observed during GS-2.1a, while values are lower during the LG, GS-1, and in the EH (Fig. 4). δ 44/40 Ca values in Galia and Laura are within uncertainty of each other. The δ 44/40 Ca values of the three bedrock samples are consistent, suggesting a homogeneous source of Ca for the three stalagmites (Fig. 4). This allows us to calculate f Ca and quantitatively estimate the amount of PCP for the stalagmites. By their nature, f Ca values mirror the δ 44/40 Ca and suggest that Galia was subject to PCP to a much higher degree than Candela and Laura, where f Ca is comparable. As for the δ 44/40 Ca, f Ca values in all three stalagmites indicate no major changes over the deglaciation, suggesting minimal changes in PCP.
Comparing the three proxies to temperature reconstructions from the Iberian Margin (Darfeuil et al., 2016), using linear interpolation to roughly match the different records, confirms a negative correlation between δ 13 C spel and temperature (−0.63 ‰ • C −1 to −0.9 ‰ • C −1 , r 2 = 0.67-0.96), while the relationship between δ 44/40 Ca and DCF to temperature is weak and/or inconsistent (Fig. 5). Since this comparison is meant to simply illustrate the concept, we refrain from using more advanced statistical methods to determine correlations between the records, and also do not consider any chronological and measurement uncertainties associated with either reconstruction.

Modelling
Each combination of soil gas pCO 2 and δ 13 C produced 363 model solutions, resulting in 13 068 solutions for the LG, 11 979 for the DEG, and 10 890 for the EH (the total number of solutions varies due to extrapolation to lower atmospheric pCO 2 during the LG and DEG). However, only a fraction of the simulations resulted in carbonate precipitation (37 % for LG, 40 % for DEG, and 44 % for EH), while for the rest, precipitation was inhibited by the solution not reaching supersaturation with respect to calcium carbonate. Supersaturation was not reached where low soil gas pCO 2 or closed-system conditions reduced the amount of carbonate being dissolved or where the difference between cave air pCO 2 and solution pCO 2 was very small or negative. Thus, there is no need to further prescribe the cave air pCO 2 as a fraction of the soil gas pCO 2 , as simulations with unrealistic parameter combinations (i.e. higher cave air pCO 2 than soil gas pCO 2 ) are automatically discarded.
Simulations from all three mixing lines produce results that match the stalagmite DCF, δ 44/40 Ca, and δ 13 C spel within measurement uncertainty (Fig. 6). Thus, the initial parameter selection was sufficient to constrain the system and the estimate of the soil respired end-member composition is accurate. Test simulations extending the mixing line to pCO 2 higher than 8000 ppmv consistently lead to overestimation of stalagmite δ 44/40 Ca values, further validating the initial parameter selection.
The matching solutions from all three mixing lines show an increasing trend in median soil gas pCO 2 values over the deglaciation (Supplement Fig. S1). Soil gas pCO 2 values  Supplement Table S2), compared to regional temperature reconstructions (TEX 86 -derived sea surface temperatures from the Iberian Margin; Darfeuil et al., 2016), a Greenland ice core δ 18 O record (NGRIP on GICC05 timescale; Wolff et al., 2010), and global CO 2 (ice core composite from Antarctica; Bereiter et al., 2015). The high-resolution δ 13 C spel record from Candela (thin green line) is shown for reference and was originally published in Moreno et al. (2010). The DCF values were calculated from 14 C measurements paired with U-Th ages. The time periods (LG, DEG, EH) at the top of the figure indicate the intervals used for the modelling to define temperature and atmospheric pCO 2 . are consistently lower during colder time periods (LG, GS-2.1a, and GS-1) and increase during warmer periods (GI-1 and EH), maximising with the onset of the Holocene. This also holds true when considering the results from all mixing lines combined (Fig. 6). Mixing lines 1 and 3 result in few matching solutions for the LG, GS-2.1a, and GS-1, a consequence of the more negative respired end-member δ 13 C used (−24.5 ‰ VPDB and −27.5 ‰ VPDB, respectively), compared to mixing line 2 (respired end-member −21.5 ‰ VPDB). Moreover, sensitivity tests using mixing lines with much higher/lower soil gas pCO 2 (10 000 and 4000 ppmv) but keeping the δ 13 C as in mixing line 1 again results in increasing soil gas pCO 2 values over the deglaciation.
The model is not very sensitive to the choice of DCF threshold. Tests using a higher DCF confidence interval (±3 %) did not lead to any meaningful change in the results. Changes in δ 44 Ca, however, are more important, and we had to increase the confidence interval from the uncertainty from the proxy measurement, as lower uncertainty led to the model not finding matching solutions for all three proxies.
The sensitivity analysis allows more degrees of freedom in the model, where soil gas pCO 2 , soil gas F 14 C, cave air pCO 2 , and gas volume are allowed to freely vary but soil gas δ 13 C is kept constant at −18 ‰. While solutions matching DCF and δ 44/40 Ca are easily found with this set of parameters, the deglacial trend in δ 13 C spel cannot be reproduced (Fig. 6). Only ∼ 2 ‰ of the ∼ 6 ‰ decrease in δ 13 C spel (between 24.9 and 8.5 ka) can be explained through processes other than changes in the soil gas δ 13 C (Fig. 7). It should be noted that the absolute value of the residual calculated from the measured and modelled δ 13 C spel is tied to the initial parameter selection and would vary if we chose differently. The relative differences, however, would remain the same, as long as the initial soil gas δ 13 C is not allowed to vary. We have chosen a relatively high initial soil gas δ 13 C as more negative values result in very few solutions matching the proxy data. This illustrates how the δ 13 C spel trend over the deglaciation requires a change in the initial soil gas δ 13 C. Holding soil gas pCO 2 constant and letting soil gas δ 13 C vary would lead to the entire 6 ‰ change in δ 13 C spel being driven by changes in the respired end-member δ 13 C. This is unrealistic, as biomelevel values of respired δ 13 C typically show little variation (e.g. Pataki et al., 2003;Fig. 3b), and therefore even a substantial deglacial transition from boreal to forested landscape would likely not lead to such a large shift in δ 13 C.

Estimation of the soil respired end-member from cave air
Combined multi-proxy analysis on three stalagmites and geochemical modelling provide strong evidence that changes in initial soil gas δ 13 C are necessary to explain the deglacial trend in δ 13 C spel observed in the NW Iberian Peninsula. Here we show that this trend is best explained by variations in soil respiration and in the relative proportion of respired vs. atmospheric CO 2 in soil gas. Soil gas CO 2 is a mixture of CO 2 from respiration and atmospheric air (Amundson et al., 1998). Therefore, the pCO 2 and isotopic composition of soil gas over depth can be modelled by a mixing line between the atmospheric and soil respired end-member (Pataki et al., 2003). While more recent research has pointed out that this approach neglects spatio-temporal fluctuations in the isotopic Figure 7. Residual δ 13 C spel calculated as the difference between measured and modelled δ 13 C spel over time.
signature of soil CO 2 sources (Goffin et al., 2014), as well as soil storage capacity and the possibility of turbulent transport (Maier et al., 2010), it still provides a valid model with which we can test the overall effects of bulk variations in soil respiration on the drip water solution. In the absence of a full soil monitoring campaign, samples collected from the cave in summer months represent a reasonable approach to estimate the isotopic value of the respired end-member contributing to soil/epikarst gas. This is because the cave, like the soil, is defined by a two end-member mixing system, which is driven by the physical ventilation of the cave.
The main fluxes of carbon in a system like El Pindal and La Vallina caves are from soil gas (mainly seeping through the host rock and into the cave) and atmospheric air (through ventilation). Like many mid-latitude and high-latitude cave systems, there is a seasonal reversal in the airflow direction in La Vallina cave (Stoll et al., 2012). In the summer, when cave air is colder than exterior air, cave air flows out the entrance and is replaced by inflow and diffusion of soil/epikarst gas. In this season, the cave air has the highest CO 2 concentrations, with an isotopic composition that falls close to the soil respired end-member on the Keeling plot (Fig. 3a). In the winter, when cave air is warmer than exterior air, exterior air flows in through the cave entrance, bringing the cave closer to the atmospheric end-member. The data from the monitoring of the cave primarily reflects CO 2 from the soil that is drawn through the karst network into the cave. It therefore likely reflects soil column-integrated conditions and the full contribution of respired CO 2 in the soil and epikarst unsaturated zone (below the soil, "ground air"). We do not consider the effect from increasing canopy cover over the deglacial transition, which may affect the soil gas δ 13 C through the canopy effect (more depleted CO 2 close to the ground; Buchmann et al., 2002), as well as atmospheric turbulence and advection effects, as these effects lead to small shifts in soil gas δ 13 C compared to the recorded signal. Any contribution of carbon from bedrock dissolution does not significantly affect our estimation of the respired endmember because the intercept defining the respired endmember is most influenced by the summer season cave air pCO 2 data (Fig. 3a). During the summer season drip flow rates are more than an order of magnitude lower than in the winter, and degassing from this drip is suppressed by the high cave air pCO 2 . In winter, when drip rates are higher and cave air pCO 2 is lower, degassing may contribute to carbon in cave air, as seen in other systems (Waring et al., 2017). However, winter monitoring data corresponding to ventilated periods are near the global atmospheric composition, suggesting an insignificant impact of degassing of dissolved limestone on our calculated mixing line. Furthermore we do not find evidence for a different Keeling intercept in winter and summer, unlike monitoring studies which infer a strong effect of degassing of a carbon source from limestone dissolution (Waring et al., 2017).
The seasonality of the modern cave air has the advantage of helping to define the modern respired end-member. From our monitoring, we do not find evidence for a different isotopic value of the respired end-member in different seasons. Thus, exploiting the seasonal ventilation of the cave to define the mixing line and respired end-member does not preclude using this respired end-member to interpret records from speleothems in which deposition is dominant in one season. We therefore argue that summer cave air can be used to estimate the isotopic composition of the respired end-member of soil CO 2 , when the competing fluxes are minimised.
5.2 Temperature sensitivity of soil respiration as the main driver for δ 13 C spel Our modelling results show a consistent pattern of increasing soil gas pCO 2 over the last deglaciation, with absolute values ranging between ∼ 410-600 ppmv during the LG (at 24.9 ka), and ∼ 1000-6000 ppmv during the EH (at 8.6 ka), depending on the mixing line used. An increase in soil respiration rates coinciding with deglacial warming is likely, as higher temperatures promote more rapid soil carbon turnover (Vaughn and Torn, 2019) and the establishment of denser forests (Vargas and Allen, 2008). Climate model simulations confirm that net primary productivity in the NW Iberian Peninsula was lower during the LG than at present (Scheff et al., 2017). Pollen studies from the NW Iberian Peninsula show significant and rapid changes in vegetation type and cover over the Pleistocene-Holocene transition (Moreno et al., 2014). While LG pollen reconstructions suggest a landscape dominated by open grassland (30 %-35 % Poaceae) with significant steppe taxa and low arboreal pollen (30 %-50 % primarily Pinus sylvestris and Betula), the EH pollen assemblage is dominated by arboreal pollen (70 %-90 %; Moreno et al., 2011). It is likely that the rapid response of pollen assemblages to climate warming is due to the region's proximity to documented tree refugia in the Mediterranean region (Fletcher et al., 2010). Assuming a temperature change of roughly 7 • C between the LG and EH, in line with TEX 86 -based temperature reconstructions from the Iberian Margin (Darfeuil et al., 2016), the sensitivity of soil respiration to temperature change (Q 10 , i.e. factor by which soil respiration increases with a 10 • C rise in temperature) derived by our modelling experiments lies between 3.6 and 14.7, depending on the initial conditions of the models. This is higher than the mean global Q 10 values of 3.0 ± 1.1 found by the soil respiration database (Bond-Lamberty and Thomson, 2010), and may be exaggerated by chronological uncertainty in the marine and speleothem records and by uncertainties in the temperature proxies (> 4 • C, Darfeuil et al., 2016;Tierney and Tingley, 2015). Changes in seasonality between EH and LG could also lead to bias in the temperature reconstruction due to shifts in nutrient availability and marine productivity (Darfeuil et al., 2016). Furthermore, a change in respired substrate over time, leading to a shift in the soil respired endmember δ 13 C (Boström et al., 2007), could lead to higher Q 10 than the global mean. We can exclude changes in vegetation assemblage from C4 to C3 plants, as there is no evidence for the widespread presence of C4 plants during the glacial in the NW Iberian Peninsula (Moreno et al., 2010) or elsewhere at temperate western European sites (Denniston et al., 2018;Genty et al., 2006Genty et al., , 2003. A change in the balance between heterotrophic and autotrophic respiration is another possibility that would influence the soil gas δ 13 C. Changes in temperature affect root and microbial respiration differently (Wang et al., 2014), as do changes in other environmental variables, e.g. precipitation regimes and nutrient cycling (Li et al., 2018). Microorganisms are typically enriched by 2 ‰-4 ‰ compared to plants (Gleixner et al., 1993), and vertical enrichment by ∼ 2.5 ‰ in soil profiles has been attributed to an increasing contribution of soil microbially derived material with depth to the overall soil carbon turnover (Boström et al., 2007). The release of older and enriched carbon from soils and long-lived plant material through respiration could provide an additional mechanism with which the soil gas δ 13 C could be shifted regardless of changes in soil respiration (Fung et al., 1997). Very high Q 10 values have also been found in winter at a Danish beech forest site, suggesting a stronger temperature (and potentially soil moisture) control during the cold season, potentially resulting in the activation of dormant microbial communities, the alteration of diffusion processes of organic molecules and of cell metabolism (Janssens and Pilegaard, 2003). It is possible that the significant environmental changes occurring over the last deglaciation resulted in similar responses from the soil microbial community, leading to the observed high Q 10 values. A higher respired δ 13 C during the LG is also suggested by the model results, where mixing line 3 with the lowest respired δ 13 C (−25.9 ‰) fails to produce solutions matching the speleothem data (Fig. S1). Given the small variation in DCF values in Candela over the deglaciation, we can exclude the possibility that changes in the fraction of bedrock carbon from changing dissolution conditions constitute an important driver of the deglacial signal.
Another intriguing possibility is that the carbon isotopic fractionation of C3 vegetation is controlled by atmospheric pCO 2 (Schubert and Jahren, 2015). A recent global compilation of speleothem records shows that, after correcting for the expected effect of precipitation and temperature on δ 13 C of C3 biomass and the temperature-dependent fractionation between CO 2 and calcite, the global average δ 13 C spel closely tracks atmospheric pCO 2 over the last 90 ka (Breecker, 2017). The magnitude of the deglacial shift in C3 plant δ 13 C has been proposed to lie around 2.1 ‰ (Schubert and Jahren, 2015). The deglacial δ 13 C spel record from the NW Iberian Peninsula, however, shows clear millennialscale variations that coincide with temperature variations, but are not driven by atmospheric CO 2 (Fig. 4). Therefore, while it is possible that a CO 2 fertilisation effect contributed to the overall decrease in δ 13 C spel over the deglaciation, this effect is likely not dominant.
Our findings will likely apply more broadly for caves in settings where soil gas pCO 2 is temperature limited. Sites where soil gas pCO 2 is moisture limited, e.g. further south on the Iberian Peninsula, will likely exhibit very different trends in δ 13 C spel over glacial-interglacial cycles, as hydroclimate and temperature may have different phasings.

Other processes affecting δ 13 C spel
While a change in soil respiration and consequently in the proportion of respired vs. atmospheric CO 2 in the soil gas can explain the deglacial trend in δ 13 C spel , a number of other, cave-specific processes could also contribute to changes in δ 13 C spel . The direct effect of the glacial-interglacial temperature change on carbonate equilibria and fractionation factors is small and taken into consideration by running the simulations with EH, DEG, and LG parameters. It is more difficult to assess whether kinetic fractionation effects affected the stalagmite at different times, potentially amplifying the δ 13 C spel signal. CaveCalc uses standard kinetic fractionation factors for the CO 2 -DIC-carbonate system (Romanek et al., 1992;Zhang et al., 1995), where DIC stands for dissolved inorganic carbon, and therefore such variations are not considered by the model. However, the high degree of coherence between δ 13 C spel records from the entire temperate western European region suggests that localised, cave-specific kinetic fractionation effects likely played a minor role in driving the deglacial trend (Fig. 1).
Changes in the amount of PCP the drip water experiences en route to the speleothem can lead to significant variability in δ 13 C spel records (Fohlmeister et al., 2020), and they are tightly coupled to changes in cave air pCO 2 and cave ventilation dynamics. Higher cave air pCO 2 and a reduced CO 2 gradient between the supersaturated drip water solution and the cave air result in less PCP and vice versa for lower cave air pCO 2 . It is likely that cave air pCO 2 was lower during the last glacial at the study sites, and indeed this is also suggested by our model results (Supplement Fig. S2). Cave air pCO 2 is coupled to soil gas pCO 2 , which provides its upper limit, and model results automatically filter out unrealistic scenarios, as no speleothem precipitation occurs when cave air pCO 2 is equal to or higher than soil gas pCO 2 . Our multiproxy dataset allows us to evaluate the importance of PCP for δ 13 C spel quantitatively, as δ 44/40 Ca can provide quantitative PCP reconstructions over time (Owen et al., 2016). Mg/Ca ratios are also often used as a proxy for PCP; however, caution is required in their interpretation in El Pindal Cave because in the Holocene, Mg/Ca is also affected by increasing surf-zone marine aerosol contributions as rising sea level brought the coastline to the foot of the sea cliff in which the cave has its entrance (Supplement Figs. S3 and S4). Over the last deglaciation, δ 44/40 Ca and f Ca varied only minimally in both Candela and Galia (Fig. 4), suggesting that changes in PCP were small. This is also reflected in the sensitivity analysis, where changes in δ 13 C spel cannot be reproduced while also fitting the δ 44/40 Ca curve (Fig. 6). CaveCalc uses cave air pCO 2 to match the degree to which drip water has lost its initial Ca due to calcite precipitation, giving us a measure for PCP. A solution equilibrated with a high soil gas pCO 2 would lose the majority of its carbonate in a simulation where cave air pCO 2 is atmospheric, due to the high degree of oversaturation of the drip water solution compared to cave air. If δ 44/40 Ca provides evidence that only a small portion of Ca has been precipitated, then the simulation must match the data by prescribing a higher cave air pCO 2 . In reality, the fraction of Ca precipitated from drip waters depends not only on the oversaturation of the solution, but also on the time the water is present as a thin film on the cave ceiling and stalagmite surface before being replaced by a new water parcel (i.e. drip interval; Fohlmeister et al., 2020;Stoll et al., 2012). When the drip interval is short, none of the water parcels will have enough time to fully degas CO 2 and equilibrate with the cave atmosphere, and PCP is lower than what would be possible given the cave air pCO 2 . CaveCalc does not model drip interval, and therefore the cave air pCO 2 inferred from the simulations might be overestimated. We test the effect of drip interval length changes on f Ca and PCP using the forward model ISTAL (Stoll et al., 2012), which explicitly models this parameter. Two model scenarios simulate full glacial and Holocene conditions, including changes in temperature, cave air pCO 2 , and soil gas pCO 2 for "winter" (i.e. atmospheric) and "summer" (i.e. elevated) cave air pCO 2 (Supplement Fig. S4). The effect of the glacial-interglacial temperature change is only significant for high drip intervals during the cold season, where PCP is slightly higher during interglacial conditions. At high drip intervals, the temperature increase leads to a change in f Ca of ∼ −0.1, which translates to a ∼ 0.04 ‰ change in δ 44/40 Ca and a −0.7 ‰ VPDB change in δ 13 C spel . This corroborates our expectation from the δ 44/40 Ca record and CaveCalc model results, suggesting that only a small part of the shift in Candela δ 13 C spel over the last deglaciation was due to changes in PCP.
While it is likely that some or all of these processes affected the deglacial δ 13 C spel to some extent, their magnitude is not large enough to explain the measured ∼ 6 ‰ shift, suggesting that changes in soil gas pCO 2 played a significant role.

Insights into regional hydroclimate over the last deglaciation
Our new multi-proxy record from stalagmites from the NW Iberian Peninsula also offers nuanced insights into local hydroclimate conditions over the last deglaciation. While DCF mainly responds to changes in carbonate dissolution conditions and therefore is sensitive to changes in infiltration, δ 44/40 Ca is driven by both infiltration dynamics (determining the initial oversaturation of drip water and the degassing timescale) and cave atmospheric pCO 2 (determining the amount of PCP occurring). The Candela record suggests no substantial shift in infiltration dynamics or PCP occurring between LG and EH (Fig. 6), as both proxies fluctuate around a mean value without long-term trends. This result suggests that the glacial hydroclimate in the NW Iberian Peninsula was not significantly different from the Holocene and stands at odds with previous mainly pollen-based studies that often point towards a drier glacial but with considerable variability over millennial timescales (Fletcher et al., 2010). Recent modelling results have challenged the interpretation of the glacial being cold and dry, suggesting instead that, while precipitation was lower during the LG, topsoil moisture was actually higher than at present (Scheff et al., 2017). Our new stalagmite data support this interpretation, suggesting that temperature, and not hydroclimate conditions, was the main driver of ecosystem productivity over the deglaciation.

Conclusions
We have combined multi-proxy (δ 13 C, δ 44/40 Ca, and DCF) data from three speleothems and quantitative geochemical modelling to show that the temperature sensitivity of δ 13 C spel over the last deglaciation in western Europe is best explained by increasing soil respiration. Generating a large ensemble of forward models of processes in soil, karst, and cave allows the estimation of their likely importance and variability over time. Speleothem geochemical proxies that are sensitive to different components of the soil-karst-cave system can be employed to extract the most likely model solutions from the ensembles, thus quantifying the system's initial conditions, particularly soil gas pCO 2 . Our approach involved the coupling of soil gas pCO 2 and δ 13 C values, as expected when following a mixing line between a soil respired and an atmospheric end-member, thus allowing us to model changes in soil respiration. While uncertainties remain, in particular with respect to possible changes in the soil respired endmember δ 13 C value over time, we find that an increase in soil respiration is necessary to explain the large shifts in δ 13 C spel over the last deglaciation in the NW Iberian Peninsula. Given the exceptional regional coherency of δ 13 C spel records over temperate western Europe, it is likely that this effect is of broader regional significance. Our study is the first to quantitatively model environmental processes in karst systems using a multi-proxy approach and paves the way towards more nuanced interpretations of δ 13 C spel records. Moreover, our multi-proxy records support recent climate model results that reject the long-standing "drier and colder glacial" notion in western Europe, pointing instead toward a dominant forcing of temperature, rather than hydroclimate, on ecosystem productivity.
Author contributions. FAL, HS, and GMH designed the study and acquired funding for the project. FAL, NH, and CCD performed the geochemical analysis on speleothem samples. OK collected and measured cave air samples from La Vallina Cave and acquired funding for the monitoring work. FAL and MW performed the modelling experiments in CaveCalc and wrote the R code for the data-model evaluation. FAL wrote the paper and generated the figures. HMS and CCD provided additional input to the text. All authors provided feedback on the paper and approved it before submission.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
thank Yu-Te (Alan) Hsieh for assistance with the δ 44/40 Ca measurements at the University of Oxford, Madalina Jaggi at ETH Zurich for measurements of δ 13 C and trace elements, and Saul Gonzalez-Lemos for cave air sampling. We thank two anonymous reviewers and handling editor Marie-France Loutre for their fair and competent assessment of our paper.
Financial support. This research has been supported by the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no. P400P2_180789), by ETH core, and from ETH doctoral fellowship ETH-13 18-1.
Review statement. This paper was edited by Marie-France Loutre and reviewed by two anonymous referees.