Thermal impact of Heinrich stadials in cave temperature and speleothem oxygen isotope records

Abstract

During each Heinrich stadial (HS), temperatures in southern Europe typically dropped several degrees during several hundred to few thousand years. We have developed a one-dimensional thermal conduction model that transfers the typical surface temperature anomaly of a HS to a series of hypothetical underlying caves. The results show that with increasing depth, the thermal anomaly is attenuated, the lag time increases, and the signal structure experiences larger modifications. The model suggests that in most cases, it is not acceptable to assume a synchronous thermal variability and similar average temperature values between the surface atmosphere and the cave interior at millennial timescales. We also simulated the thermal impact of the modeled HS on speleothem δ¹⁸O records. The outputs of most model scenarios suggest that temperature changes associated with the HS produce δ¹⁸O anomalies capable of contributing significantly or even decisively to the speleothem isotope variability. Therefore, despite controls other than temperature often being considered more important when interpreting Pleistocene speleothem δ¹⁸O records in temperate climates, this research suggests that temperature is expected to be one of the major controls of δ¹⁸O values in most cave sites outside the tropics and should be included as a significant parameter affecting Pleistocene speleothem δ¹⁸O records.

INTRODUCTION

Ocean sediments from the midlatitudes of the North Atlantic record the occurrence of “armadas” of icebergs in proxies such as the ice-rafted debris (IRD). The occurrence of these events that displaced the discharge of detrital material from icebergs to unusually south latitudes are known as Heinrich events, named after the scientist who first identified them and related them to the climate changes associated with them (Heinrich, 1988). Heinrich events occurred during glacial periods of the Quaternary and had durations of several hundred years to a few millennia, during which the climate in most regions of the Northern Hemisphere typically shifted to significantly colder and/or drier conditions (e.g., McManus et al., 1999; Wang et al., 2001; North Greenland Ice Core Project Members, 2004). The origin of these climate changes involves complex interactions between continental ice and ocean dynamics, with icebergs being a consequence rather than a cause (Barker et al., 2015). Additionally, the duration of these exceptional cold and arid climate conditions is not limited to the periods when IRD is recorded in ocean sediments (Hemming, 2004; Barker et al., 2015). Therefore, when focusing on the climate changes related to these events rather than the oceanographic phenomenon in particular, it is more accurate to use the term “Heinrich stadials,” hereafter referred as HSs or HS when referring a singular stadial. Not all sites record a drop in temperature in relation to a HS, and in those sites where negative thermal anomalies were identified, their magnitude is not homogeneous. Different proxies from ocean sediments in mid- and high-latitude sites often record thermal anomalies of several degrees (e.g., Martrat et al., 2007). Paleotemperature estimations from continental proxies are not calculated as straightforwardly as estimations from ocean sediments, although a conservative approximation is to extend the sea-surface thermal anomalies to land as supported by climate models (e.g., Pedro et al., 2018), keeping in mind at all times the limitations of this assumption.

Temperature in caves is a critical parameter for the cave seasonal dynamics, affecting processes such as cave air ventilation or speleothem growth rates (Spötl et al., 2005; Banner et al., 2007; Kowalczk and Froelich, 2010). Potentially, changes in temperature also have significant impacts on the underground biota (Mammola et al., 2019) and affect the preservation of cave art (Bourges et al., 2014). During the precipitation of speleothems, temperature is a key control on the fractionation of isotopes or the partition of trace elements (Epstein et al., 1953; Huang and Fairchild, 2001). Consequently, changes in cave air temperature are potentially a significant factor affecting the interpretation of the geochemical records in terms of climate changes (McDermott, 2004; Fairchild and Treble, 2009). Cave temperature is especially important in areas where permafrost conditions above the cave affect the flow of seepage water (Genty et al., 2003) or cause the formation of cryogenic carbonates (Žák et al., 2004). Also, temperature is of key importance in caves with ice or snow accumulations (Luetscher et al., 2008). Despite the critical role that temperature has in caves, there are only a limited number of studies on the transfer of climate changes to caves or tunnels (Perrier et al., 2005; Domínguez-Villar et al., 2015). Beyond cave entrances, where external air has a strong influence due to enhanced air advection (Wigley and Brown, 1976; Smithson, 1991), most cave interiors record a limited thermal variability, with temperature values very close to the annual average temperature recorded on the surface over the cave (Moore 1964; Moore and Sullivan, 1964). For those caves that record little thermal variability, the physical mechanism that transfers heat from the surface to the cave is in most cases thermal conduction (Cropley, 1965; Moore and Sullivan, 1978). Exceptions to the general pattern are found in caves where dynamic air currents circulate throughout the cavity or there are significant water streams (Kranjc and Opara 2002; Covington et al., 2011). In these exceptions, advection and radiation are more efficient mechanisms than conduction for transferring heat. Seasonal cave ventilation affects the air temperature in most caves, although the amplitude of the seasonal thermal anomalies is often very limited (e.g., Milanolo and Gabrovšek, 2009). In these cases, the background temperature and the thermal amplitude of the seasonal anomalies are still controlled by thermal conduction, despite advection partially contributes to the variability of the thermal signal (Domínguez-Villar, 2012; Domínguez-Villar et al., 2015).

