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

HAL
open science

# Thermal modeling of caves ventilated by chimney effect

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. Thermal modeling of caves ventilated by chimney effect. International Journal of Thermal Sciences, 2025, 212, pp.109757. 10.1016/j.ijthermalsci.2025.109757. hal-04934720

HAL Id: hal-04934720

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

Submitted on 7 Feb 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermal modeling of caves ventilated by chimney effect

Amir Sedaghatkish<sup>a,b,*</sup>, Claudio Pastore<sup>a,b,●</sup>, Frédéric Doumenc<sup>c,d</sup>, Pierre-Yves Jeannin<sup>a,b</sup>,
Marc Luetscher<sup>a,●</sup>

## ABSTRACT

Cave ventilation significantly increases the depth of natural thermal oscillations and decreases the time of propagation compared to heat conduction in the rock mass. This makes it necessary to develop and test thermal models for the prediction of temperature fields in ventilated karst massifs. Here, we develop a thermal model of a single conduit ventilated by chimney effect. The model is based on the diffusion equation in the rock mass coupled to the conservation of energy and water vapor mass in the airflow. The effect of the latent heat of evaporation and condensation is considered. In parallel, the main conduit of a ventilated cave has been equipped with a flowmeter and several temperature sensors. The model is tested against field data collected during a complete year. The relevance of the model assumptions (geometry simplification, initial and boundary conditions, use of transfer coefficients to couple the air and the conduit wall) is thoroughly analyzed. The model correctly predicts the temperature fluctuations at daily and yearly scale, but underestimates the annual mean temperatures inside the cave. A biased assessment of the ground temperature seems to explain this discrepancy. The effect of condensation and evaporation on the cave climate turns out to be low on cave temperature, but significant on air humidity with consequences for ecology or paleoclimatology. This study is a first step towards the elaboration and validation of models providing a quantitative assessment of caves' thermal response at any location and time scale.

## 1. Introduction

Extensive cave networks in karst massifs are traversed by airflows that considerably enhance heat and mass transfer between the bedrock and the external environment [1]. The impact of airflow must be considered in a wide range of areas including paleoclimatology and ecological issues. Indeed, isotope fractionation that occurs during calcite precipitation is temperature dependent. Airflow also controls evaporation and condensation of water, the transport of carbon dioxide, and their consequences on precipitation and dissolution rates of calcite [2, 3]. A good understanding of heat and mass transfer associated with airflow is thus requested for paleoclimate reconstructions [4,5]. It is also relevant for the study of subterranean fauna and flora [6]. Indeed, caves are home to many temperature sensitive species, including bats, beetles, or thysanurans [7]. These subterranean ecosystems and their biodiversity are particularly vulnerable to climate change [8]. For instance, Rizzo et al. [9] observed that the survival rate of Pyrenean beetles over 7 days was 100 % at 6–20 °C, and collapsed to 0 % at 23–25 °C. This example illustrates the need for reliable thermal model for predicting the consequences of global warming on the temperature field in karst massifs.

The fluctuations of the external pressure can play a role in cave ventilation [10]. Large caves with a single narrow entrance are ventilated by barometric effect [11], and wind pressure may be the dominant driving force in caves with multiple entrances at similar heights [12]. In caves with multiple entrances at different heights, the most significant cause of ventilation is the chimney effect resulting from the density contrast between the air inside and outside the cave [1]. In temperate climates, this density contrast is mainly due to temperature [13]. Because of the high thermal inertia of the massif, the airflow direction is upward during most of the winter, when the temperature inside the cave is warmer than the atmospheric temperature. Conversely, the air flows downward during most of the summer. A noticeable consequence of this seasonal inversion of the airflow is the thermal anomalies observed in the entrance regions [1]. In lower entrance areas, the annual mean temperature is shifted to a lower value than outside because this part of the cave receives cold air at the atmospheric temperature in winter and air cooled down after circulating in the massif in summer. The same mechanism produces a warm anomaly in upper entrance areas.

Airflow modifies the rock temperature in the cave, which in turn changes the air temperature profile inside the cave. Therefore, both the air and the rock mass must be considered in a model aiming at predicting

<footnote><sup>a</sup> Swiss Institute for Speleology and Karst Studies (SISKA), CH-2300, La Chaux-de-Fonds, Switzerland</footnote>

<footnote><sup>b</sup> Center for Hydrogeology and Geothermics (CHYN), University of Neuchâtel, 2000, Neuchâtel, Switzerland</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, CNRS, FAST, 91405, Orsay, France</footnote>

<footnote><sup>d</sup> Sorbonne Université, UFR 919, 4 place Jussieu, F-75252, Paris Cedex 05, France</footnote>

<footnote>* Corresponding author. Swiss Institute for Speleology and Karst Studies (SISKA), CH-2300, La Chaux-de-Fonds, Switzerland. E-mail address: amir.sedaghatkish@isska.ch (A. Sedaghatkish).</footnote>
the climate of a ventilated cave. Lismonde [1] and Gabrovšek [13] developed such models based on the 1D energy balance in the airflow and the heat conduction equation in the rock mass. Both equations were coupled by the Newton's law of cooling, which states that the local heat flux through the cave wall is proportional to the difference between the wall temperature and the mixing temperature of the air. The coefficient of proportionality, the so-called heat transfer coefficient, depends on the conduit geometry, the fluid properties, the wall roughness and the fluid velocity. It was assessed using standard correlations for fully developed forced convection in pipes (see for instance Ref. [14]). This approach allowed Lismonde [1] and Gabrovšek [13] to reproduce qualitatively in numerical simulations some field observations as for instance the seasonal hysteresis of the flow rate or the thermal anomaly in the entrance areas.

Recently, we used a similar model to investigate the propagation of thermal perturbations by the airflow inside a karst massif in a quantitative way [15]. We considered the simple geometry of a straight horizontal conduit of constant circular cross-section. Numerical simulations showed that turbulent airflow, very common in ventilated caves, significantly modifies the rock temperature. Numerical simulations supported the partition of the cave into three regions as first proposed by Cropley [16]: (a) a short diffusive region, where heat mainly propagates from the 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. The length of the diffusive region is of the order of a few meters, whatever the values of the air flowrate and conduit diameter. The length of the convective region, defined as the convection length, corresponds to both the extent of the thermal anomaly and the distance of propagation of the annual temperature fluctuations. Numerical simulations showed that the convection length can commonly reach a few tens to a few hundreds of meters. It is approximately proportional to the amplitude of the flowrate annual fluctuations divided by the square root of the cave radius [15]. The orders of magnitude predicted by the model were consistent with field data from a mine tunnel and two caves.

In the present paper, we want to go beyond orders of magnitude by carrying out detailed and systematic comparisons between numerical simulations and field data. Indeed, numerical simulation of ventilated cave climate is still a challenge. To our knowledge, direct and quantitative comparisons between models and field data are still lacking. The present article provides such a comparison, with the aim of testing the validity of the physical assumptions included in the model. We focused on the following issues.

* A cave is a succession of bends, conduit contractions or enlargements, and obstacles of any kind. The first of the problems posed by this complex geometry is the level of accuracy required for its description. More specifically, is it necessary to consider a full 3D problem, or can we define an equivalent 2D geometry to save computational resources?

* The irregular geometry of caves also has important consequences on the assessment of the heat transfer coefficient. Correlations commonly used in thermal models are valid for straight rough pipes of uniform cross-section [14]. Geometrical singularities are expected to increase the turbulence level, which in turn can significantly enhance the heat transfer coefficient. Some corrections exist to take into account the effect of the well-defined singularities commonly found in heat exchangers [17]. However, they can hardly be used for the tortuous geometries typically encountered in caves.

