# Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves

Amir Sedaghatkish, Claudio Pastore, Frédéric Doumenc, Pierre-yves Jeannin, Marc Luetscher

► To cite this version:

Amir Sedaghatkish, Claudio Pastore, Frédéric Doumenc, Pierre-yves Jeannin, Marc Luetscher. Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves. Journal of Geophysical Research: Earth Surface, 2024, 129 (6), pp.e2024JF007646. 10.1029/2024JF007646. hal-04592106

HAL Id: hal-04592106

https://hal.science/hal-04592106

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves

Amir Sedaghatkish, Claudio Pastore, Frédéric Doumenc, Pierre-yves Jeannin, Marc Luetscher

► To cite this version:

Amir Sedaghatkish, Claudio Pastore, Frédéric Doumenc, Pierre-yves Jeannin, Marc Luetscher. Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves. Journal of Geophysical Research: Earth Surface, 2024, 129 (6), pp.e2024JF007646. 10.1029/2024JF007646. hal-04592106

HAL Id: hal-04592106

https://hal.science/hal-04592106v1

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Earth Surface

## RESEARCH ARTICLE

10.1029/2024JF007646

**Key Points:**

* A numerical model was developed to investigate heat transfer inside ventilated caves driven by chimney effect

* The airflow can propagate annual temperature fluctuations over a few tens to a few hundreds of meters from the cave entrance (over kilometers in extreme cases)

* The distance of propagation (the convection length) is approximately proportional to the airflow rate and inversely proportional to the square root of the cave diameter

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Sedaghatkish,
amir.sedaghatkish@isska.ch

**Citation:**

Sedaghatkish, A., Pastore, C., Doumenc, F., Jeannin, P.-Y., & Luetscher, M. (2024). Modeling heat transfer for assessing the convection length in ventilated caves. *Journal of Geophysical Research: Earth Surface, 129*, e2024JF007646. https://doi.org/10.1029/2024JF007646

Received 16 JAN 2024
Accepted 13 MAY 2024

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves

Amir Sedaghatkish¹,² ID, Claudio Pastore¹,² ID, Frédéric Doumenc³,⁴ ID, Pierre-Yves Jeannin¹, and Marc Luetscher¹ ID

¹Swiss Institute for Speleology and Karst Studies (SISKA), La Chaux-de-Fonds, Switzerland, ²Center for Hydrogeology and Geothermics (CHYN), University of Neuchâtel, Neuchâtel, Switzerland, ³Université Paris-Saclay, CNRS, FAST, Orsay, France, ⁴Sorbonne Université, UFR 919, Paris Cedex 05, France

**Abstract** The present study focuses on heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. We use a numerical model that couples the convective heat transfer due to the airflow in a single karst conduit with the conductive heat transfer in the rock mass. Assuming dry air and a simplified geometry, we investigate the propagation of thermal perturbations inside the karst massif. We perform a parametric study to identify general trends regarding the effect of the air flowrate and conduit size on the amplitude and spatial extent of thermal perturbations. Numerical results support the partition of a cave into three regions: (a) a short (few meters) diffusive region, where heat mainly propagates from the external atmosphere by conduction in the rock mass; (b) a convective region where heat is mainly transported by the air flow; (c) a deep karst region characterized by quasi-constant temperatures throughout the year. Numerical simulations show that the length of the convective region is approximately proportional to the amplitude of the flowrate annual fluctuations divided by the square root of the cave radius. This result is tested against field data from a mine tunnel and two caves. Our study provides first estimates to identify climate sensitive regions for speleothem science and/or ecosystemic studies.

**Plain Language Summary** Karsts are landscapes formed from the dissolution of soluble rocks. The chemical erosion due to rainwater results in the formation of an extensive network of caves traversed by air and water flows. Understanding heat transfer in karst is a key issue in many fields, as diverse as study of underground biota or paleoclimate reconstruction. The present work focuses on heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. In winter, the air inside the massif is hotter (and thus lighter) than the air outside. This situation results in an upward airflow as in a chimney. The airflow is reversed in summer (colder and thus heavier air inside the massif). Assuming a simplified geometry, we investigate by numerical simulations the propagation of thermal perturbations inside a single ventilated conduit. The results of the numerical simulations are tested against field data from a mine and two caves. We show that in ventilated caves, the airflow can propagate the annual temperature fluctuations over distances from the entrance ranging from a few tens to a few hundreds of meters. In extreme cases, it could go up to kilometers. Our study provides first estimates to identify climate sensitive regions for speleothem science and/or ecosystemic studies.

## 1. Introduction

Understanding heat transfer in karst systems is a key issue for underground biota (Mammola et al., 2019), preservation of cave art (Bourges et al., 2014), speleothem growth rates (Banner et al., 2007; Spötl et al., 2005), or paleoclimate reconstruction (Borsato et al., 2016; Casteel & Banner, 2015; Domínguez-Villar et al., 2021). However, heat transfer in karst results from an intricate coupling between several mechanisms, including heat conduction in the rock mass (Quindos et al., 1987), convection due to water or air flow in caves (Cropley, 1965), or radiative transfer between cave walls (Guerrier et al., 2019). Latent heat exchanges are also present because of evaporation and condensation on cave walls (Dreybrodt et al., 2005) or ice formation and permafrost (Luetscher et al., 2008). A key issue for the understanding of heat transfer in karst is to determine the processes that dominate in a given configuration, and those that can be neglected.

The simplest heat transfer model of a vadose karst assumes 1D conduction heat transfer in a rock mass of infinite depth with periodic temperature fluctuations at the ground surface (Domínguez-Villar et al., 2023; Salmon et al., 2023; Villar et al., 1983). This model predicts that the depth corresponding to 99% attenuation of the
# JGR Earth Surface

## RESEARCH ARTICLE

10.1029/2024JF007646

**Key Points:**

* A numerical model was developed to investigate heat transfer inside ventilated caves driven by chimney effect

* The airflow can propagate annual temperature fluctuations over a few tens to a few hundreds of meters from the cave entrance (over kilometers in extreme cases)

* The distance of propagation (the convection length) is approximately proportional to the airflow rate and inversely proportional to the square root of the cave diameter

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Sedaghatkish,
amir.sedaghatkish@isska.ch

**Citation:**

Sedaghatkish, A., Pastore, C., Doumenc, F., Jeannin, P.-Y., & Luetscher, M. (2024). Modeling heat transfer for assessing the convection length in ventilated caves. *Journal of Geophysical Research: Earth Surface, 129*, e2024JF007646. https://doi.org/10.1029/2024JF007646

Received 16 JAN 2024
Accepted 13 MAY 2024

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves

Amir Sedaghatkish¹,² ID, Claudio Pastore¹,² ID, Frédéric Doumenc³,⁴ ID, Pierre-Yves Jeannin¹, and Marc Luetscher¹ ID

¹Swiss Institute for Speleology and Karst Studies (SISKA), La Chaux-de-Fonds, Switzerland, ²Center for Hydrogeology and Geothermics (CHYN), University of Neuchâtel, Neuchâtel, Switzerland, ³Université Paris-Saclay, CNRS, FAST, Orsay, France, ⁴Sorbonne Université, UFR 919, Paris Cedex 05, France

**Abstract** The present study focuses on heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. We use a numerical model that couples the convective heat transfer due to the airflow in a single karst conduit with the conductive heat transfer in the rock mass. Assuming dry air and a simplified geometry, we investigate the propagation of thermal perturbations inside the karst massif. We perform a parametric study to identify general trends regarding the effect of the air flowrate and conduit size on the amplitude and spatial extent of thermal perturbations. Numerical results support the partition of a cave into three regions: (a) a short (few meters) diffusive region, where heat mainly propagates from the external atmosphere by conduction in the rock mass; (b) a convective region where heat is mainly transported by the air flow; (c) a deep karst region characterized by quasi-constant temperatures throughout the year. Numerical simulations show that the length of the convective region is approximately proportional to the amplitude of the flowrate annual fluctuations divided by the square root of the cave radius. This result is tested against field data from a mine tunnel and two caves. Our study provides first estimates to identify climate sensitive regions for speleothem science and/or ecosystemic studies.

**Plain Language Summary** Karsts are landscapes formed from the dissolution of soluble rocks. The chemical erosion due to rainwater results in the formation of an extensive network of caves traversed by air and water flows. Understanding heat transfer in karst is a key issue in many fields, as diverse as study of underground biota or paleoclimate reconstruction. The present work focuses on heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. In winter, the air inside the massif is hotter (and thus lighter) than the air outside. This situation results in an upward airflow as in a chimney. The airflow is reversed in summer (colder and thus heavier air inside the massif). Assuming a simplified geometry, we investigate by numerical simulations the propagation of thermal perturbations inside a single ventilated conduit. The results of the numerical simulations are tested against field data from a mine and two caves. We show that in ventilated caves, the airflow can propagate the annual temperature fluctuations over distances from the entrance ranging from a few tens to a few hundreds of meters. In extreme cases, it could go up to kilometers. Our study provides first estimates to identify climate sensitive regions for speleothem science and/or ecosystemic studies.

## 1. Introduction

Understanding heat transfer in karst systems is a key issue for underground biota (Mammola et al., 2019), preservation of cave art (Bourges et al., 2014), speleothem growth rates (Banner et al., 2007; Spötl et al., 2005), or paleoclimate reconstruction (Borsato et al., 2016; Casteel & Banner, 2015; Domínguez-Villar et al., 2021). However, heat transfer in karst results from an intricate coupling between several mechanisms, including heat conduction in the rock mass (Quindos et al., 1987), convection due to water or air flow in caves (Cropley, 1965), or radiative transfer between cave walls (Guerrier et al., 2019). Latent heat exchanges are also present because of evaporation and condensation on cave walls (Dreybrodt et al., 2005) or ice formation and permafrost (Luetscher et al., 2008). A key issue for the understanding of heat transfer in karst is to determine the processes that dominate in a given configuration, and those that can be neglected.

The simplest heat transfer model of a vadose karst assumes 1D conduction heat transfer in a rock mass of infinite depth with periodic temperature fluctuations at the ground surface (Domínguez-Villar et al., 2023; Salmon et al., 2023; Villar et al., 1983). This model predicts that the depth corresponding to 99% attenuation of the
temperature fluctuations at the ground surface is a few tens of centimeters for the daily fluctuations and approximately 10 m for the annual fluctuations. Quindos and co-workers (Quindos et al., 1987; Villar et al., 1983) measured the monthly mean temperatures in Altamira Cave during a year. They found that the periodic 1D conduction model correctly predicts the amplitude and the phase shift of the annual fluctuations at the cave roof whose depth varied from 3.5 to 17.5 m. In contrast, the floor temperature was much closer to the roof temperature than expected from a 1D conduction model. The authors attributed this effect to the homogenization of the temperature field by radiative transfer between roofs and floors. Guerrier et al. (2019) and Qaddah et al. (2023) confirmed by numerical simulations the significance of radiative transfer inside weakly ventilated shallow caves, and showed that turbulent free convection also had to be considered. Both heat transfer mechanisms contribute to the homogenization of the wall temperature and the deformation of the temperature field in the rock matrix around the cave.

The situation is different when forced convection is present. It has long been known that, in large cave systems, the convective transport of heat by air currents and streams can propagate thermal perturbations over much larger distances than predicted by a pure conduction model. Cropley (1965) found that temperature measurements performed in two caves in West Virginia supported the definition of three zones: zone I immediately adjacent to the cave entrance, where the cave temperature follows the surface temperature fluctuations; zone II where the cave temperature is driven by streams and air currents; zone III where the temperature is approximately constant and close to the annual mean temperature (AMT) at the surface. Areas where the amplitude of temperature fluctuations was lower than 1°F (ca. 0.56°C) around the AMT were located at more than 1,500 m from the cave entrances. This long distance was mainly due to cold water carried into the caves by winter water recharge. Luetscher and Jeannin (2004) defined the heterothermic vadose zone as the surficial part of a karst massif showing seasonal inversions of the temperature gradient, in contrast with the homothermic vadose zone characterized by a temperature gradient close to the lapse rate of humid air (ca. 0.5°C 100 m⁻¹). These authors estimated the depth of the heterothermic zone at about 100 m in most cases.

Considerable work has been undertaken to simulate heat transfer by forced convection due to water flow in the vadose and phreatic zones of karst (Gong et al., 2019; Long & Gilcrease, 2009; Sinokrot & Stefan, 1993). Covington et al. (2011) showed that the relative significance of different heat fluxes including convection and conduction are timescale dependent. Conduction through the rock surrounding a conduit determines heat flux at times of the order of weeks and longer. Coupling convection in water with conduction in rock is thus necessary to get a realistic model valid at all time scales.

Forced convection due to airflow must be considered when a network of karst conduits has at least two entrances. In this configuration, several mechanisms can contribute to the production of airflow through the cave, resulting in significant convective heat transfer between the cave and the external atmosphere. Potential driving forces are barometric fluctuations (Gomell et al., 2021), dynamic pressure effect due to external winds (Kukuljan et al., 2021), diphasic flow due to water circulation (Atkinson, 1977), or the buoyancy due to the density contrast between the air inside and outside of the cave (Gabrovšek, 2023). The latter mechanism is likely the most significant in most cases. The air density depends on moisture content, CO₂ concentration and temperature. In temperate climates, temperature fluctuations are the main cause of density variations (Gabrovšek, 2023). Therefore, the flow direction is upward when the temperature inside the cave is warmer (or downward when inside the cave is colder) than the external atmosphere. Because of the thermal inertia of the massif, a lower entrance operates as an inlet of cold external air during most of the winter, whereas it acts as an outlet during most of the summer. In the latter case, the air cooled down when crossing the massif before reaching the lower entrance. The lower entrance thus receives in summer an airflow colder than the external temperature. Because of this seasonal asymmetry, the AMT is shifted to colder values at lower entrances (cold thermal anomaly). The same mechanism produces a shift of the AMT to hotter values at upper entrances (hot thermal anomaly) (Lismonde, 2002). In alpine ice caves, where this ventilation pattern represents the normal unless the cave is clogged with sediments, these thermal anomalies usually extend over a few hundred meters from the cave entrances (Luetscher et al., 2008).

Wigley and Brown (1971) calculated the temperature and moisture content profiles in the air as a function of the distance from the entrance. They used a 1D model based on the energy and water mass conservation in the airflow. The convective heat transfer coefficient between the air and the wall was estimated using an empirical correlation valid for forced convection in smooth pipes. The model assumes prescribed uniform wall temperature, a major
temperature fluctuations at the ground surface is a few tens of centimeters for the daily fluctuations and approximately 10 m for the annual fluctuations. Quindos and co-workers (Quindos et al., 1987; Villar et al., 1983) measured the monthly mean temperatures in Altamira Cave during a year. They found that the periodic 1D conduction model correctly predicts the amplitude and the phase shift of the annual fluctuations at the cave roof whose depth varied from 3.5 to 17.5 m. In contrast, the floor temperature was much closer to the roof temperature than expected from a 1D conduction model. The authors attributed this effect to the homogenization of the temperature field by radiative transfer between roofs and floors. Guerrier et al. (2019) and Qaddah et al. (2023) confirmed by numerical simulations the significance of radiative transfer inside weakly ventilated shallow caves, and showed that turbulent free convection also had to be considered. Both heat transfer mechanisms contribute to the homogenization of the wall temperature and the deformation of the temperature field in the rock matrix around the cave.

The situation is different when forced convection is present. It has long been known that, in large cave systems, the convective transport of heat by air currents and streams can propagate thermal perturbations over much larger distances than predicted by a pure conduction model. Cropley (1965) found that temperature measurements performed in two caves in West Virginia supported the definition of three zones: zone I immediately adjacent to the cave entrance, where the cave temperature follows the surface temperature fluctuations; zone II where the cave temperature is driven by streams and air currents; zone III where the temperature is approximately constant and close to the annual mean temperature (AMT) at the surface. Areas where the amplitude of temperature fluctuations was lower than 1°F (ca. 0.56°C) around the AMT were located at more than 1,500 m from the cave entrances. This long distance was mainly due to cold water carried into the caves by winter water recharge. Luetscher and Jeannin (2004) defined the heterothermic vadose zone as the surficial part of a karst massif showing seasonal inversions of the temperature gradient, in contrast with the homothermic vadose zone characterized by a temperature gradient close to the lapse rate of humid air (ca. 0.5°C 100 m⁻¹). These authors estimated the depth of the heterothermic zone at about 100 m in most cases.

Considerable work has been undertaken to simulate heat transfer by forced convection due to water flow in the vadose and phreatic zones of karst (Gong et al., 2019; Long & Gilcrease, 2009; Sinokrot & Stefan, 1993). Covington et al. (2011) showed that the relative significance of different heat fluxes including convection and conduction are timescale dependent. Conduction through the rock surrounding a conduit determines heat flux at times of the order of weeks and longer. Coupling convection in water with conduction in rock is thus necessary to get a realistic model valid at all time scales.

Forced convection due to airflow must be considered when a network of karst conduits has at least two entrances. In this configuration, several mechanisms can contribute to the production of airflow through the cave, resulting in significant convective heat transfer between the cave and the external atmosphere. Potential driving forces are barometric fluctuations (Gomell et al., 2021), dynamic pressure effect due to external winds (Kukuljan et al., 2021), diphasic flow due to water circulation (Atkinson, 1977), or the buoyancy due to the density contrast between the air inside and outside of the cave (Gabrovšek, 2023). The latter mechanism is likely the most significant in most cases. The air density depends on moisture content, CO₂ concentration and temperature. In temperate climates, temperature fluctuations are the main cause of density variations (Gabrovšek, 2023). Therefore, the flow direction is upward when the temperature inside the cave is warmer (or downward when inside the cave is colder) than the external atmosphere. Because of the thermal inertia of the massif, a lower entrance operates as an inlet of cold external air during most of the winter, whereas it acts as an outlet during most of the summer. In the latter case, the air cooled down when crossing the massif before reaching the lower entrance. The lower entrance thus receives in summer an airflow colder than the external temperature. Because of this seasonal asymmetry, the AMT is shifted to colder values at lower entrances (cold thermal anomaly). The same mechanism produces a shift of the AMT to hotter values at upper entrances (hot thermal anomaly) (Lismonde, 2002). In alpine ice caves, where this ventilation pattern represents the normal unless the cave is clogged with sediments, these thermal anomalies usually extend over a few hundred meters from the cave entrances (Luetscher et al., 2008).

Wigley and Brown (1971) calculated the temperature and moisture content profiles in the air as a function of the distance from the entrance. They used a 1D model based on the energy and water mass conservation in the airflow. The convective heat transfer coefficient between the air and the wall was estimated using an empirical correlation valid for forced convection in smooth pipes. The model assumes prescribed uniform wall temperature, a major
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Cave geometry (the computational domain is displayed in red). The rock massif is assumed infinite in both vertical directions.</figcaption></figure>

simplification that allows to obtain simple closed-form expressions for the temperature and moisture profiles. However, this approach implicitly assumes that the heat flux through the cave wall is convection-limited, and that the rock mass does not play any role. Lismonde (2002) pointed out that the airflow modifies the rock temperature, which in turn changes the air temperature profile in the cave. This author developed a model to predict heat transfer in a straight inclined duct of constant diameter included in a rock massif. The model coupled 1D radial conduction in the rock mass with convective transfer in the gas, using a sinusoidal function of time as inlet temperature. The airflow rate was predicted by considering the interaction between the temperature field in the gas and the airflow through the buoyancy term in the momentum balance equation of the gas. This model reproduced qualitatively some field observations, for example, the hysteresis of the flow rate or the thermal anomaly in the entrance areas. Gabrovšek (2023) used a similar model to investigate by numerical simulations the effect of conduit shape on the airflow pattern in ventilated caves.

The present study focuses on the numerical simulation of heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. Our objective is to go beyond the qualitative results obtained from previous numerical studies available in the literature. We aim at assessing orders of magnitudes of the spatial extent of the karst region where the airflow induces significant thermal perturbations. We use for that a numerical model close to that already developed by Lismonde (2002) or Gabrovšek (2023), and we apply it to the simplest possible configuration in order to quantify general trends common to most ventilated caves. The convective heat transfer due to dry airflow in a single horizontal karst conduit of constant diameter is coupled with the conductive heat transfer in the impermeable rock mass.

