Modeling Heat Transfer for Assessing the Convection Length in Ventilated Caves

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


Introduction
Understanding heat transfer in karst systems is a key issue for underground biota (Mammola et al., 2019), preservation of cave art (Bourges et al., 2014), speleothem growth rates (Banner et al., 2007;Spötl et al., 2005), or paleoclimate reconstruction (Borsato et al., 2016;Casteel & Banner, 2015;Domínguez-Villar et al., 2021).However, heat transfer in karst results from an intricate coupling between several mechanisms, including heat conduction in the rock mass (Quindos et al., 1987), convection due to water or air flow in caves (Cropley, 1965), or radiative transfer between cave walls (Guerrier et al., 2019).Latent heat exchanges are also present because of evaporation and condensation on cave walls (Dreybrodt et al., 2005) or ice formation and permafrost (Luetscher et al., 2008).A key issue for the understanding of heat transfer in karst is to determine the processes that dominate in a given configuration, and those that can be neglected.
The simplest heat transfer model of a vadose karst assumes 1D conduction heat transfer in a rock mass of infinite depth with periodic temperature fluctuations at the ground surface (Domínguez-Villar et al., 2023;Salmon et al., 2023;Villar et al., 1983).This model predicts that the depth corresponding to 99% attenuation of the • The distance of propagation (the convection length) is approximately proportional to the airflow rate and inversely proportional to the square root of the cave diameter

Supporting Information:
Supporting Information may be found in the online version of this article.temperature fluctuations at the ground surface is a few tens of centimeters for the daily fluctuations and approximately 10 m for the annual fluctuations.Quindos and co-workers (Quindos et al., 1987;Villar et al., 1983) measured the monthly mean temperatures in Altamira Cave during a year.They found that the periodic 1D conduction model correctly predicts the amplitude and the phase shift of the annual fluctuations at the cave roof whose depth varied from 3.5 to 17.5 m.In contrast, the floor temperature was much closer to the roof temperature than expected from a 1D conduction model.The authors attributed this effect to the homogenization of the temperature field by radiative transfer between roofs and floors.Guerrier et al. (2019) and Qaddah et al. (2023) confirmed by numerical simulations the significance of radiative transfer inside weakly ventilated shallow caves, and showed that turbulent free convection also had to be considered.Both heat transfer mechanisms contribute to the homogenization of the wall temperature and the deformation of the temperature field in the rock matrix around the cave.
The situation is different when forced convection is present.It has long been known that, in large cave systems, the convective transport of heat by air currents and streams can propagate thermal perturbations over much larger distances than predicted by a pure conduction model.Cropley (1965) found that temperature measurements performed in two caves in West Virginia supported the definition of three zones: zone I immediately adjacent to the cave entrance, where the cave temperature follows the surface temperature fluctuations; zone II where the cave temperature is driven by streams and air currents; zone III where the temperature is approximately constant and close to the annual mean temperature (AMT) at the surface.Areas where the amplitude of temperature fluctuations was lower than 1°F (ca.0.56°C) around the AMT were located at more than 1,500 m from the cave entrances.This long distance was mainly due to cold water carried into the caves by winter water recharge.Luetscher and Jeannin (2004) defined the heterothermic vadose zone as the surficial part of a karst massif showing seasonal inversions of the temperature gradient, in contrast with the homothermic vadose zone characterized by a temperature gradient close to the lapse rate of humid air (ca.0.5°C 100 m 1 ).These authors estimated the depth of the heterothermic zone at about 100 m in most cases.
Considerable work has been undertaken to simulate heat transfer by forced convection due to water flow in the vadose and phreatic zones of karst (Gong et al., 2019;Long & Gilcrease, 2009;Sinokrot & Stefan, 1993).Covington et al. (2011) showed that the relative significance of different heat fluxes including convection and conduction are timescale dependent.Conduction through the rock surrounding a conduit determines heat flux at times of the order of weeks and longer.Coupling convection in water with conduction in rock is thus necessary to get a realistic model valid at all time scales.Forced convection due to airflow must be considered when a network of karst conduits has at least two entrances.In this configuration, several mechanisms can contribute to the production of airflow through the cave, resulting in significant convective heat transfer between the cave and the external atmosphere.Potential driving forces are barometric fluctuations (Gomell et al., 2021), dynamic pressure effect due to external winds (Kukuljan et al., 2021), diphasic flow due to water circulation (Atkinson, 1977), or the buoyancy due to the density contrast between the air inside and outside of the cave (Gabrovšek, 2023).The latter mechanism is likely the most significant in most cases.The air density depends on moisture content, CO 2 concentration and temperature.In temperate climates, temperature fluctuations are the main cause of density variations (Gabrovšek, 2023).Therefore, the flow direction is upward when the temperature inside the cave is warmer (or downward when inside the cave is colder) than the external atmosphere.Because of the thermal inertia of the massif, a lower entrance operates as an inlet of cold external air during most of the winter, whereas it acts as an outlet during most of the summer.In the latter case, the air cooled down when crossing the massif before reaching the lower entrance.The lower entrance thus receives in summer an airflow colder than the external temperature.Because of this seasonal asymmetry, the AMT is shifted to colder values at lower entrances (cold thermal anomaly).The same mechanism produces a shift of the AMT to hotter values at upper entrances (hot thermal anomaly) (Lismonde, 2002).In alpine ice caves, where this ventilation pattern represents the normal unless the cave is clogged with sediments, these thermal anomalies usually extend over a few hundred meters from the cave entrances (Luetscher et al., 2008).Wigley and Brown (1971) calculated the temperature and moisture content profiles in the air as a function of the distance from the entrance.They used a 1D model based on the energy and water mass conservation in the airflow.The convective heat transfer coefficient between the air and the wall was estimated using an empirical correlation valid for forced convection in smooth pipes.The model assumes prescribed uniform wall temperature, a major simplification that allows to obtain simple closed-form expressions for the temperature and moisture profiles.However, this approach implicitly assumes that the heat flux through the cave wall is convection-limited, and that the rock mass does not play any role.Lismonde (2002) pointed out that the airflow modifies the rock temperature, which in turn changes the air temperature profile in the cave.This author developed a model to predict heat transfer in a straight inclined duct of constant diameter included in a rock massif.The model coupled 1D radial conduction in the rock mass with convective transfer in the gas, using a sinusoidal function of time as inlet temperature.The airflow rate was predicted by considering the interaction between the temperature field in the gas and the airflow through the buoyancy term in the momentum balance equation of the gas.This model reproduced qualitatively some field observations, for example, the hysteresis of the flow rate or the thermal anomaly in the entrance areas.Gabrovšek (2023) used a similar model to investigate by numerical simulations the effect of conduit shape on the airflow pattern in ventilated caves.
The present study focuses on the numerical simulation of heat transfer in ventilated caves for which the airflow is driven by the temperature contrast between the cave and the external atmosphere.Our objective is to go beyond the qualitative results obtained from previous numerical studies available in the literature.We aim at assessing orders of magnitudes of the spatial extent of the karst region where the airflow induces significant thermal perturbations.We use for that a numerical model close to that already developed by Lismonde (2002) or Gabrovšek (2023), and we apply it to the simplest possible configuration in order to quantify general trends common to most ventilated caves.The convective heat transfer due to dry airflow in a single horizontal karst conduit of constant diameter is coupled with the conductive heat transfer in the impermeable rock mass.This article contains a first theoretical part (Sections 2-5) dedicated to numerical simulations based on a model including many simplifying assumptions.The relevance of these assumptions to achieve our objective, which consists in the assessment of general orders of magnitude, is evaluated in a second part (Section 6) by comparison with the field data obtained from three sites, a mine tunnel and two caves.