* Numerical simulations are based on mathematical models that require initial and boundary conditions. Ideally, the temperature field at any point of the massif should be known at the initial time, which is never the case. Another issue is the boundary condition at the external surface (i.e., the interface between the massif and the atmosphere). The simplest choice is to assume that the external surface temperature is equal to the atmospheric temperature [13,

[15]. However, the earth surface does not only interact with the surrounding air. Its temperature results from an energy balance including convection with the air, solar irradiation (short wavelength), radiative exchange with the surrounding environment and the sky (long wavelength), latent heat due to evaporation or condensation, and conduction inside the rock [18]. This raises the question of the most relevant boundary condition at the external surface.

* Condensation and evaporation of water are commonly observed on the walls of ventilated caves [19–21]. Since the latent heat of evaporation or condensation induces cooling or warming, its contribution to cave climate must be evaluated.

With the aim of addressing these issues, flowrate and temperatures in the ventilated cave of Longeaigne (Switzerland) have been recorded during a complete year. In parallel, we developed a thermal model of the cave, for which results have been tested against the field data collected in the cave.

## 2. Field data

### 2.1. Study site

The "Baume de Longeaigne" displayed in Fig. 1 is located in the Swiss Jura Mountains (Val-de-Travers; Upper entrance: 46°52'22" N 6°31'09" E, 917 m a.s.l.; Lower entrance: 46°52'18" N 6°31'09" E, 820 m a.s.l.). The difference of elevation between both entrances is 97 m. This induces an intense buoyancy-driven airflow inside the conduit that connects both entrances. A secondary conduit (blue arrow in Fig. 1), likely connecting the cave to the atmosphere through a network of fractures, intersects the main conduit at approximately 16 m from the upper entrance. The possibility of a leak through this secondary conduit cannot be excluded.

The total length of the main conduit of Longeaigne cave is 311 m. Its shape and size are highly variable. The cross-sectional area A and perimeter P are displayed in Fig. 2 as a function of x, the distance from the upper entrance following the conduit axis. Both values were obtained using a laser-scanner with a $10^{-2}$ m resolution every 10 m along the cave. Fig. 3 displays the angle $\theta(x)$ between the upward vertical direction and the velocity vector of the air when flowing from the upper to the lower entrance.

The lower entrance is a temporary spring acting as an overflow of the "Raies" hydrogeological system, whose catchment area is approximately 60–80 km². The maximum water discharge-rate at Longeaigne entrance does probably reach approximately 5 m³/s [22]. A perennial lake has formed in the lower part of the main conduit (see Fig. 1). During flood events, the rise of the water level is sufficient to close the conduit, interrupting the airflow. This is an important feature of Longeaigne Cave that must be considered in the thermal model.

Since August 2020, several stations monitor the cave temperature as well as the airflow along the main conduit of the cave (red symbols in Fig. 1).

### 2.2. Temperature measurements

The temperature data logging stations indicated by H in Fig. 1 were equipped with Hobo Water Temperature Pro v2 (sensor type U22-01, accuracy ±0.21 °C, resolution 0.02 °C). Reefnet Sensus Ultra temperature (accuracy ±0.3 °C, resolution 0.01 °C) were deployed at the stations indicated by R. Two additional Reefnet sensors were used to monitor the atmospheric temperature close to the upper and lower entrances. Another temperature sensor was installed in Anemo station close to the upper entrance (sensor type MS8607 (Adafruit), accuracy ±1 °C, resolution 0.01 °C). The distance of each station from the upper (and lower) entrance are listed in Table 1. To make the conduit open enough for the commuting of cavers, all the sensors were located close to
<figure><img src="image_3.png" /><figcaption>Fig. 1. Longeaglique cave topography (developed vertical profile; main conduit in black, sensor locations in red). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Variation of the cross-sectional area A and perimeter P as a function of the distance x from the upper entrance.</figcaption></figure>

the wall, at a distance ranging from 2 to 4 cm.

With the exception of H1, the temperatures used in this paper were recorded from August 13, 2020 to August 13, 2021, which is exactly one year. The Hobo water probe installed in H1 at a few centimeters from the wall was complemented with a Pt100 sensor (TMC1-HD from Hobo Company, resolution: 0.02 °C, accuracy: 0.1) suspended approximately at the center of the cross-section. The temperatures of these two probes were recorded for approximately one month, from October 11, 2022 to November 06, 2022 (sampling rate for all temperature measurements: one measuring point per hour).

External upper and lower temperatures are displayed in Fig. 4a. The dashed lines represent the function:

$$T_{\text{atm}}(t) = T_m - \Delta T \sin\left[\frac{2\pi}{\tau}(t - t_0)\right] \quad (1)$$

where τ = 1 year, t₀ is a parameter that accounts for the origin of time, Tm is the annual mean temperature (AMT) and ΔT the amplitude of the annual temperature fluctuations (ATF). t₀, Tm and ΔT are deduced from the Fast Fourier Transform (FFT) of the atmospheric temperatures displayed in Fig. 4b, which also provides the amplitude of all other modes,
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Angle $\theta$ between the gravity vector and the velocity vector of the air when flowing from the upper to lower entrance as a function of the distance $x$ from the upper entrance ($\theta$ is larger than 90 °C in the conduit segment between P and R4 stations, see Fig. 1).</figcaption></figure>

<table><caption>**Table 1**
Location of measuring stations.</caption><thead><tr><td>Station</td><td>Distance from the upper (lower) entrance</td></tr></thead><tbody><tr><td>Up_ext</td><td>Outside, near the upper entrance</td></tr><tr><td>Anemo</td><td>4 m (307 m)</td></tr><tr><td>H1</td><td>16 m (295 m)</td></tr><tr><td>H3</td><td>55 m (256 m)</td></tr><tr><td>H4</td><td>97 m (214 m)</td></tr><tr><td>R2</td><td>140 m (171 m)</td></tr><tr><td>P</td><td>216 m (95 m)</td></tr><tr><td>H6</td><td>269 m (42 m)</td></tr><tr><td>R4</td><td>294 m (17 m)</td></tr><tr><td>Low_ext</td><td>Outside, near the lower entrance</td></tr></tbody></table>

including that of the daily temperature fluctuations (DTF). As expected, the modes showing the highest amplitudes are the ATF and the DTF. The ATF amplitudes at the upper and lower entrances (respectively 8.34 and 8.28 °C) are very close to each other. The amplitudes of the DTF differ a little more (3.4 °C and 2.4 °C at the upper and lower entrances, respectively). The AMT is higher at the upper entrance (7.38 °C) than at the lower entrance (6.75 °C). The mean vertical thermal gradient between both entrances is thus 6.5 °C/km.

## 2.3. Airflow measurements

The Sensirion SFM3003-300-CE (from now on SFM) is a digital bidirectional mass flowmeter installed at Anemo station, at 4 m from the upper entrance (Fig. 1). The device is compact and robust. It shows a high data resolution (up to 16 bit) and low power requirement (typically 3.3 V and 3.8 mA). Furthermore, the bidirectional cave airflow is well recognized by the device which records the flow direction as negative or positive values. The device measures the flowrate that passes through a tube with an inner diameter of 20 mm. This flowrate is then converted in local air velocity through appropriate calibration with an accuracy of ±5 % for speeds higher than 0.30 m/s and ±26 % for speeds in the range from 0.15 to 0.30 m/s [23]. The flowrate through the conduit is obtained by multiplying this local velocity by the cross-sectional area. The main cause of uncertainty is the non-uniform velocity field over the cross-section. Complementary measurements performed by other means (including manual measurements with an anemometer, CO₂ gauging, comparison with other SFM) suggest that: (a) the flowrate measured at Anemo station is certainly underestimated, (b) the ratio of the real flowrate over the measured flowrate cannot be larger than 1.5.

The air mass flowrate was recorded simultaneously with the temperatures using the same sampling rate (see Fig. 5a). It varies approximately in the range from -1.5 kg s⁻¹ in winter to 1 kg s⁻¹ in summer

(positive values correspond to airflow direction from the upper to the lower entrance). The sporadic interruptions of airflow when the lake completely fills the lower part of the conduit are clearly observed in Fig. 5a. A pressure sensor has been installed in point P (see Fig. 1 and Table 1) to measure the lake level. It provides the time ranges during which the main conduit is blocked by the lake (called “cave closures” in the following). Fig. 5a shows that the flowrate measured at Anemo station during these cave closures is not always zero. This can be due to the existence of a leakage rate through the secondary conduit (see Section. 2.1 and Fig. 1), or to sporadic free convection cells developing from the upper entrance, inducing local velocities recorded by the SFM. Conversely, it can happen that no airflow is measured whereas the pressure sensor indicates that the cave is open (e.g., for time ranging from 196 to 213 days). This happens when the flowrate is too small to be detected by the SFM, or because of an uncertainty on the lake level corresponding to cave closures. During the monitoring year, the cave was closed for 203 days (56 % of the time). The airflow was directed from the upper to the lower entrance for 66 days and in the reverse direction for 95 days (respectively 18 % and 26 % of the time).

The FFT of the mass flowrate is displayed in Fig. 5b. In contrast with the temperature spectrum that was dominated by the annual and daily fluctuations, many other modes emerge with amplitudes comparable or even higher than the annual flowrate fluctuations (AFF) or daily flow-rate fluctuations (DFF). These modes with intermediate frequencies originate from the intermittent lake closures. Moreover, the annual mean flowrate (AMF) is strictly negative and its absolute value is equal to half the AFF (see Fig. 5b). This seasonal asymmetry is due to the highest flowrates reached in winter. The dashed lines in Fig. 5a represent the same sine function as in Eq. (1) after substituting temperatures for mass flowrates.

# 3. Numerical simulations

## 3.1. Simplifying assumptions

Only the main conduit connecting the upper to the lower entrance is considered in the numerical simulations (see Fig. 1). For the sake of simplification, it has been “unfolded” to yield a rectilinear conduit of the same length $L_{cave}$ as the real cave (see Fig. 6). A circular cross-section of variable diameter $D_p(x)$ is assumed. $D_p(x)$ is set to get the same perimeter $P(x)$ as the true conduit (i.e., $D_p(x) = P(x)/\pi$, where the perimeter $P(x)$ is known from the cave survey, see Fig. 2). The sensitivity study presented in Appendix A suggests that this choice, which preserves the exchange surface between the air and the rock, is a good approximation in most cases. The conduit is located in an impermeable rock domain of outer radius $R_{dom}$ where conduction is the only heat transfer process. The outer radius must be large enough to ensure that the corresponding domain boundary is not significantly influenced by the airflow in the conduit and is thus approximately adiabatic. These simplifications allow to consider a 2D axisymmetric conduction problem in the rock mass, resulting in a significant reduction in the computational resources compared to the real 3D geometry. Moreover, we assume that the temperature of the external surface is equal to the measured atmospheric temperature. This boundary condition offers the advantage of simplicity, but neglects the effect of radiative transfer and evaporation or condensation on the temperature of the external surface.

The need for an aeraulic model is avoided by using the air flowrate measured by the SFM. During cave closure periods, we assume that the air only circulates in the part of the conduit between the upper entrance and the junction with the secondary conduit (at 16 m from the upper entrance). We impose zero flowrate in the rest of the conduit (i.e., between the junction with the secondary conduit and the lower entrance). Although the conduit is assumed rectilinear, the variation of potential energy with the altitude will be considered in the air energy balance using the conduit tilt angle $\theta(x)$ displayed in Fig. 3.

The effect of the latent heat of evaporation or condensation inside
<figure><img src="image_6.png" /><figcaption>Fig. 4. Atmospheric temperature close to the upper and lower entrances. a) Time series; b) Amplitude as a function of frequency deduced from FFT (due to the logarithmic scale the AMT, which corresponds to zero frequency, is located at an arbitrary origin).</figcaption></figure>

the cave is an important issue that must be evaluated. A difficulty is that the cave walls might be covered with a water film or completely dry, depending on the season and the weather conditions. To avoid the complicated task of defining a model predicting the amount of water lying on the walls, we define two simple limiting cases representing a lower and an upper bound of the latent heat effects.

- Case 1: evaporation and condensation are disregarded (no latent heat effect).

- Case 2: a permanent thin film of water covers all the cave walls, so that evaporation cannot be interrupted by the lack of water on the walls. In addition, the water film is assumed thin enough to neglect its thermal resistance compared to that due to conduction in the rock or convection. Condensation in the gas phase is disregarded. We assume that it only takes place on the conduit wall (same assumptions as in Qaddah et al. [24]).

The advection of heat by water flowing down the conduit is disregarded. This assumption seems reasonable when the water flow reduces to a thin liquid film due to condensation or percolation through the rock porosity. In contrast, it is questionable during cave flooding

induced by intense water recharge. This point is discussed in Section 4.2.

The heat flux through the conduit wall, which couples the energy balance in the air to the conduction equation in the rock mass, is estimated using the Newton's law of cooling based on a heat transfer coefficient. When condensation and evaporation are considered, the mass balance of the water vapor transported by the humid air is also included in the model. This equation contains a mass transfer coefficient for the assessment of the evaporation/condensation fluxes at the cave wall. The Lewis analogy [25] states that the heat and mass transfer coefficients follow the same physical laws. To overcome the difficulty that arises from the complex geometry of a cave, we adopt the same method as in our previous work [15]. We perform two distinct set of simulations to assess a lower bound and a higher bound of the heat and vapor fluxes through the cave wall. In a first set, the transfer coefficients are estimated from standard correlations valid for forced convection and fully developed flow in pipes, which is expected to yield a lower bound. In a second set of simulations, an upper bound of these fluxes is obtained assuming infinite transfer coefficients. Combining dry or humid air along with finite or infinite transfer coefficients yields the four limiting cases listed in Table 2.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Air mass flowrate measured at Anemo station; a) Time series (dashed line corresponds to Eq. (1)) where $T_m$ and $\Delta T$ are replaced by the AMF $\dot{m}_m$ and the AFF $\Delta \dot{m}$, respectively). The cave closure periods are shown by blue solid line determining if the cave is open (one) or closed (zero). b) Amplitudes versus frequency deduced from FFT (due to the logarithmic scale, the AMF, which corresponds to zero frequency, is located at an arbitrary origin). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