This article contains a first theoretical part (Sections 2–5) dedicated to numerical simulations based on a model including many simplifying assumptions. The relevance of these assumptions to achieve our objective, which consists in the assessment of general orders of magnitude, is evaluated in a second part (Section 6) by comparison with the field data obtained from three sites, a mine tunnel and two caves.

## 2. Model Definition

### 2.1. Cave Geometry and Computational Domain

We consider the idealized ventilated cave displayed in Figure 1, located in a rock massif of infinite extent in both vertical directions. The cave consists of a vertical conduit connected by two horizontal conduits to two entrances at different elevations. The rock mass is assumed impermeable. The horizontal conduits are supposed long
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Cave geometry (the computational domain is displayed in red). The rock massif is assumed infinite in both vertical directions.</figcaption></figure>

simplification that allows to obtain simple closed-form expressions for the temperature and moisture profiles. However, this approach implicitly assumes that the heat flux through the cave wall is convection-limited, and that the rock mass does not play any role. Lismonde (2002) pointed out that the airflow modifies the rock temperature, which in turn changes the air temperature profile in the cave. This author developed a model to predict heat transfer in a straight inclined duct of constant diameter included in a rock massif. The model coupled 1D radial conduction in the rock mass with convective transfer in the gas, using a sinusoidal function of time as inlet temperature. The airflow rate was predicted by considering the interaction between the temperature field in the gas and the airflow through the buoyancy term in the momentum balance equation of the gas. This model reproduced qualitatively some field observations, for example, the hysteresis of the flow rate or the thermal anomaly in the entrance areas. Gabrovšek (2023) used a similar model to investigate by numerical simulations the effect of conduit shape on the airflow pattern in ventilated caves.

The present study focuses on the numerical simulation of heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere. Our objective is to go beyond the qualitative results obtained from previous numerical studies available in the literature. We aim at assessing orders of magnitudes of the spatial extent of the karst region where the airflow induces significant thermal perturbations. We use for that a numerical model close to that already developed by Lismonde (2002) or Gabrovšek (2023), and we apply it to the simplest possible configuration in order to quantify general trends common to most ventilated caves. The convective heat transfer due to dry airflow in a single horizontal karst conduit of constant diameter is coupled with the conductive heat transfer in the impermeable rock mass.

This article contains a first theoretical part (Sections 2–5) dedicated to numerical simulations based on a model including many simplifying assumptions. The relevance of these assumptions to achieve our objective, which consists in the assessment of general orders of magnitude, is evaluated in a second part (Section 6) by comparison with the field data obtained from three sites, a mine tunnel and two caves.

## 2. Model Definition

### 2.1. Cave Geometry and Computational Domain

We consider the idealized ventilated cave displayed in Figure 1, located in a rock massif of infinite extent in both vertical directions. The cave consists of a vertical conduit connected by two horizontal conduits to two entrances at different elevations. The rock mass is assumed impermeable. The horizontal conduits are supposed long
enough for the vertical conduit to be fully included in the homothermic area where the temperature gradient reduces to the adiabatic lapse rate. The aeraulic and thermal problems are uncoupled in this simplified configuration, since the temperature field in the horizontal conduits have no effect on buoyancy. The air flowrate thus only depends on the temperature contrast between the temperature in the vertical conduit, independent of time, and the temperature of the atmosphere outside the cave, a known function of time.

Another source of simplification is that the thermal problems in the regions of the upper and lower entrances are uncoupled. They can thus be treated separately. We arbitrarily focus on the upper entrance, but all the results can be easily translated to a lower entrance (the only difference is that the air flows through the upper entrance inward during summer and outward during winter, and vice-versa through the lower entrance).

The computational domain is a cylinder of length $L_{dom}$ and outer radius $R_{dom}$ displayed in red in Figure 1. It contains a conduit of same axis and length and constant radius $R \ll R_{dom}$, inside which the air circulates. The length $L_{dom}$ and outer radius $R_{dom}$ are set so that the computational domain includes the whole thermal perturbation induced in the rock mass by the airflow. Practically, in a given configuration, $L_{dom}$ and $R_{dom}$ are increased until they no longer influence the results.

## 2.2. Atmosphere Temperature and Mass Flow Rate

For the sake of simplification, only the annual fluctuation of the atmosphere temperature $T_{atm}$ is considered. It is thus assumed that $T_{atm}(t)$ follows the simple law:

$$T_{atm}(t) = T_m + \Delta T_y \sin\left(2\pi \frac{t}{\tau_y}\right) \quad (1)$$

where $T_m$ and $\Delta T_y$ are the AMT and the amplitude of the annual temperature fluctuation (ATF), respectively, in the external atmosphere at the elevation of the upper entrance. The corresponding period is $\tau_y = 1$ year.

The air mass flow rate $\dot{m}$ is deduced from the momentum balance applied to the air inside the cave. Assuming negligible inertia, the balance between friction and buoyancy yields (Lismonde, 2002):

$$K \dot{m} |\dot{m}| = -(\bar{\rho}_{atm} - \bar{\rho}_m) gH, \quad (2)$$

where $H$ is the length of the vertical conduit and $g$ is the gravitational acceleration. The cave aeraulic resistance $K$ can be assumed to be constant in the turbulent regime. $\bar{\rho}_{atm}$ and $\bar{\rho}_m$ are the densities of the air outside the cave and inside the homothermic zone, respectively, averaged over the cave height. Assuming that: (a) the difference of densities depends only on the temperature drop between the external atmosphere and the homothermic zone, and (b) the temperature field in the vertical conduit follows the AMT in the external atmosphere, we get after linearization:

$$\bar{\rho}_{atm}(t) - \bar{\rho}_m = -\frac{M_a P_{atm}}{R_g T_m^2} (T_{atm}(t) - T_m) \quad (3)$$

where $P_{atm}$ is the atmospheric pressure, $M_a$ the molar mass of air and $R_g$ the ideal gas constant. Equations 2 and 3 show that the air flowrate is proportional to the square root of the temperature drop ($T_{atm} - T_m$). Injecting Equations 1 and 3 in Equation 2 yields

$$\dot{m}(t) = S \Delta \dot{m} \sqrt{\left| \sin\left(\frac{2\pi}{\tau_y} t\right) \right|} \quad (4)$$

where $S = -1$ for $T_{atm} < T_m$ (wintertime, upward flow) and $S = 1$ for $T_{atm} > T_m$ (summertime, downward flow). The positive constant $\Delta \dot{m}$ is the amplitude of the annual fluctuation of the flow rate. It is an input of the model.

## 2.3. Governing Equations

The air temperature verifies the energy balance
enough for the vertical conduit to be fully included in the homothermic area where the temperature gradient reduces to the adiabatic lapse rate. The aeraulic and thermal problems are uncoupled in this simplified configuration, since the temperature field in the horizontal conduits have no effect on buoyancy. The air flowrate thus only depends on the temperature contrast between the temperature in the vertical conduit, independent of time, and the temperature of the atmosphere outside the cave, a known function of time.

Another source of simplification is that the thermal problems in the regions of the upper and lower entrances are uncoupled. They can thus be treated separately. We arbitrarily focus on the upper entrance, but all the results can be easily translated to a lower entrance (the only difference is that the air flows through the upper entrance inward during summer and outward during winter, and vice-versa through the lower entrance).

The computational domain is a cylinder of length $L_{dom}$ and outer radius $R_{dom}$ displayed in red in Figure 1. It contains a conduit of same axis and length and constant radius $R \ll R_{dom}$, inside which the air circulates. The length $L_{dom}$ and outer radius $R_{dom}$ are set so that the computational domain includes the whole thermal perturbation induced in the rock mass by the airflow. Practically, in a given configuration, $L_{dom}$ and $R_{dom}$ are increased until they no longer influence the results.

## 2.2. Atmosphere Temperature and Mass Flow Rate

For the sake of simplification, only the annual fluctuation of the atmosphere temperature $T_{atm}$ is considered. It is thus assumed that $T_{atm}(t)$ follows the simple law:

$$T_{atm}(t) = T_m + \Delta T_y \sin\left(2\pi \frac{t}{\tau_y}\right) \quad (1)$$

where $T_m$ and $\Delta T_y$ are the AMT and the amplitude of the annual temperature fluctuation (ATF), respectively, in the external atmosphere at the elevation of the upper entrance. The corresponding period is $\tau_y = 1$ year.

The air mass flow rate $\dot{m}$ is deduced from the momentum balance applied to the air inside the cave. Assuming negligible inertia, the balance between friction and buoyancy yields (Lismonde, 2002):

$$K \dot{m} |\dot{m}| = -(\bar{\rho}_{atm} - \bar{\rho}_m) gH, \quad (2)$$

where $H$ is the length of the vertical conduit and $g$ is the gravitational acceleration. The cave aeraulic resistance $K$ can be assumed to be constant in the turbulent regime. $\bar{\rho}_{atm}$ and $\bar{\rho}_m$ are the densities of the air outside the cave and inside the homothermic zone, respectively, averaged over the cave height. Assuming that: (a) the difference of densities depends only on the temperature drop between the external atmosphere and the homothermic zone, and (b) the temperature field in the vertical conduit follows the AMT in the external atmosphere, we get after linearization:

$$\bar{\rho}_{atm}(t) - \bar{\rho}_m = -\frac{M_a P_{atm}}{R_g T_m^2} (T_{atm}(t) - T_m) \quad (3)$$

where $P_{atm}$ is the atmospheric pressure, $M_a$ the molar mass of air and $R_g$ the ideal gas constant. Equations 2 and 3 show that the air flowrate is proportional to the square root of the temperature drop ($T_{atm} - T_m$). Injecting Equations 1 and 3 in Equation 2 yields

$$\dot{m}(t) = S \Delta \dot{m} \sqrt{\left| \sin\left(\frac{2\pi}{\tau_y} t\right) \right|} \quad (4)$$

where $S = -1$ for $T_{atm} < T_m$ (wintertime, upward flow) and $S = 1$ for $T_{atm} > T_m$ (summertime, downward flow). The positive constant $\Delta \dot{m}$ is the amplitude of the annual fluctuation of the flow rate. It is an input of the model.

## 2.3. Governing Equations

The air temperature verifies the energy balance
$$ \dot{m} c_{p,a} \frac{\partial T_a}{\partial x} = \varphi_w P \qquad (5) $$

where $c_{p,a}$ is the air heat capacity at constant pressure, $T_a$ is the mixing temperature of the air inside the conduit, $x$ the distance from the external rock surface, $P$ the conduit perimeter ($P = 2\pi R$ for a circular cross-section) and $\varphi_w$ the conductive flux at the conduit wall, positive when going from the rock to the air. Equation 5 is a balance between the energy advected by the air flow (LHS) and the conduction flux at the conduit wall (RHS). Air thermal inertia has been neglected compared to advection (quasi-steady approximation). Equation 5 requires a single boundary condition at the conduit inlet located at $x = 0$ for downward flow ($S = 1$) or at $x = L_{dom}$ for upward flow ($S = -1$). We thus impose the boundary condition

$$ T_a(x = 0, t) = T_{atm}(t) \text{ for } S = 1 \text{ or } T_a(x = L_{dom}, t) = T_m \text{ for } S = -1 \qquad (6) $$

The conductive flux $\varphi_w$ is required to solve Equation 5. Therefore, the conduction equation must be solved in the impermeable rock matrix around the cave. The transient 2D axisymmetric conduction equation reads

$$ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial T_r}{\partial r} \right) + \frac{\partial^2 T_r}{\partial x^2} = \frac{1}{\alpha_r} \frac{\partial T_r}{\partial t} \qquad (7) $$

where $T_r(x, r, t)$ is the rock temperature, $\alpha_r$ the rock thermal diffusivity, and $r$ the distance from the cave axis. The boundary conditions are as follows. The atmosphere temperature is imposed on the external rock surface:

$$ T_r(x = 0, r, t) = T_{atm}(t) \text{ for } R < r < R_{dom} \qquad (8) $$

The boundary at $r = R_{dom}$ can be assumed adiabatic because the temperature field only depends on $x$ far from the cave (i.e., for a large value of $R_{dom}$), and the radial component of the temperature gradient is thus close to zero. The boundary condition at $x = L_{dom}$ is also adiabatic because this boundary is located in the homothermic zone. We get

$$ \frac{\partial T_r}{\partial r}(x, R_{dom}, t) = 0 \text{ for } 0 < x < L_{dom} \text{ and } \frac{\partial T_r}{\partial x}(L_{dom}, r, t) = 0 \text{ for } R < r < R_{dom} \qquad (9) $$

The last boundary condition is at the conduit wall. It is given by Newton's law of cooling and the heat flux continuity at the conduit wall:

$$ \varphi_w = \lambda_r \frac{\partial T_r}{\partial r}(x, R, t) = h_{th}(t) (T_r(x, R, t) - T_a(x, t)) \qquad (10) $$

where $\lambda_r$ is the rock thermal conductivity, $T_r(x, R, t) = T_w(x, t)$ is the temperature of the conduit wall and $h_{th}$ the heat transfer coefficient. The latter is time-dependent since it depends on the air flowrate $\dot{m}(t)$. All the conditions at the boundaries of the rock domain are displayed in Figure 2. No initial condition is required because we are looking for the periodic regime, that is, the solution asymptotically reached by the model at infinite time.

## 2.4. Determination of the Heat Transfer Coefficient

The heat transfer coefficient $h_{th}$ in Equation 10 is a key parameter. The estimation of $h_{th}$ must therefore be performed with great care to get a reliable model. Assuming forced convection in conduits, Wigley and Brown (1971) and Gabrovšek (2023) used a correlation close to the Colburn correlation (Bergman et al., 2017), valid for fully developed turbulent flow in smooth pipes. Lismonde (2002) pointed out that a cave cannot be considered as a smooth pipe, and multiplied by two the numerical prefactor of the Colburn correlation to account for the effect of wall roughness. Covington et al. (2011) used the Gnielinski correlation (Bergman et al., 2017), which explicitly considers the effect of wall roughness.

However, it is unlikely that a cave can be considered as a rough pipe. Indeed, caves can have complicated shapes, ranging from sub-circular conduits with several meters in diameter to complex cave passages associated with the
$$ \dot{m} c_{p,a} \frac{\partial T_a}{\partial x} = \varphi_w P \qquad (5) $$

where $c_{p,a}$ is the air heat capacity at constant pressure, $T_a$ is the mixing temperature of the air inside the conduit, $x$ the distance from the external rock surface, $P$ the conduit perimeter ($P = 2\pi R$ for a circular cross-section) and $\varphi_w$ the conductive flux at the conduit wall, positive when going from the rock to the air. Equation 5 is a balance between the energy advected by the air flow (LHS) and the conduction flux at the conduit wall (RHS). Air thermal inertia has been neglected compared to advection (quasi-steady approximation). Equation 5 requires a single boundary condition at the conduit inlet located at $x = 0$ for downward flow ($S = 1$) or at $x = L_{dom}$ for upward flow ($S = -1$). We thus impose the boundary condition

$$ T_a(x = 0, t) = T_{atm}(t) \text{ for } S = 1 \text{ or } T_a(x = L_{dom}, t) = T_m \text{ for } S = -1 \qquad (6) $$

The conductive flux $\varphi_w$ is required to solve Equation 5. Therefore, the conduction equation must be solved in the impermeable rock matrix around the cave. The transient 2D axisymmetric conduction equation reads

$$ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial T_r}{\partial r} \right) + \frac{\partial^2 T_r}{\partial x^2} = \frac{1}{\alpha_r} \frac{\partial T_r}{\partial t} \qquad (7) $$

where $T_r(x, r, t)$ is the rock temperature, $\alpha_r$ the rock thermal diffusivity, and $r$ the distance from the cave axis. The boundary conditions are as follows. The atmosphere temperature is imposed on the external rock surface:

$$ T_r(x = 0, r, t) = T_{atm}(t) \text{ for } R < r < R_{dom} \qquad (8) $$

The boundary at $r = R_{dom}$ can be assumed adiabatic because the temperature field only depends on $x$ far from the cave (i.e., for a large value of $R_{dom}$), and the radial component of the temperature gradient is thus close to zero. The boundary condition at $x = L_{dom}$ is also adiabatic because this boundary is located in the homothermic zone. We get

$$ \frac{\partial T_r}{\partial r}(x, R_{dom}, t) = 0 \text{ for } 0 < x < L_{dom} \text{ and } \frac{\partial T_r}{\partial x}(L_{dom}, r, t) = 0 \text{ for } R < r < R_{dom} \qquad (9) $$

The last boundary condition is at the conduit wall. It is given by Newton's law of cooling and the heat flux continuity at the conduit wall:

$$ \varphi_w = \lambda_r \frac{\partial T_r}{\partial r}(x, R, t) = h_{th}(t) (T_r(x, R, t) - T_a(x, t)) \qquad (10) $$

where $\lambda_r$ is the rock thermal conductivity, $T_r(x, R, t) = T_w(x, t)$ is the temperature of the conduit wall and $h_{th}$ the heat transfer coefficient. The latter is time-dependent since it depends on the air flowrate $\dot{m}(t)$. All the conditions at the boundaries of the rock domain are displayed in Figure 2. No initial condition is required because we are looking for the periodic regime, that is, the solution asymptotically reached by the model at infinite time.

## 2.4. Determination of the Heat Transfer Coefficient

The heat transfer coefficient $h_{th}$ in Equation 10 is a key parameter. The estimation of $h_{th}$ must therefore be performed with great care to get a reliable model. Assuming forced convection in conduits, Wigley and Brown (1971) and Gabrovšek (2023) used a correlation close to the Colburn correlation (Bergman et al., 2017), valid for fully developed turbulent flow in smooth pipes. Lismonde (2002) pointed out that a cave cannot be considered as a smooth pipe, and multiplied by two the numerical prefactor of the Colburn correlation to account for the effect of wall roughness. Covington et al. (2011) used the Gnielinski correlation (Bergman et al., 2017), which explicitly considers the effect of wall roughness.

However, it is unlikely that a cave can be considered as a rough pipe. Indeed, caves can have complicated shapes, ranging from sub-circular conduits with several meters in diameter to complex cave passages associated with the
<figure><img src="image_4.png" /><figcaption>Figure 2. Rock domain with the boundary conditions of the heat conduction Equation 7.</figcaption></figure>

collapse of the host-rock. A succession of bends, conduit contractions or enlargements, or obstacles of any kind are expected to increase the transfer of heat between the air and the wall. Indeed, all of these singularities not only increase the heat exchange surface between the wall and the fluid but also enhance the heat transfer coefficient by generating secondary flows and by increasing the turbulence level. It is well known that heat exchanger performance can be improved by the insertion of twisted tapes, longitudinal fins, or helical ribs inside the tubes (Bergman et al., 2017). For instance, spirally corrugated tubes can increase the heat transfer coefficient by a factor 3 (Pethkool et al., 2011; Promthaisong et al., 2016). However, although empirical correlations exist to predict the heat transfer coefficient for a wide range of well-defined geometries used in heat exchangers, no correlation is available for the irregular and tortuous geometries commonly encountered in caves.

To get around this difficulty, two distinct cases will be considered. In the first case (case A in the following), correlations for forced convection and fully developed flow in pipes will be used to estimate the heat transfer coefficient. This approach is expected to give a lower bound of the heat transfer coefficient, and thus a lower bound of the heat flux at the conduit wall. In the second case (case B in the following), a higher bound of the wall heat flux will be obtained assuming an infinite heat transfer coefficient. Reality must lie between these two limiting cases.

### Case A: Forced convection and fully developed flow

The Nusselt number $Nu = \frac{h_h D_h}{k_a}$ used for the estimation of $h_{th}$ is displayed in Figure 3 as a function of the Reynolds number $Re_t = \frac{\nu(t) D_h}{\nu_a}$, for three values of the relative roughness $\epsilon$ (ratio of the wall roughness over $D_h$). $k_a$ and $\nu_a$ are the conductivity and the kinematic viscosity of the air, $V$ is the mean air velocity at time $t$, $D_h = 4 A/P$ is the hydraulic diameter with $A$ and $P$ the cross-sectional area and the perimeter of the conduit. $Nu$ is computed as follows:

