Assessing Landscape and Seasonal Controls on Soil CO₂ Fluxes in a Karst Sinkhole

Exchange of carbon dioxide (CO₂) from the land surface to the atmosphere is an important process in the global carbon cycle, yet the controls and magnitudes of CO₂ emissions from the soil surface remain poorly constrained (Jian et al., 2022). Carbon dynamics are of particular interest in karst regions, which are characterized by relatively soluble carbonate rock, usually dolomite or limestone rock, and which occupy as much as 15 percent of the global land surface (Ford and Williams, 2007). Over timescales of thousands to hundreds of thousands of years, karst regions absorb CO₂ during chemical weathering of carbonate rocks. Much of the CO₂ involved in this weathering process comes from overlying soils or the atmosphere, with gas diffusion found to be a more important process than off-gassing from aqueous solutions (Breecker et al., 2012). Thus, understanding the magnitudes and directions of CO₂ fluxes in karst soils is necessary to perform accurate regional and global carbon budgeting.

Subsurface transport of gases such as CO₂ occurs via two processes, diffusion and advection. Advection occurs as a mass flow in which gases are typically transported via flowing wind or water. Diffusion occurs due to concentration gradients driving gas movement from areas of high to low concentration, and it is the more important process for gas exchange in most soils (Pumpanen et al., 2003; Laemmel et al., 2017). The rate of gas diffusion depends on two factors: (1) the diffusion coefficient, D_s, and (2) the concentration gradient, (i.e., ΔC/Δx for one-dimensional transport). The diffusion coefficient is a measure of the rate of movement of a gas under constant temperature and pressure conditions through a uniform medium. In soils, D_s is affected by several factors, with the volume of air-filled pore space representing a particularly important term (Jin and Jury, 1996; Lafond et al., 2011).

Concentration gradients can be quantified by measuring partial pressures, or equivalent molar concentrations, at multiple depths. In typical soils, CO₂ concentrations tend to increase with depth due to respiration and absence of any large-scale sinks (Tang et al., 2003). Conversely, studies conducted in karst areas have identified conditions under which soil CO₂ concentrations near the surface are greater than those at depth (Chen, 2019). One possibility for these observations is CO₂ consumption during dissolution of carbonate minerals (Chen, 2019; Zhao et al., 2021). A second possibility is that venting of cave air removes CO₂ from deeper soil and epikarst; for example, Benavente et al. (2010) found that 30–50 percent of CO₂ in a Mediterranean cave originated from the overlying soil profile, suggesting downward diffusion of CO₂ from soils into the cave. This venting process typically occurs during cooler seasons of the year, when air inside the cave system is warmer and less dense than the surrounding air (Mattey et al., 2016; Krajnc et al., 2017; Riechelmann et al., 2019). These observations imply that a zero-flux CO₂ plane can exist in the subsurface, from which gas diffuses both upward toward the soil surface and downward through the subsurface, similar in concept to hydraulic zero-flux planes that form in soils (Brutsaert, 2014). However, limited research has been conducted to assess the frequency with which near-surface zero-flux planes exist, and the corresponding implications for gas diffusion remain poorly understood.

Variation in wetness, soil depth, and other environmental factors associated with topography can also affect CO₂ fluxes. Previous work has shown that riparian areas often have greater soil respiration than surrounding uplands. This result is likely due to wetter conditions encouraging greater ecosystem productivity (Webster et al., 2008; Pacific et al., 2009, 2011), which can supersede smaller D_s values found in riparian compared to upland soils (Pacific et al., 2008). Other work focused on hillslope profiles has shown that ridge tops often have greater soil respiration rates than lower-lying positions (Tian et al., 2019), particularly in wet years (Kopp et al., 2022). Slope angle can also affect soil respiration and carbon allocation, with flatter or concave slopes found to have greater respiration than steeper or convex slopes (Saggar et al., 1999; Fissore et al., 2017). It is likewise probable that topography affects soil respiration and gas diffusion in karst areas, but the controls on CO₂ gas movement in such systems have not yet been critically examined.