## 3.2. Governing equations

We first detail the complete set of equations for the HA case (humid air and finite transfer coefficients), which is the most general. Then we give the simplifications for the simulation of dry air or infinite transfer coefficients.

### 3.2.1. Rock mass

The temperature field in the impermeable rock mass follows the heat conduction equation in the cylindrical coordinates system:

$$ \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} \quad (2) $$

$T_r(x, r)$ is the rock temperature at distances $x$ from the upper entrance and $r$ from the conduit axis. $\alpha_r$ is the thermal diffusivity of the rock.

Defining right boundary condition at the external ground surface is challenging. For the sake of simplicity, the atmospheric temperature is imposed at the external surfaces located at $x = 0$ and $x = L_{cave}$:

$$ T_r(0, r, t) = T_{atmU}(t) \text{ for } R_p(0) \le r \le R_{dom}, \quad (3) $$

$$ T_r(L_{cave}, r, t) = T_{atmL}(t) \quad \text{for } R_p(L_{cave}) \le r \le R_{dom}, \quad (4) $$

where $R_p(x) = D_p(x)/2$ is the cave radius that yields the perimeter obtained from the cave survey (Fig. 2). $T_{atmL}$ and $T_{atmU}$ are the atmospheric temperatures measured in the vicinity of the lower and upper entrances, respectively (Fig. 4). The impact of the boundary conditions defined in Eqs. (3) and (4) will be analyzed in detail in section 5.1.

Assuming an adiabatic boundary at $r = R_{dom}$ yields

$$ \frac{\partial T_r}{\partial r}(x, R_{dom}, t) = 0 \quad \text{for} \quad 0 \le x \le L_{cave} \quad (5) $$

The energy conservation at the conduit wall reads:

$$ k_r \frac{\partial T_r}{\partial r} (x, R_p, t) = \varphi_c + L_v J_w \quad (6) $$

where $k_r$ is the thermal conductivity of the rock, $\varphi_c$ is the convective heat flux at the conduit wall (positive when directed from the rock to the air), $L_v$ is the molar latent heat of evaporation (J.mol⁻¹) and $J_w$ is the molar flux at the conduit wall (mol.m⁻².s⁻¹), positive for evaporation, negative for condensation. Eq. (6) states that the sum of the convective and
<figure><img src="image_8.png" /><figcaption>Fig. 6. The simplified geometry of the ventilated cave with a single conduit and two entrances at different elevations.</figcaption></figure>

**Table 2**

Abbreviated names of the four cases considered in the numerical simulations. The estimation of the transfer coefficients required for DA and HA models is detailed in Appendix B.

<table><thead><tr><th rowspan="2"></th><th colspan="2">Transfer coefficients</th></tr><tr><th>Finite</th><th>Infinite</th></tr></thead><tbody><tr><td>1- Dry air</td><td>DA</td><td>DA<sub>∞</sub></td></tr><tr><td>2- Humid air</td><td>HA</td><td>HA<sub>∞</sub></td></tr></tbody></table>

latent heat fluxes is equal to the conductive heat flux leaving the rock mass. This equation couples the rock mass with the airflow.

### 3.2.2. Airflow

Since the water vapor is always dilute (its mass fraction in the humid air never exceeds 2 %), the heat capacity of the water vapor can be neglected against that of the air. With this assumption, the energy balance in the air reduces to:

$$ \rho_a c_{v,a} A(x) \frac{\partial \bar{T}_a}{\partial t} + \dot{m} \left( c_{p,a} \frac{\partial \bar{T}_a}{\partial x} - g \cos(\theta(x)) \right) = P(x) \phi_c, \quad (7) $$

where $\bar{T}_a(x)$ is the mixing temperature of the air, i.e., the air temperature averaged over the conduit cross-section [14]:

$$ \bar{T}_a(x) = \frac{1}{\dot{m}} \iint T_a(x, r) \rho_a u(x, r) dA \quad (8) $$

where $T_a(x, r)$ is the local air temperature and $u(x, r)$ the axial component of the local velocity vector. $A(x)$ and $P(x)$ are the cross-sectional area and perimeter obtained from the cave survey and displayed in Fig. 2. $\dot{m}$ is the air mass flowrate measured by the SFM-device and displayed in Fig. 5 (positive when the air enters through the upper entrance). $\rho_a$, $c_{v,a}$ and $c_{p,a}$ are the air density, specific heat at constant volume and specific heat at constant pressure.

The first term in the left-hand side (LHS) of Eq. (7) is an approximate expression of the rate of variation of the air internal energy. This term is necessary to regularize the solution when rapid fluctuations occur during flow reversals. Most of the time, it is negligible compared to the other terms. Therefore, Eq. (7) mainly expresses a balance between the energy transported by the fluid (thermal and gravitational potential energy in the second term of the LHS) and the heat transferred to the conduit wall, in the right-hand side (RHS).

The inlet temperature is equal to the atmospheric temperature $T_{atmL}$ or $T_{atmU}$ according to the direction of the flow:

$$ \bar{T}_a(0, t) = T_{atmU}(t) \text{ for } \dot{m} > 0 \quad (9) $$

$$ \bar{T}_a(L_{cave}, t) = T_{atmL}(t) \text{ for } \dot{m} < 0 \quad (10) $$

The water vapor balance in the humid air and its associated boundary conditions are similar to Eqs. 7–10:

$$ A(x) \frac{\partial \bar{c}_v}{\partial t} + \frac{\dot{m}}{\rho_a} \frac{\partial \bar{c}_v}{\partial x} = P(x) J_w \quad (11) $$

$$ \bar{c}_v(0, t) = \phi_{ext} c_{s, atmU} \text{ for } \dot{m} > 0 \quad (12) $$

$$ \bar{c}_v(L_{cave}, t) = \phi_{ext} c_{s, atmL} \text{ for } \dot{m} < 0 \quad (13) $$

$\bar{c}_v$ is the water vapor molar concentration (mol.m$^{-3}$) averaged over the cross-section similarly to the temperature in Eq. (8). Index “s” indicates saturated condition defined as:

$$ c_s = \frac{p_{sat}(T)}{RT} \quad (14) $$

where $p_{sat}$ is the saturated water vapor pressure at temperature T [26]. R is the ideal gas constant. Thus, $c_{s, atmU}$ and $c_{s, atmL}$ are the saturated water vapor concentrations based on the atmospheric temperature close to the upper and lower entrances, respectively. $\phi_{ext}$ is the atmospheric relative humidity. The same constant value of 75 % is assumed at the lower and
upper entrances.

### 3.2.3. Heat and vapor fluxes at the conduit wall

The Newton's law of cooling is a general law that applies to any convection process [14]. This law states that the heat flux transferred from the wall to a flowing fluid is proportional to the temperature difference between the wall and the fluid:

$$q_c(x, t) = h_t(T_w(x, t) - \bar{T}_a(x, t)) \quad \text{for } 0 \le x \le L_{\text{cave}}, \qquad (15)$$

where $T_w(x, t) = T_r(x, R_p, t)$ is the wall temperature and $\bar{T}_a(x, t)$ the air mixing temperature defined in Eq. (8). The heat and mass transfer analogy provide the same kind of relation for the water vapor flux:

$$J_w(x, t) = h_m(c_{sw}(x, t) - \bar{c}_v(x, t)) \quad \text{for } 0 \le x \le L_{\text{cave}}, \qquad (16)$$

where $c_{sw}$ is the saturated water vapor concentration obtained by substituting $T$ for $T_w$ in Eq. (14). $h_t$ and $h_m$ are the heat and mass transfer coefficients, respectively. In a pipe flow, they depend on the physical properties of air (Table 3) and the conduit geometry defined by the cross-section $A(x)$ and the perimeter $P(x)$ (Fig. 2). In the turbulent regime, they also depend on the air flowrate $\dot{m}$ (Fig. 5) and the wall relative roughness $\epsilon$, defined as the ratio of the roughness over the hydraulic diameter $D_h(x) = 4A(x)/P(x)$. We set $\epsilon = 0.01$ in all simulations.

Assessing the heat and mass transfer coefficients requires the evaluation of the Reynolds number $Re$ that characterizes the ratio of inertia to viscous friction in the airflow:

$$Re(x) = \frac{\rho_a |\bar{u}(x)| D_h(x)}{\mu_a} = \frac{4 |\dot{m}|}{P(x) \mu_a} \qquad (17)$$

where $\bar{u} = \frac{\dot{m}}{(\rho_a A)}$ is the mean velocity over the cross-section and $\mu_a$ the dynamic viscosity of the air. At fixed air flowrate, the Reynolds number is inversely proportional to the conduit perimeter $P(x)$. At the maximum mass flowrate (of the order of 1 kg s⁻¹, see Fig. 5) $Re$ ranges from $10^4$ (in

<caption>**Table 3**</caption>

Physical properties used in the numerical simulations (reference temperature: 12 °C).

<table><thead><tr><th>Properties</th><th>Value and unit</th><th>Reference</th></tr></thead><tbody><tr><td>Rock density</td><td>ρ<sub>r</sub> = 2325 kg/m<sup>3</sup></td><td>[49]</td></tr><tr><td>Rock heat capacity</td><td>c<sub>p,r</sub> = 841.09 J/(kg K)</td><td>[49]</td></tr><tr><td>Rock thermal conductivity</td><td>k<sub>r</sub> = 2.302 W/(m·K)</td><td>[49]</td></tr><tr><td>Rock thermal diffusivity</td><td>α<sub>r</sub> = 1.177 × 10<sup>-6</sup> m<sup>2</sup>/s</td><td>[49]</td></tr><tr><td>Air dynamic viscosity</td><td>μ<sub>a</sub> = 1.77 × 10<sup>-5</sup> Pa·s</td><td>[14]</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub> = 1.23 kg/m<sup>3</sup></td><td>[14]</td></tr><tr><td>Air kinematic viscosity</td><td>ν = μ<sub>a</sub>/ρ<sub>a</sub> = 1.44 × 10<sup>-5</sup> m<sup>2</sup>/s</td><td>-</td></tr><tr><td>Air heat capacity at constant pressure</td><td>c<sub>p,a</sub> = 1007 J/(kg K)</td><td>[14]</td></tr><tr><td>Air heat capacity at constant volume</td><td>c<sub>v,a</sub> = 719.29 J/(kg K)</td><td>-</td></tr><tr><td>Air thermal conductivity</td><td>k<sub>a</sub> = 0.0251 W/(m·K)</td><td>[14]</td></tr><tr><td>Air thermal diffusivity</td><td>α<sub>a</sub> = 2.02 × 10<sup>-5</sup> m<sup>2</sup>/s</td><td>-</td></tr><tr><td>Ideal gas constant</td><td>R = 8.314 J/mol·K</td><td>-</td></tr><tr><td>Molar mass of water vapor</td><td>M<sub>w</sub> = 18.015 g/mol</td><td>-</td></tr><tr><td>Diffusion coefficient of water vapor in air</td><td>D<sub>w</sub> = 2.43 × 10<sup>-5</sup> m<sup>2</sup>/s</td><td>[25]</td></tr><tr><td>Molar specific vaporization latent heat</td><td>L<sub>v</sub> = 44.55 kJ/mol</td><td>[50]</td></tr><tr><td>Prandtl number</td><td>Pr = v/α<sub>a</sub> = 0.71</td><td>-</td></tr><tr><td>Schmidt number</td><td>Sc = v/D<sub>w</sub> = 0.59</td><td>-</td></tr></tbody></table>

the largest parts of the conduit) to $10^5$ (in the narrowest parts) denoting turbulent flow throughout the cave. The method used to estimate the transfer coefficients from standard correlations for fully developed flow in pipes in the laminar and turbulent regimes is detailed in Appendix B.