Cave Geometry and Computational Domain
We consider the idealized ventilated cave displayed in Figure 1, located in a rock massif of infinite extent in both vertical directions.The cave consists of a vertical conduit connected by two horizontal conduits to two entrances at different elevations.The rock mass is assumed impermeable.The horizontal conduits are supposed long enough for the vertical conduit to be fully included in the homothermic area where the temperature gradient reduces to the adiabatic lapse rate.The aeraulic and thermal problems are uncoupled in this simplified configuration, since the temperature field in the horizontal conduits have no effect on buoyancy.The air flowrate thus only depends on the temperature contrast between the temperature in the vertical conduit, independent of time, and the temperature of the atmosphere outside the cave, a known function of time.
Another source of simplification is that the thermal problems in the regions of the upper and lower entrances are uncoupled.They can thus be treated separately.We arbitrarily focus on the upper entrance, but all the results can be easily translated to a lower entrance (the only difference is that the air flows through the upper entrance inward during summer and outward during winter, and vice-versa through the lower entrance).
The computational domain is a cylinder of length L dom and outer radius R dom displayed in red in Figure 1.It contains a conduit of same axis and length and constant radius R ≪ R dom , inside which the air circulates.The length L dom and outer radius R dom are set so that the computational domain includes the whole thermal perturbation induced in the rock mass by the airflow.Practically, in a given configuration, L dom and R dom are increased until they no longer influence the results.

Atmosphere Temperature and Mass Flow Rate
For the sake of simplification, only the annual fluctuation of the atmosphere temperature T atm is considered.It is thus assumed that T atm (t) follows the simple law: where T m and ΔT y are the AMT and the amplitude of the annual temperature fluctuation (ATF), respectively, in the external atmosphere at the elevation of the upper entrance.The corresponding period is τ y = 1 year.
The air mass flow rate ṁ is deduced from the momentum balance applied to the air inside the cave.Assuming negligible inertia, the balance between friction and buoyancy yields (Lismonde, 2002): where H is the length of the vertical conduit and g is the gravitational acceleration.The cave aeraulic resistance K can be assumed to be constant in the turbulent regime.ρ atm and ρ m are the densities of the air outside the cave and inside the homothermic zone, respectively, averaged over the cave height.Assuming that: (a) the difference of densities depends only on the temperature drop between the external atmosphere and the homothermic zone, and (b) the temperature field in the vertical conduit follows the AMT in the external atmosphere, we get after linearization: where P atm is the atmospheric pressure, M a the molar mass of air and R g the ideal gas constant.Equations 2 and 3 show that the air flowrate is proportional to the square root of the temperature drop (T atm T m ).Injecting Equations 1 and 3 in Equation 2 yields where S = 1 for T atm < T m (wintertime, upward flow) and S = 1 for T atm > T m (summertime, downward flow).The positive constant Δ ṁ is the amplitude of the annual fluctuation of the flow rate.It is an input of the model.

Governing Equations
The air temperature verifies the energy balance where c p,a is the air heat capacity at constant pressure, T a is the mixing temperature of the air inside the conduit, x the distance from the external rock surface, P the conduit perimeter (P = 2πR for a circular cross-section) and φ w the conductive flux at the conduit wall, positive when going from the rock to the air.Equation 5 is a balance between the energy advected by the air flow (LHS) and the conduction flux at the conduit wall (RHS).Air thermal inertia has been neglected compared to advection (quasi-steady approximation).Equation 5 requires a single boundary condition at the conduit inlet located at x = 0 for downward flow (S = 1) or at x = L dom for upward flow (S = 1).We thus impose the boundary condition The conductive flux φ w is required to solve Equation 5. Therefore, the conduction equation must be solved in the impermeable rock matrix around the cave.The transient 2D axisymmetric conduction equation reads where T r (x, r, t) is the rock temperature, α r the rock thermal diffusivity, and r the distance from the cave axis.The boundary conditions are as follows.The atmosphere temperature is imposed on the external rock surface: The boundary at r = R dom can be assumed adiabatic because the temperature field only depends on x far from the cave (i.e., for a large value of R dom ), and the radial component of the temperature gradient is thus close to zero.
The boundary condition at x = L dom is also adiabatic because this boundary is located in the homothermic zone.We get The last boundary condition is at the conduit wall.It is given by Newton's law of cooling and the heat flux continuity at the conduit wall: where λ r is the rock thermal conductivity, T r (x, R, t) = T w (x, t) is the temperature of the conduit wall and h th the heat transfer coefficient.The latter is time-dependent since it depends on the air flowrate ṁ(t) .All the conditions at the boundaries of the rock domain are displayed in Figure 2. No initial condition is required because we are looking for the periodic regime, that is, the solution asymptotically reached by the model at infinite time.