$$Nu = \max(Nu_L, Nu_T) \qquad (11)$$

where $Nu_L = 3.66$ is the Nusselt number in the laminar regime (independent of the roughness) and $Nu_T$ is given by the Gnielinski correlation (Bergman et al., 2017), which takes into account the effect of the wall roughness expected in the turbulent regime:

$$Nu_T = \frac{\left(\frac{f_d}{8}\right) (Re_t - 1000) Pr}{1 + 12.7 \left(\frac{f_d}{8}\right)^{0.5} (Pr^{2/3} - 1)} \qquad (12)$$

where $Pr$ is the air Prandtl number ($Pr = 0.71$). $f_d$ is the Darcy friction factor which depends on the Reynolds number and the wall relative roughness $\epsilon$. $f_d$ was estimated using the Haaland correlation (Haaland, 1983):

<figure><img src="image_5.png" /><figcaption>Figure 3. Nusselt number as a function of the Reynolds number (Prandtl number: $Pr = 0.71$). The red curve was used for the simulations of case A.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2. Rock domain with the boundary conditions of the heat conduction Equation 7.</figcaption></figure>

collapse of the host-rock. A succession of bends, conduit contractions or enlargements, or obstacles of any kind are expected to increase the transfer of heat between the air and the wall. Indeed, all of these singularities not only increase the heat exchange surface between the wall and the fluid but also enhance the heat transfer coefficient by generating secondary flows and by increasing the turbulence level. It is well known that heat exchanger performance can be improved by the insertion of twisted tapes, longitudinal fins, or helical ribs inside the tubes (Bergman et al., 2017). For instance, spirally corrugated tubes can increase the heat transfer coefficient by a factor 3 (Pethkool et al., 2011; Promthaisong et al., 2016). However, although empirical correlations exist to predict the heat transfer coefficient for a wide range of well-defined geometries used in heat exchangers, no correlation is available for the irregular and tortuous geometries commonly encountered in caves.

To get around this difficulty, two distinct cases will be considered. In the first case (case A in the following), correlations for forced convection and fully developed flow in pipes will be used to estimate the heat transfer coefficient. This approach is expected to give a lower bound of the heat transfer coefficient, and thus a lower bound of the heat flux at the conduit wall. In the second case (case B in the following), a higher bound of the wall heat flux will be obtained assuming an infinite heat transfer coefficient. Reality must lie between these two limiting cases.

### Case A: Forced convection and fully developed flow

The Nusselt number $Nu = \frac{h_h D_h}{k_a}$ used for the estimation of $h_{th}$ is displayed in Figure 3 as a function of the Reynolds number $Re_t = \frac{\nu(t) D_h}{\nu_a}$, for three values of the relative roughness $\epsilon$ (ratio of the wall roughness over $D_h$). $k_a$ and $\nu_a$ are the conductivity and the kinematic viscosity of the air, $V$ is the mean air velocity at time $t$, $D_h = 4 A/P$ is the hydraulic diameter with $A$ and $P$ the cross-sectional area and the perimeter of the conduit. $Nu$ is computed as follows:

$$Nu = \max(Nu_L, Nu_T) \quad (11)$$

where $Nu_L = 3.66$ is the Nusselt number in the laminar regime (independent of the roughness) and $Nu_T$ is given by the Gnielinski correlation (Bergman et al., 2017), which takes into account the effect of the wall roughness expected in the turbulent regime:

$$Nu_T = \frac{\left(\frac{f_d}{8}\right) (Re_t - 1000) Pr}{1 + 12.7 \left(\frac{f_d}{8}\right)^{0.5} (Pr^{2/3} - 1)} \quad (12)$$

where $Pr$ is the air Prandtl number ($Pr = 0.71$). $f_d$ is the Darcy friction factor which depends on the Reynolds number and the wall relative roughness $\epsilon$. $f_d$ was estimated using the Haaland correlation (Haaland, 1983):

<figure><img src="image_5.png" /><figcaption>Figure 3. Nusselt number as a function of the Reynolds number (Prandtl number: $Pr = 0.71$). The red curve was used for the simulations of case A.</figcaption></figure>
<table><caption>**Table 1**

Thermophysical Properties (Estimated at Temperature $T_m$ = 285 K = 12°C and Atmospheric Pressure $p_{ref}$ = 101,325 Pa)</caption><thead><tr><th>Properties</th><th>Value and unit</th><th>Reference</th></tr></thead><tbody><tr><td>Conduit relative roughness</td><td>ε = 0.01</td><td>(see Section 2.4)</td></tr><tr><td>Rock density</td><td>ρ<sub>r</sub> = 2320 kg/m<sup>3</sup></td><td>Covington et al. (2011)</td></tr><tr><td>Rock heat capacity</td><td>c<sub>p,r</sub> = 810 J/(kg K)</td><td>Covington et al. (2011)</td></tr><tr><td>Rock thermal conductivity</td><td>k<sub>r</sub> = 1.656 W/(m·K)</td><td>Guerrier et al. (2019)</td></tr><tr><td>Rock thermal diffusivity</td><td>α<sub>r</sub> = <sup>k<sub>r</sub></sup>/<sub>ρ<sub>r</sub>·c<sub>p,r</sub></sub> = 8.81 × 10<sup>-7</sup> m<sup>2</sup>/s</td><td>–</td></tr><tr><td>Air dynamic viscosity</td><td>μ<sub>a</sub> = 1.77 × 10<sup>-5</sup> Pa·s</td><td>Bergman et al. (2017)</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub> = 1.23 kg/m<sup>3</sup></td><td>Bergman et al. (2017)</td></tr><tr><td>Air kinematic viscosity</td><td>ν<sub>a</sub> = <sup>μ<sub>a</sub></sup>/<sub>ρ<sub>a</sub></sub> = 1.43 × 10<sup>-5</sup> m<sup>2</sup>/s</td><td>–</td></tr><tr><td>Air thermal conductivity</td><td>k<sub>a</sub> = 0.0251 W/(m·K)</td><td>Bergman et al. (2017)</td></tr><tr><td>Air Prandtl number</td><td>Pr = <sup>ν<sub>a</sub></sup>/<sub>α<sub>a</sub></sub> = 0.71</td><td>–</td></tr><tr><td>Molar mass of air</td><td>M<sub>a</sub> = 28.97 g/mol</td><td>–</td></tr><tr><td>Ideal gas constant</td><td>R<sub>g</sub> = 8.314 J/mol·K</td><td>–</td></tr></tbody></table>

$$ \frac{1}{\sqrt{f_d}} = -1.8 \log \left[ \left( \frac{\varepsilon}{3.7} \right)^{1.11} + \frac{6.9}{R_{e_i}} \right] \quad (13) $$

All the simulations of case A were done using $\epsilon = 0.01$ as a lower bound of the relative roughness in a cave (red curve in Figure 3).

**Case B:** Infinite heat transfer coefficient

In this case, the conduit wall temperature is equal to the bulk air temperature and Equation 10 thus reduces to

$$ \varphi_w = \lambda_r \frac{\partial T_r}{\partial r}(x, R, t) \text{ and } T_r(x, R, t) = T_a(x, t) \quad (14) $$

## 2.5. Physical Properties

All the physical properties are gathered in Table 1. They are assumed constant and estimated at the temperature $T_m = 12°C$.

## 2.6. Dimensionless Equations

We define the dimensionless temperatures $\theta_a = (T_a - T_m)/\Delta T_y$ and $\theta_r = (T_r - T_m)/\Delta T_y$ in the air and the rock, respectively. The dimensionless temperature 0 thus corresponds to the external AMT and 1 is the amplitude of the external annual fluctuations. We take the period $\tau_y = 1$ year as the unit of time and the diffusion length $L_{dif} = \sqrt{\alpha_r \tau_y} \approx 5.3$ m as the unit of length. The dimensionless atmospheric temperature $\theta_{atm}$ and air flowrate $\mu$ are deduced from Equations 1 and 4:

$$ \theta_{atm}(\tilde{t}) = \sin(2\pi\tilde{t}) \quad (15) $$

$$ \mu(\tilde{t}) = \frac{\dot{m}(t)}{\Delta\dot{m}} = S\sqrt{|\sin(2\pi\tilde{t})|} \quad (16) $$

The dimensionless counterparts of Equations 5–10 read

$$ \mu(\tilde{t}) \text{Re} \frac{\partial \theta_a}{\partial \tilde{x}} = \left( 4 \frac{k_r}{k_a} Pr^{-1} \right) \tilde{\phi}_w, \quad (17) $$
<table><caption>**Table 1**

Thermophysical Properties (Estimated at Temperature $T_m$ = 285 K = 12°C and Atmospheric Pressure $p_{ref}$ = 101,325 Pa)</caption><thead><tr><th>Properties</th><th>Value and unit</th><th>Reference</th></tr></thead><tbody><tr><td>Conduit relative roughness</td><td>ε = 0.01</td><td>(see Section 2.4)</td></tr><tr><td>Rock density</td><td>ρ<sub>r</sub> = 2320 kg/m<sup>3</sup></td><td>Covington et al. (2011)</td></tr><tr><td>Rock heat capacity</td><td>c<sub>p,r</sub> = 810 J/(kg K)</td><td>Covington et al. (2011)</td></tr><tr><td>Rock thermal conductivity</td><td>k<sub>r</sub> = 1.656 W/(m K)</td><td>Guerrier et al. (2019)</td></tr><tr><td>Rock thermal diffusivity</td><td>α<sub>r</sub> = <sup>k<sub>r</sub></sup>/<sub>ρ<sub>r</sub>ε<sub>p,r</sub></sub> = 8.81 × 10<sup>-7</sup> m<sup>2</sup>/s</td><td>–</td></tr><tr><td>Air dynamic viscosity</td><td>μ<sub>a</sub> = 1.77 × 10<sup>-5</sup> Pa.s</td><td>Bergman et al. (2017)</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub> = 1.23 kg/m<sup>3</sup></td><td>Bergman et al. (2017)</td></tr><tr><td>Air kinematic viscosity</td><td>ν<sub>a</sub> = <sup>μ<sub>a</sub></sup>/<sub>ρ<sub>a</sub></sub> = 1.43 × 10<sup>-5</sup> m<sup>2</sup>/s</td><td>–</td></tr><tr><td>Air thermal conductivity</td><td>k<sub>a</sub> = 0.0251 W/(m K)</td><td>Bergman et al. (2017)</td></tr><tr><td>Air Prandtl number</td><td>Pr = <sup>ν<sub>a</sub></sup>/<sub>α<sub>a</sub></sub> = 0.71</td><td>–</td></tr><tr><td>Molar mass of air</td><td>M<sub>a</sub> = 28.97 g/mol</td><td>–</td></tr><tr><td>Ideal gas constant</td><td>R<sub>g</sub> = 8.314 J/mol K</td><td>–</td></tr></tbody></table>

$$ \frac{1}{\sqrt{f_d}} = -1.8 \log \left[ \left( \frac{\epsilon}{3.7} \right)^{1.11} + \frac{6.9}{Re_i} \right] \quad (13) $$

All the simulations of case A were done using $\epsilon = 0.01$ as a lower bound of the relative roughness in a cave (red curve in Figure 3).

### Case B: Infinite heat transfer coefficient

In this case, the conduit wall temperature is equal to the bulk air temperature and Equation 10 thus reduces to

$$ \varphi_w = \lambda_r \frac{\partial T_r}{\partial r}(x, R, t) \text{ and } T_r(x, R, t) = T_a(x, t) \quad (14) $$

## 2.5. Physical Properties

All the physical properties are gathered in Table 1. They are assumed constant and estimated at the temperature $T_m = 12°C$.

## 2.6. Dimensionless Equations

We define the dimensionless temperatures $\theta_a = (T_a - T_m)/\Delta T_y$ and $\theta_r = (T_r - T_m)/\Delta T_y$ in the air and the rock, respectively. The dimensionless temperature 0 thus corresponds to the external AMT and 1 is the amplitude of the external annual fluctuations. We take the period $\tau_y = 1$ year as the unit of time and the diffusion length $L_{dif} = \sqrt{\alpha_r \tau_y} \approx 5.3$ m as the unit of length. The dimensionless atmospheric temperature $\theta_{atm}$ and air flowrate $\mu$ are deduced from Equations 1 and 4:

$$ \theta_{atm}(\tilde{t}) = \sin(2\pi\tilde{t}) \quad (15) $$

$$ \mu(\tilde{t}) = \frac{\dot{m}(t)}{\Delta\dot{m}} = S\sqrt{|\sin(2\pi\tilde{t})|} \quad (16) $$

The dimensionless counterparts of Equations 5–10 read

$$ \mu(\tilde{t}) \text{Re} \frac{\partial \theta_a}{\partial \tilde{x}} = \left( 4 \frac{k_r}{k_a} Pr^{-1} \right) \tilde{\phi}_w, \quad (17) $$
$$ \theta_a(\tilde{x} = 0, t) = \sin(2\pi \tilde{t}) \text{ for } S = 1 \text{ or } \theta_a(\tilde{x} = \tilde{L}_{dom}, \tilde{t}) = 0 \text{ for } S = -1, \quad (18) $$

$$ \frac{1}{\tilde{r}} \frac{\partial}{\partial \tilde{r}} \left( \tilde{r} \frac{\partial \theta_r}{\partial \tilde{r}} \right) + \frac{\partial^2 \theta_r}{\partial \tilde{x}^2} = \frac{\partial \theta_r}{\partial \tilde{t}}, \quad (19) $$

$$ \theta_r(\tilde{x} = 0, \tilde{r}, \tilde{t}) = \sin(2\pi \tilde{t}) \text{ for } \tilde{R} < \tilde{r} < \tilde{R}_{dom}, \quad (20) $$

$$ \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}_{dom}, \tilde{t}) = 0 \text{ for } 0 < \tilde{x} < \tilde{L}_{dom} \text{ and } \frac{\partial \theta_r}{\partial \tilde{x}}(\tilde{L}_{dom}, \tilde{r}, \tilde{t}) = 0 \text{ for } \tilde{R} < \tilde{r} < \tilde{R}_{dom}, \quad (21) $$

$$ \tilde{\phi}_w = \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}, \tilde{t}) = \eta(\tilde{t}) Bi (\theta_r(\tilde{x}, \tilde{R}, \tilde{t}) - \theta_a(\tilde{x}, \tilde{t})), \quad (22) $$

where the tilde (~) indicates dimensionless time or length. $\eta(\tilde{t}) = \frac{h_{th}(t)}{H_{th}}$ is the heat transfer coefficient scaled by $H_{th}$, its maximum value reached for $\dot{m} = \Delta\dot{m}$. The Reynolds number $Re = \frac{4 \Delta \dot{m}}{\rho \mu_a}$ (with $\mu_a$ the air dynamic viscosity) and the Biot number $Bi = \frac{H_{th}}{\frac{\sqrt{a_r + \tilde{r}}}{k_r}}$ are based on the maximum values $\Delta\dot{m}$ and $H_{th}$, respectively. $\mu(\tilde{t}) Re$ and $\eta(\tilde{t}) Bi$ are the Reynolds and Biot numbers at a given time $\tilde{t}$. In a given configuration, $\tilde{R}_{dom}$ and $\tilde{L}_{dom}$ that define the size of the computational domain are increased until they do not modify the results, as stated in Section 2.1.

With the constant physical properties defined in Table 1, a configuration is well-defined if only two dimensionless numbers are known: the dimensionless conduit radius $\tilde{R}$ and the Reynolds number $Re$. In case A, that is, when fully developed forced convection is assumed for the estimation of the heat transfer coefficient, the Biot number in Equation 22 is related to the Nusselt number by the simple relation $\eta(\tilde{t}) Bi = \frac{Nu}{\tilde{R}} \left(\frac{k_a}{2 k_r}\right)$, where $Nu$ depends on $Re_{\tilde{t}} = \mu(\tilde{t}) Re$ through Equations 11 and 12. In case B, the Biot number is infinite, that is, $\eta(\tilde{t}) Bi \to \infty$. Equation 22 thus reduces to

$$ \tilde{\phi}_w = \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}, \tilde{t}) \quad \text{and} \quad \theta_r(\tilde{x}, \tilde{R}, \tilde{t}) = \theta_a(\tilde{x}, \tilde{t}) \quad (23) $$

### 3. Method of Resolution

We are looking for the periodic solution of the mathematical model defined in Section 2.6. The periodic regime could be obtained by time integration starting from an arbitrary initial condition, after simulating a number of cycles large enough to get a good approximation of the solution at infinite time. However, the slow convergence of this method can generate large computational times. We thus prefer to use another method based on the Fourier series. All time-varying variables, including the temperatures in the rock and in the air, are approximated by truncated Fourier series and inserted in the mathematical model of Section 2.6. This yields a set of coupled partial differential equations (PDE), which does not include the time. The numerical resolution by the Galerkin method of this set of PDE gives the spatial distribution of the amplitude and phase of the Fourier modes. The temporal evolution of the dependent variables is then recovered from Fourier series. This method is detailed in Supporting Information S1. The validation by comparison with a standard time integration method is presented in the same appendix. All the numerical simulations were performed with the commercial software Comsol Multiphysics, version 6.1.

The method of resolution based on the Fourier series used in this work is particularly relevant for the simplified time variation of the external temperature considered in Equation 1. Considering all the Fourier modes included in the real time series of the external temperature would clearly increase both the complexity of the method implementation and the computational resources required for the simulations. In that case, only the most significant modes of the input temperature should be considered. In contrast, using Fourier series to simulate the complex geometry of real caves would not cause any particular problem.
$$ \theta_a(\tilde{x} = 0, t) = \sin(2\pi \tilde{t}) \text{ for } S = 1 \text{ or } \theta_a(\tilde{x} = \tilde{L}_{dom}, \tilde{t}) = 0 \text{ for } S = -1, \quad (18) $$

$$ \frac{1}{\tilde{r}} \frac{\partial}{\partial \tilde{r}} \left( \tilde{r} \frac{\partial \theta_r}{\partial \tilde{r}} \right) + \frac{\partial^2 \theta_r}{\partial \tilde{x}^2} = \frac{\partial \theta_r}{\partial \tilde{t}}, \quad (19) $$

$$ \theta_r(\tilde{x} = 0, \tilde{r}, \tilde{t}) = \sin(2\pi \tilde{t}) \text{ for } \tilde{R} < \tilde{r} < \tilde{R}_{dom}, \quad (20) $$

$$ \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}_{dom}, \tilde{t}) = 0 \text{ for } 0 < \tilde{x} < \tilde{L}_{dom} \text{ and } \frac{\partial \theta_r}{\partial \tilde{x}}(\tilde{L}_{dom}, \tilde{r}, \tilde{t}) = 0 \text{ for } \tilde{R} < \tilde{r} < \tilde{R}_{dom}, \quad (21) $$

$$ \tilde{\phi}_w = \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}, \tilde{t}) = \eta(\tilde{t}) Bi (\theta_r(\tilde{x}, \tilde{R}, \tilde{t}) - \theta_a(\tilde{x}, \tilde{t})), \quad (22) $$

where the tilde (~) indicates dimensionless time or length. $\eta(\tilde{t}) = \frac{h_{th}(t)}{H_{th}}$ is the heat transfer coefficient scaled by $H_{th}$, its maximum value reached for $\dot{m} = \Delta\dot{m}$. The Reynolds number $Re = \frac{4\Delta\dot{m}}{\rho\mu_a}$ (with $\mu_a$ the air dynamic viscosity) and the Biot number $Bi = \frac{H_{th}}{\frac{\sqrt{a_r+\tilde{r}}}{k_r}}$ are based on the maximum values $\Delta\dot{m}$ and $H_{th}$, respectively. $\mu(\tilde{t})Re$ and $\eta(\tilde{t})Bi$ are the Reynolds and Biot numbers at a given time $\tilde{t}$. In a given configuration, $\tilde{R}_{dom}$ and $\tilde{L}_{dom}$ that define the size of the computational domain are increased until they do not modify the results, as stated in Section 2.1.