### 3.2.4. Infinite transfer coefficient and/or dry air

Eqs. 2–16 assume humid air and finite transfer coefficients (HA). In the limiting case of infinite heat transfer coefficients (HA<sub>∞</sub>), the conditions at the conduit wall Eqs. 15 and 16 are replaced by

$$T_a(x, t) = T_w(x, t) \quad \text{for } 0 \le x \le L_{\text{cave}}, \qquad (18)$$

$$c_v(x, t) = c_{sw}(x, t) \quad \text{for } 0 \le x \le L_{\text{cave}}. \qquad (19)$$

Injecting Eqs. 18 and 19 in Eqs. (7) and (11) provides the relations for $\varphi_w$ and $J_w$ required by Eq. (6). A noticeable consequence of Eqs. 18 and 19 is that the humid air inside the cave is saturated at the wall temperature all over the conduit.

In the case of dry air with finite transfer coefficients (DA), Eqs. (11)–(14) and (16) are ignored and we impose $J_w = 0$ in Eq. (6). When, in addition, the heat transfer coefficient is infinite (DA<sub>∞</sub>), Eq. (15) is replaced by Eq. (18).

### 3.2.5. Initial condition

The initial condition required for solving Eq. (2) is the temperature field at any point of the rock mass at the beginning of the cave monitoring (i.e., on August 13, 2020). Obviously, this information is not directly available, since the temperature is only known at a few locations where sensors have been installed. Because of the thermal inertia of the massif, the initial temperature field depends in a non-trivial way on the history of the atmospheric temperature and the air flowrate.

Interestingly enough, the upper entrance of Longeague was enlarged by cavers in 1982. This artificial modification likely induced a significant increase of the air flowrate through the main conduit. Although it is unlikely that the air flowrate before the opening was strictly equal to zero, it was certainly much lower than after the opening. We thus do the approximation that the main conduit was closed before this date, and the air flowrate insignificant. Considering this specific feature of Longeague Cave, an approximate initial temperature field was assessed assuming a simplified history consisting of two steps. We first assume no airflow in the conduit, resulting in an adiabatic conduit wall. More precisely, we solved the pure conduction problem defined by Eqs. (2)–(6), with $\varphi_c = 0$ and $J_w = 0$ in Eq. (6). The atmospheric temperatures in Eqs. (3) and (4) were approximated by Eq. (1) (displayed by the dashed lines in Fig. 4a), which only considers the annual fluctuations. Starting from an arbitrary uniform initial temperature field, the simulations of this 1D conduction problem was carried out over 200 years, a time long enough to reach the periodic regime. In a second step, we used this result as the initial temperature field to simulate the 38 years between the opening of the conduit in 1982 and the beginning of the monitoring in 2020. The full HA model defined in Section 3.2 was used with the simplified atmospheric temperatures of Eq. (1) as in the first step, along with the same kind of simplified function for the air flowrate (dashed line in Fig. 5a). The resulting temperature field was used as the initial condition for the numerical simulation of the year during which the cave was monitored.

### 3.3. Physical properties

The physical properties used in the numerical simulations are assumed constant. They are listed in Table 3. The saturated vapor pressure of water as a function of temperature was interpolated from data in Ref. [26].
## 3.4. Numerical methods

All the numerical simulations were performed with the commercial software Comsol Multiphysics, version 6.1. This software solves partial differential equations by finite elements (Galerkin method). Time discretization was based on implicit backward differentiation formula. The mesh was made of Lagrangian quadratic elements, including 12600 2D elements in the rock mass (mapped mesh) and 200 1D elements in the air. The hyperbolic Eqs. (7) and (11) required the implementation of the stabilization techniques detailed in Appendix C.

# 4. Results

## 4.1. Field data

The AMT, ATF and DTF deduced from the FFT of the field data are displayed by black symbols in Fig. 7a to c. Thermal anomalies are clearly visible in Fig. 7a. Indeed, the AMT passes through a local maximum at ~9.5 °C at H3 station located at 55 m from the upper entrance (warm anomaly). Conversely, a local minimum at ~5.7 °C is observed at R4 station at 20 m from the lower entrance (cold anomaly). The warm anomaly is more pronounced than the cold one (approximately +2 °C and -1 °C relative to the outside temperatures at the upper and lower entrances, respectively). Between these two extrema (i.e., between H3 and R4 stations), the AMT continuously decreases with increasing x. The absence of plateau in the temperature profile suggests that the convection length is larger than half the cave length, and that there is no deep karst region in Longeaigne Cave.

Fig. 7b displays the ATF. A sharp decrease is observed close to the entrances, as expected in diffusive regions [15]. The ATF then gradually declines in the convective region to reach a minimum nearby the middle of the conduit.

The DTF is displayed in Fig. 7c. A sharp decrease is observed close to the entrances, similarly to the ATF. However, in the convective region, the DTF decreases over a much shorter distance than the ATF. It becomes hardly measurable beyond a few tens of meters from the entrances. This suggests that the distance of propagation of a thermal perturbation advected by the airflow should depend on the frequency.

The time series of field data from three selected stations are displayed by solid black lines in Fig. 8. Station H3 (Fig. 8a) is close to the upper entrance, station R2 (Fig. 8b) is approximately in the middle of the cave, and station R4 (Fig. 8c) is close to the lower entrance. Comparisons with the atmospheric temperatures displayed in Fig. 4a confirms the strong damping of the daily temperature fluctuations. In contrast, the impact of cave closures on the temperatures is clearly observed at all stations displayed in Fig. 8 (see for instance the time range from 120 to 196 days).

## 4.2. Numerical simulations

Fig. 7a displays the AMT of the conduit wall and the air as a function of the distance from the upper entrance, for the four models defined in Table 2 (notice that in $DA_\infty$ and $HA_\infty$ models, wall and air temperatures are equal). All models significantly underestimate the AMT computed from the field data (approximately by 1 or 2 °C). $DA_\infty$ and $HA_\infty$ predict larger thermal anomalies compared to DA and HA. $HA_\infty$ seems to better reproduce the shape of the field data close to the upper entrance.

Regarding the ATF (Fig. 7b), the four models reproduce qualitatively the main trends observed in the field data. All of them predict a sharp decrease over the first few meters from the entrances (i.e., in the diffusive regions) followed by a more gradual decline. $DA_\infty$ globally underestimates the field ATF. $HA_\infty$ is in good agreement with the field data in the lower half of the cave, but overestimates the ATF at $x = 55$ m (by approximately 1 °C). The air and wall temperatures predicted by DA or HA models significantly differ from each other. The simulated air temperatures of both models overestimate the field data by approximately

1 °C all over the cave. A significantly better agreement is obtained with the simulated wall temperatures. HA model shows a maximum error of approximately half a degree at R4 station, at 17 m from the lower entrance.

Similar comments apply to the DTF displayed in Fig. 7c. $DA_\infty$ and $HA_\infty$ slightly overestimate the field data close to the upper entrance. The best fit is obtained for the wall temperature with DA and HA, whereas the air DTF predicted by these models is significantly above the field data. This is confirmed by the time series at R4 station, at 20 m from the lower entrance (Fig. 9). This station is close enough to the lower entrance to get a significant DTF when the lower entrance operates as an inlet. Compared to the field data, the wall temperature predicted by the HA model is shifted to lower values, with comparable (slightly lower) amplitudes. In contrast, the model predicts that the air temperature fluctuates with a much larger amplitude than the field temperatures.

The comparison between the simulated wall temperature and the field temperatures is completed by the time series displayed in Fig. 8a to c. At H3 ($x = 55$ m) and R2 ($x = 140$ m) stations, the simulated wall temperatures follow most of the time the field temperatures with a quasi-constant shift to lower values, in agreement with the difference between the simulated and field AMT observed in Fig. 7a. An exception is the first 30 days of the monitoring period. During this time range, the field temperatures globally decrease, whereas all models predict increasing temperatures. This discrepancy is due to the approximated initial condition used in the simulations. The actual history of the atmospheric temperature before the initial time, different from the periodic function assumed in the simulations (see Section 3.2.5), results in significant errors during the first month. This time lapse is short compared to the duration of the monitoring (1 year).

The behavior observed in Fig. 8c at R4 station ($x = 294$ m), located in the lower part of the cave, is different from H3 and R2. The evolution of simulated and measured temperatures not only differs during the first month of the simulated time range, but the field temperature also shows several peaks not predicted by the simulations. They are indicated by red arrows in Fig. 8c. The most significant event is not localized at a specific time, but extends over approximately two months, from time $t = 310$ days to the end of the monitoring period (see Fig. 8c). The same behavior has been observed at station H6 (figure not shown). These events are correlated with the lake level. During floods, H6 and R4 stations, both located in the lower part of the cave, are submerged, or very close to the water stream which modifies the temperature field in the cave. Because of thermal anomalies, a downward water flow induces heat transfer from the warmer upper part of the cave to the colder lower part (compare for instance the time series of the field temperatures at stations H3 ($x = 55$ m) in Fig. 8a and R4 ($x = 294$ m) in Fig. 8c). This explains that, in the lower part of the cave, floods always result in temperature increases. In Longeaigne cave, the rise of the atmospheric temperature with the altitude (see Section 2.2) might also play a role.

As stated in Section 3.1, heat advection by water flow is not considered in our models. Taking into account this effect would require the development of a hydraulic model, an intricate task beyond the scope of this article. However, comparisons between field and simulated temperatures reveal that the effect of water flow is only perceptible in the lower part of the conduit. At higher elevation, the climate is controlled by the airflow coupled to heat conduction in the rock mass. This is likely a general pattern in case of diffuse recharge through the epikarst. At the top of the massif, water infiltration is distributed all over the catchment area inducing low water velocity and thus negligible advection compared to heat conduction in the rock mass. When flowing downward, water concentrates in specific conduits until it reaches the spring where the water velocity, and thus the advection of heat, are maximum. The situation would be different in case of concentrated recharge (swallowing stream), as described by Ref. [27–29]. In this case, the effect of water advection on the temperature field would be significant throughout the conduit, including its upper part.