The transport of heat from the external atmosphere to the cave by conduction has two main effects on cave air temperature: (1) The magnitude of thermal anomalies recorded in the external atmosphere is reduced with depth. As a result of the progressive attenuation of thermal anomalies with depth, seasonal cycles transferred by conduction are completely filtered out at 20 m below the surface (Pollack and Huang, 2000), although the thermal anomalies of more durable thermal cycles continue their propagation to deeper sections of karst. (2) The thermal signal from the surface takes some time to be transferred to a certain depth; consequently, cave temperature signals lag those from the exterior (Cropley, 1965; Moore and Sullivan, 1978; Badino, 2004). Thermal signals are transferred underground as waves characterized by particular periods, amplitudes, and phase shifts (Smerdom and Stieglitz, 2006). Although the record of surface air temperature (SAT) in nature is more complex than simple waves, every signal can be decomposed into a series of independent waves or harmonics, and those simple waves are transferred underground. The propagation velocities underground depend on the period of the harmonic (Villar et al., 1983); for example, thermal anomalies associated with annual cycles are transferred to a certain depth faster than those anomalies associated with 11 yr cycles. The complexity of thermal conduction has limited the study of underground thermal conduction to recent decades or centuries (e.g., Perrier et al., 2005; Pollack et al., 2005; Beltrami et al., 2006), and studies that consider climate changes with cycles of thousands of years are rare (Beltrami et al., 2014). However, when studying the long-term evolution of cave microclimate in detail, or when interpreting proxy records from speleothems that cover periods of thousands of years, the impact of long-term changes in cave temperature should not be ignored.

Although the stable air temperature of many cave interiors is similar to the mean annual temperature measured outside the cave, the coupling between the annual SAT and the cave air temperature is often not perfect due to multiple causes (Beltrami and Kellman, 2003; Domínguez-Villar et al., 2013b; Cuthbert et al., 2014). Additionally, because a transit time is required for thermal anomalies recorded at the surface to be transferred underground by conduction, lags are often significant at interannual timescales, enhancing the thermal decoupling. This thermal decoupling has significant impact on speleothem δ¹⁸O records (Domínguez-Villar et al., 2013b), one of the most commonly used proxies in speleothem paleoclimate research. Fractionation of oxygen isotopes that are involved in the formation of speleothems is deeply impacted by temperature (Lachniet, 2009). In midlatitudes, the empirical fractionation factors of water condensation and equilibration with the atmosphere often have a similar magnitude and different sign than the fractionation under equilibrium conditions that takes place during the precipitation of calcite in caves. For this reason, in midlatitude caves, temperature is often neglected as a significant control on the speleothem δ¹⁸O record, and most of its variability is assumed to originate from changes in the isotope composition of drip water (e.g., Domínguez-Villar et al., 2008; Rossi et al., 2018). However, in other midlatitude cave sites, these fractionation factors do not fully counteract each other, and large temperature changes such as those from the Pleistocene are expected to have a significant impact on the speleothem δ¹⁸O records (Wainer et al., 2011). Therefore, cave temperature control should be carefully evaluated when studying speleothem δ¹⁸O records from midlatitude sites.

Here we provide a theoretical model that transfers the thermal anomaly of a typical HS to a series of caves located at different depths. The model is based on the current knowledge of heat conduction in karst regions, which was previously validated with observations in real caves (Domínguez-Villar et al., 2013b, 2015). Because conduction is the main mechanism to transfer thermal anomalies to the inner part of most caves, the conclusions of this study are transferable to a wide range of real caves. The differences in amplitude and lag time between external and cave thermal signals are also used to evaluate the impact in one of the most commonly used proxies in speleothem paleoclimate research: the ratio of oxygen isotopes (δ¹⁸O). The findings of this research have important implications for multiple fields of research such as paleoclimate based on speleothems, cave art conservation, ecological evolution in caves, human and/or faunal evolution, long-term changes in cave microclimate, long-term history of mountain karst aquifers, and the evolution of cave processes (such as modifications in the formation of speleothems, preservation of ice or snow in caves, etc.).

METHODS

Surface air temperature

