Carbon dioxide (CO2) gas diffusion is an important component of carbon cycling in soils. This process is particularly relevant in karst landscapes, which contain easily weathered rock, subsurface fractures, and cave networks. We instrumented three soil profiles—the shoulder, back slope, and toe slope of a sinkhole—above a karst cave in Virginia. Each profile had solid-state CO2 sensors and soil water content/temperature sensors at 20 and 60 cm depth that collected hourly measurements from 2017 to 2019. We calculated CO2 fluxes using Fick’s first law along with measured soil and assumed atmospheric CO2 concentrations. With this approach, we identified occasional near-surface zero-flux planes, in which CO2 likely diffused both upward and downward. All profiles had upward CO2 fluxes during warm-season months, with maximum fluxes of 1.2 μmol CO2 m−2 s−1 in the shoulder and back slope versus 2.0 μmol CO2 m−2 s−1 in the toe slope. During cool-season months, upward CO2 fluxes were smaller (0–0.3 μmol CO2 m−2 s−1) and were often counteracted by downward fluxes in the toe slope, possibly driven by ventilation into the underlying cave. The toe slope had a cumulative annual efflux of 14.5 mol CO2 m−2, which was >3 times greater than the other profiles. Fluxes were sensitive to soil porosity, with an order-of-magnitude difference when porosity was assumed to be 0.40 versus 0.56 cm3 cm−3. The results of this study offer new insight into short-term and seasonal variations in diffusive CO2 gas transport in karst soils, and they may inform other investigations of non-uniform diffusion processes.

Exchange of carbon dioxide (CO2) from the land surface to the atmosphere is an important process in the global carbon cycle, yet the controls and magnitudes of CO2 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 CO2 during chemical weathering of carbonate rocks. Much of the CO2 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 CO2 fluxes in karst soils is necessary to perform accurate regional and global carbon budgeting.