Determination of the Heat Transfer Coefficient
The heat transfer coefficient h th in Equation 10 is a key parameter.The estimation of h th must therefore be performed with great care to get a reliable model.Assuming forced convection in conduits, Wigley and Brown (1971) and Gabrovšek (2023) used a correlation close to the Colburn correlation (Bergman et al., 2017), valid for fully developed turbulent flow in smooth pipes.Lismonde (2002) pointed out that a cave cannot be considered as a smooth pipe, and multiplied by two the numerical prefactor of the Colburn correlation to account for the effect of wall roughness.Covington et al. (2011) used the Gnielinski correlation (Bergman et al., 2017), which explicitly considers the effect of wall roughness.
However, it is unlikely that a cave can be considered as a rough pipe.Indeed, caves can have complicated shapes, ranging from sub-circular conduits with several meters in diameter to complex cave passages associated with the collapse of the host-rock.A succession of bends, conduit contractions or enlargements, or obstacles of any kind are expected to increase the transfer of heat between the air and the wall.Indeed, all of these singularities not only increase the heat exchange surface between the wall and the fluid but also enhance the heat transfer coefficient by generating secondary flows and by increasing the turbulence level.It is well known that heat exchanger performance can be improved by the insertion of twisted tapes, longitudinal fins, or helical ribs inside the tubes (Bergman et al., 2017).For instance, spirally corrugated tubes can increase the heat transfer coefficient by a factor 3 (Pethkool et al., 2011;Promthaisong et al., 2016).However, although empirical correlations exist to predict the heat transfer coefficient for a wide range of well-defined geometries used in heat exchangers, no correlation is available for the irregular and tortuous geometries commonly encountered in caves.
To get around this difficulty, two distinct cases will be considered.In the first case (case A in the following), correlations for forced convection and fully developed flow in pipes will be used to estimate the heat transfer coefficient.This approach is expected to give a lower bound of the heat transfer coefficient, and thus a lower bound of the heat flux at the conduit wall.In the second case (case B in the following), a higher bound of the wall heat flux will be obtained assuming an infinite heat transfer coefficient.Reality must lie between these two limiting cases.
Case A: Forced convection and fully developed flow The Nusselt number Nu = h th D h k a used for the estimation of h th is displayed in Figure 3 as a function of the Reynolds number Re t = V(t) D h ν a , for three values of the relative roughness ε (ratio of the wall roughness over D h ).k a and ν a are the conductivity and the kinematic viscosity of the air, V is the mean air velocity at time t, D h = 4 A/P is the hydraulic diameter with A and P the cross-sectional area and the perimeter of the conduit.Nu is computed as follows: where Nu L = 3.66 is the Nusselt number in the laminar regime (independent of the roughness) and Nu T is given by the Gnielinski correlation (Bergman et al., 2017), which takes into account the effect of the wall roughness expected in the turbulent regime: where Pr is the air Prandtl number (Pr = 0.71).f d is the Darcy friction factor which depends on the Reynolds number and the wall relative roughness ε. f d was estimated using the Haaland correlation (Haaland, 1983):  1 All the simulations of case A were done using ε = 0.01 as a lower bound of the relative roughness in a cave (red curve in Figure 3).

Case B: Infinite heat transfer coefficient
In this case, the conduit wall temperature is equal to the bulk air temperature and Equation 10 thus reduces to x,R,t) and T r (x,R,t) = T a (x,t) (14)

Physical Properties
All the physical properties are gathered in Table 1.They are assumed constant and estimated at the temperature T m = 12°C.

Dimensionless Equations
We define the dimensionless temperatures θ a = (T a T m )/ΔT y and θ r = (T r T m )/ΔT y in the air and the rock, respectively.The dimensionless temperature 0 thus corresponds to the external AMT and 1 is the amplitude of the external annual fluctuations.We take the period τ y = 1 year as the unit of time and the diffusion length ≈ 5.3 m as the unit of length.The dimensionless atmospheric temperature θ atm and air flowrate μ are deduced from Equations 1 and 4: The dimensionless counterparts of Equations 5-10 read where the tilde (∼) indicates dimensionless time or length.η t) = h th (t) H th is the heat transfer coefficient scaled by H th , its maximum value reached for ṁ = Δ ṁ.The Reynolds number Re = 4 Δ m Pμ a (with μ a the air dynamic viscosity) and the Biot number Bi = are based on the maximum values Δ ṁ and H th , respectively.μ t) Re and η t) Bi are the Reynolds and Biot numbers at a given time t.In a given configuration, Rdom and Ldom that define the size of the computational domain are increased until they do not modify the results, as stated in Section 2.1.
With the constant physical properties defined in Table 1, a configuration is well-defined if only two dimensionless numbers are known: the dimensionless conduit radius R and the Reynolds number Re.In case A, that is, when fully developed forced convection is assumed for the estimation of the heat transfer coefficient, the Biot number in Equation 22is related to the Nusselt number by the simple relation η , where Nu depends on Re t = μ t) Re through Equations 11 and 12.In case B, the Biot number is infinite, that is, η t) Bi→ ∞.Equation 22 thus reduces to

Method of Resolution
We are looking for the periodic solution of the mathematical model defined in Section 2.6.The periodic regime could be obtained by time integration starting from an arbitrary initial condition, after simulating a number of cycles large enough to get a good approximation of the solution at infinite time.However, the slow convergence of this method can generate large computational times.We thus prefer to use another method based on the Fourier series.All time-varying variables, including the temperatures in the rock and in the air, are approximated by truncated Fourier series and inserted in the mathematical model of Section 2.6.This yields a set of coupled partial differential equations (PDE), which does not include the time.The numerical resolution by the Galerkin method of this set of PDE gives the spatial distribution of the amplitude and phase of the Fourier modes.The temporal evolution of the dependent variables is then recovered from Fourier series.This method is detailed in Supporting Information S1.The validation by comparison with a standard time integration method is presented in the same appendix.All the numerical simulations were performed with the commercial software Comsol Multiphysics, version 6.1.
The method of resolution based on the Fourier series used in this work is particularly relevant for the simplified time variation of the external temperature considered in Equation 1. Considering all the Fourier modes included in the real time series of the external temperature would clearly increase both the complexity of the method implementation and the computational resources required for the simulations.In that case, only the most significant modes of the input temperature should be considered.In contrast, using Fourier series to simulate the complex geometry of real caves would not cause any particular problem.