To evaluate the magnitudes and directions of CO₂ fluxes in a karst landscape, we instrumented three locations along a sinkhole slope—shoulder, back slope, and toe slope—and measured hourly CO₂ concentrations at 20 and 60 cm depths from winter 2017 to summer 2019. We used those measurements, along with assumed atmospheric CO₂ concentration and data for soil physical properties (e.g., porosity, water content, temperature), to calculate CO₂ fluxes. Our first objective was to compare upward and downward CO₂ fluxes for the three sinkhole positions. We hypothesized that the upward fluxes would be similar in magnitude between the three soil profiles, but the downward fluxes would be largest in the toe slope and smallest in the shoulder profile, based on differences in vertical distance from the underlying cave system. Our second objective was to examine seasonal trends in CO₂ fluxes. Here, we hypothesized that CO₂ fluxes would be smaller during cooler months of the year, due to limited soil respiration along with wet conditions limiting gas exchange, and fluxes would also include a sizeable downward component due to convective venting from the underlying cave system. Conversely, we hypothesized that CO₂ fluxes would be larger and primarily upward during warmer months of the year, because respiration typically increases with temperature, and gas diffusion rates typically increase with lower soil water contents (Fang and Moncrieff, 2005; Ali et al., 2018). Our third objective was to evaluate the sensitivity of gas fluxes to model parameters used to estimate gas diffusivity. We expected that the estimated flux would be sensitive to available pore space in the soil profiles (i.e., the difference between soil porosity and soil water content), and that the effects would be strongest in the non-growing season, when water contents are typically high, compared to the growing season.