Although we aim to provide realistic parameters to reproduce climate changes related to HSs, this is a theoretical study and does not focus on a particular location or HS. The magnitude and duration of the modeled HS considered data from the North Atlantic off the Iberian Peninsula and the Mediterranean during the Late Pleistocene (Cacho et al., 1999; Pérez-Folgado, et al., 2003; Essallami et al., 2007; Martrat et al., 2007; Castaneda et al., 2010; Voelker and de Abreu, 2011). The hypothetical HS that we reproduce does not consider HSs that occurred during terminations, as their duration, evolution, and climate conditions often differ from HSs that occur during full-glacial climate conditions (e.g., Denton et al., 2006; Domínguez-Villar et al., 2020). We assume that the continents and sea surface recorded similar thermal anomalies. Alkenone derived sea-surface temperature reconstructions were preferred as paleotemperature records, as thermal histories based on fossil assemblages provided inconsistent results in the different regions considered. The modeled HS lasts 1000 yr and has a thermal anomaly of 4°C. The background temperature outside the HS is set in our model to 8°C. We selected this value because we assumed that the cave would be several degrees warmer during interglacial conditions. Thus, during interglacial conditions the temperature of our hypothetical caves would be within the current range of annual average temperatures observed in wide regions of western and southern Europe. Additionally, the selected background temperature does not cause underground freezing during the modeled HS.

The SAT is constructed from a basic wave and up to five of its harmonics to produce a square-like wave. Every wave is characterized by its period, amplitude, and wave phase. The computation of waves as well as the thermal model (see following section) was implemented in a standard spreadsheet program using basic trigonometric, logical, and arithmetic formulas. So, no programming codes were required to perform any calculation in this study. The basic temporal unit of the model is accounted in days, and every time step lasts a month. Every 4 yr, one additional day is added to February to account for leap years. The period of the harmonic with the shorter duration is 181 yr. Thermal anomalies with smaller periods, including the seasonality, are filtered out in this study to emphasize long-term thermal changes. We produce a SAT signal of 5000 yr divided in three different intervals. The SAT signal during the interval from 1000 to 3000 modeled years is the sum of six sinusoidal waves that eventually reproduce the HS from 1500 to 2500 modeled years (Fig. 1). Note that the time scale in all graphs is in modeled years and not in years before present. Consequently, earlier events occur at the left of the graphs and subsequent events to the right. We also provide a background thermal variability before and after the interval from 1000 to 3000 modeled years. In these earlier and later intervals, the SAT signal is composed of only two of the harmonics used to construct the HS signal. This background thermal variability is introduced for the thermal history to be more realistic and to visualize how the HS anomaly interacts underground with earlier and subsequent thermal changes.

Figure 1. Surface atmosphere temperature (SAT) signal considered in the model. The 5000 yr modeled are divided into three intervals. Between 1000 and 3000 modeled years, the record results of the sum of six sinusoidal signals that reproduce a common thermal anomaly during Heinrich stadials. During the previous and subsequent intervals, the record is composed of the sum of two sinusoidal signals that provide background variability. Notice that the timescale is reported as modeled years and not as years ago.

Underground transfer of the SAT

We used a one-dimensional model in which the SAT signal is transferred underground by thermal conduction exclusively. The model assumes a homogeneous infinite half-space with uniform thermal diffusivity. It applies only to caves and karst terrains in the vadose zone, where the phreatic level is deep enough not to interfere with the thermal regime of considered depths. The model was applied to four hypothetical caves located at different depths. The term “depth” is used in this paper for simplicity, although in every case we refer more specifically to the thickness of bedrock above the cave. The selected caves are located 10, 50, 100, and 500 m underground, which covers most cave settings in vadose karsts and is enough to illustrate the importance of the cave depth for the transfer of thermal anomalies. Some alpine caves (i.e., located under high-elevation ranges) can exceed the range of depths here considered. However, at millennial timescales, climate or environmental changes at the surface have a negligible impact on the temperature transferred by conduction to such deep caves, and the thermal consequences would be similar to the ones described in our deeper scenario. Once the thermal signal is conducted from the surface to the ceiling and walls of the cave, the air temperature in the gallery assimilates the temperature of the host rock mainly by radiation in a quick process (Guerrier et al., 2019). We assume that our modeled caves have no significant water streams and that they have a limited air advection that is not enough to disturb the underground host rock temperature signal. The geothermal gradient and the thermal gradients found in the vadose zone of karst as a result of adiabatic gradients, vertical water flow, and latent heat exchanges are not considered here (Badino, 1995; Luetscher and Jeannin, 2004). In the case of shallow caves, these are negligible controls, and in deep alpine caves, they only produce limited and systematic shifts in absolute temperature (Domínguez-Villar, 2012).

The thermal signal of every harmonic is transferred underground, and the anomalies are added together according to Eq. 1, where T_z is the underground air temperature at depth z (measured in °C), A_z is the amplitude of each signal recorded at the depth z for every wave (measured in °C), t is the time (measured in days), P is the period of each wave (measured in days), ɛ is the wave phase of each wave (measured in radians from 0 to 2π), and Φ_z is the lag time of each wave at the depth z.

(1)$$T_z = \sum {A_z}\,\,\ast\,\, {\sin } \Big( \displaystyle{{2\,\,\ast\, \pi \,\,\ast\,\, t} \over P}\Big) + \varepsilon -\Phi _z$$

The value of A_z is calculated from the attenuation of the amplitude (a_A), a parameter obtained for every wave at particular depths from Eq. 3, and the amplitude that every wave has at the surface (A_s).