With the constant physical properties defined in Table 1, a configuration is well-defined if only two dimensionless numbers are known: the dimensionless conduit radius $\tilde{R}$ and the Reynolds number $Re$. In case A, that is, when fully developed forced convection is assumed for the estimation of the heat transfer coefficient, the Biot number in Equation 22 is related to the Nusselt number by the simple relation $\eta(\tilde{t})Bi = \frac{Nu}{\tilde{R}}(\frac{k_a}{2k_r})$, where $Nu$ depends on $Re_{\tilde{t}} = \mu(\tilde{t})Re$ through Equations 11 and 12. In case B, the Biot number is infinite, that is, $\eta(\tilde{t})Bi \to \infty$. Equation 22 thus reduces to

$$ \tilde{\phi}_w = \frac{\partial \theta_r}{\partial \tilde{r}}(\tilde{x}, \tilde{R}, \tilde{t}) \quad \text{and} \quad \theta_r(\tilde{x}, \tilde{R}, \tilde{t}) = \theta_a(\tilde{x}, \tilde{t}) \quad (23) $$

### 3. Method of Resolution

We are looking for the periodic solution of the mathematical model defined in Section 2.6. The periodic regime could be obtained by time integration starting from an arbitrary initial condition, after simulating a number of cycles large enough to get a good approximation of the solution at infinite time. However, the slow convergence of this method can generate large computational times. We thus prefer to use another method based on the Fourier series. All time-varying variables, including the temperatures in the rock and in the air, are approximated by truncated Fourier series and inserted in the mathematical model of Section 2.6. This yields a set of coupled partial differential equations (PDE), which does not include the time. The numerical resolution by the Galerkin method of this set of PDE gives the spatial distribution of the amplitude and phase of the Fourier modes. The temporal evolution of the dependent variables is then recovered from Fourier series. This method is detailed in Supporting Information S1. The validation by comparison with a standard time integration method is presented in the same appendix. All the numerical simulations were performed with the commercial software Comsol Multiphysics, version 6.1.