The field site was located at James Cave in southwest Virginia (Supplemental Material Figure S1, https://www.aegweb.org/e-eg-supplements). James Cave is a doline karst cave system with substantial amounts of limestone and dolomite rock (Eagle et al., 2015; Schreiber et al., 2015). Based on the U.S. Department of Agriculture Soil Survey (Soil Survey Staff, 2019), the soils surrounding James Cave are classified as Slabtown silt loam, Lowell silt loam, and Wurno-Newbern-Faywood silt loam, with soil thickness ranging from 0.5 to 2.0 m (Eagle et al., 2015). The field site was managed as a beef cattle pasture, similar to the predominant land use of the surrounding area. The vegetation was a mixed stand of cool-season grasses, including tall fescue (Festuca arundinacea).

Sensors were installed in December 2016 in the sinkhole that surrounds the entrance to James Cave, with vertical profiles instrumented at the shoulder, back slope, and toe slope of the sinkhole slope (Figure 1). Each profile contained two CO₂ sensors (Vaisala GMM221, 5 percent maximum concentration, Vaisala Corporation, Vantaa, Finland) and two soil water content and temperature sensors (CS655, Campbell Scientific, Logan, UT) located at 20 and 60 cm below the surface. The 60 cm soil water content/temperature sensor did not properly function at the back-slope location, so an additional sensor was installed at that depth (5TM, Decagon Devices, Pullman, WA). The data were measured every hour using a CR1000 logger (Campbell Scientific, Logan, UT).

All installed sensors functioned from 7 February 2017 through 8 July 2018, when the 20-cm-deep CO₂ sensor failed at the back slope. The 60-cm-deep CO₂ sensor failed at the shoulder location on 24 September 2018, and the 20-cm-deep CO₂ sensors failed at the toe slope on 19 September 2019. All three profiles therefore had a common 1.5 year period in which both CO₂ sensors functioned, and the toe slope had an additional year of data. The barometric pressure and hourly precipitation data were downloaded through the National Oceanic and Atmospheric Administration’s National Centers for Environmental Information from the Integrated Surface Database for the Blacksburg Airport station, which was located approximately 17 km from the study site. Temperature readings came from the adjacent soil temperature sensors.

Volumetric soil cores were collected from the site in order to determine water retention properties, which we used to infer total porosity and macro- versus micro-porosity. The cores were collected from depths of 15, 20, 35, and 40 cm for the shoulder location and from 15, 20, and 40 cm depths for the back-slope and toe-slope locations. The cores were first saturated with 0.01 M CaCl₂ solution and were then analyzed for water retention using a positive pressure system (SoilMoisture, Inc., Santa Barbara, CA). Soil water content was determined at pressure heads of 0 cm (saturation), −340 cm (assumed to represent field capacity), −1,020 cm, and −3,060 cm. Porosity was calculated as the water content at saturation, macro-porosity was calculated as the difference in water content between saturation and −340 cm, and micro-porosity was calculated as the water content at −340 cm.

For every hourly time step, we calculated two fluxes: one flux representing the upper 20 cm (J_0–20) that used the measured C_corrected value at 20 cm and an assumed atmospheric CO₂ concentration of 400 μmol mol⁻¹, and a second flux representing the 20–60 cm depth interval (J_20–60) that used C_corrected values from those two depths. To calculate the gas diffusivity constant for J_0–20, we used measured T and θ values from 20 cm, and for J_20–60, we used the mean T and θ values from the 20 and 60 cm sensors. The measured C_corrected values at 20 cm were always greater than the assumed atmospheric concentration. Therefore, whenever C_corrected at 20 cm was greater than C_corrected at 60 cm, we assumed that a near-surface zero-flux existed for CO₂ diffusion, and we calculated an upward flux in the upper 20 cm of the soil profile (i.e., J_upward = J_0–20) and a downward flux from 20 to 60 cm (i.e., J_downward = J_20–60). Whenever C_corrected at 20 cm was less than or equal to C_corrected at 60 cm, we calculated only an upward flux using a depth-weighted average of the two fluxes (i.e., J_upward = 1/3 × J_0–20 + 2/3 × J_20–60). This approach assumed isothermal conditions within the soil profile (Jones, 2014). There were some brief periods in which the CO₂ concentrations at one or more locations were out of sensor range (i.e., >55,000 ppm_v or 5.5 percent); fluxes were not estimated for those times.

We also quantified seasonal and annual CO₂ fluxes from each profile, with the upward and downward fluxes assessed separately for each profile. Annual fluxes were calculated by summing hourly fluxes between 1 April 2017 and 31 March 2018. For the seasonal fluxes, we divided the data into the warm, growing season from 1 April 2017 to 30 September 2017, and the cool, non-growing season from 1 October 2017 to 31 March 2018. We used linear interpolation to fill gaps in the hourly flux records, using the flux values that were calculated immediately before and after the missing records. This approach gave similar values as using only records with measured fluxes (data not shown), indicating that the linear interpolation process did not overly influence the results.

We performed a brief sensitivity analysis to test how the porosity (ϕ) term in Eq. 3 and Eq. 4 affected calculations of D_s and J_upward. First, we calculated D_s through time for each profile using five assumed porosities: ϕ = 0.40, 0.44, 0.48, 0.52, and 0.56 cm³ cm⁻³. These porosities reflected the range of values measured in the soil samples used for soil characterization. The D_s term was assumed to equal 0 any time that the measured volumetric water content was equal to or greater than the assumed porosity. We next used Eq. 2 to Eq. 5 to determine upward CO₂ fluxes for each profile and then summed the cumulative CO₂ efflux for the period between 1 April 2017 and 31 March 2018.

The water retention curves were similar between most soil cores (Supplemental Material Figure S2, https://www.aegweb.org/e-eg-supplements). For example, at field capacity (i.e., pressure head = −340 cm), volumetric water content values were between 0.32 and 0.41 cm³ cm⁻³ for the toe-slope location, 0.33 and 0.39 cm³ cm⁻³ for the back-slope location, and 0.33 and 0.38 cm³ cm⁻³ for the shoulder location. Similar ranges of water contents were seen for the −1,020 cm and −3,060 cm pressure heads.

We used these water retention values to calculate total, macro-, and micro-porosities (Supplemental Material Table S1, https://www.aegweb.org/e-eg-supplements). Porosity varied between 0.41 and 0.55 cm³ cm⁻³ for the toe-slope location, 0.41 and 0.49 cm³ cm⁻³ for the back-slope location, and 0.40 and 0.47 cm³ cm⁻³ for the shoulder location. We therefore used mean values of ϕ = 0.48 cm³ cm⁻³ for the toe-slope profile, ϕ = 0.45 cm³ cm⁻³ for the back-slope profile, and ϕ = 0.44 cm³ cm⁻³ for the shoulder when calculating D_s. The macro-porosities varied from 0.06 to 0.15 cm³ cm⁻³, with the toe-slope location having greater macro-porosity at 15–20 cm compared to the other two profiles. Micro-porosity varied from 0.30 to 0.40 cm³ cm⁻³, with similar sample-to-sample variation in all three profiles.

The largest precipitation events occurred in the summer and early fall, with maximum rates of 1.0 to 2.5 cm hr⁻¹ (Supplemental Material Figure S3a, https://www.aegweb.org/e-eg-supplements). Volumetric water contents were lower in the summer and fall months, with a range of 0.10 to 0.35 cm³ cm⁻³. During the cool season (i.e., the months of October through March), the volumetric water contents were higher, varying between 0.30 and 0.45 cm³ cm⁻³. The shoulder and toe-slope locations tended to experience the greatest seasonal differences in volumetric water contents, while the back slope had smaller-magnitude changes throughout the study period.

Soil temperatures also fluctuated seasonally, varying between 15°C and 26°C for the warmest months (June–September) and between 1°C and 13°C in the coolest months (December–March; Supplemental Material Figure S3b, https://www.aegweb.org/e-eg-supplements). Soil temperatures showed marked transitions between the colder and warmer extremes during early spring (i.e., March through May) and early fall (i.e., September through November).

Calculated values for gas diffusivity, D_s, varied through time due to differences in soil temperatures and water content, and also varied among the three soil profiles due to those factors and different assumed values for ϕ (Figure 2). In the first year of monitoring, the toe slope had peak D_s values that were nearly three times the values of the other two profiles. In the second and third years, the shoulder and toe slope profiles had similar D_s values, while the back slope continued to have smaller values. In all three profiles, the D_s values were larger during the warm season and were often near 0 during the cool seasons.

The shoulder profile had relatively high CO₂ concentrations early in the warm season, at times exceeding 55,000 ppm_v (5.5 percent), and had relatively low CO₂ concentrations during the cool, non-growing season (Figure 3a). The 60 cm sensor tended to have higher concentrations than the 20 cm sensor, though there were many brief periods in which the concentration gradient became inversed (i.e., greater CO₂ concentration at 20 cm than 60 cm). Those periods were indicative of a near-surface zero-flux plane, and they translated to relatively large upward CO₂ fluxes and much smaller downward fluxes (Figure 3b). During the cool seasons, the profile had negligible fluxes in either direction. Upward CO₂ fluxes peaked between 0.5 and 1.2 μmol CO₂ m⁻² s⁻¹, while the magnitude of downward fluxes never exceeded 0.2 μmol CO₂ m⁻² s⁻¹.

The back-slope profile had similar magnitudes and seasonal trends in CO₂ concentrations as the shoulder, with CO₂ concentrations peaking above 55,000 ppm_v during the beginning of the warm season and then gradually declining (Figure 4a). The CO₂ fluxes for the back-slope location had similar magnitudes of upward fluxes as the shoulder location (Figure 4b). However, during the first warm season, the downward fluxes were consistently larger in magnitude than those detected in the shoulder profile. Furthermore, the downward and upward fluxes often had similar magnitudes in this profile, with larger fluxes in both directions during warm seasons and near-zero fluxes during the cool seasons.

The toe-slope location had more consistently elevated CO₂ concentrations compared to the other two profiles (Figure 5a). Here, the 20 cm sensor often recorded greater CO₂ concentrations than the 60 cm sensor, particularly during the cool, non-growing season. The toe-slope profile also tended to have larger CO₂ fluxes than the other two profiles (Figure 5b). The CO₂ fluxes reached peak values of 1.5–2.0 μmol CO₂ m⁻² s⁻¹ (J_upward) and −0.5 μmol CO₂ m⁻² s⁻¹ (J_downward). Downward fluxes were often present during the cool season, with larger magnitudes during the winter of 2017–2018 as compared to 2018–2019.

Annual CO₂ efflux for each profile was calculated by integrating the areas underneath the CO₂ flux curves for each profile from 1 April 2017 through 31 March 2018. The annual efflux values ranged from 3.5 mol CO₂ m⁻² in the shoulder to 14.5 mol CO₂ m⁻² in the toe slope (Table 1). Most of the efflux occurred during the growing season, with warm-season efflux representing between 76 percent (in the toe slope) and 84 percent (in the back slope) of total CO₂ emissions. The cumulative downward flux was much smaller, with the toe slope having a cumulative total of −0.36 mol CO₂ m⁻², the back slope having a cumulative total of −0.13 mol CO₂ m⁻², and the shoulder having a cumulative total of −0.04 mol CO₂ m⁻².

The sensitivity analysis showed that the magnitude of D_s varied widely based on the assumed value of ϕ (Figure 6). Here, we use the toe-slope location as an example. Under the lowest tested ϕ value (0.40 cm³ cm⁻³), the maximum D_s values, based on mean θ and T values from the 20 and 60 cm sensors, were between 1.0 and 1.2 mm² s⁻¹. The largest tested ϕ value (0.56 cm³ cm⁻³) gave a maximum D_s value of 2.8 mm² s⁻¹. In the cool season, when the soil water content was high, the variation in porosity drove the D_s value closer to zero for the minimum porosity value versus 0.1–0.2 mm² s⁻¹ when larger porosity values were used.

The differences in calculated D_s values translated to sizeable differences in the estimated amount of CO₂ that was emitted from each soil profile over a 1 year period (Figure 7). The toe-slope location had the least variation, with total CO₂ efflux varying by 7× between the lowest (ϕ = 0.40) and highest (ϕ = 0.56) porosities. Cumulative effluxes between lowest and highest porosities differed by 11× for the shoulder profile and by 14× in the back-slope location. These results show the importance of accurately constraining porosity when estimating CO₂ fluxes using solid-state sensors embedded in the soil.

We initially hypothesized that the upward CO₂ fluxes would be similar in magnitude for the three profiles, whereas the shoulder would have the smallest downward CO₂ fluxes, and the toe slope would have the largest downward fluxes. We also hypothesized that upward CO₂ fluxes would be much greater in the warm season than the cold season, whereas downward fluxes would be dominant in the cold season due to convective venting from the cave releasing CO₂ back to the atmosphere. The study results provided support for these hypotheses, though some discrepancies and nuances emerged during our analysis.

The back-slope and shoulder profiles had similar upward CO₂ fluxes, partially supporting our hypothesis. The toe-slope profile, however, had more sustained and often greater upward fluxes, with peak fluxes near twice those of the other two profiles (i.e., ∼2.0 versus ∼1.0–1.2 μmol m⁻² s⁻¹). These differences primarily occurred during the warm growing season of the first study year, when the shoulder and back-slope profiles had smaller D_s estimates than the toe slope. In the second warm season, by comparison, the shoulder and toe slope had similar D_s values, and the two profiles had similar peak effluxes (1.3–1.6 μmol m⁻² s⁻¹). These comparisons indicate that CO₂ production was likely similar between locations, and the differences in efflux were attributable to variations in soil moisture conditions and resulting estimates for air-filled porosity ε.

All three profiles also had periods in which the shallow (20 cm) sensors recorded higher CO₂ concentrations compared to their deeper (60 cm) counterparts, implying the existence of a zero-flux CO₂ plane in the shallow subsurface. Downward diffusion likely occurred whenever that condition existed, and annual total downward fluxes were estimated to be 0.04–0.36 mol CO₂ m⁻². While considerably smaller in magnitude than cumulative upward fluxes, which ranged from 3.6 to 14.5 mol CO₂ m⁻², these results show that downward CO₂ diffusion may be a non-negligible process. Comparing sinkhole locations, the toe-slope profile had the largest peak and cumulative downward fluxes, as we hypothesized. These fluxes occurred primarily in the cool, non-growing season, which also supported our hypothesis. However, downward fluxes in the shoulder and back-slope locations were detected primarily during warm months, which we had not expected. These contrasting results imply that the degree and timing of downward gas diffusion in karst soils may be governed by complex interactions among CO₂ production, CO₂ consumption during carbonate dissolution, subsurface properties, and landscape position.

In temperate climates, cave ventilation increases during the cool season, as temperature-driven differences in air density allow warmer cave air to escape via cave openings and become replaced by atmospheric air that is lower in CO₂ (Meyer et al., 2014; Mattey et al., 2016). Under such conditions, elevated CO₂ concentrations in the soil would result in downward diffusion, but only so long as D_s is sufficiently large. The toe-slope location, being closest in elevation to the cave system and likely to be underlain by fractures and fissures, could have maintained air-filled pathways between the soil profile and the deeper epikarst and cave layers (Weisbrod et al., 2009; Eagle et al., 2015). The water content measurements indicated that the toe-slope profile was rarely saturated, with sustained water contents of <0.37 cm³ cm⁻³ during the winter period compared to measured porosities of 0.41–0.55 cm³ cm⁻³. Concurrent work examining cave drip patterns within James Cave also showed a prolonged dry period with limited cave recharge throughout much of the 2017–2018 cool season (Groce-Wright et al., 2022), suggesting that larger fissures and pores above the cave may have remained unsaturated. However, these air-filled pathways may not have been connected to the shoulder and back-slope profiles, which were located at a greater distance above the cave system and had near-saturated soil water contents throughout the cool season. The shoulder and back-slope profiles were instead found to only experience downward CO₂ diffusion during dry periods of the year, when ε may have been sufficient to enable rapid gas diffusion rates (Cuezva et al., 2011; Sanchez‐Cañete et al., 2011).

Besides cave ventilation, there are other potential mechanisms by which the underlying karst layers could serve as CO₂ sinks. For one, carbon dioxide can dissolve into percolating water and become converted into carbonic acid (H₂CO₃) and then bicarbonate (HCO₃⁻) during carbonate rock weathering (Ford and Williams, 2007). This process could reduce CO₂ concentrations at depth as carbonates undergo dissolution. Previous work in the James Cave system found carbonate weathering to be highly seasonal, with increased CO₂ advection identified during December to March, when aquifer recharge occurred via cave drips (Eagle et al., 2015). This period overlaps with the period of the most consistent downward fluxes in the toe slope, suggesting a possible linkage between these processes.

The CO₂ fluxes calculated in this study were sensitive to the porosity (ϕ) value used to calculate D_s, with the cumulative efflux calculated over a 1 year period being an order of magnitude less when using the smallest versus largest ϕ values. The differences were most pronounced during the colder, non-growing season, when water contents were typically high. During those times, the calculated fluxes were negligible (i.e., near zero) for the lowest porosity values versus 0.1–0.3 μmol m⁻² s⁻¹ when the highest porosity of ϕ = 0.56 cm³ cm⁻³ was used. Previous work has shown that ε is one of the most important terms for D_s (Jin and Jury, 1996; Lafond et al., 2011), as captured in the diffusion model used here (Eqs. 3–5).

We based our assumed ϕ values on measurements taken from soil cores, yet those measurements revealed substantial variability between cores (ϕ ranging from 0.40 to 0.55 cm³ cm⁻³). The actual ϕ in the field also likely differed from the values measured from the cores, due to possible errors when collecting and analyzing the cores (e.g., incomplete or over-saturation) and due to any naturally occurring differences in porosity throughout the soil profiles, such as the near-surface compaction that can result from animal traffic (Daniel et al., 2002). Discrepancies in assumed versus actual porosity values may have the greater effect on flux estimates in wet conditions, when the presence of connected air-filled pathways becomes particularly important (Jin and Jury, 1996). Our approach also did not consider the possibility of temporarily reduced gas diffusion due to near-surface wetting events, which can cause over-estimates of CO₂ efflux (Fan and Jones, 2014). These factors imply that the calculated CO₂ fluxes in this study had considerable uncertainty, and they emphasize the importance of accurately constraining porosity when estimating CO₂ fluxes using solid-state sensors embedded in the soil.

Other measurements and models for D_s have identified the importance of larger macro-pores for gas diffusion (McCourty et al., 2018; Jayarathne et al., 2020). The soil samples had estimated macro-porosities of 0.06–0.15 cm³ cm⁻³ (Supplemental Material Table S1, https://www.aegweb.org/e-eg-supplements), with the toe-slope location having a greater proportion of macro-pores compared to the other profiles. The larger pores also likely extended below the soil mantle, as highly weathered epikarst layers often contain many fractures and large pores (Eagle et al., 2015; Mattey et al., 2016; White, 2016; Cao et al., 2020). Therefore, future work should endeavor to better understand how spatial differences in vadose zone properties—including ϕ, soil macro-porosity, and epikarst fracture connectivity—influence D_s values and gas diffusion rates.

The study had several limitations. For one, many of the solid-state CO₂ sensors reached maximum readings (∼55,000 ppm_v) multiple times during the study, which resulted in data gaps. The sensors then started to fail after 1.5–2.5 years, affecting the duration of the observational record. In addition, the overall spacing between sensors in this study was large compared to other studies (e.g., 8 cm spacing in Tang et al., 2003), which increased the likelihood that the concentration gradients were not uniform between sensors. The shallowest sensor depth was 20 cm, which is substantially deeper than that used in many other studies, and this could have affected the accuracy of the calculated CO₂ efflux to the atmosphere.

We also did not directly measure CO₂ concentrations in the cave, which limited our ability to fully understand the magnitude of downward CO₂ diffusion. Other studies in similar cave systems have measured typical concentrations of <12,000 ppm_v (Breecker et al., 2012; Meyer et al., 2014; Kukuljan et al., 2021), albeit with isolated observations of CO₂ concentrations up to 30,000 ppm_v (Cowan et al., 2013). These values are, for the most part, much lower than the typical concentrations measured during the growing season in our soil profiles, which were often in excess of 25,000 ppm_v. It is therefore possible that a zero-flux CO₂ plane existed below our deepest sensor. If so, the instantaneous and cumulative downward fluxes calculated in this study likely under-estimated the true magnitude of downward diffusion. Future work could better understand this possibility by installing deeper CO₂ sensors (Benavente et al., 2010), simultaneously instrumenting soil profiles and underlying caves (Meyer et al., 2014), or using isotope discrimination techniques to trace carbon movement throughout the subsurface (Breecker et al., 2012; Mattey et al., 2016).

Due to limitations in materials and budget, the study was limited to a single site, with only one set of sensors per sinkhole slope position. The lack of replication prevented statistical comparisons between the different soil profiles or with other sinkholes. Nonetheless, the similar range of concentrations and seasonal patterns measured in all three profiles suggests that the sensors were adequately capturing CO₂ dynamics within the sinkhole. Measured CO₂ concentrations and calculated fluxes also tended to be similar between years, which indicates that sensor drift was not a major issue during the study.

Results of this study support previous work and also offer new insights into CO₂ gas movement within karst soils. For example, similar to previous work, we detected increasing soil CO₂ concentrations during the growing season, leading to peak upward fluxes in the middle of the growing season, followed by decreased concentrations during the dormant season. These seasonal patterns have been noted in many investigations in both karst and non-karst soils (e.g., Yang et al., 2012). However, our study provided additional unique findings. In particular, the use of multiple sensors allowed us to analyze both the magnitudes and directions of CO₂ diffusion in karst soils. Our data set included many periods in which measured CO₂ concentrations were lower at 60 cm compared to 20 cm depth. This condition implies the presence of a near-surface zero-flux plane for CO₂, similar to the zero-flux planes that can be seen for subsurface hydrologic gradients (e.g., Brutsaert, 2014). It is probable that downward gas diffusion occurred during these periods (assuming non-zero gas diffusivity), which is a process that has not been carefully considered prior to this study.

We posited that the underlying cave system acts as a primary CO₂ sink in the subsurface, especially during the cool season, when convective venting likely occurs. This process allows relatively rapid exchange between air in the cave and the external atmosphere (Mattey et al., 2016; Krajnc et al., 2017). The relatively low CO₂ concentrations that occur during these periods (Meyer et al., 2014; Kukuljan et al., 2021) would encourage downward diffusion through the soil profile. Future work could incorporate observations such as the ones collected in this study into biogeochemical models, and thereby better understand the relative importance of each of these mechanisms.

Finally, doline caves like James Cave occur throughout southwest Virginia and the Appalachian karst system. While our study was limited to a single site, it is likely that the findings of this study are generalizable to other karst areas with similar geologic and environmental factors. In particular, our study results suggest that CO₂ gas exchange within soil profiles above cave systems is highly dynamic over space and time, influenced by soil properties such as temperature, depth, water content, and organic matter, but also properties of the underlying cave. The magnitude of upward fluxes was determined to have similar seasonal trends between soil profiles, with most of the CO₂ efflux occurring during the warm growing season. The toe-slope profile had greater cumulative efflux than the back slope and shoulder, which was likely caused by connected air-filled pore space. Air-filled porosity was also important for downward gas diffusion, which was a much smaller but not negligible flux in the toe slope. Our study also emphasized the importance of quantifying soil hydraulic properties, such as air-filled porosity (ε), when determining subsurface CO₂ fluxes. This type of characterization should be used in future work to constrain terrestrial carbon budgets and quantify carbon storage within karst areas.

Funding for this work was provided the Global Change Center at Virginia Tech and by the Virginia Agricultural Experiment Station and the Hatch Program of the U.S. Department of Agriculture National Institute of Food and Agriculture (1026126). Support was also provided by the George Washington Carver program at Virginia Tech. We acknowledge Wil Orndorff, Daniel Doctor, and Benjamin Schwartz for discussions during the development of this project, and we thank Brandon Lester, Ayush Gyawali, Jesse Radolinski, Jinshi Jian, Caroline Wolcott, and Jonathan Grubb for their help with fieldwork and laboratory analyses.

Supplemental Material associated with this article can be found online at https://doi.org/10.2113/EEG-D-22-00082.