(2)$$A_z = a_A\,\,\ast \,\,A_s$$

(3)$$a_A = e^{( -\Phi _z) }$$

The Φ_z values are calculated from the wave vector (k) measured in radians per meter and the depth measured in meters according to Eq. 4. The k values are calculated according to Eq. 5 from the period of every wave and the thermal diffusivity (κ), a parameter taken from the karst thermal conduction literature (Domínguez-Villar et al., 2013b).

(4)$$\Phi _2 = k\ast z$$

(5)$$k = \sqrt {\displaystyle{\pi \over {P + \kappa}}} $$

We have implemented our thermal model in two different scenarios. The first scenario assumes thermal coupling between SAT and ground temperature and represents caves with a forest cover above the cave that did not change significantly as a result of the climate changes associated with the HS. Good examples of this scenario are sites where the geographical settings act as refugia for trees regardless of the ecological changes taking place in response to the HS in areas proximal to the site. The second scenario assumes a variable thermal coupling between SAT and ground temperature related to changes in the vegetal canopy over the cave site. The changes in thermal coupling conditions are the ecological response to the climate change associated with the HS, because the microclimate under different vegetal canopies impacts the coupling of SAT and ground temperature (Domínguez-Villar et al., 2013b). In this second scenario, we assume an annual thermal decoupling between SAT and the ground temperature of 2°C that is applied only during the HS interval. This implies that before and after the HS, surface air and ground temperatures are coupled and no thermal shift is applied. Thus, under this scenario, the change to a more open landscape over the hypothetical caves during the HS (Combourieu-Nebout et al., 2002; Brauer et al., 2007; Fletcher and Sanchez-Goñi, 2008; Camuera et al., 2019) results in a mean annual ground temperature 2°C warmer than the mean annual SAT, limiting the impact of the HS on the underground thermal anomaly during the HS.

Oxygen isotope fractionation

We use SAT and T_z signals to calculate an isotope model that evaluates the impact of temperature on the speleothem δ¹⁸O records. This model considers two fractionation events: first at the moment of occurrence of atmospheric precipitation, and second, when calcite is formed from the drip-water solution. The so-called temperature effect that takes place during the formation of rainfall (Dansgaard, 1964) is based on empirical relationships between rainfall δ¹⁸O values and SAT. Isotope fractionation that is temperature related occurs during condensation of rain droplets, as well as during the equilibration of rainwater with the atmospheric moisture (Gat, 2010). The ratio between rainwater δ¹⁸O (δ¹⁸O_rw) and SAT in midlatitude sites has a positive relationship, although its magnitude is slightly variable from site to site (Rozanski et al., 1993). We have selected four relationships between δ¹⁸O_rw and SAT that cover the range of the observed variability in most midlatitude Mediterranean sites (International Atomic Energy Agency, 2009). These δ¹⁸O_rw/SAT ratios are 0.23, 0.29, 0.35, and 0.41‰/°C.

