A minimal model for predicting ventilation rates of subterranean caves

Subterranean caves contain calcium-carbonate ($CaCO_{3}$) mineral deposits, known as speleothems, which record important clues about past climate. Much like tree-rings and ice cores, these formations can be used to infer past conditions and identify major climatic shifts [8, 14, 73]. Speleothem records, however, can be difficult to interpret due to a number of competing factors that influence chemical and isotopic composition, mineral fabric, and timing of deposition [30, 24, 47, 22, 77]. Early cave researchers realized the potential importance of internal cave conditions [10, 37, 36], but only within the past twenty years has in situ monitoring technology become sufficiently advanced to record microclimate data with which to calibrate speleothem records [13, 55, 64, 6, 7, 52, 46, 49, 4, 69]. These studies highlighted the fact that in order to fully understand a speleothem paleoclimate record, it is necessary to have a holistic understanding of the site-specific parameters that control stalagmite growth. These include but are not limited to rainfall amount, hydrologic saturation, vegetation and soil productivity, dripwater residence time and water-rock interaction, drip rates, and ventilation [22].

Speleothems are the end product of a series of chemical reactions: (1) rainwater combines with soil-zone carbon dioxide to form carbonic acid: $H_{2}O+CO_{2(g)}\to H_{2}CO_{3}$; (2) carbonic acid percolates downward and slowly dissolves carbonate bedrock, creating a high pCO₂ solution of calcium and bicarbonate: $H_{2}CO3+CaCO_{3}\to Ca^{2+}_{(aq)}+2HCO_{3(aq)}^{-}$: and (3) the drip reaches a cave void, where the gradient between low-pCO₂ cave air and high-pCO₂ dripwater causes the drip to degas CO₂, driving the pH of the solution up and causing re-precipitation of calcium carbonate: $Ca^{2+}_{(aq)}+2HCO_{3(aq)}^{-}\to CO_{2(g)}+H_{2}O+CaCO_{3}$ [57]. Cave ventilation modulates cave air CO₂ concentration and thus the gradient between air and drip, driving both the timing and the vigor of speleothem growth. Without ventilation, cave air CO₂ would come to chemical equilibrium with dripwater CO₂, and speleothem formation would cease.

In addition to controlling the rate and timing of deposition, ventilation has been shown in several long-term cave monitoring studies to have significant impacts on the isotopic composition of dripwater and subsequent speleothem calcite [64, 31, 48, 71, 27]. Results from other field, theoretical, and laboratory studies support the hypothesis that under strong ventilation regimes where CO₂ is quickly removed from dripwater, rapid precipitation may induce kinetic isotope fractionation, causing a shift in oxygen and carbon isotopes within the speleothem and making paleoclimatic interpretation of the record much more difficult [74, 62, 54, 58, 18, 65].

Because ventilation controls the rate of calcite precipitation, it not only influences when a speleothem forms (which seasons), but also where the speleothem forms; i.e. on or within the cave ceiling (stalactite) or on the cave floor (stalagmite). Wherever dripwater encounters low-CO₂ cave air, there is potential for calcite precipitation up-stream of the stalagmite. This phenomenon is known as Prior-Calcite-Precipitation, or PCP [23]. When PCP occurs, dripwater cation-to-calcium (X/Ca) ratios change as Ca is preferentially removed from solution and minor elements become more concentrated [22]. While PCP is often associated with climate-driven water balance above the cave, where more/less hydrologic saturation and rainfall amount above the cave results in higher/lower drip rates, lower/higher PCP, and lower/higher speleothem X/Ca ratios [42, 43, 53, 72, 70, 45, 66], strong ventilation and low cave-air CO₂ also promotes PCP and exerts control on variations in speleothem x/Ca ratios [25, 63, 76, 68].

Given that the process of cave ventilation can have significant impacts on the isotopic and trace element composition of modern calcite, we must endeavor to understand the ventilation of each cave that contains speleothems of paleoclimatic importance. This is especially true in light of recent work that recognized how cave ventilation varies seasonally with well-defined global patterns [41]. While past studies have quantified cave ventilation in terms of internal temperature measurements [17, 12, 33], such data would not be available for the task of climate reconstruction. Instead, predictive models must be developed.