The effect of the uncertainty on the airflow rate has been investigated
<figure><img src="image_9.png" /><figcaption>Fig. 7. AMT(a), ATF (b) and DTF (c) from the field data and the numerical simulations. Data at x = 0 and x = L_cave = 311 m correspond to the atmospheric temperature close to the upper and the lower entrances, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Temperature time series in three stations inside the cave a) H3 (x = 55 m); b) R2 (x = 140 m); c) R4 (x = 294 m), red arrows indicate temperature peaks in the field data not predicted by the model. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Temperature time series from field data and HA model at R4 station, at 17 m from the lower entrance. Positive airflow indicates that the air enters through the upper entrance (origin of time: August 13, 2020).</figcaption></figure>

by recalculating the outputs of DA and HA models after multiplying the airflow rate by 1.5 (see Appendix D). Although non-negligible, the impact on the results during the monitoring year is limited, and does not modify the observations drawn from the comparisons between the simulated and field temperatures. However, these observations raise issues requiring in-depth analyses developed in the next section:

1) Models systematically underestimate the annual mean temperatures obtained from field data;

2) Although the temperature sensors are located in the air, the amplitude measured in the field are much closer to temperature simulated for the cave walls than for the air;

3) Models with finite or infinite transfer coefficients yields significantly different results. What is the most relevant assumption?

# 5. Discussion

## 5.1. Annual mean temperatures

As pointed out in Section 4.2, the numerical simulations underestimate by approximately 1 or 2 °C the AMT obtained from the field data (see Fig. 7a). The predominant cause of this discrepancy, which is the main weakness of the model, must be sought in the boundary conditions at the external surfaces (i.e., the interfaces between the rock mass and the atmosphere, see Fig. 6). Indeed, Eqs. (3) and (4) impose that the external surface temperature is equal to the atmospheric temperature. However, several effects may induce a significant gap between these temperatures. The energy balance at an external surface reads [18]:

$$(-k_r \cdot \nabla T_r) \cdot \vec{n} = \varphi_{\text{conv}} + \varphi_{\text{rad,lw}} - \varphi_{\text{rad,sw}} + \varphi_{\text{evap}} \quad (20)$$

This equation states that the conduction flux leaving the rock (LHS) is the sum of all the thermal fluxes transferred from the external surface to the external environment (RHS). $\vec{n}$ is the normal unit vector pointing to the atmosphere, $\varphi_{\text{conv}}$ the convective heat flux from the rock to the atmosphere, $\varphi_{\text{rad,lw}}$ the net radiative flux lost by the rock in the long wave-length range, $\varphi_{\text{rad,sw}}$ the sun irradiation (short wave-length range) and $\varphi_{\text{evap}}$ the latent heat flux (positive for evaporation). Some of these terms tend to increase the soil temperature (e.g., the sun irradiation), others to decrease it (e.g., the latent heat of evaporation). Molnar [30] analyzed the atmosphere and soil AMT from 212 sites throughout the world. He concluded that the soil temperature is generally warmer than the atmosphere. The difference between the AMT of the soil and the atmosphere mainly depends on the land surface cover. It approximately ranges from 1 °C in wetlands and forests to 3–5 °C or more in arid or cold

regions. Fig. 10a displays the values measured at different weather stations located in Switzerland. Most data are included in the range from 0 to 2 °C, with significant variations from year to year.

The sensitivity of the conduit temperature to the temperature of the external surfaces was assessed from additional numerical simulations. Fig. 10b displays the AMT of the conduit wall obtained from the HA model in two cases: (a) equal atmospheric and wall temperatures, as specified by Eqs. (3) and (4) (blue curves), (b) external surface temperatures increased by $\Delta T = 2$ °C compared to the atmospheric temperatures (red curves). In the latter case, the boundary conditions

$$T_r(0, r, t) = T_{\text{atmU}}(t) + \Delta T \quad \text{for } R_p(0) \le r \le R_{\text{dom}}, \quad (21)$$

$$T_r(L_{\text{cave}}, r, t) = T_{\text{atml}}(t) + \Delta T \quad \text{for } R_p(L_{\text{cave}}) \le r \le R_{\text{dom}}, \quad (22)$$

have been used instead of Eqs. (3) and (4). Setting $\Delta T = 2$ °C is a rough approximation since the actual difference between the atmospheric and external surface temperatures is unknown and may vary with time and space. However, this is relevant for a sensitivity analysis. To simulate a large time range (200 years) with a reasonable computational time, we used the simplified atmospheric temperatures of Eq. (1) and the same kind of simplified function for the air flowrate (dashed line in Fig. 5a). The initial condition was determined using the method detailed in Section 3.2.5. Fig. 10b shows that increasing the external surface temperatures by 2 °C results in a similar temperature raise throughout the conduit wall, over the entire time range considered in the simulations (200 years). This result strongly suggests that the main source of discrepancy between the field and simulated AMT at the conduit wall may be the boundary conditions at the external surfaces.

Improving the accuracy of the model thus necessitates to replace Eqs. (3) and (4) with more realistic boundary conditions, based for instance on the energy balance Eq. (20). This would require appropriate models for assessing the various terms included in this equation (see for instance Ref. [31]), and thus the monitoring of additional input data (e.g., the wind velocity for the assessment of $\varphi_{\text{conv}}$ or the sky temperature for $\varphi_{\text{rad,lw}}$). An alternative could be to disseminate temperature sensors at a few tens of centimeters below the external surface and use the measured temperatures in Eqs. (3) and (4) instead of the atmospheric temperature.

## 5.2. Temperature fluctuations

The ATF and DTF of the air predicted by DA and HA models show significant discrepancies with field data (see Fig. 7b and c). $DA_{\infty}$ and $HA_{\infty}$ better fit the measured values, but are slightly less accurate than the wall temperatures yielded by DA and HA models. This is unexpected
<figure><img src="image_12.png" /><figcaption>Fig. 10. a) Annual mean temperature difference between the soil and the atmosphere ($T_{soil} - T_{am}$) for some different weather stations in Switzerland (from Federal Office for Meteorology and Climatology of Switzerland [48]). The soil temperature is measured at 5 cm depth, the atmospheric temperature at 2 m above the ground. b) Effect of the external surface temperatures on the AMT of the conduit wall (numerical simulations with simplified time functions for the air flowrate and the atmospheric temperature, HA model). Blue curves: external surface temperature equal to the atmospheric temperature (Eqs. (3) and (4)). Red curves: external surface temperatures increased by $\Delta T = 2$ °C (Eqs. 21 and 22). Dashed lines: initial conditions obtained by solving the diffusion problem. The arrows point to increasing times. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

since the temperature sensors are located in the airflow, not in the rock. There are two possible explanations for this: (1) DA and HA models fail in predicting the air temperature, $DA_{\infty}$ and $HA_{\infty}$ models based on infinite transfer coefficients and equal temperatures in the air and the wall are closer to reality; (2) The wall temperature from DA and HA models fit well the field temperatures because the temperature sensors measure the wall temperature rather than the air temperature. Although counterintuitive, the latter explanation is the most likely. A first point is that the local air temperature $T_a$ varies within the cross-section, and reaches the rock temperature at the conduit wall. Incidentally, the mixing temperature of the air $\bar{T}_a$ provided by the numerical simulations is a cross-sectional averaged temperature. Considering its definition in Eq. (8), we expect $\bar{T}_a$ to be quite close to the local temperature in the center of the conduit, where the local velocity is maximum. A second point to be considered is that a temperature sensor always measures its own temperature [32], which might differ from the temperature of the surrounding fluid.

Fig. 11 displays the temperatures measured in the air at H1 station

(at 16 m from the upper entrance) during approximately a month. Two temperature sensors were installed, a Hobo Water Pro v2 at a few centimeters from the wall (similar to the other probes used in this study) and a Pt100 sensor suspended approximately in the center of the conduit. Fig. 11 shows that the amplitude of the daily oscillations measured by the Pt100 is two or three times larger than that measured by the Hobo Water sensor.

As pointed out by Lundström and Mattsson [32], a sensor immersed in a transparent fluid (the air in the present case) receives a convective flux from the fluid and a net radiative flux from the conduit wall. Therefore, its temperature is a weighted average of the air and wall temperatures. A simple model based on the energy conservation yields the following expression for the sensor temperature [32]:

$$T_{\text{sensor}} = \frac{h_{\text{conv}}}{h_{\text{conv}} + h_{\text{rad}}} T_a + \frac{h_{\text{rad}}}{h_{\text{conv}} + h_{\text{rad}}} T_w \quad (21)$$

where $h_{\text{conv}}$ and $h_{\text{rad}}$ are the convective and radiative transfer coefficients at the sensor surface. Eq. (21) qualitatively explains the results displayed
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Temperature measured by a Hobo Water sensor at a few centimeters from the conduit wall and a Pt100 sensor in the center of the conduit. Test section H1 located at 16 m from the upper entrance. Positive flowrate indicates airflow entering through the upper entrance.</figcaption></figure>

in Fig. 11. A temperature sensor measures the air temperature if convection prevails over radiation (i.e., if $h_{conv} \gg h_{rad}$). Otherwise, it measures a temperature in between the wall and air temperatures, or even the wall temperature if $h_{conv} \ll h_{rad}$. $h_{conv}$ is an increasing function of the air velocity and a decreasing function of the sensor size. $h_{rad}$ is proportional to the emissivity of the sensor surface. Compared to the Pt100 sensor, the Hobo Water sensor is larger (30 mm against 5 mm) and has a higher emissivity (polypropylene sheath against metallic sheath). In addition, the Hobo Water sensor is located closer to the wall, where the air flows with a smaller local velocity and a local temperature $T_a$ closer to the wall temperature $T_w$. All of these effects bring the Hobo Water sensor nearer the wall temperature. This sensor thus yields a better estimate of the wall temperature $T_w$ than the air mixing temperature $\bar{T}_a$. In contrast, the temperature of the Pt100 sensor yields a better estimate of the air temperature in the center of the conduit. This interpretation, consistent with the results displayed in Fig. 11, explains the good agreement between the field data and the wall temperature predicted by DA and HA models, and the discrepancy with the air temperature. Indeed, the configuration of all sensors used in this work is similar to that of the Hobo Water sensor used in the comparison displayed in Fig. 11.

### 5.3. Finite versus infinite transfer coefficients

As pointed out in the introduction, assessing relevant values of the transfer coefficients in the complex geometry of a cave is a difficult task. However, the impact of transfer coefficients on the prediction of cave climate can be more or less significant depending on whether heat transfer is convection-limited or conduction-limited. In a previous work [15], we showed from numerical simulations that no accurate estimation of the heat transfer coefficient is required when the Reynolds number corresponding to the maximum mass flowrate is larger than $3 \times 10^5$ (convection-limited regime). Conversely, the heat transfer coefficient is a key information when the maximum Reynolds number is lower than $10^4$ (convection-limited regime).