Our isotope model also considers a second fractionation event that occurs at the moment of speleothem formation. The hypothetical speleothems are composed of calcite and precipitate under equilibrium conditions. The model assumes that no significant evaporation takes place in the karst system and that the drip water is representative of the annual weighted mean of the δ¹⁸O_rw values. These assumptions were proved right in temperate sites (Genty et al., 2014; Domínguez-Villar et al., 2018), although they might not be realistic in some arid environments (Ayalon et al., 1998; Cuthbert et al., 2014). Although the model considers an effective mixing of the rainwater in the epikarst above the cave during all months of the year, to limit the complexity of the model, we assume a direct transmission of the mixed rainwater to the drip site. We are aware that a negligible residence time (<1 yr) is an unrealistic assumption for most caves. Typical residence times of drip water measured in caves range from one or several years (e.g., Kluge et al., 2010) to decades (e.g., Yamada, et al., 2008), even when the caves are hundreds of meters under the surface (Chapman et al., 1992). To remove the high-frequency variability in the thermal records, this study considers only climate signals with periods >180 yr. Therefore, at the scale of this study, even drip waters with residence times of several decades would have limited impact on the variability of the speleothem δ¹⁸O records (i.e., <0.2‰). Drip-water temperature is assumed to be mostly influenced by cave air temperature which is controlled by bedrock temperature (Guo et al., 2019). The precipitation of calcite is set to have a −0.23‰/°C δ¹⁸O–cave temperature relationship for the temperature range explored (Kim and O'Neil, 1997; Kim et al., 2007). The chosen thermal relationship agrees within uncertainty with the values found in natural and cave-analog experiments (e.g., Johnston et al., 2013; Hansen et al., 2019). The isotope model reports the isotope anomalies related to temperature changes alone. This implies that changes in δ¹⁸O_rw values related to controls other than temperature are not considered in the modeled isotope anomalies. However, these speleothem δ¹⁸O anomalies highlight and quantify the potential impact of temperature on speleothem δ¹⁸O records and are useful for evaluating the role of temperature in speleothem δ¹⁸O record interpretations.

RESULTS

Cave temperature models

The temperature signal presented in Figure 1 was transmitted by thermal conduction to the caves located at 10, 50, 100, and 500 m underground assuming a thermal coupling between SAT and ground temperature. The original thermal signal was decomposed in six different waves with their particular periods and amplitudes. Each of these signals is recorded at every modeled cave with a different amplitude and lag time (Fig. 2). The amplitude of the thermal anomalies is exponentially muted with depth. The amplitude of the signals with longer periods is more gradually attenuated, which makes the period, and not so much the thermal anomaly, the key control to record significant thermal anomalies in deep caves. The model results show that the lag time for the SAT to be transferred to the cave also depends on the period of every signal. This is important, because the original anomalies in the SAT signal are recorded in the same cave at different times depending on the periods of the different waves that compose the SAT signal. Another obvious implication of the model is that the depth of the cave is a critical parameter controlling both the amplitude of the thermal anomaly and the lag time in relation to the SAT signal.

Figure 2. (color online) Underground attenuation of surface atmosphere thermal anomalies and their lag times at selected depths. (A) Underground attenuation of thermal amplitudes of the anomalies for the six signals considered in this study for the atmosphere–soil thermally coupled scenario. P: signal period. The thermal amplitude of the anomaly for every signal is half the difference between the maximum and minimum temperature recorded within a cycle. (B) Lag times modeled for simulated caves at 10, 50, 100, and 500 m in depth (z) as a function of the period of each signal.

The results of transferring the SAT shown in Figure 1 to the hypothetical caves under the first scenario (i.e., thermal coupling between SAT and ground temperature) are shown in Figure 3. The model shows that the temperature signal of the cave located 10 m underground is very similar to the SAT signal, with only some decades of delay. Thus, the background temperature of 8°C cools down on average to 4°C temperature during the HS in both the SAT and 10-m-deep cave temperature records. However, caves located at 50 and 100 m underground already show important attenuations of the maximum thermal amplitude recorded in the caves. Thus, the cave located at a depth of 100 m records a minimum temperature around 6°C in response to the HS. Not just the amplitude of the anomaly, but also the structure of the HS is modified in the cave thermal signal. The different structure of the thermal signal in the cave is caused by the variable periods of the waves that form the SAT signal, which reflect different transit times until they are recorded in the cave. Because the temperature recorded at the cave is the sum of all the waves at a particular moment, the result is an asymmetrical thermal signal with minimum temperature toward the end of the HS. The cave located at a depth of 500 m records a very flat temperature signal despite the existence of a 4°C thermal anomaly on the surface over the cave during the HS. The inset graph within Figure 3 shows the cave air temperature at 500 m in a different scale to display the amplitude and structure of the thermal anomaly related to the HS. Only the wave with the period of 2000 yr lasts long enough to be recorded at this depth, while all other waves are completely muted before reaching the cave, which explains why this anomaly looks like a simple sinusoidal signal. Despite the 4°C of anomaly outside the cave in relation to the HS, less than 0.1°C change is recorded in the cave located at 500 m underground. In addition, the lag time of the 2000 yr period anomaly exceeds 1000 yr at this depth, and consequently, the very small thermal anomaly recorded in the cave occurs when the HS outside the cave has already ended.

Figure 3. (color online) Temperature recorded underground in karst terrains at different depths in response to the surface atmosphere temperature (SAT) signal in an atmosphere–soil thermally coupled scenario. This scenario considers no change in vegetation cover over the cave during the Heinrich stadial. The inset graph shows the temperature record at the depth of 500 m with a different thermal scale to allow visualization of the lag time and structure of the signal.

The results of transferring the SAT shown in Figure 1 to the hypothetical caves under the second scenario (i.e., there is thermal decoupling of 2°C between SAT and ground temperature limited only to the HS) are shown in Figure 4. The temperature transferred underground is the ground-surface temperature and not the SAT. Consequently, changing environmental conditions such as the vegetal cover over the cave modifies the coupling conditions up to several degrees. This scenario assumes that a more open landscape develops over the cave during the HS interval compared with the earlier and later intervals defined in Figure 1. The change in landscape over the cave sites reduces the thermal anomaly to be transferred underground to half the value in the SAT. Apart from the amplitude, the structure of the signals is exactly the same as in the coupled scenario.

Figure 4. (color online) Temperature recorded underground in karst terrains at different depths in response to the surface atmosphere temperature (SAT) signal in an atmosphere–soil thermally decoupled scenario. This scenario considers a change in vegetation cover over the cave site during the Heinrich stadial. The inset graph shows the temperature record at the depth of 500 m with a different thermal scale to allow visualization of the lag time and structure of the signal.

Speleothem δ¹⁸O models

We use the results of the cave temperature records of both scenarios and the SAT record to evaluate the thermal impact on speleothem δ¹⁸O records as a result of the HS. A first event of isotope fractionation occurs in the atmosphere above the cave. The isotope models use four different δ¹⁸O_rw/SAT ratios that account for the impact of temperature during fractionation in the atmosphere, covering most of the variability found in temperate Mediterranean sites. Because the ratios are always positive, there is a direct correlation between δ¹⁸O_rw and SAT. The second fractionation event considered in this model is also temperature dependent and occurs at the moment of speleothem calcite precipitation. During the precipitation of calcite, the oxygen isotopes have a negative relationship with temperature. We considered this value to be constant due to the limited thermal variability in caves.

The impact of temperature in the cave after these two fractionation events is shown in Figure 5 for the first scenario (i.e., thermal coupling between SAT and ground temperature). The results show that the thermal changes related to HS under the boundary conditions constrained in this model have variable impact on caves at different depths. Both the cave depth and the δ¹⁸O_rw/SAT ratios are important parameters that determine the significance of the impact that temperature changes during HS have on speleothem δ¹⁸O records. In caves located 10 m from the surface and having a δ¹⁸O_rw/SAT ratio of +0.23‰/°C, the thermal impact on speleothem δ¹⁸O records is mostly indistinguishable from analytical uncertainties. At the other extreme, caves located 500 m underground and having a δ¹⁸O_rw/SAT ratio of +0.41‰/°C are capable of recording a thermal impact in their speleothem δ¹⁸O records of >1.6‰ during the HS. The attenuation of the temperature anomalies with depth causes thermal gradients with SAT that have an impact on the speleothem δ¹⁸O records. This impact is enhanced in deeper caves, as their temperature is less affected by changes in SAT. When significant speleothem δ¹⁸O anomalies related to the HS result from the model, the speleothem δ¹⁸O signals are generally shifted toward more negative values. The lag times related to the transfer of thermal signal underground also impact the structure of the isotope anomalies recorded in the speleothems, which is more obvious in the modeled caves located 50 and 100 m underground. During the HS, the speleothem δ¹⁸O anomaly is not constant in response to the progressive change in cave temperature. Additionally, at the end of the HS, the caves 50 and 100 m underground record a shift in their modeled speleothem δ¹⁸O records toward less-negative δ¹⁸O values. This positive δ¹⁸O anomaly is a product of the lag times and can be >0.5‰ at its peak, lasting several decades.

Figure 5. (color online) Impact of surface and underground temperature changes in calcite speleothem δ¹⁸O records (reported as δ¹⁸O anomalies) for caves at different depths (z) in an atmosphere–ground surface thermally coupled scenario. This scenario implies no changes in the vegetation cover over the cave as a result of the Heinrich stadial (HS). Isotope calculations consider that during the precipitation of calcite, fractionation occurs under equilibrium conditions. The model assumes that drip-water isotope composition only changes due to changes in surface atmosphere temperature (SAT). The model is implemented for four cases in which the relationship between rainwater δ¹⁸O values and the SAT (δ¹⁸O_rw/SAT ratio) during fractionation of oxygen isotopes in the atmosphere differs. (A) Case of δ¹⁸O_rw/SAT ratio = +0.23‰/°C; (B) case of δ¹⁸O_rw/SAT ratio = +0.29‰/°C; (C) case of δ¹⁸O_rw/SAT ratio = +0.35‰/°C; and (D) case of δ¹⁸O_rw/SAT ratio = +0.41‰/°C. The variability of the input signal of δ¹⁸O_rw anomaly in each of the four cases coincides with the variability of the δ¹⁸O_cc signal at the depth of 500 m, because there is no significant change in temperature at this depth as result of the changes in SAT during the HS.

The speleothem δ¹⁸O model under the second scenario (i.e., thermal decoupling of 2°C between SAT and ground temperature limited only to the HS) shows either remarkable differences or nearly identical results compared with the first scenario, depending on the depth of the cave (Fig. 6). The temperature of the cave located 500 m underground is almost unaffected by SAT or the changes in vegetation, and consequently, the outputs of the speleothem δ¹⁸O model are nearly identical to those of the first scenario. On the other hand, the temperature of the cave located 10 m underground changes the most as a result of the variable thermal decoupling, and its speleothem δ¹⁸O record does too. Under the variable decoupling scenario, the cave 10 m deep has negative speleothem δ¹⁸O anomalies ranging from 0.4 to 1.2‰ during the HS depending on the δ¹⁸O_rw/SAT ratio considered. The caves located 50 and 100 m underground also record slightly more-negative speleothem δ¹⁸O anomalies during the HS under this second scenario. The impact of lag times on the speleothem δ¹⁸O records is less obvious compared with the first scenario, because thermal anomalies are reduced. However, all the speleothem δ¹⁸O signal structures reported in the first scenario are also recorded here, despite their differences in amplitude.

Figure 6. (color online) Impact of surface and underground temperature changes in calcite speleothems δ¹⁸O records (reported as δ¹⁸O anomalies) for caves at different depths (z) in a variable atmosphere–ground surface thermally decoupled scenario. This scenario implies changes in the vegetation cover over the cave as a result of the Heinrich stadial (HS). Isotope calculations consider that during the precipitation of calcite, fractionation occurs under equilibrium conditions. The model assumes that drip-water isotope composition only changes due to changes in surface atmosphere temperature (SAT). The model is implemented for four cases in which the relationship between rainwater δ¹⁸O values and the SAT (δ¹⁸O_rw/SAT ratio) during fractionation of oxygen isotopes in the atmosphere differs. (A) Case of δ¹⁸O_rw/SAT ratio = +0.23‰/°C; (B) case of δ¹⁸O_rw/SAT ratio = +0.29‰/°C; (C) case of δ¹⁸O_rw/SAT ratio = +0.35‰/°C; and (D) case of δ¹⁸O_rw/SAT ratio = +0.41‰/°C. The variability of the input signal of δ¹⁸O_rw anomaly in each of the four cases coincides with the variability of the δ¹⁸O_cc signal at the depth of 500 m, because at this depth there is no significant change in temperature as result of the changes in SAT or vegetation cover during the HS.

DISCUSSION

Impact of millennial-scale SAT changes in caves

We have applied our thermal conduction model to two scenarios that consider coupling and variable decoupling of the SAT and ground temperature. There is a possible third scenario that we have not considered yet, in which the ground records a systematic temperature shift compared with SAT. This scenario implies that there is a temperature difference between the ground and the atmosphere, but the long-term mean of this value does not change. The thermal variability under this constantly decoupled scenario will be identical to that of the coupled scenario, although the absolute values will differ as much as the magnitude of the decoupling. The same absolute differences will be transferred underground to the cave, whereas the cave temperature variability will be exactly the same as in the coupled scenario. This third scenario is likely to occur, for example, in caves under a rocky terrain in which trees or shrubs are unlikely to grow independent of the climate conditions. In the second scenario, we have assumed a variable thermal decoupling to account for changes in the landscape over the cave. Vegetation cover is a major control on the coupling conditions between SAT and ground-surface temperature (Lewis and Wang, 1998). However, other controls such as snow cover, irradiation, or evaporation are potentially also significant (Beltrami and Kellman, 2003; Smerdon, et al., 2006; Yazaki et al., 2013). We have selected a decoupling value of 2°C, which is less than the annually averaged thermal difference between SAT and the ground temperature under Mediterranean forest and grassland (Domínguez-Villar et al., 2013b). This muted decoupling value was selected to account for the additional impact of other potential controls such as the increase of snow cover and the decreased evaporation during the HS or the diminished irradiance during glacial times compared with the present.

The model shows that the temperature of cave interiors can differ by several degrees from the long-term SAT. These long-term thermal gradients are more pronounced in deeper caves. This creates a thermal gradient between the cave and the external atmosphere that affects the cave ventilation by changing the gradient of air densities (de Freitas et al., 1982; Kowalczk and Froelich, 2010). This gradient is not constant but fluctuates, as we have to consider that our thermal model has filtered out the seasonality. Therefore, every year the SAT is expected to oscillate above and below the cave temperature, forcing different ventilation modes. Nevertheless, this long-term thermal gradient is expected to affect the intensity and duration of seasonal ventilation modes. The seasonal period of enhanced ventilation that occurs when SAT is below the cave interior temperature is expected to last longer and be more vigorous as a consequence of cave temperature being warmer than the mean annual SAT. On the other hand, the season of diminished ventilation is expected to be shortened. The extended seasonal period of enhanced ventilation caused by the long-term warmer cave temperature is expected to favor the growth rate of speleothems (Banner et al., 2007; Boch et al., 2011). However, other key controls on speleothem growth rate are drip rate and calcium activity in the drip solution (Baker et al., 1998), which under dry conditions are likely to counteract or even surpass the impact of enhanced ventilation on speleothem growth rate. So, despite the expected changes in ventilation in some caves in relation to HS, no general pattern on the changes in speleothem growth rate can be predicted, as local conditions may determine the dominant control.

Cave microclimate and underground ecological environments can also be potentially affected by the extended and more dynamic seasonal ventilation during HSs. A significantly drier climate may impact the water availability in the cave, which together with the enhanced advection due to the more dynamic ventilation, is potentially capable of decreasing seasonally the relative humidity of the cave interior. When cave atmosphere is not close to saturation in water (i.e., relative humidity clearly below 100%), clay and other particles are easier to move as aerosols, and dust can be observed in the cave atmosphere (Dredge et al., 2013) and potentially recorded in the trace elements of speleothems. The occurrence of moist–dry cycles on cave walls as a result of unsaturated cave air in water can favor the conditions for weathering of carbonates (Karbowska-Berent, 2003; Dreybrodt et al., 2005; Krklec et al., 2016) and impact the cave ecological niches, having important consequences for the preservation of cave art (Saiz-Jimenez et al., 2011). Evaporation can cause significant changes in growth rate by affecting the saturation index that impacts calcite fabrics and trace elements, as well as controlling the speleothem δ¹⁸O and δ¹³C values (Frisia et al., 2000; Mickler et al., 2004; Fairchild and Treble, 2009). Although the water availability in caves has a critical local component, caves under arid and semiarid climates are more likely candidates to be affected by evaporation as a result of HS.

Impact of SAT changes in speleothem isotope records

We used four different constant δ¹⁸O_rw/SAT ratios in our isotope models. The assumption that a site will have a constant δ¹⁸O_rw/SAT ratio through time is likely unrealistic. However, little is known yet concerning the long-term changes in the δ¹⁸O_rw/SAT relationship (Fricke and O'Neil, 1999). Our assumption is still a good approximation to evaluate potential impacts of HS on speleothem δ¹⁸O values, keeping in mind the boundary conditions of the model and its limitations when doing interpretations.

We mentioned in the previous section that enhanced ventilation could favor evaporation and consequently affect the speleothem δ¹⁸O and δ¹³C values. However, in the isotope model we present, we are assuming that the caves maintain a saturated atmosphere, preventing the occurrence of kinetic fractionation during calcite precipitation that would affect the speleothem δ¹⁸O record. Yet the long-term thermal gradient between the cave and the external atmosphere is expected to affect the ventilation modes, potentially affecting the speleothem δ¹³C record. The degassing of CO₂ from the dissolved inorganic carbon in the solution can be a major control on the speleothem δ¹³C values for some speleothems (Frisia et al., 2011). In those cases, a shift in the speleothem δ¹³C record is expected as a result of the thermal changes of the HS independent of the changes in the vegetation cover, soil properties on top of the cave, or the system's carbonate dissolution style. However, the controls affecting speleothem δ¹³C are site specific, and degassing is not always a major contributor to the speleothem δ¹³C variability (Mattey et al., 2008; Dreybrodt and Scholz, 2011; Tremaine et al., 2011; Griffiths, et al., 2012). Therefore, although the long-term thermal gradient caused by the HS has the potential to affect speleothem δ¹³C records, this effect is not expected to have a significant contribution to the δ¹³C variability of every speleothem.

The results of our isotope model suggest that under equilibrium conditions, the isolated impact of temperature on speleothem δ¹⁸O values can be >1.6‰. As many factors affect the δ¹⁸O_rw values (Lachniet, 2009), during HS the changes in the hydrological cycle that impact δ¹⁸O_rw values are not limited to temperature changes. However, our model is not designed to represent the real anomalies that speleothems would record, but the anomalies that could only be attributed to the thermal impact. So the impact of temperature on speleothem δ¹⁸O records will be significant only when other controls taken all together will not counteract the thermal effect. Because the δ¹⁸O variability associated with HSs and similar cold episodes in temperate Mediterranean speleothems is often <2‰ (Drysdale et al., 2007; Fleitmann et al., 2009; Moreno et al., 2010; Domínguez-Villar et al., 2013a), it is clear that the temperature effect is likely to be a significant control on the speleothem δ¹⁸O records at most sites that share the boundary conditions used in our model.

CONCLUSION

This research has developed a one-dimensional thermal conduction model to show how millennial-scale temperature drops, such as those recorded during HSs in Atlantic and Mediterranean sites, are transferred to caves. The model shows that the depth of the cave and the duration of the anomaly are important controls for the thermal anomaly to be recorded underground. Caves hundreds of meters deep can record little evidence of the HSs temperature changes recorded at the surface. On the other hand, shallow caves are more likely to record significant thermal changes during HSs. A thermal anomaly is progressively attenuated with the increasing depth of the cave, but the recorded thermal signal is also affected by increased delays and more drastic modifications in its structure as a function of depth. The results provided here support that not all caves are in thermal equilibrium with the exterior atmosphere during HSs.

We have also discussed the impact that this thermal disequilibrium may have on different cave dynamics and speleothem proxies. The results support the occurrence of long-term thermal gradients between the caves and the exterior as a consequence of HSs. These gradients are more pronounced in deep caves and less significant in shallow caves. Thermal gradients are likely to cause modifications on the ventilation regime in caves and affect multiple processes, some of them being recorded in speleothems. Although changes in ventilation can affect multiple speleothem proxies, local controls are too important to expect a unitary response in cave systems. However, we can evaluate the impact of temperature during HSs on speleothem δ¹⁸O records. The isotope models here reported show that for most temperate and Mediterranean cave sites, the thermal impact of HSs cannot be neglected when interpreting speleothem δ¹⁸O records. Therefore, despite the anticipated important modifications of the hydrological cycle affecting δ¹⁸O_rw values during HS, Pleistocene speleothem δ¹⁸O values still record significant changes in relation to temperature alone.

This research has put the focus once more on temperature as a significant control of speleothem isotope records during Pleistocene climate changes. Also, we have demonstrated how the depth of the cave can impact its dynamics and speleothem records. Therefore, we encourage scientists working in paleoclimate research to always report the depths of the caves they are studying.