Subsurface transport of gases such as CO2 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, Ds, and (2) the concentration gradient, (i.e., ΔCx 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, Ds 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, CO2 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 CO2 concentrations near the surface are greater than those at depth (Chen, 2019). One possibility for these observations is CO2 consumption during dissolution of carbonate minerals (Chen, 2019; Zhao et al., 2021). A second possibility is that venting of cave air removes CO2 from deeper soil and epikarst; for example, Benavente et al. (2010) found that 30–50 percent of CO2 in a Mediterranean cave originated from the overlying soil profile, suggesting downward diffusion of CO2 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 CO2 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 CO2 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 Ds 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 CO2 gas movement in such systems have not yet been critically examined.

To evaluate the magnitudes and directions of CO2 fluxes in a karst landscape, we instrumented three locations along a sinkhole slope—shoulder, back slope, and toe slope—and measured hourly CO2 concentrations at 20 and 60 cm depths from winter 2017 to summer 2019. We used those measurements, along with assumed atmospheric CO2 concentration and data for soil physical properties (e.g., porosity, water content, temperature), to calculate CO2 fluxes. Our first objective was to compare upward and downward CO2 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 CO2 fluxes. Here, we hypothesized that CO2 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 CO2 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.

Site Description

The field site was located at James Cave in southwest Virginia (Supplemental Material Figure S1, 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 CO2 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 CO2 sensor failed at the back slope. The 60-cm-deep CO2 sensor failed at the shoulder location on 24 September 2018, and the 20-cm-deep CO2 sensors failed at the toe slope on 19 September 2019. All three profiles therefore had a common 1.5 year period in which both CO2 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.

Soil Characterization

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 CaCl2 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.

CO2 Flux Calculations

The measured CO2 concentrations were corrected for temperature and pressure based on the ideal gas law (Vaisala, 2017):
where Ccorrected is the CO2 concentration corrected for temperature and pressure [N N−1], Cmeasured is the measured CO2 concentration [N N−1], Pstd is the standard barometric pressure used in the instrument calibration (101.3 kPa), Tstd is the standard temperature (293 K), and R is the ideal gas law constant (8.31 × 10−3 kPa m3 mol−1 K−1).
We next calculated CO2 fluxes, J [N L−2 t−1], using Fick’s first law:
where ΔC is the difference in CO2 concentrations [N N−1] for two locations separated by a vertical distance Δx [L], and Ds is the diffusion coefficient [L2 t−1].
Ds was calculated following Tang et al. (2003):
where ε is the air-filled porosity [L3 L−3], ϕ is the soil porosity [L3 L−3], θ is the soil water content [L3 L−3], T is temperature [K], P is air pressure [M L−1 t−2], and Da0 is the gas diffusivity in pure air under standard temperature and pressure conditions [L2 t−1], assumed here as Da0 = 1.4 × 10−5 m2 s−1 (Pritchard and Currie, 1961). Porosity, ϕ, was determined based on the mean volumetric saturated water content (θs) from the soil cores collected in each profile.

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

We also quantified seasonal and annual CO2 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.

Sensitivity Analysis

We performed a brief sensitivity analysis to test how the porosity (ϕ) term in Eq. 3 and Eq. 4 affected calculations of Ds and Jupward. First, we calculated Ds through time for each profile using five assumed porosities: ϕ = 0.40, 0.44, 0.48, 0.52, and 0.56 cm3 cm−3. These porosities reflected the range of values measured in the soil samples used for soil characterization. The Ds 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 CO2 fluxes for each profile and then summed the cumulative CO2 efflux for the period between 1 April 2017 and 31 March 2018.

Soil Properties

The water retention curves were similar between most soil cores (Supplemental Material Figure S2, For example, at field capacity (i.e., pressure head = −340 cm), volumetric water content values were between 0.32 and 0.41 cm3 cm−3 for the toe-slope location, 0.33 and 0.39 cm3 cm−3 for the back-slope location, and 0.33 and 0.38 cm3 cm−3 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, Porosity varied between 0.41 and 0.55 cm3 cm−3 for the toe-slope location, 0.41 and 0.49 cm3 cm−3 for the back-slope location, and 0.40 and 0.47 cm3 cm−3 for the shoulder location. We therefore used mean values of ϕ = 0.48 cm3 cm−3 for the toe-slope profile, ϕ = 0.45 cm3 cm−3 for the back-slope profile, and ϕ = 0.44 cm3 cm−3 for the shoulder when calculating Ds. The macro-porosities varied from 0.06 to 0.15 cm3 cm−3, 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 cm3 cm−3, with similar sample-to-sample variation in all three profiles.

Environmental Conditions

The largest precipitation events occurred in the summer and early fall, with maximum rates of 1.0 to 2.5 cm hr−1 (Supplemental Material Figure S3a, Volumetric water contents were lower in the summer and fall months, with a range of 0.10 to 0.35 cm3 cm−3. 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 cm3 cm−3. 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, 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).

Gas Diffusivity Constant and CO2 Fluxes

Calculated values for gas diffusivity, Ds, 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 Ds 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 Ds values, while the back slope continued to have smaller values. In all three profiles, the Ds values were larger during the warm season and were often near 0 during the cool seasons.

The shoulder profile had relatively high CO2 concentrations early in the warm season, at times exceeding 55,000 ppmv (5.5 percent), and had relatively low CO2 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 CO2 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 CO2 fluxes and much smaller downward fluxes (Figure 3b). During the cool seasons, the profile had negligible fluxes in either direction. Upward CO2 fluxes peaked between 0.5 and 1.2 μmol CO2 m−2 s−1, while the magnitude of downward fluxes never exceeded 0.2 μmol CO2 m−2 s−1.

The back-slope profile had similar magnitudes and seasonal trends in CO2 concentrations as the shoulder, with CO2 concentrations peaking above 55,000 ppmv during the beginning of the warm season and then gradually declining (Figure 4a). The CO2 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 CO2 concentrations compared to the other two profiles (Figure 5a). Here, the 20 cm sensor often recorded greater CO2 concentrations than the 60 cm sensor, particularly during the cool, non-growing season. The toe-slope profile also tended to have larger CO2 fluxes than the other two profiles (Figure 5b). The CO2 fluxes reached peak values of 1.5–2.0 μmol CO2 m−2 s−1 (Jupward) and −0.5 μmol CO2 m−2 s−1 (Jdownward). Downward fluxes were often present during the cool season, with larger magnitudes during the winter of 2017–2018 as compared to 2018–2019.

Annual CO2 efflux for each profile was calculated by integrating the areas underneath the CO2 flux curves for each profile from 1 April 2017 through 31 March 2018. The annual efflux values ranged from 3.5 mol CO2 m−2 in the shoulder to 14.5 mol CO2 m−2 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 CO2 emissions. The cumulative downward flux was much smaller, with the toe slope having a cumulative total of −0.36 mol CO2 m−2, the back slope having a cumulative total of −0.13 mol CO2 m−2, and the shoulder having a cumulative total of −0.04 mol CO2 m−2.

Sensitivity Analysis

The sensitivity analysis showed that the magnitude of Ds 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 cm3 cm−3), the maximum Ds values, based on mean θ and T values from the 20 and 60 cm sensors, were between 1.0 and 1.2 mm2 s−1. The largest tested ϕ value (0.56 cm3 cm−3) gave a maximum Ds value of 2.8 mm2 s−1. In the cool season, when the soil water content was high, the variation in porosity drove the Ds value closer to zero for the minimum porosity value versus 0.1–0.2 mm2 s−1 when larger porosity values were used.

The differences in calculated Ds values translated to sizeable differences in the estimated amount of CO2 that was emitted from each soil profile over a 1 year period (Figure 7). The toe-slope location had the least variation, with total CO2 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 CO2 fluxes using solid-state sensors embedded in the soil.

CO2 Fluxes Increased during the Growing Season and Were Larger in the Sinkhole Toe Slope

We initially hypothesized that the upward CO2 fluxes would be similar in magnitude for the three profiles, whereas the shoulder would have the smallest downward CO2 fluxes, and the toe slope would have the largest downward fluxes. We also hypothesized that upward CO2 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 CO2 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 CO2 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−2 s−1). These differences primarily occurred during the warm growing season of the first study year, when the shoulder and back-slope profiles had smaller Ds estimates than the toe slope. In the second warm season, by comparison, the shoulder and toe slope had similar Ds values, and the two profiles had similar peak effluxes (1.3–1.6 μmol m−2 s−1). These comparisons indicate that CO2 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 CO2 concentrations compared to their deeper (60 cm) counterparts, implying the existence of a zero-flux CO2 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 CO2 m−2. While considerably smaller in magnitude than cumulative upward fluxes, which ranged from 3.6 to 14.5 mol CO2 m−2, these results show that downward CO2 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 CO2 production, CO2 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 CO2 (Meyer et al., 2014; Mattey et al., 2016). Under such conditions, elevated CO2 concentrations in the soil would result in downward diffusion, but only so long as Ds 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 cm3 cm−3 during the winter period compared to measured porosities of 0.41–0.55 cm3 cm−3. 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 CO2 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 CO2 sinks. For one, carbon dioxide can dissolve into percolating water and become converted into carbonic acid (H2CO3) and then bicarbonate (HCO3) during carbonate rock weathering (Ford and Williams, 2007). This process could reduce CO2 concentrations at depth as carbonates undergo dissolution. Previous work in the James Cave system found carbonate weathering to be highly seasonal, with increased CO2 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.

Sensitivity Analysis

The CO2 fluxes calculated in this study were sensitive to the porosity (ϕ) value used to calculate Ds, 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−2 s−1 when the highest porosity of ϕ = 0.56 cm3 cm−3 was used. Previous work has shown that ε is one of the most important terms for Ds (Jin and Jury, 1996; Lafond et al., 2011), as captured in the diffusion model used here (Eqs. 35).

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 cm3 cm−3). 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 CO2 efflux (Fan and Jones, 2014). These factors imply that the calculated CO2 fluxes in this study had considerable uncertainty, and they emphasize the importance of accurately constraining porosity when estimating CO2 fluxes using solid-state sensors embedded in the soil.

Other measurements and models for Ds 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 cm3 cm−3 (Supplemental Material Table S1,, 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 Ds values and gas diffusion rates.

Limitations of the Study

The study had several limitations. For one, many of the solid-state CO2 sensors reached maximum readings (∼55,000 ppmv) 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 CO2 efflux to the atmosphere.

We also did not directly measure CO2 concentrations in the cave, which limited our ability to fully understand the magnitude of downward CO2 diffusion. Other studies in similar cave systems have measured typical concentrations of <12,000 ppmv (Breecker et al., 2012; Meyer et al., 2014; Kukuljan et al., 2021), albeit with isolated observations of CO2 concentrations up to 30,000 ppmv (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 ppmv. It is therefore possible that a zero-flux CO2 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 CO2 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 CO2 dynamics within the sinkhole. Measured CO2 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 CO2 gas movement within karst soils. For example, similar to previous work, we detected increasing soil CO2 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 CO2 diffusion in karst soils. Our data set included many periods in which measured CO2 concentrations were lower at 60 cm compared to 20 cm depth. This condition implies the presence of a near-surface zero-flux plane for CO2, 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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

Supplementary data