The maximum Reynolds number in Longeague Cave ranges from $10^4$ to $10^5$, depending on the cave section (see Section 3.2.3). Longeague Cave is thus in the intermediate regime, in which both convection in the air and conduction in the rock mass must be considered. Consequently, the temperature fields predicted by DA and HA (finite transfer coefficients) or by $DA_\infty$ and $HA_\infty$ (infinite transfer coefficients) are significantly different. If we assume that:

(1) the temperature sensors better estimate the wall temperatures than the air temperatures (see Section 5.2),

(2) the general underestimation of the AMT by all models is mainly due to the initial and boundary conditions (see Section 5.1),

then DA and HA models yield more accurate results compared to $DA_\infty$ and $HA_\infty$ (see Fig. 7). With this in mind, the fact that $DA_\infty$ and $HA_\infty$ better reproduce the measured AMT in the upper part of the cave (see Fig. 7a) is likely a coincidence resulting from error compensation. In line with this, Fig. 7b shows that, in the same region of the cave, the ATF is underestimated by $DA_\infty$ and overestimated by $HA_\infty$.

As pointed out in Section 3.1, DA and HA models are two limiting cases (no latent heat effect in the former case, permanent water film on the cave walls in the latter case). Comparisons between the outputs of both models suggest that the effect of the latent heat on the cave climate is rather weak (see Figs. 7 and 8). Considering the uncertainties on the airflow rate and the temperature measurements, the accuracy of DA or HA models for the prediction of the temperature fields are approximately equivalent.

## 6. Condensation and evaporation

### 6.1. Vapor transfer rate and consequences on paleoclimatology

Fig. 12a displays the mass fluxes of water vapor at the conduit wall (positive and negative for evaporation and condensation, respectively) at two selected times corresponding to high air flowrates (HA model). Both mass fluxes are close to the maximum values reached during the year, one in summer ($t = 32.6$ days, $\dot{m} = 1.26 \text{ kg.s}^{-1}$, air intake from the upper entrance), the other in winter ($t = 220.3$ days, $\dot{m} = -1.92 \text{ kg.s}^{-1}$, air intake from the lower entrance).

In summer, evaporation is only observed in areas very close to the entrances: less than a meter from the upper entrance, a few meters from the lower entrance. Evaporation at the upper entrance arises from the need to increase the relative humidity from 75 % (the atmospheric humidity) to 100 % before condensation starts. At the lower entrance, evaporation results from the increase in the wall temperature in the diffusive region. With the exception of these two short regions, condensation takes place everywhere over the conduit wall. As expected, the condensation rate is maximum close to the inlet (the upper entrance), where the temperature contrast between the air and the wall is greatest. However, the condensation rate does not decrease monotonically with increasing distance from the upper entrance. Oscillations
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** a) Vapor mass flux at the conduit wall corresponding to high air flowrate events in summer (t = 32.6 days) and winter (t = 220.3 days); insets: zooms in the first 10 m from the entrances b) Evaporation flux integrated over the conduit wall as a function of time. In both figures, positive and negative values denote evaporation and condensation, respectively (HA model).</figcaption></figure>
are correlated with conduit size variations (compare Fig. 12a with the conduit perimeter displayed in Fig. 2). We demonstrate in Appendix B that, in the turbulent regime and for a given air flowrate, the transfer coefficients are inversely proportional to the product of the conduit perimeter by the hydraulic diameter. This quadratic dependence results in a strong effect of the conduit size on the vapor transfer rate, with highest values in the narrowest parts of the conduit.

In winter, condensation is limited to a few meters from the upper entrance (the air outlet in winter). Condensation in this short region results from the decrease in the conduit wall temperature in the diffusive region. Evaporation takes place anywhere else. As expected, the evaporation rate is maximum close to the inlet (the lower entrance). Despite the higher air flowrate in winter, the evaporation rate is lower than the condensation rate in summer, for two reasons: (a) lower temperatures in winter (especially near the inlets) resulting in lower saturated vapor pressure and thus lower vapor concentration in the air; (b) larger conduit size in the lower part of the cave compared to the upper part (see Fig. 2), resulting in lower transfer coefficients. The total amount of condensation and evaporation is assessed and discussed in Section 6.2.

A quantitative estimate of a cave's thermal response to the outside environment is essential for the interpretation of paleoclimate records from speleothems [33]. Although variations in the drip water composition largely control the proxy partitioning in speleothems [34], trace elements and oxygen isotope fractionation during CaCO₃ precipitation are also temperature-dependent [35,36] and become significant for temperature amplitudes exceeding ±0.5 °C. The model developed in this study is capable of calculating this amplitude on daily and yearly scale at any location within the cave and illustrates the temperature sensitivity in the convection zone. Moreover, Fig. 12a (red line) reveals that, in winter, evaporation is present all along the cave at significant rates but, in particular, close to the lower cave entrance. There, the evaporation flux may reach $2 \times 10^{-3} \text{ g m}^{-2} \text{ s}^{-1}$ during intense upward air flow ($t = 220.3$ days; $\dot{m} = -1.92 \text{ kg s}^{-1}$). Whilst this rate is still significantly lower than the potential evaporation rate considered in geochemical models [37], it may become locally important at very low drip rates feeding speleothems. Results from the highly ventilated Longeigue cave thus suggest that the role of evaporation on the isotope proxy record can often be neglected in hydrologically active caves but must be considered in semi-arid environments with only sporadic drips

## 6.2. Role of condensation in water production

Evaporation and condensation consume and produce water inside the cave. Fig. 12b displays the evaporation rate integrated throughout the conduit wall as a function of time (HA model, negative values refer to condensation). Condensation takes place during approximately two months in summer. The total amounts of condensed and evaporated water during the year were $-8.8 \times 10^3$ and $16.7 \times 10^3$ kg respectively. The mass of evaporated water throughout the year is twice the mass of condensed water. These values are low (<0.1 %) compared to the water flow-through in the groundwater catchment of Raies/Longeigue system. Nonetheless, our model also shows that condensation may happen, preferentially at the upper entrance, during downward ventilation regimes. Whilst the total amount of condensed water during the 2020-21 annual cycle reached 8800 kg, it is still negligible with respect to the hydrological mass balance of Longeigue cave. In arid environments, this amount may nonetheless be sufficient to maintain a moist atmosphere in the upper entrance zone of a cave system and thus support the local ecosystem. Because cave-adapted species may be sensitive to even small changes in the subsurface environment [38,39], prediction of the amplitude and frequency of temperature changes might be crucial for cave conservation issues.

## 7. Conclusion

The main conduit of Longeigue cave is ventilated by a strong airflow driven by chimney effect. In order to predict the temperature field in this conduit, we developed a thermal model based on the diffusion equation in the rock mass along with the conservation of energy and vapor mass in the airflow. The conduit wall and the air are coupled by transfer coefficients assessed from standard correlations for fully developed forced convection in pipes. The model has been tested against the field data collected in the cave during a complete year.

The complex geometry of a cave is a serious difficulty for thermal modeling. A major simplification was done by assuming a rectilinear conduit of circular cross-section with variable diameter. In this context, the conduction problem in the rock mass is 2D axisymmetric, and only two independent parameters are required to characterize the conduit cross-section: the effective diameter $D_p$ based on the conduit perimeter and the hydraulic diameter $D_h$. Using $D_p$ in the conservation equations preserves the exchange surface between the air and the rock, while $D_h$ is the relevant geometrical characteristic to be used in the estimation of the transfer coefficients.

The numerical results underestimate by 1 or 2 °C the annual mean temperatures obtained from field data. This is the most significant weakness of the model. This discrepancy is mainly due to the assumption of equal rock and atmospheric temperatures at the external surface of the massif. More realistic boundary conditions are proposed.

In contrast, the model accurately predicts the temperature fluctuations from daily to yearly time scales, which includes in the specific case of Longeigue Cave the temperature fluctuations generated by intermittent cave closures. The impact of the initial condition assessed from a simplified history is limited to the first month after the beginning of the simulation. Moreover, comparisons between field and simulated temperatures reveal that the effect of water flow due to diffuse water recharge is only perceptible in the lower part of the cave. At higher elevation, the temperature field is controlled by the airflow coupled to heat conduction in the rock mass. The latent heat effect of evaporation and condensation seems to play a minor role on the cave climate. However, predicting the order of magnitude of condensation or evaporation rate is valuable for applications related to ecology or paleoclimatology. In arid environments, maintaining moist atmosphere in the upper part of caves might be crucial for the ecosystems. The order of magnitude of evaporation rates also show that, in temperate European caves, the role of evaporation on the isotope proxy record can often be neglected.

Compared to heat diffusion in a rock mass, cave ventilation significantly increases the depth of natural thermal oscillations and decreases the time of propagation. This makes it necessary to develop and test thermal models for the prediction of temperature fields in ventilated karst massifs. We provided and discussed simplifying assumptions that allow the accurate prediction of the temperature field in a ventilated cave from daily to yearly time scale. However, extrapolating the results of this study to larger time horizons (centuries or more) must be undertaken with great caution. The long-term impact of geometry simplifications, or initial and boundary conditions, will require further research. This study is a first step towards the elaboration and validation of models capable of tackling these issues. Our model opens the way for a quantitative assessment of the cave's thermal response at any location, providing a well-known cave geometry.

The thermal model presented in this article is applied to caves naturally ventilated by chimney effect. However, some issues raised in this article are relevant to engineering applications with artificial ventilation, as mines [40,41] and tunnels [42]. Defining the right boundary conditions or convection model also matters in these configurations. More specifically, there is a growing interest for the recovery of geothermal energy from tunnels [43]. Taking into account the advection of heat by the airflow and the coupling between convection and conduction in the rock mass should improve performance predictions for
these complex geotechnical structures.

### CRediT authorship contribution statement

**Amir Sedaghatkish:** Writing – original draft, Visualization, Software, Methodology, Investigation, Conceptualization. **Claudio Pastore:** Methodology, Investigation, Data curation. **Frédéric Doumenc:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. **Pierre-Yves Jeannin:** Writing – review & editing, Supervision, Project administration. **Marc Luetscher:** Writing – review & editing, Supervision, Conceptualization.

### Data availability

The Comsol Multiphysics file for HA scenario (humid air with finite heat transfer) is available in https://doi.org/10.5281/zenodo.14535340

associated with [47] but a license of the software is necessary to use it.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This project was supported by the Swiss National Science Foundation (project no. 200021_188636). Open access funding is provided by university of Neuchâtel. We would like to thank Sven Friedel for his support and providing the COMSOL Multiphysics license for the Thermokarst project.

## Appendix A. Effect of the cross-sectional shape on wall temperature and heat flux