Here, we construct a theoretical framework to predict ventilation rates from a minimal set of external information. Building on previous work [75, 17, 12, 46], we model ventilation as arising from buoyancy-driven flows created by the internal-external temperature difference. [21] showed that the temperature explained more than 99% of variations in cave air density, thus, we used temperature difference to model ventilation rates. Unlike previous studies, though, our model does not rely on any internal measurements. It only requires basic physical parameters of the cave itself and easily obtainable outside temperature data. To close the model, we require the flow to be critical via the composite Froude number — a condition borrowed from the field of oceanography [2, 3, 15, 35, 59]. Analysis of the relevant scales, along with a few key simplifications, produces a system that can be solved through perturbation methods, ultimately giving explicit formulas for how ventilation depends on the system’s physical parameters.

To test our ventilation model, we use two tracers common to all limestone caves; carbon dioxide (CO₂) and radon-222 (²²²Rn). As previously discussed, cave air CO₂ is derived from microbial decay and plant-root respiration in the soil zone, as well as dissolved limestone bedrock. Cave air CO₂ concentrations are modulated by a combination of input from both dripwaters and gaseous transport through cracks and fissures above the cave, import of atmospheric CO₂ via ventilation, and removal via condensation corrosion [61, 32, 5, 51]. ²²²Rn is a daughter isotope in a radioactive decay chain that begins with ²³⁸U and ends with ²⁰⁶Pb. Limestone contains an average 1.3-2.5ppm ²³⁸U, and therefore acts as a continuous source of radon production and advection into cave chambers through dripwaters, surfaces, and cracks [28, 34]. ²²²Rn is an ideal tracer for cave ventilation because atmospheric concentrations of radon are negligible and it is an inert gas with a 3.82-day half-life, which is sufficiently long when compared with most cave-air exchange rates [28, 55, 60, 46]. We report a set of time-resolved measurements of cave air gas concentrations from a subterranean cave located in Florida Caverns State Park [44]. The theory accounts for qualitative trends in the measurements throughout the year, and, over certain periods, even shows quantitative agreement with weekly-scale fluctuations. This comparison also allows the production rates of CO₂ and ²²²Rn to be estimated, which could complement other, more invasive, estimation techniques [20, 46].

Dragon’s Tooth Cave (DTC) is a pristine, karst-hosted dry cave located in Florida Caverns State Park, Marianna, northwest Florida (Fig. 1a-b). The cave formed within a carbonate bluff along the nearby Chipola River, and has two primary passages that trend northwest to southeast in alignment with regional bedrock fractures. The two large passages are intersected by a third east-west trending corridor that dips at approximately 45 degrees from west to east, serving as an entrance to both sections. The northern section was first surveyed and mapped in 1978. The floor of the northern section lies at 18 meters above sea level (ASL) and is composed of mud-clay, which is typical of caves in Marianna that are subjected to seasonal groundwater inundation when Chipola River levels exceed 3-4 meters above the riverbed [69]. Passages within DTC are 2-3 meters tall, thus floodwaters tend to fill the passages and remove any stalactites or soda straws that may have formed, minimizing the potential of flowstone decoration.

In 1984, spelunkers observed airflow through a pile of breakdown near the entrance and dug a narrow tunnel making a connection to the southern section of the cave, finding a previously unexplored room that was highly decorated with beautiful speleothem formations. This room later came to be known as the Dragon’s Belly (DB) [50]. The floor of the Belly lies at 21 meters ASL; groundwater slowly fills the Belly from below only during Chipola flooding events that exceed 5 meters above riverbed. Topography above the cave suggests that the Belly is the remains of a large diameter NW-to-SE conduit that has collapsed at both ends of the room. The porous nature of the breakdown materials promotes significant volumetric airflow (ventilation) and therefore the Belly is highly decorated with pristine white stalagmites, flowstones, and moon-milk atop several large sections of breakdown blocks. In 2012, our group mapped the DTC using Leica ScanStationP20 to obtain highly accurate measurements of its dimensions. The main chamber of the DB is nearly rectangular – 41 m long, 12 m wide, and 6 m high with 45 degree slopes at each end of the room that pinch off at the ceiling, and it connects to the entrance of the cave system via a single, narrow passageway (0.6 m high, 0.3 m wide and 3 m long). These geometrically simple features make the DB an ideal site for testing basic ventilation models. Once a firm theoretical foundation is established, future studies could examine more complex cave systems (e.g. complex geometries and/or multiple openings).

Beginning in November 2011, we established a micrometeorological monitoring program inside the DB cave to continuously collect spatial and temporal variations of ²²²Rn, CO₂, and temperature. ²²²Rn and CO₂ levels were measured via two cave stations deployed 15 meters apart (see Fig. 1b). These stations were equipped with Durridge RAD7 detectors to measure radon, LiCor-820 gas analyzers to measure CO₂, as well as Vaisala HMP45C temperature and relative humidity sensors. Four additional Onset HOBO Pro v2 External Temperature Data Loggers were also deployed throughout the Belly. All data was recorded hourly, with the exception of occasional outages due to flooding events or power failures. Meanwhile, we obtained external temperature data from daily measurements taken at Marianna Municipal Airport, located 8 km from the study site.