Characterization of the Thermal Perturbation
Figure 4 displays the wall and air temperatures inside the cave at different distances from the entrance, for case A with R = 0.189 and Re = 1.8 × 10 5 .Figure 4a shows the amplitude of the Fourier modes that stems directly from the numerical simulation.The mode k refers to the dimensionless period 1/k.The amplitudes of the modes k = 0 and k = 1 are thus the AMT and the amplitude of the ATF, respectively.The modes k = 2, k = 3, and so on refer to the periods 6 months, 4 months, etc.Although the external temperature imposed at x = 0 contains the single mode k = 1 (since it reduces to Equation 15), the coupling between modes due to convection results in the existence of other modes than k = 1 in cave temperatures (modes k = 0 to 18 were considered in the simulations, see Supporting Information S1). Figure 4b displays the periodic time series built from the Fourier series.The amplitude of the temperature fluctuations decreases when the distance from the entrance increases, as expected.A hot thermal anomaly expected at the vicinity of an upper entrance is also observed.Indeed, for x = 10 to 10 3 , the dimensionless temperature is positive throughout the year, which means that it is always higher than the external AMT, even in winter.
Figure 5a displays the temperature profiles in the air and at the cave wall at times t = 0, 0.25, 0.5, and 0.75, corresponding to the external temperatures 0, 1, 0 and 1, respectively.Following an approach similar to that of Cropley (1965), we split the cave into three zones.
1.A first zone, from x = 0 to x ≈ 1 (the diffusion length), where the cave temperature follows the fluctuations of the external temperature.More precisely, the sign of the temperature gradient changes in summer and winter (see the inset in Figure 5a).This first zone is dominated by the temperature fluctuations of the external wall, propagating inside the rock mass by heat conduction.It will be called the diffusive region in the following.2. A second zone characterized by a positive dimensionless temperature and reduced but significant fluctuations.
In the specific case of Figure 5a, it approximately extends from x ≈ 1 to 10 3 .In this zone, heat diffusion from the external wall is negligible compared to convection.It will be called convective region in the following.
Because of the seasonal flow reversals, the inlet temperature of the air flow is 1 in summer (inward flow) and 0 in winter (outward flow, coming from the deep karst).As a consequence, this zone does not "feel" the negative winter temperatures.This results in: (a) a shift of the AMT to positive values, the so-called thermal anomaly, (b) a reduction of the amplitude of the fluctuations, because the temperature oscillates between 0 and 1 whereas the external temperature fluctuates between 1 and 1. 3. A third zone, called deep karst in the following, where the dimensionless temperature is approximately constant and close to zero (the external AMT).
These three regions appear clearly on the amplitude of the Fourier modes plotted in Figure 5b as a function of the distance from the entrance x.The amplitude of the mode k = 1 (i.e., the ATF) decreases in two steps.It first drops down from 1 to 0.3 in the diffusive region (from x = 0 to x ≈ 1), then slowly decreases in the convective region, over a distance much longer than 1 (of the order of 10 3 in the specific case of Figure 5).The amplitude of the mode k = 0 (i.e., the AMT) increases in the diffusive region until it reaches a quasi-plateau.Then it slowly decreases in the convective region, over a distance of the same order as for the mode k = 1.The modes such that k ≳ 2 follow qualitatively the behavior of mode k = 0, with a lower amplitude and a shorter attenuation length.In the following, we will focus on both modes with the highest amplitude, k = 0 and k = 1, as representative of the AMT (the thermal anomaly) and the ATF, respectively.
In all the simulations we did, the transition between the diffusive and convective regions always take place at a distance from the entrance of the order of 1, as expected.Conversely, the characteristics of the convective domain (amplitude of the thermal perturbation and its spatial extent) depend on the model parameters: the Reynolds number Re (containing the air flow rate) and the cave radius R.This dependence will systematically investigated in the parametric study of the next section.In order to facilitate the interpretation of the results, we define a small number of relevant parameters characterizing the perturbation of the cave temperature field by the air flow.The magnitude of the perturbation is characterized by the AMT and the ATF at the beginning of the convective region, where the impact of convection is maximum, and the effect of thermal conduction from the external wall is negligible.We define Δθ w as the maximum amplitude of the mode k = 1 in the convective region, for the wall temperature.Practically, it is evaluated at the change of slope noticed on the solid blue line in the inset of Figure 5b.Similarly, θ w is the maximum amplitude reached by the mode k = 0 for the wall temperature (maximum of the red solid line in the inset of Figure 5b).The corresponding amplitudes for the air temperature, Δθ a and θ a , are estimated at the same location as Δθ w and θ w , respectively.For example, in the specific case of Figure 5, we get Δθ w ≈ 0.30 and Δθ a ≈ 0.34 (maximum ATF), θ w ≈ 0.425 and θ a ≈ 0.44 (maximum AMT).
The four parameters θ w , θ a , Δθ w , and Δθ a represent the maximum amplitudes of the perturbation in the convective region.These amplitudes decrease with the distance from the entrance, as the heat transported by the air flow is gradually dissipated by conduction in the rock mass.We arbitrarily assume that a Fourier mode has been damped when its amplitude has been divided by 10. Figure 6a shows the thermal perturbation in a rock massif with the black solid lines corresponding to an AMT equal to θ w / 10 = 0.0425.The rock domain thus defined can be characterized by the length L0 along the cave axis and the maximum width W0 along the radial direction.Note that this maximum width is observed at a certain distance from the external wall because of the competition between radial diffusion from the cave which transports the thermal anomaly through the massif, and axial diffusion from the external wall which reduces it.Similarly, the black solid lines in Figure 6b correspond to an ATF equal to Δθ w /10 = 0.030.The size of the corresponding domain is characterized by its length L1 and its width W1 .The latter is taken close to the entrance, but out of the diffusive region.The lengths L0 and L1 can also be estimated from the air temperature (length such that the amplitude of AMT and ATF in the air are equal to θ a / 10 and Δθ a /10, respectively).