The method of resolution based on the Fourier series used in this work is particularly relevant for the simplified time variation of the external temperature considered in Equation 1. Considering all the Fourier modes included in the real time series of the external temperature would clearly increase both the complexity of the method implementation and the computational resources required for the simulations. In that case, only the most significant modes of the input temperature should be considered. In contrast, using Fourier series to simulate the complex geometry of real caves would not cause any particular problem.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wall and air temperatures at different distances $\tilde{x}$ from the entrance (case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Amplitude of the Fourier modes. (b) Time series.</figcaption></figure>

**4. Characterization of the Thermal Perturbation**

Figure 4 displays the wall and air temperatures inside the cave at different distances from the entrance, for case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$. Figure 4a shows the amplitude of the Fourier modes that stems directly from the numerical simulation. The mode $k$ refers to the dimensionless period $1/k$. The amplitudes of the modes $k = 0$ and $k = 1$ are thus the AMT and the amplitude of the ATF, respectively. The modes $k = 2, k = 3$, and so on refer to the periods 6 months, 4 months, etc. Although the external temperature imposed at $\tilde{x} = 0$ contains the single mode $k = 1$ (since it reduces to Equation 15), the coupling between modes due to convection results in the existence of other modes than $k = 1$ in cave temperatures (modes $k = 0$ to 18 were considered in the simulations, see Supporting Information S1). Figure 4b displays the periodic time series built from the Fourier series. The amplitude of the temperature fluctuations decreases when the distance from the entrance increases, as expected. A hot thermal anomaly expected at the vicinity of an upper entrance is also observed. Indeed, for $\tilde{x} = 10$ to $10^3$, the dimensionless temperature is positive throughout the year, which means that it is always higher than the external AMT, even in winter.

Figure 5a displays the temperature profiles in the air and at the cave wall at times $t = 0, 0.25, 0.5$, and $0.75$, corresponding to the external temperatures 0, 1, 0 and $-1$, respectively. Following an approach similar to that of Cropley (1965), we split the cave into three zones.

1. A first zone, from $\tilde{x} = 0$ to $\tilde{x} \approx 1$ (the diffusion length), where the cave temperature follows the fluctuations of the external temperature. More precisely, the sign of the temperature gradient changes in summer and winter (see the inset in Figure 5a). This first zone is dominated by the temperature fluctuations of the external wall, propagating inside the rock mass by heat conduction. It will be called the diffusive region in the following.

2. A second zone characterized by a positive dimensionless temperature and reduced but significant fluctuations. In the specific case of Figure 5a, it approximately extends from $\tilde{x} \approx 1$ to $10^3$. In this zone, heat diffusion from the external wall is negligible compared to convection. It will be called convective region in the following. Because of the seasonal flow reversals, the inlet temperature of the air flow is 1 in summer (inward flow) and 0 in winter (outward flow, coming from the deep karst). As a consequence, this zone does not “feel” the negative winter temperatures. This results in: (a) a shift of the AMT to positive values, the so-called thermal anomaly, (b) a reduction of the amplitude of the fluctuations, because the temperature oscillates between 0 and 1 whereas the external temperature fluctuates between $-1$ and $1$.

3. A third zone, called deep karst in the following, where the dimensionless temperature is approximately constant and close to zero (the external AMT).

These three regions appear clearly on the amplitude of the Fourier modes plotted in Figure 5b as a function of the distance from the entrance $\tilde{x}$. The amplitude of the mode $k = 1$ (i.e., the ATF) decreases in two steps. It first drops
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wall and air temperatures at different distances $\tilde{x}$ from the entrance (case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Amplitude of the Fourier modes. (b) Time series.</figcaption></figure>

**4. Characterization of the Thermal Perturbation**

Figure 4 displays the wall and air temperatures inside the cave at different distances from the entrance, for case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$. Figure 4a shows the amplitude of the Fourier modes that stems directly from the numerical simulation. The mode $k$ refers to the dimensionless period $1/k$. The amplitudes of the modes $k = 0$ and $k = 1$ are thus the AMT and the amplitude of the ATF, respectively. The modes $k = 2, k = 3$, and so on refer to the periods 6 months, 4 months, etc. Although the external temperature imposed at $\tilde{x} = 0$ contains the single mode $k = 1$ (since it reduces to Equation 15), the coupling between modes due to convection results in the existence of other modes than $k = 1$ in cave temperatures (modes $k = 0$ to 18 were considered in the simulations, see Supporting Information S1). Figure 4b displays the periodic time series built from the Fourier series. The amplitude of the temperature fluctuations decreases when the distance from the entrance increases, as expected. A hot thermal anomaly expected at the vicinity of an upper entrance is also observed. Indeed, for $\tilde{x} = 10$ to $10^3$, the dimensionless temperature is positive throughout the year, which means that it is always higher than the external AMT, even in winter.

Figure 5a displays the temperature profiles in the air and at the cave wall at times $t = 0, 0.25, 0.5$, and $0.75$, corresponding to the external temperatures 0, 1, 0 and $-1$, respectively. Following an approach similar to that of Cropley (1965), we split the cave into three zones.

1. A first zone, from $\tilde{x} = 0$ to $\tilde{x} \approx 1$ (the diffusion length), where the cave temperature follows the fluctuations of the external temperature. More precisely, the sign of the temperature gradient changes in summer and winter (see the inset in Figure 5a). This first zone is dominated by the temperature fluctuations of the external wall, propagating inside the rock mass by heat conduction. It will be called the diffusive region in the following.

2. A second zone characterized by a positive dimensionless temperature and reduced but significant fluctuations. In the specific case of Figure 5a, it approximately extends from $\tilde{x} \approx 1$ to $10^3$. In this zone, heat diffusion from the external wall is negligible compared to convection. It will be called convective region in the following. Because of the seasonal flow reversals, the inlet temperature of the air flow is 1 in summer (inward flow) and 0 in winter (outward flow, coming from the deep karst). As a consequence, this zone does not “feel” the negative winter temperatures. This results in: (a) a shift of the AMT to positive values, the so-called thermal anomaly, (b) a reduction of the amplitude of the fluctuations, because the temperature oscillates between 0 and 1 whereas the external temperature fluctuates between $-1$ and $1$.

3. A third zone, called deep karst in the following, where the dimensionless temperature is approximately constant and close to zero (the external AMT).

These three regions appear clearly on the amplitude of the Fourier modes plotted in Figure 5b as a function of the distance from the entrance $\tilde{x}$. The amplitude of the mode $k = 1$ (i.e., the ATF) decreases in two steps. It first drops
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Wall (solid lines) and air (dashed lines) temperatures as a function of the distance $\tilde{x}$ from the entrance (case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Temperature profiles at various times (Inset: zoom in the diffusive region). (b) Amplitude of the Fourier modes $k = 0$ to 3 (Inset: log scale on the $\tilde{x}$-axis).</figcaption></figure>

down from 1 to 0.3 in the diffusive region (from $\tilde{x} = 0$ to $\tilde{x} \approx 1$), then slowly decreases in the convective region, over a distance much longer than 1 (of the order of $10^3$ in the specific case of Figure 5). The amplitude of the mode $k = 0$ (i.e., the AMT) increases in the diffusive region until it reaches a quasi-plateau. Then it slowly decreases in the convective region, over a distance of the same order as for the mode $k = 1$. The modes such that $k \gtrsim 2$ follow qualitatively the behavior of mode $k = 0$, with a lower amplitude and a shorter attenuation length. In the following, we will focus on both modes with the highest amplitude, $k = 0$ and $k = 1$, as representative of the AMT (the thermal anomaly) and the ATF, respectively.

In all the simulations we did, the transition between the diffusive and convective regions always take place at a distance from the entrance of the order of 1, as expected. Conversely, the characteristics of the convective domain (amplitude of the thermal perturbation and its spatial extent) depend on the model parameters: the Reynolds number $Re$ (containing the air flow rate) and the cave radius $\tilde{R}$. This dependence will be systematically investigated in the parametric study of the next section. In order to facilitate the interpretation of the results, we define a small number of relevant parameters characterizing the perturbation of the cave temperature field by the air flow. The magnitude of the perturbation is characterized by the AMT and the ATF at the beginning of the convective region, where the impact of convection is maximum, and the effect of thermal conduction from the external wall is negligible. We define $\Delta\theta_w$ as the maximum amplitude of the mode $k = 1$ in the convective region, for the wall temperature. Practically, it is evaluated at the change of slope noticed on the solid blue line in the inset of Figure 5b. Similarly, $\bar{\theta}_w$ is the maximum amplitude reached by the mode $k = 0$ for the wall temperature (maximum of the red solid line in the inset of Figure 5b). The corresponding amplitudes for the air temperature, $\Delta\theta_a$ and $\bar{\theta}_a$, are estimated at the same location as $\Delta\theta_w$ and $\bar{\theta}_w$, respectively. For example, in the specific case of Figure 5, we get $\Delta\theta_w \approx 0.30$ and $\Delta\theta_a \approx 0.34$ (maximum ATF), $\bar{\theta}_w \approx 0.425$ and $\bar{\theta}_a \approx 0.44$ (maximum AMT).

The four parameters $\bar{\theta}_w$, $\bar{\theta}_a$, $\Delta\theta_w$, and $\Delta\theta_a$ represent the maximum amplitudes of the perturbation in the convective region. These amplitudes decrease with the distance from the entrance, as the heat transported by the air flow is gradually dissipated by conduction in the rock mass. We arbitrarily assume that a Fourier mode has been damped when its amplitude has been divided by 10. Figure 6a shows the thermal perturbation in a rock massif with the black solid lines corresponding to an AMT equal to $\bar{\theta}_w/10 = 0.0425$. The rock domain thus defined can be characterized by the length $\tilde{L}_0$ along the cave axis and the maximum width $\tilde{W}_0$ along the radial direction. Note that this maximum width is observed at a certain distance from the external wall because of the competition between radial diffusion from the cave which transports the thermal anomaly through the massif, and axial diffusion from the external wall which reduces it. Similarly, the black solid lines in Figure 6b correspond to an ATF equal to $\Delta\theta_w/10 = 0.030$. The size of the corresponding domain is characterized by its length $\tilde{L}_1$ and its width $\tilde{W}_1$. The latter is taken close to the entrance, but out of the diffusive region. The lengths $\tilde{L}_0$ and $\tilde{L}_1$ can also
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Wall (solid lines) and air (dashed lines) temperatures as a function of the distance $\tilde{x}$ from the entrance (case A with $\tilde{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Temperature profiles at various times (Inset: zoom in the diffusive region). (b) Amplitude of the Fourier modes $k = 0$ to 3 (Inset: log scale on the $\tilde{x}$-axis).</figcaption></figure>

down from 1 to 0.3 in the diffusive region (from $\tilde{x} = 0$ to $\tilde{x} \approx 1$), then slowly decreases in the convective region, over a distance much longer than 1 (of the order of $10^3$ in the specific case of Figure 5). The amplitude of the mode $k = 0$ (i.e., the AMT) increases in the diffusive region until it reaches a quasi-plateau. Then it slowly decreases in the convective region, over a distance of the same order as for the mode $k = 1$. The modes such that $k \gtrsim 2$ follow qualitatively the behavior of mode $k = 0$, with a lower amplitude and a shorter attenuation length. In the following, we will focus on both modes with the highest amplitude, $k = 0$ and $k = 1$, as representative of the AMT (the thermal anomaly) and the ATF, respectively.

In all the simulations we did, the transition between the diffusive and convective regions always take place at a distance from the entrance of the order of 1, as expected. Conversely, the characteristics of the convective domain (amplitude of the thermal perturbation and its spatial extent) depend on the model parameters: the Reynolds number $Re$ (containing the air flow rate) and the cave radius $\tilde{R}$. This dependence will be systematically investigated in the parametric study of the next section. In order to facilitate the interpretation of the results, we define a small number of relevant parameters characterizing the perturbation of the cave temperature field by the air flow. The magnitude of the perturbation is characterized by the AMT and the ATF at the beginning of the convective region, where the impact of convection is maximum, and the effect of thermal conduction from the external wall is negligible. We define $\Delta\theta_w$ as the maximum amplitude of the mode $k = 1$ in the convective region, for the wall temperature. Practically, it is evaluated at the change of slope noticed on the solid blue line in the inset of Figure 5b. Similarly, $\bar{\theta}_w$ is the maximum amplitude reached by the mode $k = 0$ for the wall temperature (maximum of the red solid line in the inset of Figure 5b). The corresponding amplitudes for the air temperature, $\Delta\theta_a$ and $\bar{\theta}_a$, are estimated at the same location as $\Delta\theta_w$ and $\bar{\theta}_w$, respectively. For example, in the specific case of Figure 5, we get $\Delta\theta_w \approx 0.30$ and $\Delta\theta_a \approx 0.34$ (maximum ATF), $\bar{\theta}_w \approx 0.425$ and $\bar{\theta}_a \approx 0.44$ (maximum AMT).

The four parameters $\bar{\theta}_w$, $\bar{\theta}_a$, $\Delta\theta_w$, and $\Delta\theta_a$ represent the maximum amplitudes of the perturbation in the convective region. These amplitudes decrease with the distance from the entrance, as the heat transported by the air flow is gradually dissipated by conduction in the rock mass. We arbitrarily assume that a Fourier mode has been damped when its amplitude has been divided by 10. Figure 6a shows the thermal perturbation in a rock massif with the black solid lines corresponding to an AMT equal to $\bar{\theta}_w/10 = 0.0425$. The rock domain thus defined can be characterized by the length $\tilde{L}_0$ along the cave axis and the maximum width $\tilde{W}_0$ along the radial direction. Note that this maximum width is observed at a certain distance from the external wall because of the competition between radial diffusion from the cave which transports the thermal anomaly through the massif, and axial diffusion from the external wall which reduces it. Similarly, the black solid lines in Figure 6b correspond to an ATF equal to $\Delta\theta_w/10 = 0.030$. The size of the corresponding domain is characterized by its length $\tilde{L}_1$ and its width $\tilde{W}_1$. The latter is taken close to the entrance, but out of the diffusive region. The lengths $\tilde{L}_0$ and $\tilde{L}_1$ can also
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Amplitude contour of the thermal perturbation in a symmetry plane of the rock massif (case A with $\bar{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Mode $k = 0$; the black solid lines correspond to the amplitude $\bar{\theta}_w/10 = 0.0425$ (n.b.: at the scale of the figure, the cave reduces to a line). (b) Mode $k = 1$; the black solid lines correspond to the amplitude $\Delta\bar{\theta}_w/10 = 0.030$ (n.b.: the top view focuses on the entrance region; the bottom view has been plotted with different scales on $\bar{x}$ and $\bar{r}$ to display the whole perturbed region).</figcaption></figure>

be estimated from the air temperature (length such that the amplitude of AMT and ATF in the air are equal to $\bar{\theta}_a/10$ and $\Delta\theta_a/10$, respectively).

# 5. Parametric Study

## 5.1. Amplitude of the Thermal Perturbation Induced by the Air Flow

Figure 7 displays the AMT ($\bar{\theta}_w$ and $\bar{\theta}_a$) and the ATF ($\Delta\theta_w$ and $\Delta\theta_a$) at the beginning of the convective region. Cases A and B yield close values of $\bar{\theta}_a$, which is also little sensitive to the model parameters $Re$ and $\bar{R}$. Indeed, $\bar{\theta}_a$ varies only from 0.25 to 0.6 when $Re$ and $\bar{R}$ are increased by three and two orders of magnitude, respectively. With

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Maximum thermal perturbation at the beginning of the convective region for cases A (Biot number estimated from forced convection correlations) and B (infinite Biot number). (a) Maximum annual mean temperature ($\bar{\theta}_w$ and $\bar{\theta}_a$ for the wall and the air, resp.). (b) Maximum annual temperature fluctuation ($\Delta\bar{\theta}_w$ and $\Delta\bar{\theta}_a$ for the wall and the air, resp.). The vertical dashed lines correspond to the limits of the laminar and turbulent regimes.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Amplitude contour of the thermal perturbation in a symmetry plane of the rock massif (case A with $\bar{R} = 0.189$ and $Re = 1.8 \times 10^5$). (a) Mode $k = 0$; the black solid lines correspond to the amplitude $\bar{\theta}_w/10 = 0.0425$ (n.b.: at the scale of the figure, the cave reduces to a line). (b) Mode $k = 1$; the black solid lines correspond to the amplitude $\Delta\bar{\theta}_w/10 = 0.030$ (n.b.: the top view focuses on the entrance region; the bottom view has been plotted with different scales on $\bar{x}$ and $\bar{r}$ to display the whole perturbed region).</figcaption></figure>

be estimated from the air temperature (length such that the amplitude of AMT and ATF in the air are equal to $\bar{\theta}_a/10$ and $\Delta\theta_a/10$, respectively).

# 5. Parametric Study

## 5.1. Amplitude of the Thermal Perturbation Induced by the Air Flow

Figure 7 displays the AMT ($\bar{\theta}_w$ and $\bar{\theta}_a$) and the ATF ($\Delta\theta_w$ and $\Delta\theta_a$) at the beginning of the convective region. Cases A and B yield close values of $\bar{\theta}_a$, which is also little sensitive to the model parameters $Re$ and $\bar{R}$. Indeed, $\bar{\theta}_a$ varies only from 0.25 to 0.6 when $Re$ and $\bar{R}$ are increased by three and two orders of magnitude, respectively. With

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Maximum thermal perturbation at the beginning of the convective region for cases A (Biot number estimated from forced convection correlations) and B (infinite Biot number). (a) Maximum annual mean temperature ($\bar{\theta}_w$ and $\bar{\theta}_a$ for the wall and the air, resp.). (b) Maximum annual temperature fluctuation ($\Delta\bar{\theta}_w$ and $\Delta\bar{\theta}_a$ for the wall and the air, resp.). The vertical dashed lines correspond to the limits of the laminar and turbulent regimes.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8. Characteristic lengths $\tilde{L}_0$ (a) and $\tilde{L}_1$ (b) as a function of the Reynolds number $Re$ for various radius $R$, computed from wall or gas temperatures. The dashed lines correspond to Equation 24. Symbols: the same legend as in Figure 7.</figcaption></figure>

the exception of small conduits ($R \approx 0.2$) at small Reynolds number ($Re \lesssim 3 \times 10^3$), a similar remark can be made on $\Delta\theta_a$.

The evolution of the wall temperature is more intricate. In case B, the assumption of infinite Biot number imposes $\bar{\theta}_w = \bar{\theta}_a$ and $\Delta\theta_w = \Delta\theta_a$ (since the air and wall temperatures are equal). Assuming that Case A provides a lower bound of the actual Biot number in a cave, Equation 14 allows us to estimate the consequences of an error in the determination of the Biot number on the estimation of the wall temperature at the beginning of the convective region. At high Reynolds number ($Re \gtrsim 3 \times 10^5$), no accurate estimation of the relation of the function $Nu(Re)$ is required because the heat transfer is completely driven by heat conduction in the rock massif. Conversely, the relation $Nu(Re)$ is a key information for convection-limited heat transfer at small Reynolds number ($Re \lesssim 10^4$). At intermediate $Re$, an accurate estimation of $Nu$ is required to get the ATF of the wall in a cave with a large radius ($R \gtrsim 1$).

## 5.2. Spatial Extent of the Thermal Perturbation Induced by the Air Flow

The characteristic lengths $\tilde{L}_0$ and $\tilde{L}_1$ for the damping of the AMT and ATF along the cave axis are displayed in Figure 8. The values of $\tilde{L}_0$ computed in case B (infinite Biot number) support the simple relation obtained from curve fitting:

$$ \tilde{L}_0 \approx 3.1 \times 10^{-3} Re^{1.06} R^{0.61} \quad (24) $$

In case A, estimating $\tilde{L}_0$ from the air or the wall temperatures give the same results (the lengths required to get $\bar{\theta}_a$ or $\bar{\theta}_w$ divided by 10 are the same, even when $\bar{\theta}_a$ and $\bar{\theta}_w$ are different). At high Reynolds number, $\tilde{L}_0$ estimated from case A follows Equation 24 derived for case B, as expected when heat transfer is conduction-limited. When $Re$ is decreased, using the forced convection correlations of case A reduces the heat flux at the cave wall compared to the infinite Biot number of case B, leading to larger values of $\tilde{L}_0$ than predicted by Equation 24.

$\tilde{L}_1$, the distance corresponding to the damping of the ATF, is nearly equal to $\tilde{L}_0$, with the exception of slight differences at low Reynolds number for $R \approx 2$ (see Figure 8b). The length of the convective region (called convection length in the following) can thus be indifferently determined from the damping of the AMT or ATF, based on the wall or the air temperatures. A higher bound of the convection length is provided by case A, while case B provides a lower bound. We can see in Figure 8 that the uncertainty on the convection length increases when the Reynolds number $Re$ is decreased and when the cave diameter $R$ is increased.

The width $\tilde{W}_0$ corresponding to the damping of the AMT in radial directions inside the rock mass is displayed in Figure 9. It is approximately proportional to $\tilde{L}_0$ and an order of magnitude smaller.
<figure><img src="image_10.png" /><figcaption>Figure 8. Characteristic lengths $\tilde{L}_0$ (a) and $\tilde{L}_1$ (b) as a function of the Reynolds number $Re$ for various radius $R$, computed from wall or gas temperatures. The dashed lines correspond to Equation 24. Symbols: the same legend as in Figure 7.</figcaption></figure>

the exception of small conduits ($R \approx 0.2$) at small Reynolds number ($Re \lesssim 3 \times 10^3$), a similar remark can be made on $\Delta\theta_a$.

The evolution of the wall temperature is more intricate. In case B, the assumption of infinite Biot number imposes $\bar{\theta}_w = \bar{\theta}_a$ and $\Delta\theta_w = \Delta\theta_a$ (since the air and wall temperatures are equal). Assuming that Case A provides a lower bound of the actual Biot number in a cave, Equation 14 allows us to estimate the consequences of an error in the determination of the Biot number on the estimation of the wall temperature at the beginning of the convective region. At high Reynolds number ($Re \gtrsim 3 \times 10^5$), no accurate estimation of the relation of the function $Nu(Re)$ is required because the heat transfer is completely driven by heat conduction in the rock massif. Conversely, the relation $Nu(Re)$ is a key information for convection-limited heat transfer at small Reynolds number ($Re \lesssim 10^4$). At intermediate $Re$, an accurate estimation of $Nu$ is required to get the ATF of the wall in a cave with a large radius ($R \gtrsim 1$).

## 5.2. Spatial Extent of the Thermal Perturbation Induced by the Air Flow

The characteristic lengths $\tilde{L}_0$ and $\tilde{L}_1$ for the damping of the AMT and ATF along the cave axis are displayed in Figure 8. The values of $\tilde{L}_0$ computed in case B (infinite Biot number) support the simple relation obtained from curve fitting:

$$ \tilde{L}_0 \approx 3.1 \times 10^{-3} Re^{1.06} R^{0.61} \qquad (24) $$

In case A, estimating $\tilde{L}_0$ from the air or the wall temperatures give the same results (the lengths required to get $\bar{\theta}_a$ or $\bar{\theta}_w$ divided by 10 are the same, even when $\bar{\theta}_a$ and $\bar{\theta}_w$ are different). At high Reynolds number, $\tilde{L}_0$ estimated from case A follows Equation 24 derived for case B, as expected when heat transfer is conduction-limited. When $Re$ is decreased, using the forced convection correlations of case A reduces the heat flux at the cave wall compared to the infinite Biot number of case B, leading to larger values of $\tilde{L}_0$ than predicted by Equation 24.

$\tilde{L}_1$, the distance corresponding to the damping of the ATF, is nearly equal to $\tilde{L}_0$, with the exception of slight differences at low Reynolds number for $R \approx 2$ (see Figure 8b). The length of the convective region (called convection length in the following) can thus be indifferently determined from the damping of the AMT or ATF, based on the wall or the air temperatures. A higher bound of the convection length is provided by case A, while case B provides a lower bound. We can see in Figure 8 that the uncertainty on the convection length increases when the Reynolds number $Re$ is decreased and when the cave diameter $R$ is increased.

The width $\tilde{W}_0$ corresponding to the damping of the AMT in radial directions inside the rock mass is displayed in Figure 9. It is approximately proportional to $\tilde{L}_0$ and an order of magnitude smaller.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Width $\tilde{W}_0$ corresponding to the damping of the annual mean temperature in the radial directions inside the rock mass.</figcaption></figure>

# 6. Comparison With Field Data

## 6.1. Baulmes Mine (Switzerland)

The Baulmes mine (46°47′33″N, 6°31′35″E, 655 m a.s.l.) is located in western Switzerland. This mine comprises a network of 11 sub horizontal levels with a cumulated length of about 17 km (Figure 10). The difference in elevation between the uppermost and lowest entrances is 90 m, ensuring a strong chimney effect. The sub horizontal conduits at intermediate elevations between the lower and upper entrances, connected by shafts, lead to dead ends in the horizontal directions (for more information about Baulmes mine, see Bourret et al. (2007)). The temperature delivered by Mathod weather station (46°44′13″N, 6°34′4″E, 435 m a.s.l., at 7 km SW from the study site) is used as an estimation of the external atmosphere temperature. A correction has been applied to account for the difference of elevation with the mine entrance, using the vertical thermal gradient -5.2°C/km (mean value from the network of Swiss meteorological stations for the time range 1991–2020). The resulting external temperature is displayed in Figure 11 for 1 year from 1st July 2022 at 0h00 (this date is taken as time $t = 0$ for all the data from Baulmes mine). The sampling period was 1 hr. The Fast Fourier Transform (FFT) of the external atmosphere temperature yielded $T_m = 10.4°C$ (after correction) and $\Delta T_y = 9.4°C$ for the external AMT and ATF, respectively.

The 360 m long tunnel #10 is connected horizontally to the lower entrance with no connection to other branches. This is the zone of interest of study. It is indicated by the yellow shading in Figure 10. This conduit has a typical horseshoe cross section. The geometrical characteristics of this conduit, including the mean perimeter and cross-sectional area, their standard deviation as well as its total length are summarized in Table 2. We instrumented this tunnel with seven temperature probes at different distances from the entrance. They are listed in Table 3. The first six temperature probes are Reefnet (Sensus Ultra) with ±0.3°C accuracy. The last one located at 360 m from

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Map of Baulmes mine (vertical section).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Width $\tilde{W}_0$ corresponding to the damping of the annual mean temperature in the radial directions inside the rock mass.</figcaption></figure>

# 6. Comparison With Field Data

## 6.1. Baulmes Mine (Switzerland)

The Baulmes mine (46°47′33″N, 6°31′35″E, 655 m a.s.l.) is located in western Switzerland. This mine comprises a network of 11 sub horizontal levels with a cumulated length of about 17 km (Figure 10). The difference in elevation between the uppermost and lowest entrances is 90 m, ensuring a strong chimney effect. The sub horizontal conduits at intermediate elevations between the lower and upper entrances, connected by shafts, lead to dead ends in the horizontal directions (for more information about Baulmes mine, see Bourret et al. (2007)). The temperature delivered by Mathod weather station (46°44′13″N, 6°34′4″E, 435 m a.s.l., at 7 km SW from the study site) is used as an estimation of the external atmosphere temperature. A correction has been applied to account for the difference of elevation with the mine entrance, using the vertical thermal gradient -5.2°C/km (mean value from the network of Swiss meteorological stations for the time range 1991–2020). The resulting external temperature is displayed in Figure 11 for 1 year from 1st July 2022 at 0h00 (this date is taken as time $t = 0$ for all the data from Baulmes mine). The sampling period was 1 hr. The Fast Fourier Transform (FFT) of the external atmosphere temperature yielded $T_m = 10.4°C$ (after correction) and $\Delta T_y = 9.4°C$ for the external AMT and ATF, respectively.

The 360 m long tunnel #10 is connected horizontally to the lower entrance with no connection to other branches. This is the zone of interest of study. It is indicated by the yellow shading in Figure 10. This conduit has a typical horseshoe cross section. The geometrical characteristics of this conduit, including the mean perimeter and cross-sectional area, their standard deviation as well as its total length are summarized in Table 2. We instrumented this tunnel with seven temperature probes at different distances from the entrance. They are listed in Table 3. The first six temperature probes are Reefnet (Sensus Ultra) with ±0.3°C accuracy. The last one located at 360 m from

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Map of Baulmes mine (vertical section).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** External atmosphere temperature estimated at Baulmes mine (daily averaged data). The horizontal arrow displays the time range during which the temperature inside the cave was measured by the sensors listed in Table 3. The vertical red bars correspond to the dates of the air flowrate measurements listed in Table 4. The black dashed line displays the mode *k* = 0 and *k* = 1 of the external temperature.</figcaption></figure>

entrance is HOBO (u22-001) with ±0.2°C accuracy. The temperature inside the conduit was recorded by these seven probes from 3rd February 2023 (*t* = 0.596 year) to 7th July 2023 (*t* = 1.017 year). This time range is displayed by the horizontal arrow in Figure 11. The sampling period was 1 hr.

The air mass flow rate was measured sporadically with a manual hot wire anemometer (Testo425 with an accuracy of 0.03 m/s ±5% of the measured value). The measurements were carried out at the four dates listed in Table 4 and displayed by vertical red bars in Figure 11. There are two dates in winter and two in summer. The mean flowrates in winter and summer are 2.4 and 5.0 kg/s, respectively. In the following, we consider Δ*ṁ* ~ (2.4 + 5.0)/2 ~ 4 kg/s as a rough estimation of the annual amplitude of the air flowrate. With the mean conduit perimeter *P* = 9.7 m from Table 2 and the air viscosity *μ*<sub>a</sub> ≈ 1.77 × 10<sup>-5</sup> Pa·s from Table 1, we get the order of magnitude of the Reynolds number *Re* = (4 *Δ* *ṁ*) / *μ*<sub>a*<sup>4</sup></sub> ~ 10<sup>5</sup>. The airflow is thus in the turbulent regime most of the time. The conduit radius reads *R* = *P*/(2π) ≈ 1.5 m that yields the scaled value *R̃* = *R*/*L*<sub>dif</sub> ≈ 0.3. The conduit is thus close to the configuration displayed in Figures 4 and 5 corresponding to *Re* = 1.8 × 10<sup>5</sup> and *R̃* = 0.189. In this configuration, cases A and B yield approximately the same results. The heat transfer is thus limited by conduction in the rock rather than convection. Therefore, Equation 24, for which conduction-limited heat transfer is assumed, is expected to provide a fair estimation of the convection length, and not only a lower bound. We get *L̃*<sub>0</sub> ≈ *L̃*<sub>1</sub> ~ 3 × 10<sup>2</sup>. The theory thus predicts a convection length significantly larger than the conduit length *L̃*<sub>c</sub> = *L*<sub>c</sub>/*L*<sub>dif</sub> ≈ 68.

The time series of the external atmosphere temperature and the temperatures inside the conduit have been plotted in Figure 12a. The dimensionless temperatures inside the conduit are always negative (i.e., lower than the external AMT), even in summer. This is the cold thermal anomaly expected at a lower entrance.

In the following, we assume that the temperature profiles θ<sub>w</sub>(x̄) at time t̄ = 0.611 and θ<sub>s</sub>(x̄) at time t̄ = 1.016 give an assessment of the coldest and hottest temperatures reached inside the conduit in winter and summer, respectively. These two times are displayed by vertical bars in Figure 12a, and the corresponding profiles are plotted in Figure 12b (the temperature at x̄ = 0 is the external atmosphere temperature). In winter, the dimensionless temperature continuously increases from −1.6 at the inlet (x̄ = 0) to −0.56 at the end of the conduit (x̄ = 68). Spatial temperature variations are thus significant over the entire length of the conduit. This suggests that the convection length should be larger than the conduit length. In summer, the reverse flow direction yields a

<table><caption>**Table 2**

Geometrical Characteristics of the Instrumented Conduit of Baulmes Mine (Length *L*<sub>c</sub>, Mean Value and Standard Deviation of the Perimeter *P* and the Cross-Sectional Area *A* From the Survey of 28 Cross-Sections)</caption><thead><tr><th rowspan="2">L<sub>c</sub> (m)</th><th colspan="2">P (m)</th><th colspan="2">A (m<sup>2</sup>)</th></tr><tr><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td>360</td><td>9.7</td><td>1.1</td><td>6.6</td><td>1.5</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** External atmosphere temperature estimated at Baulmes mine (daily averaged data). The horizontal arrow displays the time range during which the temperature inside the cave was measured by the sensors listed in Table 3. The vertical red bars correspond to the dates of the air flowrate measurements listed in Table 4. The black dashed line displays the mode *k* = 0 and *k* = 1 of the external temperature.</figcaption></figure>

entrance is HOBO (u22-001) with ±0.2°C accuracy. The temperature inside the conduit was recorded by these seven probes from 3rd February 2023 (*t* = 0.596 year) to 7th July 2023 (*t* = 1.017 year). This time range is displayed by the horizontal arrow in Figure 11. The sampling period was 1 hr.

The air mass flow rate was measured sporadically with a manual hot wire anemometer (Testo425 with an accuracy of 0.03 m/s ±5% of the measured value). The measurements were carried out at the four dates listed in Table 4 and displayed by vertical red bars in Figure 11. There are two dates in winter and two in summer. The mean flowrates in winter and summer are 2.4 and 5.0 kg/s, respectively. In the following, we consider Δ*ṁ* ~ (2.4 + 5.0)/2 ~ 4 kg/s as a rough estimation of the annual amplitude of the air flowrate. With the mean conduit perimeter *P* = 9.7 m from Table 2 and the air viscosity *μ*<sub>a</sub> ≈ 1.77 × 10<sup>-5</sup> Pa·s from Table 1, we get the order of magnitude of the Reynolds number *Re* = (4 *Δ* *ṁ* / *P* *μ*<sub>a</sub>) ~ 10<sup>5</sup>. The airflow is thus in the turbulent regime most of the time. The conduit radius reads *R* = *P*/(2π) ≈ 1.5 m that yields the scaled value *R̃* = *R*/*L*<sub>dif</sub> ≈ 0.3. The conduit is thus close to the configuration displayed in Figures 4 and 5 corresponding to *Re* = 1.8 × 10<sup>5</sup> and *R̃* = 0.189. In this configuration, cases A and B yield approximately the same results. The heat transfer is thus limited by conduction in the rock rather than convection. Therefore, Equation 24, for which conduction-limited heat transfer is assumed, is expected to provide a fair estimation of the convection length, and not only a lower bound. We get *L̃*<sub>0</sub> ≈ *L̃*<sub>1</sub> ~ 3 × 10<sup>2</sup>. The theory thus predicts a convection length significantly larger than the conduit length *L̃*<sub>c</sub> = *L*<sub>c</sub>/*L*<sub>dif</sub> ≈ 68.

The time series of the external atmosphere temperature and the temperatures inside the conduit have been plotted in Figure 12a. The dimensionless temperatures inside the conduit are always negative (i.e., lower than the external AMT), even in summer. This is the cold thermal anomaly expected at a lower entrance.

In the following, we assume that the temperature profiles θ<sub>w</sub>(x̄) at time t̄ = 0.611 and θ<sub>s</sub>(x̄) at time t̄ = 1.016 give an assessment of the coldest and hottest temperatures reached inside the conduit in winter and summer, respectively. These two times are displayed by vertical bars in Figure 12a, and the corresponding profiles are plotted in Figure 12b (the temperature at x̄ = 0 is the external atmosphere temperature). In winter, the dimensionless temperature continuously increases from −1.6 at the inlet (x̄ = 0) to −0.56 at the end of the conduit (x̄ = 68). Spatial temperature variations are thus significant over the entire length of the conduit. This suggests that the convection length should be larger than the conduit length. In summer, the reverse flow direction yields a

<table><caption>**Table 2**

Geometrical Characteristics of the Instrumented Conduit of Baulmes Mine (Length *L*<sub>c</sub>, Mean Value and Standard Deviation of the Perimeter *P* and the Cross-Sectional Area *A* From the Survey of 28 Cross-Sections)</caption><thead><tr><th rowspan="2">L<sub>c</sub> (m)</th><th colspan="2">P (m)</th><th colspan="2">A (m<sup>2</sup>)</th></tr><tr><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td>360</td><td>9.7</td><td>1.1</td><td>6.6</td><td>1.5</td></tr></tbody></table>
<table><caption>**Table 3**
Location of the Seven Temperature Probes Inside the Instrumented Conduit of Baulmes Mine</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th></tr></thead><tbody><tr><td>x (m)</td><td>10</td><td>26</td><td>50</td><td>90</td><td>128</td><td>220</td><td>360</td></tr><tr><td>x&#773; = x/L<sub>dif</sub></td><td>1.9</td><td>4.9</td><td>9.4</td><td>17</td><td>24</td><td>42</td><td>68</td></tr></tbody></table>

Note. x is the distance from the lower entrance.

different temperature profile. The temperature at the end of the conduit ($\bar{x}$ = 68) is -0.06, close to the external AMT. When getting closer to the entrance, the temperature decreases to -0.16 at the first probe ($\bar{x}$ = 1.9). The temperature then suddenly increases within the diffusive region, from -0.16 at the first probe ($\bar{x}$ = 1.9) to the positive value 0.78 at the entrance ($\bar{x}$ = 0).

The difference $\Delta\theta = (\theta_s - \theta_w)/2$ between summer and winter profiles yields an assessment of the amplitude of the ATF in the conduit. $\Delta\theta$ drops from 1.2 to 0.58 between the entrance ($\bar{x}$ = 0) and the first probe ($\bar{x}$ = 1.9). This sudden decrease takes place in the diffusive region. $\Delta\theta$ then smoothly decreases from

0.58 (first probe at $\bar{x}$ = 1.9) to 0.25 at the end of the conduit ($\bar{x}$ = 68). Assuming that the first sensor yields the temperature at the beginning of the convective region, the amplitude of the temperature fluctuation is thus divided by approximately 2.3 in the convective domain of the conduit. This confirms that the conduit length is significantly shorter than the convection length, defined as the distance required to divide by 10 the ATF measured at the beginning of the convective region (see Section 4).

## 6.2. D7.1 Cave (Sieben Hengste, Switzerland)

D7.1 cave is located in the Sieben Hengste karst area (46°45'43"N, 7°49'51"E, 1,750 m a.s.l.). It is part of a large ventilated cave system (>160 km of cumulated length) subject to chimney effect (see the map in Figure 13). Field observations reveal a seasonal airflow consistent most of the time with an upper entrance of the cave system (i.e., air inflow during the summer season). The ingress of warm external air during summer is marked by some condensation within the first 20 m of the meander, drained as a streamlet at the base of the meander. Otherwise, the cave passage remains remarkably dry. Unexpectedly, temperature and flow rate measurements show periods of time in winter during which the cave sporadically behaves as a lower entrance, with cold air inflow. The explanation of these unexpected flow reversals is still being searched.

External air temperatures displayed in Figure 14a are provided by an automatic weather station located at 1,850 m elevation (Swiss Federal Institute for Snow and Avalanche Research), approximately 2 km SE of our study site (sampling period: 1 hr). The external atmosphere temperature was corrected from the elevation as in Section 6.1. The FFT yielded $T_m = 5.0^{\circ}\text{C}$ (after correction) and $\Delta T_y = 7.1^{\circ}\text{C}$ for the external AMT and ATF, respectively.

The ca. 100 m-long entrance meander of D7.1 Cave has a simple geometry with an average hydraulic diameter of ca. 1 m, despite local variations in the cross-sectional area and perimeter. The perimeter $P \approx 3$ m and the scaled radius $R \approx 0.1$ will be considered in the calculations. The meander was instrumented for temperature and airflow measurements between 25/08/2021 and 10/08/2022. Five temperature probes were installed along the main conduit at distances from the cave entrance listed in Table 5. Probes #1 to #4 were Reefnets Sensus Ultra with $\pm 0.3^{\circ}\text{C}$ accuracy (sampling period: 30 min). Probe #5 was a HOBO U22-001 with $\pm 0.2^{\circ}\text{C}$ accuracy (sampling period: 2 hr). It is located at the base of a 20 m shaft following the entrance meander. Time series from the five probes are displayed in Figure 14b.

The cave airflow was measured using a low pressure drop flowmeter (Sensirion SFM 3000 with accuracy of 5% and 20% when velocity is higher and lower than 0.3 m/s, respectively; 50 m from cave entrance) and validated with independent flow measurements (CO₂ gauging) on 10/08/2022. The FFT of the air flow yields the order of magnitude $\Delta\dot{m} \sim 0.1$ kg/s rate (red dashed line in Figure 14a), which results in $Re \sim 7 \times 10^3$. The cave is thus most of the time in the transition between laminar and turbulent regimes. Equation 24 yields the theoretical convection length $\bar{L}_0 \approx \bar{L}_1 \sim 10$. It is important to note that this value must be considered as a lower bound of the convection

length. Indeed, Equation 24 stands for case B (infinite Biot number). At $Re \sim 7 \times 10^3$, case A (finite Biot number) yields a convection length significantly larger than case B (see Figure 8).

The AMT inside D7.1 Cave is displayed in Figure 17a. It decreases to reach a negative quasi-plateau in the range from -0.26 to -0.23 in between probes #3 ($\bar{x}$ = 6.2), and #5 ($\bar{x}$ = 24). This is not the behavior expected for an upper entrance, for which the AMT is supposed to be positive in the convective region before decreasing to zero at some distance from the entrance. Sporadic flow reversals in winter can contribute to this negative

<table><caption>**Table 4**
Air Mass Flow Rate Measured in Winter and Summer</caption><thead><tr><th rowspan="2">Date</th><th colspan="2">Winter (upward flow)</th><th colspan="2">Summer (downward flow)</th></tr><tr><th>02/02/2023</th><th>15/02/2023</th><th>06/06/2023</th><th>07/07/2023</th></tr></thead><tbody><tr><td>t&#773; = t/τ<sub>y</sub></td><td>0.595</td><td>0.628</td><td>0.933</td><td>1.018</td></tr><tr><td>ṁ (kg/s)</td><td>1.9</td><td>2.8</td><td>4.5</td><td>5.4</td></tr></tbody></table>
<table><caption>**Table 3**
Location of the Seven Temperature Probes Inside the Instrumented Conduit of Baulmes Mine</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th></tr></thead><tbody><tr><td>x (m)</td><td>10</td><td>26</td><td>50</td><td>90</td><td>128</td><td>220</td><td>360</td></tr><tr><td>x&#773; = x/L<sub>dif</sub></td><td>1.9</td><td>4.9</td><td>9.4</td><td>17</td><td>24</td><td>42</td><td>68</td></tr></tbody></table>

Note. x is the distance from the lower entrance.

different temperature profile. The temperature at the end of the conduit ($\bar{x}$ = 68) is -0.06, close to the external AMT. When getting closer to the entrance, the temperature decreases to -0.16 at the first probe ($\bar{x}$ = 1.9). The temperature then suddenly increases within the diffusive region, from -0.16 at the first probe ($\bar{x}$ = 1.9) to the positive value 0.78 at the entrance ($\bar{x}$ = 0).

The difference $\Delta\theta = (\theta_s - \theta_w)/2$ between summer and winter profiles yields an assessment of the amplitude of the ATF in the conduit. $\Delta\theta$ drops from 1.2 to 0.58 between the entrance ($\bar{x}$ = 0) and the first probe ($\bar{x}$ = 1.9). This sudden decrease takes place in the diffusive region. $\Delta\theta$ then smoothly decreases from

0.58 (first probe at $\bar{x}$ = 1.9) to 0.25 at the end of the conduit ($\bar{x}$ = 68). Assuming that the first sensor yields the temperature at the beginning of the convective region, the amplitude of the temperature fluctuation is thus divided by approximately 2.3 in the convective domain of the conduit. This confirms that the conduit length is significantly shorter than the convection length, defined as the distance required to divide by 10 the ATF measured at the beginning of the convective region (see Section 4).

## 6.2. D7.1 Cave (Sieben Hengste, Switzerland)

D7.1 cave is located in the Sieben Hengste karst area (46°45'43"N, 7°49'51"E, 1,750 m a.s.l.). It is part of a large ventilated cave system (>160 km of cumulated length) subject to chimney effect (see the map in Figure 13). Field observations reveal a seasonal airflow consistent most of the time with an upper entrance of the cave system (i.e., air inflow during the summer season). The ingress of warm external air during summer is marked by some condensation within the first 20 m of the meander, drained as a streamlet at the base of the meander. Otherwise, the cave passage remains remarkably dry. Unexpectedly, temperature and flow rate measurements show periods of time in winter during which the cave sporadically behaves as a lower entrance, with cold air inflow. The explanation of these unexpected flow reversals is still being searched.

External air temperatures displayed in Figure 14a are provided by an automatic weather station located at 1,850 m elevation (Swiss Federal Institute for Snow and Avalanche Research), approximately 2 km SE of our study site (sampling period: 1 hr). The external atmosphere temperature was corrected from the elevation as in Section 6.1. The FFT yielded $T_m = 5.0^{\circ}\text{C}$ (after correction) and $\Delta T_y = 7.1^{\circ}\text{C}$ for the external AMT and ATF, respectively.

The ca. 100 m-long entrance meander of D7.1 Cave has a simple geometry with an average hydraulic diameter of ca. 1 m, despite local variations in the cross-sectional area and perimeter. The perimeter $P \approx 3$ m and the scaled radius $R \approx 0.1$ will be considered in the calculations. The meander was instrumented for temperature and airflow measurements between 25/08/2021 and 10/08/2022. Five temperature probes were installed along the main conduit at distances from the cave entrance listed in Table 5. Probes #1 to #4 were Reefnets Sensus Ultra with $\pm 0.3^{\circ}\text{C}$ accuracy (sampling period: 30 min). Probe #5 was a HOBO U22-001 with $\pm 0.2^{\circ}\text{C}$ accuracy (sampling period: 2 hr). It is located at the base of a 20 m shaft following the entrance meander. Time series from the five probes are displayed in Figure 14b.

The cave airflow was measured using a low pressure drop flowmeter (Sensirion SFM 3000 with accuracy of 5% and 20% when velocity is higher and lower than 0.3 m/s, respectively; 50 m from cave entrance) and validated with independent flow measurements (CO₂ gauging) on 10/08/2022. The FFT of the air flow yields the order of magnitude $\Delta\dot{m} \sim 0.1$ kg/s rate (red dashed line in Figure 14a), which results in $Re \sim 7 \times 10^3$. The cave is thus most of the time in the transition between laminar and turbulent regimes. Equation 24 yields the theoretical convection length $L_0 \approx L_1 \sim 10$. It is important to note that this value must be considered as a lower bound of the convection

length. Indeed, Equation 24 stands for case B (infinite Biot number). At $Re \sim 7 \times 10^3$, case A (finite Biot number) yields a convection length significantly larger than case B (see Figure 8).

The AMT inside D7.1 Cave is displayed in Figure 17a. It decreases to reach a negative quasi-plateau in the range from -0.26 to -0.23 in between probes #3 ($\bar{x}$ = 6.2), and #5 ($\bar{x}$ = 24). This is not the behavior expected for an upper entrance, for which the AMT is supposed to be positive in the convective region before decreasing to zero at some distance from the entrance. Sporadic flow reversals in winter can contribute to this negative

<table><caption>**Table 4**
Air Mass Flow Rate Measured in Winter and Summer</caption><thead><tr><th rowspan="2">Date</th><th colspan="2">Winter (upward flow)</th><th colspan="2">Summer (downward flow)</th></tr><tr><th>02/02/2023</th><th>15/02/2023</th><th>06/06/2023</th><th>07/07/2023</th></tr></thead><tbody><tr><td>t&#773; = t/τ<sub>y</sub></td><td>0.595</td><td>0.628</td><td>0.933</td><td>1.018</td></tr><tr><td>ṁ (kg/s)</td><td>1.9</td><td>2.8</td><td>4.5</td><td>5.4</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Temperatures in Baulmes conduit. (a) Daily averaged temperatures as a function of time (the temperature at $\tilde{x} = 0$ is the external atmosphere temperature measured at Mathod weather station). (b) Temperature profiles $\theta_s(\tilde{x})$ and $\theta_w(\tilde{x})$ measured respectively in summer and winter, and amplitude of the annual variation $\Delta\theta = (\theta_s - \theta_w)/2$. Dimensionless temperatures defined as $\theta = (T - T_m)/\Delta T_y$.</figcaption></figure>

AMT (an instance of such an event can be observed around time $\tilde{t} \approx 0.33$ in Figure 14). Moreover, the cave entrance is in a topographic low, resulting in local cold air trapping during winter. In addition, snow is likely to accumulate near the cave entrance. Both may induce local free convection cells in the entrance zone in winter resulting in the cooling of the cave below the external AMT.

If we assess $\tilde{L}_0$ from the distance over which the AMT significantly varies inside the cave, it would range between $\tilde{x} = 2.1$ (probe #2) and $\tilde{x} = 6.2$ (probe #3). This estimation is confirmed by the variation of the ATF with the distance from the entrance (Figure 17b). The ATF at probe #1 ($\tilde{x} = 0.38$), located in the diffusive region, is close to 1. Then it decreases to 0.36 at probe #2 ($\tilde{x} = 2.1$), the first probe in the convective region, and drops down to 0.02 at probe #3 ($\tilde{x} = 6.2$). The ATF being reduced by a factor larger than 10 between probes #2 and #3, $\tilde{L}_1$ should be in between these two probes. It is worthwhile to note that this estimation is consistent with the limit of condensation, which falls in between probes #2 and #3 (20 m corresponds to $\tilde{x} \approx 4$). The convection length of D7.1 Cave is thus in the range from 2.1 to 6.2.

## 6.3. Bärenhöhle Cave (Austria)

The cave is located in Vorarlberg, Austria (47°22′17″N, 9°52′52″E, 901 m asl). The external atmosphere temperature is obtained from the weather station located at St. Gallen (47°25′31.7094″N, 9°23′54.7008″E, 776 m a.s.l.), Switzerland about 37 km West-Northwest of our study site (the closer Schopernau station (47°19′N, 10°1′E, 835 m a.s.l.; The Central Institute for Meteorology and Geodynamics of Austria 2012) is subject to thermal inversions during the winter season and thus considered as being less representative). The temperature measured at St. Gallen has been corrected from the elevation as in Section 6.1. It is displayed in Figure 16a. The FFT yielded $T_m = 8.7°C$ (after correction) and $\Delta T_y = 7.6°C$ for the external AMT and ATF, respectively.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Map of D7.1 Cave. The red circles indicate temperature sensors.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Temperatures in Baulmes conduit. (a) Daily averaged temperatures as a function of time (the temperature at $\tilde{x} = 0$ is the external atmosphere temperature measured at Mathod weather station). (b) Temperature profiles $\theta_s(\tilde{x})$ and $\theta_w(\tilde{x})$ measured respectively in summer and winter, and amplitude of the annual variation $\Delta\theta = (\theta_s - \theta_w)/2$. Dimensionless temperatures defined as $\theta = (T - T_m)/\Delta T_y$.</figcaption></figure>

AMT (an instance of such an event can be observed around time $\tilde{t} \approx 0.33$ in Figure 14). Moreover, the cave entrance is in a topographic low, resulting in local cold air trapping during winter. In addition, snow is likely to accumulate near the cave entrance. Both may induce local free convection cells in the entrance zone in winter resulting in the cooling of the cave below the external AMT.

If we assess $\tilde{L}_0$ from the distance over which the AMT significantly varies inside the cave, it would range between $\tilde{x} = 2.1$ (probe #2) and $\tilde{x} = 6.2$ (probe #3). This estimation is confirmed by the variation of the ATF with the distance from the entrance (Figure 17b). The ATF at probe #1 ($\tilde{x} = 0.38$), located in the diffusive region, is close to 1. Then it decreases to 0.36 at probe #2 ($\tilde{x} = 2.1$), the first probe in the convective region, and drops down to 0.02 at probe #3 ($\tilde{x} = 6.2$). The ATF being reduced by a factor larger than 10 between probes #2 and #3, $\tilde{L}_1$ should be in between these two probes. It is worthwhile to note that this estimation is consistent with the limit of condensation, which falls in between probes #2 and #3 (20 m corresponds to $\tilde{x} \approx 4$). The convection length of D7.1 Cave is thus in the range from 2.1 to 6.2.

## 6.3. Bärenhöhle Cave (Austria)

The cave is located in Vorarlberg, Austria (47°22′17″N, 9°52′52″E, 901 m asl). The external atmosphere temperature is obtained from the weather station located at St. Gallen (47°25′31.7094″N, 9°23′54.7008″E, 776 m a.s.l.), Switzerland about 37 km West-Northwest of our study site (the closer Schopernau station (47°19′N, 10°1′E, 835 m a.s.l.; The Central Institute for Meteorology and Geodynamics of Austria 2012) is subject to thermal inversions during the winter season and thus considered as being less representative). The temperature measured at St. Gallen has been corrected from the elevation as in Section 6.1. It is displayed in Figure 16a. The FFT yielded $T_m = 8.7°C$ (after correction) and $\Delta T_y = 7.6°C$ for the external AMT and ATF, respectively.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Map of D7.1 Cave. The red circles indicate temperature sensors.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Field data from D7.1 cave (all the data are daily averaged). (a) External temperature and mass flow rate; black and red dashed lines represent the modes *k* = 0 and *k* = 1 computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively. (b) Dimensionless temperatures outside and inside the cave ($\theta = (T - T_m)/\Delta T_y$).</figcaption></figure>

Chimney effect drives the airflow in the cave which behaves as a lower entrance (air inflow in winter, see Figure 15). The difference in elevation between the lower and upper entrances is ca. 130 m for a conduit length along the main air passage approximating 240 m (Perkmann & Luetscher, 2013). The average hydraulic diameter is assessed to about 1 m, so that $P \approx 3$ m and $\tilde{R} \approx 0.1$. Air temperature was recorded at 1h interval using three HOBO (U22-001) temperature probes at distances from the entrance listed in Table 6. The corresponding time series are displayed in Figure 16b.

Air flow was recorded between 16/04/2013 and 16/04/2014 at the lower entrance using a thermo-anemometer calibrated to $\pm 0.1$ m/s. The FFT of the air flow rate (red dashed line in Figure 16a) yields the orders of magnitude $\Delta \dot{m} \sim 0.6$ kg/s and $Re \sim 5 \times 10^4$. The cave is thus in the turbulent regime most of the time. In this configuration ($\tilde{R} \approx 0.1$ and $Re \sim 5 \times 10^4$), cases A and B predict approximately the same theoretical convection length (see Figure 8). Equation 24 yields $\tilde{L}_0 \approx \tilde{L}_1 \sim 70$.

The AMT of Bärenhöhle Cave is displayed in Figure 17a. Negative values are consistent with a lower entrance. The lowest AMT ($-0.27$) is measured at probe #2 ($\tilde{x} = 16$). It increases slightly between probe #2 and the last probe #3. The ATF displayed in Figure 17b decreases over the whole instrumented conduit, from 1 outside the cave to 0.11 at the last probe ($\tilde{x} = 23$). These two pieces of information suggest that the convection length should be longer than 23, the cave length.

## 7. Discussion

The main results from the three investigated sites are summarized in Table 7. For two sites (Baulmes mine and Bärenhöhle Cave), we do not observe any contradiction between the model and the field data since both yield a convection length larger than the instrumented conduit (a hundred meters for Bärenhöhle, a few hundreds of meters for Baulmes). In D7.1 Cave, the field data suggest a shorter convection length, in the range from 10 to 30 m, whereas the theoretical prediction yields a lower bound of 50 m. Although the theoretical order of magnitude is correct, the theory overestimates the convection length. In a cave with a small Reynolds number and

a short convection length, local effects in the entrance region, not considered in the model, could overcome the effect of the mean flow generated by the chimney effect. In the case of D7.1 Cave, cold air ingress in winter might cool down the cave passage more than considered in the model, thus producing an actual shorter thermal length.

Using a theoretical approach, we defined the characteristic lengths $\tilde{L}_0$ and $\tilde{L}_1$ from the damping in the convective region of the AMT (mode $k = 0$) and ATF (mode $k = 1$), respectively. $\tilde{L}_0$ and $\tilde{L}_1$ are close to each other, and similar results are obtained if the estimation of these lengths is based on the wall

<table><caption>**Table 5**
Location of the Five Temperature Probes Inside the Instrumented Conduit of D7.1 Cave</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th></tr></thead><tbody><tr><td>x (m)</td><td>2</td><td>11</td><td>33</td><td>93</td><td>125</td></tr><tr><td>x&#773; = x/L<sub>dif</sub></td><td>0.38</td><td>2.1</td><td>6.2</td><td>18</td><td>24</td></tr></tbody></table>

Note. $x$ is the distance from the entrance.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Field data from D7.1 cave (all the data are daily averaged). (a) External temperature and mass flow rate; black and red dashed lines represent the modes *k* = 0 and *k* = 1 computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively. (b) Dimensionless temperatures outside and inside the cave ($\theta = (T - T_m)/\Delta T_y$).</figcaption></figure>

Chimney effect drives the airflow in the cave which behaves as a lower entrance (air inflow in winter, see Figure 15). The difference in elevation between the lower and upper entrances is ca. 130 m for a conduit length along the main air passage approximating 240 m (Perkmann & Luetscher, 2013). The average hydraulic diameter is assessed to about 1 m, so that $P \approx 3$ m and $\tilde{R} \approx 0.1$. Air temperature was recorded at 1h interval using three HOBO (U22-001) temperature probes at distances from the entrance listed in Table 6. The corresponding time series are displayed in Figure 16b.

Air flow was recorded between 16/04/2013 and 16/04/2014 at the lower entrance using a thermo-anemometer calibrated to $\pm 0.1$ m/s. The FFT of the air flow rate (red dashed line in Figure 16a) yields the orders of magnitude $\Delta \dot{m} \sim 0.6$ kg/s and $Re \sim 5 \times 10^4$. The cave is thus in the turbulent regime most of the time. In this configuration ($\tilde{R} \approx 0.1$ and $Re \sim 5 \times 10^4$), cases A and B predict approximately the same theoretical convection length (see Figure 8). Equation 24 yields $\tilde{L}_0 \approx \tilde{L}_1 \sim 70$.

The AMT of Bärenhöhle Cave is displayed in Figure 17a. Negative values are consistent with a lower entrance. The lowest AMT ($-0.27$) is measured at probe #2 ($\tilde{x} = 16$). It increases slightly between probe #2 and the last probe #3. The ATF displayed in Figure 17b decreases over the whole instrumented conduit, from 1 outside the cave to 0.11 at the last probe ($\tilde{x} = 23$). These two pieces of information suggest that the convection length should be longer than 23, the cave length.

## 7. Discussion

The main results from the three investigated sites are summarized in Table 7. For two sites (Baulmes mine and Bärenhöhle Cave), we do not observe any contradiction between the model and the field data since both yield a convection length larger than the instrumented conduit (a hundred meters for Bärenhöhle, a few hundreds of meters for Baulmes). In D7.1 Cave, the field data suggest a shorter convection length, in the range from 10 to 30 m, whereas the theoretical prediction yields a lower bound of 50 m. Although the theoretical order of magnitude is correct, the theory overestimates the convection length. In a cave with a small Reynolds number and

a short convection length, local effects in the entrance region, not considered in the model, could overcome the effect of the mean flow generated by the chimney effect. In the case of D7.1 Cave, cold air ingress in winter might cool down the cave passage more than considered in the model, thus producing an actual shorter thermal length.

Using a theoretical approach, we defined the characteristic lengths $\tilde{L}_0$ and $\tilde{L}_1$ from the damping in the convective region of the AMT (mode $k = 0$) and ATF (mode $k = 1$), respectively. $\tilde{L}_0$ and $\tilde{L}_1$ are close to each other, and similar results are obtained if the estimation of these lengths is based on the wall

<table><caption>**Table 5**
Location of the Five Temperature Probes Inside the Instrumented Conduit of D7.1 Cave</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th></tr></thead><tbody><tr><td>x (m)</td><td>2</td><td>11</td><td>33</td><td>93</td><td>125</td></tr><tr><td>x&#773; = x/L<sub>dif</sub></td><td>0.38</td><td>2.1</td><td>6.2</td><td>18</td><td>24</td></tr></tbody></table>

Note. $x$ is the distance from the entrance.
<figure><img src="image_17.png" /><figcaption>Figure 15. Map of Bärenhöhle. The red circles indicate temperature sensors.</figcaption></figure>

temperature or that of the air (see Figure 8). A single convection length is thus sufficient to characterize the extent of the AMT and the ATF along the cave axis. The case B (infinite Biot number) provides a lower bound of the convection length through the simple Expression 24. If the Reynolds number is large enough (typically, if $Re > 3 \times 10^4$ for a conduit with a radius of 1 m, see Figure 8), the approximation of infinite Biot number is valid and Equation 24 is expected to give a right estimation of the convection length. Turning this relation into unscaled variables yields

$$L_0 \approx L_1 \approx 4.0 \times 10^2 \frac{\Delta \dot{m}^{1.06}}{R^{0.45}} \quad (25)$$

where $L_0$, $L_1$, and $R$ are in meters and $\Delta \dot{m}$ in kg/s. Equation 25 gives the convection length resulting from the balance between heat advection by the air flow and heat conduction in the rock mass, in the conduction-limited regime. In this case, the convection length is approximately proportional to the amplitude of the flowrate annual fluctuations $\Delta \dot{m}$ divided by the square root of the cave radius $R$. Ventilated conduits with $R$ of the order of 1 m traversed by an airflow with $\Delta \dot{m}$ ranging from 0.1 to 1 kg/s are commonly encountered in mines and caves (see the case studies of Section 6). For such conduits, Equation 25 yields a convection length along the conduit axis

<figure><img src="image_18.png" /><figcaption>Figure 16. Field data of Bärenhöhle Cave (all the data are daily averaged). (a) External temperature and mass flow rate; black and red dashed lines represent the modes $k = 0$ and $k = 1$ computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively. (b) Dimensionless temperatures outside and inside the cave ($\theta = (T - T_m)/\Delta T_y$).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15. Map of Bärenhöhle. The red circles indicate temperature sensors.</figcaption></figure>

temperature or that of the air (see Figure 8). A single convection length is thus sufficient to characterize the extent of the AMT and the ATF along the cave axis. The case B (infinite Biot number) provides a lower bound of the convection length through the simple Expression 24. If the Reynolds number is large enough (typically, if $Re > 3 \times 10^4$ for a conduit with a radius of 1 m, see Figure 8), the approximation of infinite Biot number is valid and Equation 24 is expected to give a right estimation of the convection length. Turning this relation into unscaled variables yields

$$L_0 \approx L_1 \approx 4.0 \times 10^2 \frac{\Delta \dot{m}^{1.06}}{R^{0.45}} \quad (25)$$

where $L_0$, $L_1$, and $R$ are in meters and $\Delta \dot{m}$ in kg/s. Equation 25 gives the convection length resulting from the balance between heat advection by the air flow and heat conduction in the rock mass, in the conduction-limited regime. In this case, the convection length is approximately proportional to the amplitude of the flowrate annual fluctuations $\Delta \dot{m}$ divided by the square root of the cave radius $R$. Ventilated conduits with $R$ of the order of 1 m traversed by an airflow with $\Delta \dot{m}$ ranging from 0.1 to 1 kg/s are commonly encountered in mines and caves (see the case studies of Section 6). For such conduits, Equation 25 yields a convection length along the conduit axis

<figure><img src="image_18.png" /><figcaption>Figure 16. Field data of Bärenhöhle Cave (all the data are daily averaged). (a) External temperature and mass flow rate; black and red dashed lines represent the modes $k = 0$ and $k = 1$ computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively. (b) Dimensionless temperatures outside and inside the cave ($\theta = (T - T_m)/\Delta T_y$).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 17. Annual mean temperature (a) and annual temperature fluctuation (b) in D7.1 and Bärenhöhle Caves.</figcaption></figure>

between a few tens and a few hundreds of meters. Thermal perturbations also propagate in the radial direction, that is, along the directions perpendicular to the conduit. As expected, the amplitude of the annual fluctuations vanishes over the diffusion length, that is, a few meters. Conversely, the thermal anomaly (defined as the shift of the AMT), propagates along distances of the order of $L_0/10$ (see Figure 9), that is, over distances of the order of a few meters to a few tens of meters. It is worthwhile to note that in some cases, the extent of the thermal perturbation might be much larger than mentioned above. Luetscher and Jeannin (2004) measured air flow rates around 10 kg/s in Hölloch Cave (Switzerland) and La Diau (France). Recently, airflow rates around 40 kg/s were measured in Shpella Shtares Cave, in Albania (42.31766° N, 19.85403°E, 1,427 m a.s.l., (see Pastore et al. (2019) for more details about this cave). Taking a couple of meters as a rough estimation of the conduit radius in Shpella Shtares, the thermal perturbation should extend over kilometers and hundreds of meters in the axial and radial directions, respectively.

Assessing the thermal penetration length at different time-scales is fundamental to quantify bio-geochemical processes in karst systems. Our study demonstrates that seasonal temperature fluctuations propagate from a few tens of meters up to several thousands of meters into a karst system, depending on ventilation rate and conduit geometry. By way of illustration, let us cite the large ice caves in Austria which host perennial ice formations over several hundreds of meters in the lower entrance zone of ventilated cave systems (e.g., Eisriesenwelt (Obleitner & Spötl, 2011)). However, it is the advent of speleothems as precisely and accurately dated paleoclimate archives which rises the significance of our study. Stalagmites and flowstones are commonly used to reconstruct climate fluctuations based on their isotope (oxygen and carbon) and trace element concentrations (Fairchild & Baker, 2012). In both cases, the proxy partitioning depends strongly on the cave temperature. While it is commonly assumed that the cave temperature at depth is constant, this is only true beyond the convection length which strongly depends on the time-scale of interest.

Because the temperature-dependent oxygen isotope fractionation between water and calcite ranges between -0.18 and -0.21‰/°C (Kim & O'Neil, 1997; Tremaine et al., 2011), seasonal temperature changes >0.5°C represent a significant source of uncertainty in the speleothem-proxy interpretation and it might be of interest to sample beyond the convective region. Meanwhile, it can be strategic to focus speleothem research on the convective region associated with decadal to millennial cycles, as this area will respond sensitively to outside climate changes without picking up short-term temperature fluctuations.

## 8. Conclusion

We developed a model coupling buoyancy-induced convection (chimney effect) in a single karst conduit with conduction in the rock mass. We used this model to investigate the propagation of thermal perturbations inside a karst massif. Assuming dry air and a simplified geometry, and reducing the time variations of the external temperature to the annual fluctuations, we performed a parametric study to identify general trends regarding the effect of the air flowrate and conduit size on the amplitude and relaxation length of the

<table><caption>**Table 6**
Location of the Three Temperature Probes Inside Bärenhöhle</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th></tr></thead><tbody><tr><td>x (m)</td><td>5</td><td>85</td><td>120</td></tr><tr><td>x&#772; = x/L<sub>dif</sub></td><td>0.95</td><td>16</td><td>23</td></tr></tbody></table>

Note. x is the distance from the lower entrance.
<figure><img src="image_19.png" /><figcaption>Figure 17. Annual mean temperature (a) and annual temperature fluctuation (b) in D7.1 and Bärenhöhle Caves.</figcaption></figure>

between a few tens and a few hundreds of meters. Thermal perturbations also propagate in the radial direction, that is, along the directions perpendicular to the conduit. As expected, the amplitude of the annual fluctuations vanishes over the diffusion length, that is, a few meters. Conversely, the thermal anomaly (defined as the shift of the AMT), propagates along distances of the order of $L_0/10$ (see Figure 9), that is, over distances of the order of a few meters to a few tens of meters. It is worthwhile to note that in some cases, the extent of the thermal perturbation might be much larger than mentioned above. Luetscher and Jeannin (2004) measured air flow rates around 10 kg/s in Hölloch Cave (Switzerland) and La Diau (France). Recently, airflow rates around 40 kg/s were measured in Shpella Shtares Cave, in Albania (42.31766° N, 19.85403°E, 1,427 m a.s.l., (see Pastore et al. (2019) for more details about this cave). Taking a couple of meters as a rough estimation of the conduit radius in Shpella Shtares, the thermal perturbation should extend over kilometers and hundreds of meters in the axial and radial directions, respectively.

Assessing the thermal penetration length at different time-scales is fundamental to quantify bio-geochemical processes in karst systems. Our study demonstrates that seasonal temperature fluctuations propagate from a few tens of meters up to several thousands of meters into a karst system, depending on ventilation rate and conduit geometry. By way of illustration, let us cite the large ice caves in Austria which host perennial ice formations over several hundreds of meters in the lower entrance zone of ventilated cave systems (e.g., Eisriesenwelt (Obleitner & Spötl, 2011)). However, it is the advent of speleothems as precisely and accurately dated paleoclimate archives which rises the significance of our study. Stalagmites and flowstones are commonly used to reconstruct climate fluctuations based on their isotope (oxygen and carbon) and trace element concentrations (Fairchild & Baker, 2012). In both cases, the proxy partitioning depends strongly on the cave temperature. While it is commonly assumed that the cave temperature at depth is constant, this is only true beyond the convection length which strongly depends on the time-scale of interest.

Because the temperature-dependent oxygen isotope fractionation between water and calcite ranges between -0.18 and -0.21‰/°C (Kim & O'Neil, 1997; Tremaine et al., 2011), seasonal temperature changes >0.5°C represent a significant source of uncertainty in the speleothem-proxy interpretation and it might be of interest to sample beyond the convective region. Meanwhile, it can be strategic to focus speleothem research on the convective region associated with decadal to millennial cycles, as this area will respond sensitively to outside climate changes without picking up short-term temperature fluctuations.

## 8. Conclusion

We developed a model coupling buoyancy-induced convection (chimney effect) in a single karst conduit with conduction in the rock mass. We used this model to investigate the propagation of thermal perturbations inside a karst massif. Assuming dry air and a simplified geometry, and reducing the time variations of the external temperature to the annual fluctuations, we performed a parametric study to identify general trends regarding the effect of the air flowrate and conduit size on the amplitude and relaxation length of the

<table><caption>**Table 6**
Location of the Three Temperature Probes Inside Bärenhöhle</caption><thead><tr><th>Probe</th><th>#1</th><th>#2</th><th>#3</th></tr></thead><tbody><tr><td>x (m)</td><td>5</td><td>85</td><td>120</td></tr><tr><td>x&#772; = x/L<sub>dif</sub></td><td>0.95</td><td>16</td><td>23</td></tr></tbody></table>

Note. x is the distance from the lower entrance.
<table><caption>**Table 7**
Summary Table of the Three Investigated Sites</caption><thead><tr><th></th><th></th><th>Baulmes mine</th><th>D7.1 cave</th><th>Bärenhöhle Cave</th></tr></thead><tbody><tr><td>R</td><td></td><td>1.5 m</td><td>0.5 m</td><td>0.5 m</td></tr><tr><td>R̃</td><td></td><td>0.3</td><td>0.1</td><td>0.1</td></tr><tr><td>Δṁ</td><td></td><td>4 kg/s</td><td>0.1 kg/s</td><td>0.6 kg/s</td></tr><tr><td>Re</td><td></td><td>10<sup>5</sup></td><td>7 × 10<sup>3</sup></td><td>5 × 10<sup>4</sup></td></tr><tr><td>L̃<sub>0</sub> ≈ L̃<sub>1</sub></td><td>Field data</td><td>>68</td><td>[2.1, 6.2]</td><td>>23</td></tr><tr><td>L̃<sub>0</sub> ≈ L̃<sub>1</sub></td><td>Theory</td><td>~3 × 10<sup>2</sup></td><td>>10</td><td>~70</td></tr><tr><td>L<sub>0</sub> ≈ L<sub>1</sub></td><td>Field data</td><td>>360 m</td><td>[10 m, 30 m]</td><td>>120 m</td></tr><tr><td>L<sub>0</sub> ≈ L<sub>1</sub></td><td>Theory</td><td>~2 × 10<sup>3</sup> m</td><td>>50 m</td><td>~4 × 10<sup>2</sup> m</td></tr></tbody></table>

Note. *R*: Cave radius in meters; *R̃*: dimensionless cave radius; *Δṁ*: amplitude of the air flow rate in kg/s; *Re*: Reynolds number; *L̃*<sub>0</sub> and *L̃*<sub>1</sub>: dimensionless convection lengths deduced from the AMT and the ATF, respectively; *L*<sub>0</sub> and *L*<sub>1</sub>: same variables in meters.

thermal perturbations induced by the air flow. The thermal anomaly (defined as a shift in the AMT) and the annual temperature fluctuations both propagate over the same distances from the entrance. We call this distance the convection length and we show that at high Reynolds number, it is approximately proportional to the air flow rate divided by the square root of the cave radius. The order of magnitude of the theoretical convection length seem to compare satisfactorily with field data obtained from a mine and two caves. Convection lengths of a few tens to a few hundreds of meters should be commonly encountered in ventilated caves. In extreme cases, it could go up to kilometers.

## Data Availability Statement

Text files including the numerical results plotted in Figures 4–9 are available at https://doi.org/10.5281/zenodo.10953068 cited as Sedaghatkish and Doumenc (2024). The Comsol Multiphysics file containing the data and results related to the reference case of Figures 4–6 is also available, but a license of Comsol Multiphysics software is necessary to use it. All the field data (including the three instrumented sites) can be found in https://doi.org/10.5281/zenodo.10514555 associated with Sedaghatkish et al. (2024).

## References

Atkinson, T. C. (1977). Carbon dioxide in the atmosphere of the unsaturated zone: An important control of groundwater hardness in limestones. *Journal of Hydrology*, 35(1–2), 111–123. https://doi.org/10.1016/0022-1694(77)90080-4

Banner, J. L., Guilfoyle, A., James, E. W., Stern, L. A., & Musgrove, M. (2007). Seasonal variations in modern speleothem calcite growth in central Texas, USA. *Journal of Sedimentary Research*, 77(8), 615–622. https://doi.org/10.2110/jsr.2007.065

Bergman, T. L., Lavine, A. S., Incropera, F. P., & DeWitt, D. P. (2017). *Fundamentals of heat and mass transfer*. Wiley.

Borsato, A., Johnston, V. E., Frisia, S., Miorandi, R., & Corradini, F. (2016). Temperature and altitudinal influence on karst dripwater chemistry: Implications for regional-scale palaeoclimate reconstructions from speleothems. *Geochimica et Cosmochimica Acta*, 177, 275–297. https://doi.org/10.1016/j.gca.2015.11.043

Bourges, F., Genthon, P., Genty, D., Lorblanchet, M., Mauduit, E., & d'Hulst, D. (2014). Conservation of prehistoric caves and stability of their inner climate: Lessons from Chauvet and other French caves. *Science of the Total Environment*, 493, 79–91. https://doi.org/10.1016/j.scitotenv.2014.05.137

Bourret, F., Jeannin, P.-Y., Lalou, J.-C., Lambelet, J., Pauli, C., Spring, D., & Thévoz, P.-Y. (2007). Inventaire spéléologique de la suisse, Tome v, Nord vaudois, La Commission de Spéléologie de l'Académie suisse des sciences naturelles.

Casteel, R. C., & Banner, J. L. (2015). Temperature-driven seasonal calcite growth and drip water trace element variations in a well-ventilated Texas cave: Implications for speleothem paleoclimate studies. *Chemical Geology*, 392, 43–58. https://doi.org/10.1016/j.chemgeo.2014.11.002

Covington, M. D., Luhmann, A. J., Gabrovšek, F., Saar, M. O., & Wicks, C. M. (2011). Mechanisms of heat exchange between water and rock in karst conduits. *Water Resources Research*, 47(10), W10514. https://doi.org/10.1029/2011wr010683

Cropley, J. B. (1965). Influence of surface conditions on temperatures in large cave systems. *Bulletin of the National Speleological Society*, 27, 1–10.

Domínguez-Villar, D., Krklec, K., López-Sáez, J. A., & Sierro, F. J. (2021). Thermal impact of Heinrich stadials in cave temperature and speleothem oxygen isotope records. *Quaternary Research*, 101, 1–14. https://doi.org/10.1017/qua.2020.99

Domínguez-Villar, D., Krklec, K., & Sierro, F. J. (2023). Thermal conduction in karst terrains dominating cave atmosphere temperatures: Quantification of thermal diffusivity. *International Journal of Thermal Sciences*, 189, 108282. https://doi.org/10.1016/j.ijthermalsci.2023.108282

Dreybrodt, W., Gabrovšek, F., & Perne, M. (2005). Condensation corrosion: A theoretical approach. *Acta Carsologica*, 34(2). https://doi.org/10.3986/ac.v34i2.262

## Acknowledgments

This research was funded in whole or in part by the Swiss National Science Foundation (SNSF) [200021_188636]. The data from Bärenhöhle were acquired in the frame of an APART-fellowship from the Austrian Academy of Science granted to ML. Open access funding provided by Université de Neuchâtel.
<table><caption>**Table 7**

Summary Table of the Three Investigated Sites</caption><thead><tr><th></th><th></th><th>Baulmes mine</th><th>D7.1 cave</th><th>Bärenhöhle Cave</th></tr></thead><tbody><tr><td>R</td><td></td><td>1.5 m</td><td>0.5 m</td><td>0.5 m</td></tr><tr><td>R̃</td><td></td><td>0.3</td><td>0.1</td><td>0.1</td></tr><tr><td>Δṁ</td><td></td><td>4 kg/s</td><td>0.1 kg/s</td><td>0.6 kg/s</td></tr><tr><td>Re</td><td></td><td>10<sup>5</sup></td><td>7 × 10<sup>3</sup></td><td>5 × 10<sup>4</sup></td></tr><tr><td>L̃<sub>0</sub> ≈ L̃<sub>1</sub></td><td>Field data</td><td>>68</td><td>[2.1, 6.2]</td><td>>23</td></tr><tr><td></td><td>Theory</td><td>~3 × 10<sup>2</sup></td><td>>10</td><td>~70</td></tr><tr><td>L<sub>0</sub> ≈ L<sub>1</sub></td><td>Field data</td><td>>360 m</td><td>[10 m, 30 m]</td><td>>120 m</td></tr><tr><td></td><td>Theory</td><td>~2 × 10<sup>3</sup> m</td><td>>50 m</td><td>~4 × 10<sup>2</sup> m</td></tr></tbody></table>

Note. *R*: Cave radius in meters; *R̃*: dimensionless cave radius; *Δṁ*: amplitude of the air flow rate in kg/s; *Re*: Reynolds number; *L̃*<sub>0</sub> and *L̃*<sub>1</sub>: dimensionless convection lengths deduced from the AMT and the ATF, respectively; *L*<sub>0</sub> and *L*<sub>1</sub>: same variables in meters.

thermal perturbations induced by the air flow. The thermal anomaly (defined as a shift in the AMT) and the annual temperature fluctuations both propagate over the same distances from the entrance. We call this distance the convection length and we show that at high Reynolds number, it is approximately proportional to the air flow rate divided by the square root of the cave radius. The order of magnitude of the theoretical convection length seem to compare satisfactorily with field data obtained from a mine and two caves. Convection lengths of a few tens to a few hundreds of meters should be commonly encountered in ventilated caves. In extreme cases, it could go up to kilometers.

## Data Availability Statement

Text files including the numerical results plotted in Figures 4–9 are available at https://doi.org/10.5281/zenodo.10953068 cited as Sedaghatkish and Doumenc (2024). The Comsol Multiphysics file containing the data and results related to the reference case of Figures 4–6 is also available, but a license of Comsol Multiphysics software is necessary to use it. All the field data (including the three instrumented sites) can be found in https://doi.org/10.5281/zenodo.10514555 associated with Sedaghatkish et al. (2024).

## References

Atkinson, T. C. (1977). Carbon dioxide in the atmosphere of the unsaturated zone: An important control of groundwater hardness in limestones. *Journal of Hydrology*, 35(1–2), 111–123. https://doi.org/10.1016/0022-1694(77)90080-4

Banner, J. L., Guilfoyle, A., James, E. W., Stern, L. A., & Musgrove, M. (2007). Seasonal variations in modern speleothem calcite growth in central Texas, USA. *Journal of Sedimentary Research*, 77(8), 615–622. https://doi.org/10.2110/jsr.2007.065

Bergman, T. L., Lavine, A. S., Incropera, F. P., & DeWitt, D. P. (2017). *Fundamentals of heat and mass transfer*. Wiley.

Borsato, A., Johnston, V. E., Frisia, S., Miorandi, R., & Corradini, F. (2016). Temperature and altitudinal influence on karst dripwater chemistry: Implications for regional-scale palaeoclimate reconstructions from speleothems. *Geochimica et Cosmochimica Acta*, 177, 275–297. https://doi.org/10.1016/j.gca.2015.11.043

Bourges, F., Genthon, P., Genty, D., Lorblanchet, M., Mauduit, E., & d'Hulst, D. (2014). Conservation of prehistoric caves and stability of their inner climate: Lessons from Chauvet and other French caves. *Science of the Total Environment*, 493, 79–91. https://doi.org/10.1016/j.scitotenv.2014.05.137

Bourret, F., Jeannin, P.-Y., Lalou, J.-C., Lambelet, J., Pauli, C., Spring, D., & Thévoz, P.-Y. (2007). Inventaire spéléologique de la suisse, Tome v, Nord vaudois, La Commission de Spéléologie de l'Académie suisse des sciences naturelles.

Casteel, R. C., & Banner, J. L. (2015). Temperature-driven seasonal calcite growth and drip water trace element variations in a well-ventilated Texas cave: Implications for speleothem paleoclimate studies. *Chemical Geology*, 392, 43–58. https://doi.org/10.1016/j.chemgeo.2014.11.002

Covington, M. D., Luhmann, A. J., Gabrovšek, F., Saar, M. O., & Wicks, C. M. (2011). Mechanisms of heat exchange between water and rock in karst conduits. *Water Resources Research*, 47(10), W10514. https://doi.org/10.1029/2011wr010683

Cropley, J. B. (1965). Influence of surface conditions on temperatures in large cave systems. *Bulletin of the National Speleological Society*, 27, 1–10.

Domínguez-Villar, D., Krklec, K., López-Sáez, J. A., & Sierro, F. J. (2021). Thermal impact of Heinrich stadials in cave temperature and speleothem oxygen isotope records. *Quaternary Research*, 101, 1–14. https://doi.org/10.1017/qua.2020.99

Domínguez-Villar, D., Krklec, K., & Sierro, F. J. (2023). Thermal conduction in karst terrains dominating cave atmosphere temperatures: Quantification of thermal diffusivity. *International Journal of Thermal Sciences*, 189, 108282. https://doi.org/10.1016/j.ijthermalsci.2023.108282

Dreybrodt, W., Gabrovšek, F., & Perne, M. (2005). Condensation corrosion: A theoretical approach. *Acta Carsologica*, 34(2). https://doi.org/10.3986/ac.v34i2.262

## Acknowledgments

This research was funded in whole or in part by the Swiss National Science Foundation (SNSF) [200021_188636]. The data from Bärenhöhle were acquired in the frame of an APART-fellowship from the Austrian Academy of Science granted to ML. Open access funding provided by Université de Neuchâtel.
Fairchild, I. J., & Baker, A. (2012). *Speleothem science: From process to past environments*. John Wiley & Sons.

Gabrovšek, F. (2023). How do caves breathe: The airflow patterns in karst underground. *PLoS One*, 18(4), e0283767. https://doi.org/10.1371/journal.pone.0283767

Gomell, A., Austin, D., Ohms, M., & Pflitsch, A. (2021). Air pressure propagation through Wind Cave and Jewel Cave: How do pressure waves travel through barometric caves? *International Journal of Speleology*, 50(3), 263–273. https://doi.org/10.5038/1827-806x.50.3.2393

Gong, X., Yang, X., Luo, Q., & Tang, L. (2019). Effects of convective heat transport in modelling the early evolution of conduits in limestone aquifers. *Geothermics*, 77, 383–394. https://doi.org/10.1016/j.geothermics.2018.10.010

Guerrier, B., Doumenc, F., Roux, A., Mergui, S., & Jeannin, P.-Y. (2019). Climatology in shallow caves with negligible ventilation: Heat and mass transfer. *International Journal of Thermal Sciences*, 146, 106066. https://doi.org/10.1016/j.ijthermalsci.2019.106066

Haaland, S. E. (1983). Simple and explicit formulas for the friction factor in turbulent pipe flow.

Kim, S.-T., & O'Neil, J. R. (1997). Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. *Geochimica et Cosmochimica Acta*, 61(16), 3461–3475. https://doi.org/10.1016/s0016-7037(97)00169-5

Kukuljan, L., Gabrovšek, F., & Covington, M. (2021). The relative importance of wind-driven and chimney effect cave ventilation: Observations in Postojna Cave (Slovenia). *International Journal of Speleology*, 50(3), 275–288. https://doi.org/10.5038/1827-806x.50.3.2392

Lismonde, B. (2002). Climatologie du monde souterrain; Tome 2-Aérologie des systèmes karstiques, Persée-Portail des revues scientifiques en SHS.

Long, A. J., & Gilcrease, P. C. (2009). A one-dimensional heat-transport model for conduit flow in karst aquifers. *Journal of Hydrology*, 378(3-4), 230–239. https://doi.org/10.1016/j.jhydrol.2009.09.024

Luetscher, M., & Jeannin, P. Y. (2004). Temperature distribution in karst systems: The role of air and water fluxes. *Terra Nova*, 16(6), 344–350. https://doi.org/10.1111/j.1365-3121.2004.00572.x

Luetscher, M., Lismonde, B., & Jeannin, P. Y. (2008). Heat exchanges in the heterothermic zone of a karst system: Monlesi cave, Swiss Jura Mountains. *Journal of Geophysical Research*, 113(F2), F02025. https://doi.org/10.1029/2007jf000892

Mammola, S., Piano, E., Cardoso, P., Vernon, P., Domínguez-Villar, D., Culver, D. C., et al. (2019). Climate change going deep: The effects of global climatic alterations on cave ecosystems. *Anthropological Review*, 6(1–2), 98–116. https://doi.org/10.1177/2053019619851594

Obleitner, F., & Spötl, C. (2011). The mass and energy balance of ice within the Eisriesenwelt cave, Austria. *The Cryosphere*, 5(1), 245–257. https://doi.org/10.5194/tc-5-245-2011

Pastore, C., Pastore, M., Marraffa, M., Marraffa, A., Romano, R., Crespo, S., et al. (2019). Albania—Tra le (R)Rughe albanesi. *Racconti di esplorazioni condivise*. *Speleologia*, 80, 38–45.

Perkmann, L., & Luetscher, M. (2013). Vermessung und hydrogeologische Untersuchung des Bärenhöhlensystems, Bregenzerwald, Vorarlberg. 64, 32–44.

Pethkool, S., Eiamsa-Ard, S., Kwankaomeng, S., & Promvonge, P. (2011). Turbulent heat transfer enhancement in a heat exchanger using helically corrugated tube. *International Communications in Heat and Mass Transfer*, 38(3), 340–347. https://doi.org/10.1016/j.icheatmasstransfer.2010.11.014

Promthaisong, P., Jedsadaratanachai, W., & Eiamsa-Ard, S. (2016). 3D numerical study on the flow topology and heat transfer characteristics of turbulent forced convection in spirally corrugated tube. *Numerical Heat Transfer Part A*, 69(6), 607–629. https://doi.org/10.1080/10407782.2015.1069670

Qaddah, B., Soucasse, L., Doumenc, F., Mergui, S., Rivière, P., & Soufiani, A. (2023). Coupled heat and mass transfer in shallow caves: Interactions between turbulent convection, gas radiative transfer and moisture transport. *International Journal of Thermal Sciences*, 194, 108556. https://doi.org/10.1016/j.ijthermalsci.2023.108556

Quindos, L. S., Bonet, A., Diaz-Caneja, N., Fernandez, P. L., Gutierrez, I., Solana, J. R., et al. (1987). Study of the environmental variables affecting the natural preservation of the Altamira Cave paintings located at Santillana del Mar, Spain. *Atmospheric Environment*, 21(3), 551–560. https://doi.org/10.1016/0004-6981(87)90037-0

Salmon, F., Lacanette, D., Lharti, H., & Sirieux, C. (2023). Heat transfer in rock masses: Application to the Lascaux Cave (France). *International Journal of Heat and Mass Transfer*, 207, 124029. https://doi.org/10.1016/j.ijheatmasstransfer.2023.124029

Sedaghatkish, A., & Doumenc, F. (2024). Software file and numerical results of modelling heat transfer for assessing the convection length in ventilated caves [Software]. *Zenodo*. https://doi.org/10.5281/zenodo.10953068

Sedaghatkish, A., Pastore, C., Doumenc, F., Jeannin, P.-Y., & Luetscher, M. (2024). Field data for modelling heat transfer for assessing the convection length in ventilated caves [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.10514555

Sinokrot, B. A., & Stefan, H. G. (1993). Stream temperature dynamics: Measurements and modeling. *Water Resources Research*, 29(7), 2299–2312. https://doi.org/10.1029/93wr00540

Spötl, C., Fairchild, I. J., & Tooth, A. F. (2005). Cave air control on dripwater geochemistry, Obir Caves (Austria): Implications for speleothem deposition in dynamically ventilated caves. *Geochimica et Cosmochimica Acta*, 69(10), 2451–2468. https://doi.org/10.1016/j.gca.2004.12.009

Tremaine, D. M., Froelich, P. N., & Wang, Y. (2011). Speleothem calcite farmed in situ: Modern calibration of δ¹⁸O and δ¹³C paleoclimate proxies in a continuously-monitored natural cave system. *Geochimica et Cosmochimica Acta*, 75(17), 4929–4950. https://doi.org/10.1016/j.gca.2011.06.005

Villar, E., Fernández, P. L., Quindos, L. S., Solana, J. R., & Soto, J. (1983). Temperature of rock surfaces in Altamira Cave (Spain). *Transactions of the British Cave Research Association*, 10, 165–170.

Wigley, T. M. L., & Brown, C. (1971). Geophysical applications of heat and mass transfer in turbulent pipe flow. *Boundary-Layer Meteorology*, 1(3), 300–320. https://doi.org/10.1007/bf02186034
Fairchild, I. J., & Baker, A. (2012). *Speleothem science: From process to past environments*. John Wiley & Sons.

Gabrovšek, F. (2023). How do caves breathe: The airflow patterns in karst underground. *PLoS One*, 18(4), e0283767. https://doi.org/10.1371/journal.pone.0283767

Gomell, A., Austin, D., Ohms, M., & Pflitsch, A. (2021). Air pressure propagation through Wind Cave and Jewel Cave: How do pressure waves travel through barometric caves? *International Journal of Speleology*, 50(3), 263–273. https://doi.org/10.5038/1827-806x.50.3.2393

Gong, X., Yang, X., Luo, Q., & Tang, L. (2019). Effects of convective heat transport in modelling the early evolution of conduits in limestone aquifers. *Geothermics*, 77, 383–394. https://doi.org/10.1016/j.geothermics.2018.10.010

Guerrier, B., Doumenc, F., Roux, A., Mergui, S., & Jeannin, P.-Y. (2019). Climatology in shallow caves with negligible ventilation: Heat and mass transfer. *International Journal of Thermal Sciences*, 146, 106066. https://doi.org/10.1016/j.ijthermalsci.2019.106066

Haaland, S. E. (1983). Simple and explicit formulas for the friction factor in turbulent pipe flow.

Kim, S.-T., & O'Neil, J. R. (1997). Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. *Geochimica et Cosmochimica Acta*, 61(16), 3461–3475. https://doi.org/10.1016/s0016-7037(97)00169-5

Kukuljan, L., Gabrovšek, F., & Covington, M. (2021). The relative importance of wind-driven and chimney effect cave ventilation: Observations in Postojna Cave (Slovenia). *International Journal of Speleology*, 50(3), 275–288. https://doi.org/10.5038/1827-806x.50.3.2392

Lismonde, B. (2002). Climatologie du monde souterrain; Tome 2-Aérologie des systèmes karstiques, Persée-Portail des revues scientifiques en SHS.

Long, A. J., & Gilcrease, P. C. (2009). A one-dimensional heat-transport model for conduit flow in karst aquifers. *Journal of Hydrology*, 378(3-4), 230–239. https://doi.org/10.1016/j.jhydrol.2009.09.024

Luetscher, M., & Jeannin, P. Y. (2004). Temperature distribution in karst systems: The role of air and water fluxes. *Terra Nova*, 16(6), 344–350. https://doi.org/10.1111/j.1365-3121.2004.00572.x

Luetscher, M., Lismonde, B., & Jeannin, P. Y. (2008). Heat exchanges in the heterothermic zone of a karst system: Monlesi cave, Swiss Jura Mountains. *Journal of Geophysical Research*, 113(F2), F02025. https://doi.org/10.1029/2007jf000892

Mammola, S., Piano, E., Cardoso, P., Vernon, P., Domínguez-Villar, D., Culver, D. C., et al. (2019). Climate change going deep: The effects of global climatic alterations on cave ecosystems. *Anthropological Review*, 6(1–2), 98–116. https://doi.org/10.1177/2053019619851594

Obleitner, F., & Spötl, C. (2011). The mass and energy balance of ice within the Eisriesenwelt cave, Austria. *The Cryosphere*, 5(1), 245–257. https://doi.org/10.5194/tc-5-245-2011

Pastore, C., Pastore, M., Marraffa, M., Marraffa, A., Romano, R., Crespo, S., et al. (2019). Albania—Tra le (R)Rughe albanesi. *Racconti di esplorazioni condivise*. *Speleologia*, 80, 38–45.

Perkmann, L., & Luetscher, M. (2013). Vermessung und hydrogeologische Untersuchung des Bärenhöhlensystems, Bregenzerwald, Vorarlberg. *64*, 32–44.

Pethkool, S., Eiamsa-Ard, S., Kwankaomeng, S., & Promvonge, P. (2011). Turbulent heat transfer enhancement in a heat exchanger using helically corrugated tube. *International Communications in Heat and Mass Transfer*, 38(3), 340–347. https://doi.org/10.1016/j.icheatmasstransfer.2010.11.014

Promthaisong, P., Jedsadaratanachai, W., & Eiamsa-Ard, S. (2016). 3D numerical study on the flow topology and heat transfer characteristics of turbulent forced convection in spirally corrugated tube. *Numerical Heat Transfer Part A*, 69(6), 607–629. https://doi.org/10.1080/10407782.2015.1069670

Qaddah, B., Soucasse, L., Doumenc, F., Mergui, S., Rivière, P., & Soufiani, A. (2023). Coupled heat and mass transfer in shallow caves: Interactions between turbulent convection, gas radiative transfer and moisture transport. *International Journal of Thermal Sciences*, 194, 108556. https://doi.org/10.1016/j.ijthermalsci.2023.108556

Quindos, L. S., Bonet, A., Diaz-Caneja, N., Fernandez, P. L., Gutierrez, I., Solana, J. R., et al. (1987). Study of the environmental variables affecting the natural preservation of the Altamira Cave paintings located at Santillana del Mar, Spain. *Atmospheric Environment*, 21(3), 551–560. https://doi.org/10.1016/0004-6981(87)90037-0

Salmon, F., Lacanette, D., Lharti, H., & Sirieux, C. (2023). Heat transfer in rock masses: Application to the Lascaux Cave (France). *International Journal of Heat and Mass Transfer*, 207, 124029. https://doi.org/10.1016/j.ijheatmasstransfer.2023.124029

Sedaghatkish, A., & Doumenc, F. (2024). Software file and numerical results of modelling heat transfer for assessing the convection length in ventilated caves [Software]. *Zenodo*. https://doi.org/10.5281/zenodo.10953068

Sedaghatkish, A., Pastore, C., Doumenc, F., Jeannin, P.-Y., & Luetscher, M. (2024). Field data for modelling heat transfer for assessing the convection length in ventilated caves [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.10514555

Sinokrot, B. A., & Stefan, H. G. (1993). Stream temperature dynamics: Measurements and modeling. *Water Resources Research*, 29(7), 2299–2312. https://doi.org/10.1029/93wr00540

Spötl, C., Fairchild, I. J., & Tooth, A. F. (2005). Cave air control on dripwater geochemistry, Obir Caves (Austria): Implications for speleothem deposition in dynamically ventilated caves. *Geochimica et Cosmochimica Acta*, 69(10), 2451–2468. https://doi.org/10.1016/j.gca.2004.12.009

Tremaine, D. M., Froelich, P. N., & Wang, Y. (2011). Speleothem calcite farmed in situ: Modern calibration of δ¹⁸O and δ¹³C paleoclimate proxies in a continuously-monitored natural cave system. *Geochimica et Cosmochimica Acta*, 75(17), 4929–4950. https://doi.org/10.1016/j.gca.2011.06.005

Villar, E., Fernández, P. L., Quindos, L. S., Solana, J. R., & Soto, J. (1983). Temperature of rock surfaces in Altamira Cave (Spain). *Transactions of the British Cave Research Association*, 10, 165–170.

Wigley, T. M. L., & Brown, C. (1971). Geophysical applications of heat and mass transfer in turbulent pipe flow. *Boundary-Layer Meteorology*, 1(3), 300–320. https://doi.org/10.1007/bf02186034