The numerical simulation of 3D geometries requires a large amount of computational resources. The objective of this appendix is to test whether a 2D (axisymmetric) circular cross-section can be substituted for the 3D geometry of a real conduit with an acceptable loss of accuracy. To address this issue, we computed the temperature field in the rock mass surrounding the conduits displayed in Fig. A1. Three cases are considered: (a) the circular shape, (b) the triangular shape as an instance of simple 3D geometry, (c) a more complicated shape obtained from a survey in a real cave. The conduit perimeter P is expected to be an important parameter since it imposes the exchange surface between the rock and the air (see Eqs. (7) and (11)). Therefore, all the comparisons between the different shapes displayed in Fig. A1 will be done between conduits of equal perimeter P.

The transient conduction equation is solved in the cross-sectional plane assuming uniform temperature in the direction parallel to the conduit. The temperature of the dry air inside the conduit follows a sinusoidal function of time. The heat flux at the conduit wall is deduced from the Newton's law of cooling, assuming uniform heat transfer coefficient all along the conduit circumference.

The problem variables are the coordinates x and y, the time t, the temperature in the rock mass $T_r$, and the temperature of the air $\bar{T}_a$. The corresponding dimensionless variables read:

$$X = \frac{x}{L_d}, \quad Y = \frac{y}{L_d}, \quad t^* = \frac{t}{\tau}, \quad \theta_r = \frac{T_r - T_m}{\Delta T}, \quad \theta_a = \frac{\bar{T}_a - T_m}{\Delta T} \qquad (A1)$$

where $T_m$ is the mean temperature of the air and $\Delta T$ the amplitude of the temperature fluctuations of period $\tau$. $L_d = \sqrt{\alpha_r \tau}$ is the diffusion length in the rock mass. $P^* = \frac{P}{L_d}$ and $D_p^* = \frac{D_p}{L_d} = \frac{P}{\pi L_d}$ are the dimensionless perimeter and perimeter-based diameter of the conduit, respectively.

The diffusion equation reads:

$$\frac{\partial^2 \theta_r}{\partial X^2} + \frac{\partial^2 \theta_r}{\partial Y^2} = \frac{\partial \theta_r}{\partial t} \qquad (A2)$$

The Newton's law of cooling yields the boundary condition at the conduit wall:

$$(-\nabla\theta_r) \cdot \vec{n} = Bi(\theta_r - \theta_a) \quad \text{with} \quad \theta_a = \sin(2\pi t^*) \qquad (A3)$$

where the Biot number $Bi = \frac{hL_d}{k}$ is the dimensionless heat transfer coefficient. $\vec{n}$ is the outward unit normal vector. From air mass flowrate measurements and correlations in Appendix B, the maximum heat transfer coefficient in Longeaigne is of the order of 20 W m⁻² K⁻¹. Therefore, $Bi$ can be as high as 3 for the daily fluctuations ($\tau = 1$ day and $L_d = 0.32$ m) and 50 for the annual fluctuations ($\tau = 1$ year and $L_d = 6.1$ m). The computational domain is sufficiently large to make the external boundaries adiabatic. The simulation time is increased until the periodic regime is reached.

Fig. A2 shows the amplitude of the local temperature fluctuation along the conduit circumference. The 3D shape of the triangular and real cross sections results in non-uniform fluctuations, lower in the corners and larger in the tips. However, the amplitudes averaged over the conduit circumference are close from each other (less than 4% difference). Therefore, if we just need the mean temperatures, the circular cross-section (a) with the same perimeter as the 3D shapes (b) and (c) yields satisfying results.
<figure><img src="image_15.png" /><figcaption>**Fig. A1.** Three different conduit shapes (the figures are not scaled). Cross-section (c) was surveyed at 47.7 m from the entrance of D7.1 Cave located in Sieben Hengste, Switzerland [15].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. A2.** Amplitude of the temperature fluctuations along the conduit circumference (solid lines) and corresponding averaged values (dashed-dotted lines), for $Bi = 15.9$ and $D_p^* = P^*/\pi = 0.378$. For case (b), the origin of the curvilinear coordinate is taken in a corner of the triangle. For case (c), the position of the maximum is indicated by a star in Fig. A1c.</figcaption></figure>

We performed a systematic parametric study varying the Biot number $Bi$ and the perimeter-based diameter $D_p^*$. The amplitude of the wall temperature variations is displayed in Fig. A3. The maximum difference between circular and 3D cross sections is less than 20 %. We are also interested in the heat flux through the conduit wall integrated over the conduit perimeter $\Phi = \int_p Bi(\theta_r - \theta_a)dl$ (where $l$ is the curvilinear coordinate over the conduit circumference). Fig. A4 compares the values obtained for the three cross-sectional shapes. The difference between the circular cross section and the other ones increases with Biot number. It never exceeds 25 %.

This simple study suggests that substituting a conduit with a circular cross-sectional shape to the real 3D geometry results in a reasonable loss of accuracy, if we content ourselves with the temperatures and heat fluxes averaged over the conduit circumference.
<figure><img src="image_17.png" /><figcaption>Fig. A3. Amplitude of the temperature fluctuations averaged over the conduit perimeter, for different values of the perimeter-based diameter $D_p^* = P^*/\pi$ and the Biot number $Bi$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. A4. Fluctuation amplitude of the wall heat flux integrated over the conduit circumference, for different values of the perimeter-based diameter $D_p^* = P^*/\pi$ and the Biot number $Bi$.</figcaption></figure>
## Appendix B. Estimation of heat and mass transfer coefficients at the conduit wall

In the turbulent regime, the heat transfer coefficient is obtained from correlations between the Reynolds (Eq. (17)), Prandtl and Nusselt numbers, respectively defined as:

$$ Pr = \frac{c_{p,a}\mu_a}{k_a} \quad \text{and} \quad Nu = \frac{h_a D_h}{k_a} \qquad (B1) $$

Nu and Re are based on the hydraulic diameter $D_h$ to account for the specific shape of the conduit. The Gnielinski correlation [44] has the advantage to take into account the wall roughness. This correlation reads:

$$ Nu = \frac{\left(\frac{f_d}{8}\right) (Re - 1000) Pr}{1 + 12.7 \left(\frac{f_d}{8}\right)^{0.5} (Pr^{2/3} - 1)} \quad \text{for } Re \ge 4000 \qquad (B2) $$

The Darcy friction factor $f_d$, which depends on the wall relative roughness $\epsilon$ and the Reynolds number, was estimated from the Haaland correlation [45]:

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

In the laminar regime, the Nusselt number is equal to a constant that depends on the conduit shape. For a circular cross-section, we get:

$$ Nu = 3.66 \quad Re \le 2000 \qquad (B4) $$

We did not try to adapt the value of the constant to the actual shape of the conduit. Indeed, the heat transfer coefficient in the laminar regime is small. It is thus not expected to have significant effect on the results.

In all simulations, the Prandtl number was set to $Pr = 0.71$ (from Table 3), and the wall relative roughness to $\epsilon = 0.01$. The resulting function $Nu = f(Re)$ is displayed in Fig. B1. A “smooth” transition between the laminar and turbulent regimes was achieved using the step function of Comsol Multiphysics for 2000 < Re < 4000.

<figure><img src="image_19.png" /><figcaption>Fig. B1. Nusselt number as a function of the Reynolds number for three different relative wall roughness and $Pr = 0.71$.</figcaption></figure>
Using the Lewis analogy [25], the mass transfer coefficient $h_m$ is obtained by substituting the Sherwood number $Sh = \frac{h_m D_h}{D_w}$ and the Schmidt number $Sc = \frac{\nu}{D_w}$ for Nu and Pr in Eqs. (B2) and (B4). It is worth noting that gases verify $Sc \approx Pr$ (from Table 3, $Sc \approx 0.59$ and $Pr \approx 0.71$). The Sherwood number $Sh$ used in the simulations is thus close to the Nusselt number $Nu$ displayed in Fig. B1.

A simple relation between transfer coefficients on one hand and the air flowrate and the conduit geometry on the other hand can be easily deduced from Fig. B1. This figure shows that, in the turbulent regime (i.e., for $Re \ge 4000$) and for a rough tube, the Nusselt number is approximately proportional to the Reynolds number (these results can also be deduced from Eqs. (B2-3)). The definitions of the Nusselt and Reynolds numbers yields:

$$h_t \propto \frac{\dot{m}}{P D_h} \qquad (B.5)$$

The same relation can be inferred for $h_m$ using the Lewis analogy. The heat and mass transfer coefficients are thus proportional to the air flowrate and inversely proportional to the product of the perimeter by the hydraulic diameter.

## Appendix C. Solver stabilization

The Galerkin method used in Comsol Multiphysics is numerically unstable for the hyperbolic partial differential Eqs. (7) and (11). Classically, stabilization is achieved by adding in the equations a minor amount of artificial diffusion.

The general 1D advection-diffusion equation for the arbitrary variable $S$ reads:

$$\frac{\partial S}{\partial t} + u \frac{\partial S}{\partial x} = c \frac{\partial^2 S}{\partial x^2} + F \qquad (C1)$$

$u$ and $c$ are the velocity and diffusivity, respectively. $F$ is a source term. It has been mathematically proven that the Galerkin method is numerically stable if the cell Peclet number $Pe$ is lower than a critical value $Pe_c$ [46]:

$$Pe = \frac{|u|\bar{h}}{c} < Pe_c \qquad (C2)$$

where $h$ is the mesh element size and $Pe_c = 2$. Obviously, Eqs. (7) and (11) cannot satisfy Eq. (C2) since, in these equations, $c = 0$ results in $Pe \to \infty$.

This problem is solved by adding in Eqs. (7) and (11) an artificial diffusion term (equivalent to the first term in the RHS of Eq. C1), with a variable diffusivity set to satisfy the condition in Eq. (C2):

$$c = \max \left( c_{min} \cdot \frac{|u| \bar{h}}{Pe_c} \right) \qquad (C3)$$

where $c_{min}$ is a lower bound of the virtual diffusivity, set to the air thermal diffusivity in Eq. (7) and the water vapor diffusivity in Eq. (11).

Eq. (C3) insures the numerical stability of the Galerkin method, but we must check that the artificial diffusion introduced in the equations does not significantly affect the results. The virtual diffusivity $c$ is a decreasing function of the mesh size. The final test for the regularization validity consists in reducing the mesh size until it has no significant effect on the results.

## Appendix D. Impact of the air flowrate uncertainty on the model outputs

Figs. D1a to D1c display the impact of the air flowrate uncertainty on the conduit wall and air temperatures obtained from HA and DA models. We assume that the actual air flowrate is in the range from the measured value to this value multiplied by 1.5.
<figure><img src="image_20.png" /></figure>

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