Spatial Extent of the Thermal Perturbation Induced by the Air Flow
The characteristic lengths L0 and L1 for the damping of the AMT and ATF along the cave axis are displayed in Figure 8.The values of L0 computed in case B (infinite Biot number) support the simple relation obtained from curve fitting: In case A, estimating L0 from the air or the wall temperatures give the same results (the lengths required to get θ a or θ w divided by 10 are the same, even when θ a and θ w are different).At high Reynolds number, L0 estimated from case A follows Equation 24 derived for case B, as expected when heat transfer is conduction-limited.When Re is decreased, using the forced convection correlations of case A reduces the heat flux at the cave wall compared to the infinite Biot number of case B, leading to larger values of L0 than predicted by Equation 24.
L1 , the distance corresponding to the damping of the ATF, is nearly equal to L0 , with the exception of slight differences at low Reynolds number for R = 2 (see Figure 8b).The length of the convective region (called convection length in the following) can thus be indifferently determined from the damping of the AMT or ATF, based on the wall or the air temperatures.A higher bound of the convection length is provided by case A, while case B provides a lower bound.We can see in Figure 8 that the uncertainty on the convection length increases when the Reynolds number Re is decreased and when the cave diameter R is increased.
The width W0 corresponding to the damping of the AMT in radial directions inside the rock mass is displayed in Figure 9.It is approximately proportional to L0 and an order of magnitude smaller.24.Symbols: the same legend as in Figure 7.

Baulmes Mine (Switzerland)
The Baulmes mine (46°47′33″N, 6°31′35″E, 655 m a.s.l.) is located in western Switzerland.This mine comprises a network of 11 sub horizontal levels with a cumulated length of about 17 km (Figure 10).The difference in elevation between the uppermost and lowest entrances is 90 m, ensuring a strong chimney effect.The sub horizontal conduits at intermediate elevations between the lower and upper entrances, connected by shafts, lead to dead ends in the horizontal directions (for more information about Baulmes mine, see Bourret et al. (2007)).The temperature delivered by Mathod weather station (46°44′13″N, 6°34′4″E, 435 m a.s.l., at 7 km SW from the study site) is used as an estimation of the external atmosphere temperature.A correction has been applied to account for the difference of elevation with the mine entrance, using the vertical thermal gradient 5.2°C/km (mean value from the network of Swiss meteorological stations for the time range 1991-2020).The resulting external temperature is displayed in Figure 11 for 1 year from 1st July 2022 at 0h00 (this date is taken as time t = 0 for all the data from Baulmes mine).The sampling period was 1 hr.The Fast Fourier Transform (FFT) of the external atmosphere temperature yielded T m = 10.4°C(after correction) and ΔT y = 9.4°C for the external AMT and ATF, respectively.
The 360 m long tunnel #10 is connected horizontally to the lower entrance with no connection to other branches.This is the zone of interest of study.It is indicated by the yellow shading in Figure 10.This conduit has a typical horseshoe cross section.The geometrical characteristics of this conduit, including the mean perimeter and crosssectional area, their standard deviation as well as its total length are summarized in Table 2.We instrumented this tunnel with seven temperature probes at different distances from the entrance.They are listed in Table 3.The first six temperature probes are Reefnet (Sensus Ultra) with ±0.3°C accuracy.The last one located at 360 m from  entrance is HOBO (u22-001) with ±0.2°C accuracy.The temperature inside the conduit was recorded by these seven probes from 3rd February 2023 (t = 0.596 year) to 7th July 2023 (t = 1.017 year).This time range is displayed by the horizontal arrow in Figure 11.The sampling period was 1 hr.
The air mass flow rate was measured sporadically with a manual hot wire anemometer (Testo425 with an accuracy of 0.03 m/s ±5% of the measured value).The measurements were carried out at the four dates listed in Table 4 and displayed by vertical red bars in Figure 11.There are two dates in winter and two in summer.The mean flowrates in winter and summer are 2.4 and 5.0 kg/s, respectively.In the following, we consider Δ ṁ ∼ (2.4 + 5.0)/2 ∼ 4 kg/ s as a rough estimation of the annual amplitude of the air flowrate.With the mean conduit perimeter P = 9.7 m from Table 2 and the air viscosity μ a ≈ 1.77 × 10 5 Pa.s from Table 1, we get the order of magnitude of the Reynolds number Re = 4 Δ ṁ Pμ a ∼ 10 5 .The airflow is thus in the turbulent regime most of the time.The conduit radius reads R = P/(2π) ≈ 1.5 m that yields the scaled value R = R/ L dif ≈ 0.3.The conduit is thus close to the configuration displayed in Figures 4 and 5 corresponding to Re = 1.8 × 10 5 and R = 0.189.In this configuration, cases A and B yield approximately the same results.The heat transfer is thus limited by conduction in the rock rather than convection.Therefore, Equation 24, for which conduction-limited heat transfer is assumed, is expected to provide a fair estimation of the convection length, and not only a lower bound.We get L0 ≈ L1 ∼ 3 × 10 2 .The theory thus predicts a convection length significantly larger than the conduit The time series of the external atmosphere temperature and the temperatures inside the conduit have been plotted in Figure 12a.The dimensionless temperatures inside the conduit are always negative (i.e., lower than the external AMT), even in summer.This is the cold thermal anomaly expected at a lower entrance.
In the following, we assume that the temperature profiles θ w (x) at time t = 0.611 and θ s (x) at time t = 1.016 give an assessment of the coldest and hottest temperatures reached inside the conduit in winter and summer, respectively.These two times are displayed by vertical bars in Figure 12a, and the corresponding profiles are plotted in Figure 12b (the temperature at x = 0 is the external atmosphere temperature).In winter, the dimensionless temperature continuously increases from 1.6 at the inlet (x = 0) to 0.56 at the end of the conduit (x = 68).Spatial temperature variations are thus significant over the entire length of the conduit.This suggests that the convection length should be larger than the conduit length.In summer, the reverse flow direction yields a  different temperature profile.The temperature at the end of the conduit (x = 68) is 0.06, close to the external AMT.When getting closer to the entrance, the temperature decreases to 0.16 at the first probe (x = 1.9).The temperature then suddenly increases within the diffusive region, from 0.16 at the first probe (x = 1.9) to the positive value 0.78 at the entrance (x = 0) .
The difference Δθ = (θ s θ w )/2 between summer and winter profiles yields an assessment of the amplitude of the ATF in the conduit.Δθ drops from 1.2 to 0.58 between the entrance (x = 0) and the first probe (x = 1.9).This sudden decrease takes place in the diffusive region.Δθ then smoothly decreases from 0.58 (first probe at x = 1.9) to 0.25 at the end of the conduit (x = 68).Assuming that the first sensor yields the temperature at the beginning of the convective region, the amplitude of the temperature fluctuation is thus divided by approximately 2.3 in the convective domain of the conduit.This confirms that the conduit length is significantly shorter than the convection length, defined as the distance required to divide by 10 the ATF measured at the beginning of the convective region (see Section 4).