We model the Dragon’s Belly Cave as a well-mixed compartment with a single opening that exchanges heat with the surrounding surfaces and the outside environment. Because the internal cave temperature $T_{c}$ is typically different from the external temperature $T_{e}$, a buoyancy-driven flow is created in the connecting passageway [40, 59, 11]. During the summer, warm outside air flows into the cave through the opening, is cooled by contact with the cave walls, and exits through the lower part of the same opening. The airflow reverses during the winter; see Figs. 2a–b for a schematic. As long as vertical mixing and viscous effects are weak (i.e. large Peclet and Reynolds numbers respectively), the flow in the opening is essentially that of two layers — a top layer of depth $H_{1}$ and speed $U_{1}$, and a bottom layer with values $H_{2}$ and $U_{2}$. Enforcing conservation of volume inside the cave, along with the Boussinesq approximation, results in a symmetric exchange flow; i.e. both the depths and speeds are the same in the two layers, $H_{1}=H_{2}=H$ and $U_{1}=U_{2}=U$. Clearly, $H$ must be the half-height of the passageway ($H=0.3$ m), but $U$ must be determined through other considerations, starting with a balance of the cave’s total thermal energy.

We now discuss the measurements taken at the Dragon’s Belly site, shown by the red curves in Figs. 3 a–c for CO₂ and d–f for radon. Cave air gas concentrations in the Belly were collected at 30-minute intervals over a 14-month period, with occasional interruptions due to flooding and/or power outages. Both the perturbation model and the numerical model were performed at 24-hr resolution, models and measured data are plotted together in Fig. 3. The highest concentrations of cave gases occur during fall and spring, while the lowest concentrations occur during summer. This result is counter to those in caves under predominantly chimney-driven airflow regimes that promote strong winter-time ventilation (e.g. Obir Cave [64]; Inner Space Caverns [7]; Hollow Ridge Cave [46]).

During the winter, carbon-dioxide levels in the Belly are near their minimum, while radon levels are intermediate. Qualitatively, these trends suggest ventilation to be lowest during fall and spring, when daily average inside-outside temperature differences are small, and highest during summer and winter, when the temperature differences are maximal. The observations therefore support the idea of ventilation being caused primarily by buoyancy-driven flows.

Predictions of CO₂ and radon are shown by the dark curves in Fig. 3. The gray curves represent the numerical model, with convection switched on/off for each surface individually, while the black curves show the perturbation solution given by Eq. (8). The two predictions are nearly indistinguishable, providing validation for the additional assumptions made in the perturbation model. We estimate the source term, $\Phi$, by minimizing the difference between model and measurements, giving $\Phi=3.1\times 10^{-9}$ kg/(m²s) for CO₂ and $\Phi=25\times 10^{-6}$ kBq/(m²s) for ²²²Rn. We note that our estimated source terms for radon and CO₂ is comparable to others found in the literature (e.g. Radon: $\Phi=9.9\times 10^{-6}$ kBq/(m²s) [20]; $\Phi=13~61\times 10^{-6}$ kBq/(m²s) [46]; CO₂: $\Phi~=~1.989~\times~10^{-10}$ kg/(m²s) [6]; $\Phi~=~1.56~\times~10^{-9}$ kg/(m²s), [20]). With the source term as the only fit parameter, both the CO₂ and radon predictions correspond well with measurements over the majority of the 14-month study. In particular, the model accounts for the higher levels of cave gases in the fall/spring and lower levels in the summer/winter. The only obvious discrepancy is in comparing radon levels in the summer, where the theory predicts a significantly higher level than was measured. Carbon-dioxide shows good agreement during this same period, though, suggesting the discrepancy to be specific to radon.