<figure><img src="image_22.png" /><figcaption>**Fig. D1.** Impact of the air flowrate uncertainty on the model outputs, (a) AMT, (b) ATF, (c) DTF. The solid and dash-dotted lines delineate the uncertainty ranges for the conduit wall and the air temperatures, respectively.</figcaption></figure>
# Data availability

We have attached the link to the software file used in this study.

# References

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

[2] L. Kukuljan, F. Gabrovšek, M.D. Covington, V.E. Johnston, CO2 dynamics and heterogeneity in a cave atmosphere: role of ventilation patterns and airflow pathways, Theor. Appl. Climatol. 146 (2021) 91–109.

[3] A.J. Kowalczak, P.N. Froelich, Cave air ventilation and CO2 outgassing by radon-222 modeling: how fast do caves breathe? Earth Planet Sci. Lett. 289 (2010) 209–219.

[4] M.O. Cuthbert, G.C. Rau, M.S. Andersen, H. Roshan, H. Rutlidge, C.E. Marjo, M. Markowska, C.N. Jex, P.W. Graham, G. Mariethoz, Evaporative cooling of speleothem drip water, Sci. Rep. 4 (2014) 5162.

[5] G.C. Rau, M.O. Cuthbert, M.S. Andersen, A. Baker, H. Rutlidge, M. Markowska, H. Roshan, C.E. Marjo, P.W. Graham, R.I. Acworth, Controls on cave drip water temperature and implications for speleothem-based paleoclimate reconstructions, Quat. Sci. Rev. 127 (2015) 19–36.

[6] D.C. Culver, T. Pipan, The Biology of Caves and Other Subterranean Habitats, Oxford University Press, 2019.

[7] L.E.O. Braack, Arthropod inhabitants of a tropical cave ‘island’ environment provisioned by bats, Biol. Conserv. 48 (1989) 77–84.

[8] A. Castano-Sánchez, G.C. Hose, A.S.P.S. Reboleira, Ecotoxicological effects of anthropogenic stressors in subterranean organisms: a review, Chemosphere 244 (2020) 125422.

[9] V. Rizzo, D. Sánchez-Fernández, J. Fresneda, A. Cieslak, I. Ribera, Lack of evolutionary adjustment to ambient temperature in highly specialized cave beetles, BMC Evol. Biol. 15 (2015) 1–9.

[10] A. Gomell, D. Austin, M. Ohms, A. Pflitsch, Air pressure propagation through Wind Cave and Jewel Cave: how do pressure waves travel through barometric caves? Int. J. Speleol. 50 (2021) 263–273, https://doi.org/10.5038/1827-806x.50.3.2393.

[11] A. Borsato, M. Samadelli, V. Martinucci, G. Manzi, Temperature fluctuations and ventilation dynamics induced by atmospheric pressure variations in Lamalunga Cave (Apulia, Italy) and their influences on speleothem growth, Quat. Res. 118 (2024) 100–115.

[12] L. Kukuljan, F. Gabrovsek, M. Covington, The relative importance of wind-driven and chimney effect cave ventilation: observations in Postojna Cave (Slovenia), Int. J. Speleol. 50 (2021) 275–288, https://doi.org/10.5038/1827-806x.50.3.2392.

[13] F. Gabrovsek, How do caves breathe: the airflow patterns in karst underground, PLoS One 18 (2023) e0283767.

[14] T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. DeWitt, Fundamentals of Heat and Mass Transfer, Wiley, 2017.

[15] A. Sedaghatkish, C. Pastore, F. Doumenc, P. Jeannin, M. Luetscher, Modeling heat transfer for assessing the convection length in ventilated caves, J. Geophys. Res. Earth Surf. 129 (2024) e2024JF007646.

[16] J.B. Cropley, Influence of surface conditions on temperatures in large cave systems, Bull. Natl. Speleol. Soc. 27 (1965) 1–10.

[17] W.M. Rohsenow, J.P. Hartnett, Y.I. Cho, Handbook of Heat Transfer, Mcgraw-hill, New York, 1998.

[18] F. Salmon, D. Lacanette, H. Lharti, C. Siriex, Heat transfer in rock masses: application to the lascaux cave (France), Int. J. Heat Mass Transf. 207 (2023) 124029.

[19] C.R. De Freitas, A. Schmekal, Prediction of condensation in caves, Speleogenes. Evol. Karst Aquifers 3 (2005).

[20] A. Fernandez-Cortes, D. Benavente, S. Cuezva, J.C. Cañaveras, M. Álvarez-Gallego, E. Garcia-Anton, V. Soler, S. Sanchez-Moral, Effect of water vapour condensation on the radon content in subsurface air in a hypogean inactive-volcanic environment in Galdar cave, Spain, Atmos. Environ. 75 (2013) 15–23.

[21] F. Gázquez, L. Quindós, D. Rábago, I. Fuente, S. Celaya, C. Sainz, The role of cave ventilation in the triple oxygen and hydrogen isotope composition of condensation waters in Altamira Cave, northern Spain, J. Hydrol. 606 (2022) 127416.

[22] P.-Y. Jeannin, La Baume de Longeagie, 30 ans de silence, 2018. Cavernes.

[23] C. Pastore, A. Sedaghatkish, E. Weber, N. Schmid, P.-Y. Jeannin, M. Luetscher, Monitoring air fluxes in caves using digital flow metres, Int. J. Speleol. 53 (7) (2024).

[24] B. Qaddah, L. Soucasse, F. Doumenc, S. Mergui, P. Rivière, A. Soufiani, Coupled heat and mass transfer in shallow caves: interactions between turbulent convection, gas radiative transfer and moisture transport, Int. J. Therm. Sci. 194 (2023) 108556.

[25] E.L. Cussler, E.L. Cussler, Diffusion: Mass Transfer in Fluid Systems, Cambridge university press, 2009.

[26] D.R. Lide, CRC Handbook of Chemistry and Physics, CRC press, 2004.

[27] M.D. Covington, A.J. Luhmann, F. Gabrovšek, M.O. Saar, C.M. Wicks, Mechanisms of heat exchange between water and rock in karst conduits, Water Resour. Res. 47 (2011).

[28] A.J. Luhmann, M.D. Covington, S.C. Alexander, S.Y. Chai, B.F. Schwartz, J. T. Groten, E.C. Alexander Jr, Comparing conservative and nonconservative tracers in karst and using them to estimate flow path geometry, J. Hydrol. 448 (2012) 201–211.

[29] S. Birk, R. Liedl, M. Sauter, Karst spring responses examined by process-based modeling, Groundwater 44 (2006) 832–836.

[30] P. Molnar, Differences between soil and air temperatures: implications for geological reconstructions of past climate, Geosphere 18 (2022) 800–824.

[31] B. Larwa, Heat transfer model to predict temperature distribution in the ground, Energies 12 (25) (2018).

[32] H. Lundström, M. Mattsson, Radiation influence on indoor air temperature sensors: experimental evaluation of measurement errors and improvement methods, Exp. Therm. Fluid Sci. 115 (2020) 110082.

[33] Y. Lyu, W. Luo, G. Zeng, Y. Wang, J. Chen, S. Wang, The effect of cave ventilation on carbon and oxygen isotopic fractionation between calcite and drip water, Sci. Total Environ. 915 (2024) 169967.

[34] I.J. Fairchild, A. Baker, Speleothem Science: from Process to Past Environments, John Wiley & Sons, 2012.

[35] D.M. Tremaine, P.N. Froelich, Y. Wang, Speleothem calcite farmed in situ: modern calibration of 618O and δ13C paleoclimate proxies in a continuously-monitored natural cave system, Geochim. Cosmochim. Acta 75 (2011) 4929–4950.

[36] J.A. Wassenburg, S. Riechelmann, A. Schröder-Ritzrau, D.F.C. Riechelmann, D. K. Richter, A. Immenhauser, M. Terente, S. Constantin, A. Hachenberg, M. Hansen, Calcite Mg and Sr partition coefficients in cave environments: implications for interpreting prior calcite precipitation in speleothems, Geochim. Cosmochim. Acta 269 (2020) 581–596.

[37] W. Dreybrodt, M. Deininger, The impact of evaporation to the isotope composition of DIC in calcite precipitating water films in equilibrium and kinetic fractionation models, Geochim. Cosmochim. Acta 125 (2014) 433–439.

[38] M.J. Medina, D. Antić, P.A.V. Borges, S. Borko, C. Fiser, S.-E. Lauritzen, J.L. Martín, P. Oromí, M. Pavlek, E. Premate, Temperature variation in caves and its significance for subterranean ecosystems, Sci. Rep. 13 (2023) 20735.

[39] T.C. Jegla, T.L. Poulson, Circadian rhythms-I. Reproduction in the cave crayfish, Orconectes pellucidus inermis, Comp. Biochem. Physiol. 33 (1970) 347–355.

[40] Y.-S. Yu, J.-H. Roh, J. Kim, A study on thermodynamic natural ventilation analysis by the field survey of underground mines in Korea, Tunn. Undergr. Sp. 23 (2013) 288–296.

[41] M.T. Parra, J.M. Villafruela, F. Castro, C. Mendez, Numerical and experimental analysis of different ventilation systems in deep mines, Build. Environ. 41 (2006) 87–93.

[42] Y. Lv, Y. Jiang, W. Hu, M. Cao, Y. Mao, A review of the effects of tunnel excavation on the hydrology, ecology, and environment in karst areas: current status, challenges, and perspectives, J. Hydrol. 586 (2020) 124891.

[43] S.C. Dornberger, A.F.R. Loria, M. Zhang, L. Bu, J.-L. Epard, P. Turberg, Heat exchange potential of energy tunnels for different internal airflow characteristics, Geomech. Energy Environ. 30 (2022) 100229.

[44] V. Gnielinski, Neue Gleichungen für den Wärme- und den Stoffübergang in turbulent durchströmten Rohren und Kanälen, Forsch. im Ingenieurwes. A 41 (1975) 8–16.

[45] S.E. Haaland, Simple and Explicit Formulas for the Friction Factor in Turbulent Pipe Flow, 1983.

[46] D. Kuzmin, J. Hämäläinen, Finite element methods for computational fluid dynamics: a practical guide, SIAM Rev. 57 (2015) 642.

[47] A. Sedaghatkish, Thermal Modelling of Caves Ventilated by Chimney Effect, Zenodo, 2024, https://doi.org/10.5281/zenodo.14535340 [software].

[48] Federal Office for Meteorology and Climatology of Switzerland. www.meteoswiss.admin.ch. (Accessed 1 March 2024).

[49] R. Tiwari, V. Boháč, P. Dieska, G. Götzl, Thermal properties of limestone rock by pulse transient technique using slab model accounting the heat transfer coefficient and heat capacity of heat source, AIP Conf. Proc. 2305 (2020) 020020. https://doi.org/10.1063/5.0033924.

[50] J.A. Riddick, W.B. Bunger, T.K. Sakano, Organic Solvents: Physical Properties and Methods of Purification, 1986.