D7.1 Cave (Sieben Hengste, Switzerland)
D7.1 cave is located in the Sieben Hengste karst area (46°45′43″N, 7°49′51″E, 1,750 m a.s.l).It is part of a large ventilated cave system (>160 km of cumulated length) subject to chimney effect (see the map in Figure 13).Field observations reveal a seasonal airflow consistent most of the time with an upper entrance of the cave system (i.e., air inflow during the summer season).The ingress of warm external air during summer is marked by some condensation within the first 20 m of the meander, drained as a streamlet at the base of the meander.Otherwise, the cave passage remains remarkably dry.Unexpectedly, temperature and flow rate measurements show periods of time in winter during which the cave sporadically behaves as a lower entrance, with cold air inflow.The explanation of these unexpected flow reversals is still being searched.
External air temperatures displayed in Figure 14a are provided by an automatic weather station located at 1,850 m elevation (Swiss Federal Institute for Snow and Avalanche Research), approximately 2 km SE of our study site (sampling period: 1 hr).The external atmosphere temperature was corrected from the elevation as in Section 6.1.The FFT yielded T m = 5.0°C (after correction) and ΔT y = 7.1°C for the external AMT and ATF, respectively.
The ca. 100 m-long entrance meander of D7.1 Cave has a simple geometry with an average hydraulic diameter of ca. 1 m, despite local variations in the cross-sectional area and perimeter.The perimeter P ≈ 3 m and the scaled radius R ≈ 0.1 will be considered in the calculations.The meander was instrumented for temperature and airflow measurements between 25/08/2021 and 10/08/2022.Five temperature probes were installed along the main conduit at distances from the cave entrance listed in Table 5.Probes #1 to #4 were Reefnets Sensus Ultra with ±0.3°C accuracy (sampling period: 30 min).Probe #5 was a HOBO U22-001 with ±0.2°C accuracy (sampling period: 2 hr).It is located at the base of a 20 m shaft following the entrance meander.Time series from the five probes are displayed in Figure 14b.
The cave airflow was measured using a low pressure drop flowmeter (Sensirion SFM 3000 with accuracy of 5% and 20% when velocity is higher and lower than 0.3 m/s, respectively; 50 m from cave entrance) and validated with independent flow measurements (CO 2 gauging) on 10/08/2022.The FFT of the air flow yields the order of magnitude Δ ṁ ∼ 0.1 kg/ s rate (red dashed line in Figure 14a), which results in Re ∼ 7 × 10 3 .The cave is thus most of the time in the transition between laminar and turbulent regimes.Equation 24 yields the theoretical convection length L0 ≈ L1 ∼ 10.It is important to note that this value must be considered as a lower bound of the convection length.Indeed, Equation 24 stands for case B (infinite Biot number).At Re ∼ 7 × 10 3 , case A (finite Biot number) yields a convection length significantly larger than case B (see Figure 8).
The AMT inside D7.1 Cave is displayed in Figure 17a.It decreases to reach a negative quasi-plateau in the range from 0.26 to 0.23 in between probes #3 (x = 6.2), and #5 (x = 24).This is not the behavior expected for an upper entrance, for which the AMT is supposed to be positive in the convective region before decreasing to zero at some distance from the entrance.Sporadic flow reversals in winter can contribute to this negative  AMT (an instance of such an event can be observed around time t ≈ 0.33 in Figure 14).Moreover, the cave entrance is in a topographic low, resulting in local cold air trapping during winter.In addition, snow is likely to accumulate near the cave entrance.Both may induce local free convection cells in the entrance zone in winter resulting in the cooling of the cave below the external AMT.
If we assess L0 from the distance over which the AMT significantly varies inside the cave, it would range between x = 2.1 (probe #2) and x = 6.2 (probe #3).This estimation is confirmed by the variation of the ATF with the distance from the entrance (Figure 17b).The ATF at probe #1 (x = 0.38), located in the diffusive region, is close to 1. Then it decreases to 0.36 at probe #2 (x = 2.1), the first probe in the convective region, and drops down to 0.02 at probe #3 (x = 6.2).The ATF being reduced by a factor larger than 10 between probes #2 and #3, L1 should be in between these two probes.It is worthwhile to note that this estimation is consistent with the limit of condensation, which falls in between probes #2 and #3 (20 m corresponds to x ≈ 4).The convection length of D7.1 Cave is thus in the range from 2.1 to 6.2.

Bärenhöhle Cave (Austria)
The cave is located in Vorarlberg, Austria (47°22′17″N, 9°52′52″E, 901 m asl).The external atmosphere temperature is obtained from the weather station located at St. Gallen (47°25′31.7094″N,9°23′54.7008″E,776 m a.s.Chimney effect drives the airflow in the cave which behaves as a lower entrance (air inflow in winter, see Figure 15).The difference in elevation between the lower and upper entrances is ca.130 m for a conduit length along the main air passage approximating 240 m (Perkmann & Luetscher, 2013).The average hydraulic diameter is assessed to about 1 m, so that P ≈ 3 m and R ≈ 0.1.Air temperature was recorded at 1h interval using three HOBO (U22-001) temperature probes at distances from the entrance listed in Table 6.The corresponding time series are displayed in Figure 16b.
Air flow was recorded between 16/04/2013 and 16/04/2014 at the lower entrance using a thermo-anemometer calibrated to ±0.1 m/s.The FFT of the air flow rate (red dashed line in Figure 16a) yields the orders of magnitude Δ ṁ ∼ 0.6 kg/ s and Re ∼ 5 × 10 4 .The cave is thus in the turbulent regime most of the time.In this configuration ( R ≈ 0.1 and Re ∼ 5 × 10 4 ), cases A and B predict approximately the same theoretical convection length (see Figure 8).Equation 24 yields L0 ≈ L1 ∼ 70.
The AMT of Bärenhöhle Cave is displayed in Figure 17a.Negative values are consistent with a lower entrance.The lowest AMT ( 0.27) is measured at probe #2 (x = 16).It increases slightly between probe #2 and the last probe #3.The ATF displayed in Figure 17b decreases over the whole instrumented conduit, from 1 outside the cave to 0.11 at the last probe (x = 23) .These two pieces of information suggest that the convection length should be longer than 23, the cave length.

Discussion
The main results from the three investigated sites are summarized in Table 7.For two sites (Baulmes mine and Bärenhöhle Cave), we do not observe any contradiction between the model and the field data since both yield a convection length larger than the instrumented conduit (a hundred meters for Bärenhöhle, a few hundreds of meters for Baulmes).In D7.1 Cave, the field data suggest a shorter convection length, in the range from 10 to 30 m, whereas the theoretical prediction yields a lower bound of 50 m.Although the theoretical order of magnitude is correct, the theory overestimates the convection length.In a cave with a small Reynolds number and a short convection length, local effects in the entrance region, not considered in the model, could overcome the effect of the mean flow generated by the chimney effect.In the case of D7.1 Cave, cold air ingress in winter might cool down the cave passage more than considered in the model, thus producing an actual shorter thermal length.
Using a theoretical approach, we defined the characteristic lengths L0 and L1 from the damping in the convective region of the AMT (mode k = 0) and ATF (mode k = 1), respectively.L0 and L1 are close to each other, and similar results are obtained if the estimation of these lengths is based on the wall  temperature or that of the air (see Figure 8).A single convection length is thus sufficient to characterize the extent of the AMT and the ATF along the cave axis.The case B (infinite Biot number) provides a lower bound of the convection length through the simple Expression 24.If the Reynolds number is large enough (typically, if Re > 3 × 10 4 for a conduit with a radius of 1 m, see Figure 8), the approximation of infinite Biot number is valid and Equation 24 is expected to give a right estimation of the convection length.Turning this relation into unscaled variables yields where L 0 , L 1 , and R are in meters and Δ ṁ in kg/s.Equation 25gives the convection length resulting from the balance between heat advection by the air flow and heat conduction in the rock mass, in the conduction-limited regime.In this case, the convection length is approximately proportional to the amplitude of the flowrate annual fluctuations Δ ṁ divided by the square root of the cave radius R. Ventilated conduits with R of the order of 1 m traversed by an airflow with Δ ṁ ranging from 0.1 to 1 kg/s are commonly encountered in mines and caves (see the case of Section 6).For such conduits, Equation 25 yields a convection length along the conduit axis  between a few tens and a few hundreds of meters.Thermal perturbations also propagate in the radial direction, that is, along the directions perpendicular to the conduit.As expected, the amplitude of the annual fluctuations vanishes over the diffusion length, that is, a few meters.Conversely, the thermal anomaly (defined as the shift of the AMT), propagates along distances of the order of L 0 /10 (see Figure 9), that is, over distances of the order of a few meters to a few tens of meters.It is worthwhile to note that in some cases, the extent of the thermal perturbation might be much larger than mentioned above.Luetscher and Jeannin ( 2004) measured air flow rates around 10 kg/s in Hölloch Cave (Switzerland) and La Diau (France).Recently, airflow rates around 40 kg/s were measured in Shpella Shtares Cave, in Albania (42.31766°N, 19.85403°E, 1,427 m a.s.l, (see Pastore et al. (2019) for more details about this cave).Taking a couple of meters as a rough estimation of the conduit radius in Shpella Shtares, the thermal perturbation should extend over kilometers and hundreds of meters in the axial and radial directions, respectively.
Assessing the thermal penetration length at different time-scales is fundamental to quantify bio-geochemical processes in karst systems.Our study demonstrates that seasonal temperature fluctuations propagate from a few tens of meters up to several thousands of meters into a karst system, depending on ventilation rate and conduit geometry.By way of illustration, let us cite the large ice caves in Austria which host perennial ice formations over several hundreds of meters in the lower entrance zone of ventilated cave systems (e.g., Eisriesenwelt (Obleitner & Spötl, 2011)).However, it is the advent of speleothems as precisely and accurately dated paleoclimate archives which rises the significance of our study.Stalagmites and flowstones are commonly used to reconstruct climate fluctuations based on their isotope (oxygen and carbon) and trace element concentrations (Fairchild & Baker, 2012).In both cases, the proxy partitioning depends strongly on the cave temperature.While it is commonly assumed that the cave temperature at depth is constant, this is only true beyond the convection length which strongly depends on the time-scale of interest.
Because the temperature-dependent oxygen isotope fractionation between water and calcite ranges between 0.18 and 0.21‰/°C (Kim & O'Neil, 1997;Tremaine et al., 2011), seasonal temperature changes >0.5°C represent a significant source of uncertainty in the speleothem-proxy interpretation and it might be of interest to sample beyond the convective region.Meanwhile, it can be strategic to focus speleothem research on the convective region associated with decadal to millennial cycles, as this area will respond sensitively to outside climate changes without picking up short-term temperature fluctuations.

Conclusion
We developed a model coupling buoyancy-induced convection (chimney effect) in a single karst conduit with conduction in the rock mass.We used this model to investigate the propagation of thermal perturbations inside a karst massif.Assuming dry air and a simplified geometry, and reducing the time variations of the external temperature to the annual fluctuations, we performed a parametric study to identify general trends regarding the effect of the air flowrate and conduit size on the amplitude and relaxation length of the  thermal perturbations induced by the air flow.The thermal anomaly (defined as a shift in the AMT) and the annual temperature fluctuations both propagate over the same distances from the entrance.We call this distance the convection length and we show that at high Reynolds number, it is approximately proportional to the air flow rate divided by the square root of the cave radius.The order of magnitude of the theoretical convection length seem to compare satisfactorily with field data obtained from a mine and two caves.Convection lengths of a few tens to a few hundreds of meters should be commonly encountered in ventilated caves.In extreme cases, it could go up to kilometers.

•
A numerical model was developed to investigate heat transfer inside ventilated caves driven by chimney effect • The airflow can propagate annual temperature fluctuations over a few tens to a few hundreds of meters from the cave entrance (over kilometers in extreme cases)

Figure 1 .
Figure 1.Cave geometry (the computational domain is displayed in red).The rock massif is assumed infinite in both vertical directions.

Figure 2 .
Figure 2. Rock domain with the boundary conditions of the heat conduction Equation 7.

Figure 3 .
Figure 3. Nusselt number as a function of the Reynolds number (Prandtl number: Pr = 0.71).The red curve was used for the simulations of case A.

Figure 4 .
Figure 4. Wall and air temperatures at different distances x from the entrance (case A with R = 0.189 and Re = 1.8 × 10 5 ).(a) Amplitude of the Fourier modes.(b) Time series.

Figure 5 .
Figure 5. Wall (solid lines) and air (dashed lines) temperatures as a function of the distance x from the entrance (case A with R = 0.189 and Re = 1.8 × 10 5 ).(a) Temperature profiles at various times (Inset: zoom in the diffusive region).(b) Amplitude of the Fourier modes k = 0 to 3 (Inset: log scale on the x-axis).

Figure 7 Figure 6 .
Figure 7 displays the AMT (θ w and θ a ) and the ATF (Δθ w and Δθ a ) at the beginning of the convective region.Cases A and B yield close values of θ a , which is also little sensitive to the model parameters Re and R. Indeed, θ a varies only from 0.25 to 0.6 when Re and R are increased by three and two orders of magnitude, respectively.With

Figure 7 .
Figure 7. Maximum thermal perturbation at the beginning of the convective region for cases A (Biot number estimated from forced convection correlations) and B (infinite Biot number).(a) Maximum annual mean temperature (θ w and θ a for the wall and the air, resp.).(b) Maximum annual temperature fluctuation (Δθ w and Δθ a for the wall and the air, resp.).The vertical dashed lines correspond to the limits of the laminar and turbulent regimes.

Figure 8 .
Figure 8. Characteristic lengths L0 (a) and L1 (b) as a function of the Reynolds number Re for various radius R, computed from wall or gas temperatures.The dashed lines correspond to Equation24.Symbols: the same legend as in Figure7.

Figure 9 .
Figure 9. Width W0 corresponding to the damping of the annual mean temperature in the radial directions inside the rock mass.

Figure 11 .
Figure 11.External atmosphere temperature estimated at Baulmes mine (daily averaged data).The horizontal arrow displays the time range during which the temperature inside the cave was measured by the sensors listed in Table 3.The vertical red bars correspond to the dates of the air flowrate measurements listed in Table 4.The black dashed line displays the mode k = 0 and k = 1 of the external temperature.
l.), Switzerland about 37 km West-Northwest of our study site (the closer Schopernau station (47°19′N, 10°1′E, 835 m a.s.l.;The Central Institute for Meteorology and Geodynamics of Austria 2012) is subject to thermal inversions during the winter season and thus considered as being less representative).The temperature measured at St. Gallen has been corrected from the elevation as in Section 6.1.It is displayed in Figure 16a.The FFT yielded T m = 8.7°C (after correction) and ΔT y = 7.6°C for the external AMT and ATF, respectively.

Figure 12 .Figure 13 .
Figure 12.Temperatures in Baulmes conduit.(a) Daily averaged temperatures as a function of time (the temperature at x = 0 is the external atmosphere temperature measured at Mathod weather station).(b) Temperature profiles θ s (x) and θ w (x) measured respectively in summer and winter, and amplitude of the annual variation Δθ = (θ s θ w )/2.Dimensionless temperatures defined as θ = (T T m )/∆T y .

Figure 14 .
Figure 14.Field data from D7.1 cave (all the data are daily averaged).(a) External temperature and mass flow rate; black and red dashed lines represent the modes k = 0 and k = 1 computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively.(b) Dimensionless temperatures outside and inside the cave (θ = (T T m )/∆T y ).

Figure 16 .
Figure 16.Field data of Bärenhöhle Cave (all the data are daily averaged).(a) External temperature and mass flow rate; black and red dashed lines represent the modes k = 0 and k = 1 computed from the Fast Fourier Transform of the external atmosphere temperature and the air flow rate, respectively.(b) Dimensionless temperatures outside and inside the cave (θ = (T T m )/∆T y ).

Figure 15 .
Figure 15.Map of Bärenhöhle.The red circles indicate temperature sensors.

Table 1
Thermophysical Properties (Estimated at Temperature T m = 285 K = 12°C and Atmospheric Pressure p ref = 101,325 Pa) SEDAGHATKISH ET AL.

Table 2
Geometrical Characteristics of the Instrumented Conduit of Baulmes Mine (Length L c , Mean Value and Standard Deviation of the Perimeter P and the Cross-Sectional Area A From the Survey of 28 Cross-Sections)

Table 3
Location of the Seven Temperature Probes Inside the Instrumented Conduit of Baulmes Mine

Table 4
Air Mass Flow Rate Measured in Winter and Summer

Table 5
Location of the Five Temperature Probes Inside the Instrumented Conduit of D7.1 Cave Note. x is the distance from the entrance.Journal

of Geophysical Research: Earth Surface
SEDAGHATKISH ET AL.

Table 6
Location of the Three Temperature Probes Inside Bärenhöhle Note. x is the distance from the lower entrance.