While not a focus of this study, hydrologic conditions within the rock above the Belly may have a significant influence on agreement between modeled and measured ²²²Rn over seasonal time periods. [1] demonstrated that radon emanation from limestone occurs from within a surface layer only a few microns thick, and that under dry conditions, alpha-recoil causes radon gas to be ejected from the crystal lattice into inter-crystal cracks and imperfections, eventually to find its way out of the limestone through relatively slow diffusion in air at approximately 10^-2 cm² s^-1. Conversely, under wet conditions within the epikarst, alpha-recoil can cause radon particles to be ejected directly into pore water. As dripwater percolates downward, it has a scavenging effect, driving a significantly higher flux of radon during periods when the epikarst is hydrologically saturated. An excellent example of this phenomenon was illustrated by [46] at Hollow Ridge Cave, when a tropical storm delivered 83 mm of rainfall, saturating the epikarst and causing a nearby river to flood and plug up the lower entrances of the cave, temporarily ceasing ventilation. Directly following this rainfall, ²²²Rn emanation rates increased from ’normal’ 48-222 Bq m^-2 hr^-1 to approximately 1200 Bq m^-2 hr^-1. Unpublished data from laboratory testing of limestone from both Dragon’s Tooth and Hollow Ridge Cave confirms that wet transport can increase ²²²Rn emanation by 1.6$\sim$2.4 times dry transport (William Burnett - FSU, personal communication, 2015).

Figure 4 illustrates that prior to August 2012, most of the monitored drip sites were essentially zero. During August, a 120 mm rainfall and several successive $\sim$40 mm rainfalls prompted hydroligcal saturation within the overburden. Approximately 1-2 weeks after each event, an increase in ²²²Rn was observed, followed by a decay to background levels. In September 2012, all monitored sites rapidly increased in drip rate as the epikarst came to full saturation. These saturated conditions persisted until May 2013, when lower rainfall and increased evapotranspiration conspired to reduce the hydrologic saturation, causing measured radon to fall below values predicted by our model.

In additional to the seasonal variation, the CO₂ and radon measurements show more rapid fluctuations on a few characteristic scales. Figures 3 b–c and e–f zoom on selected date ranges to highlight these. First, the measurements show very rapid, diurnal fluctuations, which result from the different day/night temperatures. In the model, we only input daily outside temperatures, and so these very fast fluctuations are not present. Presumably they could be captured with more temporal resolution though (e.g. hourly readings). Interestingly, there is a second, intermediate timescale — on the order of 5–9 days — that is also evident in the CO₂/radon data. This nearly weekly cycle is likely caused by synoptic-scale meteorology [39], though manmade activity may play a role too [16, 19]. Regardless, the model shows variations over exactly the same scale. This implies that the weekly cycle in cave-gas levels is ultimately tied to a weekly cycle in outside temperatures that drives an exchange flow as described by our theory. Over certain periods, for example Fig. 3e, the theory even shows quantitative agreement with tracer measurements.

Here, we have constructed a theoretical framework to predict subterranean-cave ventilation rates from a minimal set of external information — namely, knowledge of the outside temperature and the physical dimensions of the cave. The major advantage is that external temperatures are much more easily measured or, in the application of climate reconstruction, estimated, than are conditions inside the cave. Through scaling analysis and a few key modeling assumptions, we obtained explicit relationships for how ventilation depends on system parameters, given by Eq. (8). Comparison with time-resolved, in-situ measurements demonstrates the theory’s ability to accurately describe seasonal and synoptic-scale fluctuations of transported cave gases. As an immediate application, this theory could be used to improve estimates for production rates of cave gases such as CO₂ and ²²²Rd. Currently-used field techniques can be destructive to the cave’s fragile environment, but our methodology is completely non-invasive as it relies only on knowledge of outside temperatures.

A longer-term application is to use these results to improve speleothem interpretation. Inferring past climate from speleothems is inherently an inverse problem, typically treated by representing the climate’s influence on the proxy (e.g. speleothem, tree ring, etc.) via a transfer function, which is then inverted. While it has been recognized that ventilation plays a critical role in speleothem growth, such effects have not yet been quantitatively incorporated into paleoclimate inference models, likely due to the complexity involved. Our model, however, offers explicit formulas for how ventilation depends on cave and climatic variables. The transparency afforded by this model could therefore prove useful in obtaining tractable transfer functions for climate reconstruction. In essence, the simpler the forward model, the more feasible inference becomes.

Geophysical Fluid Dynamics Institute contribution number xxx. This research was supported by the NSF Grant AGS-1032403; additional support for K. Khazmutdinova was provided through teaching and research assistantships at Florida State University. The authors sincerely thank the management staff of Florida Caverns State Park, especially Kelly Banta, for providing an incredible opportunity to study the Dragon’s Tooth Cave. Thanks to Cameron Ridgewell for creating a 3D map of the Dragon’s Tooth Cave. We thank William Burnett for help with interpreting radon measurements, and William Dewar for inspiring suggestions and comments. The meteorological data from Marianna Airport is obtained from https://www.ncdc.noaa.gov/cdo-web/. The input files and results of the numerical calculations are available from the authors upon request (kk11m@my.fsu.edu).