Flux and stable isotope fractionation of CO₂ in a mesic prairie headwater stream

Headwater stream emission of carbon dioxide (CO₂) is a significant part of the short-term carbon cycle (Cole et al. 2007, 2010; Nadeau & Rains 2007; Battin et al. 2009; Alin et al. 2011; Butman & Raymond 2011; Striegl et al. 2012; Atkins et al. 2013; Dinsmore et al. 2013; Hotchkiss et al. 2015; Schade et al. 2016; Marx et al. 2017; Horgby et al. 2019; Li et al. 2021). Headwater streams are complex environments, accounting for more than half of all stream lengths globally (Nadeau & Rains 2007). The small-scale heterogeneity common in these systems makes it difficult to characterize and upscale CO₂ cycling and to quantify their contributions to atmospheric greenhouse gas concentrations. Further, recent research has highlighted the importance of karstic springs, and their contributions to headwater environments, as CO₂ sources to the atmosphere (Maas & Wicks 2017; Lee et al. 2021).

Carbon flux from streams depends on terrestrial production and storage of CO₂ in the soil and in shallow aquifers, and on the hydrogeologic pathways that connect them, in addition to in-stream CO₂ production (Bernal et al. 2022). Groundwater discharge is the primary source of dissolved gases in low-order streams (Hope et al. 2001; Doctor et al. 2008; Johnson et al. 2008; Sand-Jensen & Staehr 2012; Crawford et al. 2013; Hotchkiss et al. 2015; Deirmendjian et al. 2018), with shallow aquifers and interflow, not deep groundwater, the main contributors of CO₂ to streams (Hotchkiss et al. 2015). It is estimated that over 70% of this CO₂ stream flux is produced terrestrially, with the majority of the terrestrial contribution derived from infiltration through soil (Hotchkiss et al. 2015; Nosrati et al. 2020). Terrestrially derived CO₂ is produced during organic matter degradation, and root and organism respiration. Meteoric water infiltrates the soil, taking up soil CO₂ (Chapelle 2000). Although most soil CO₂ escapes upward to the atmosphere, downward water migration to groundwater dissolves and transports ∼1–2% of soil CO₂ (Hendry et al. 1993; Schlesinger & Lichter 2001; Tsypin & Macpherson 2012). Carbonic acid from the dissolution of CO₂ dissolves carbonate minerals, releasing inorganic carbon into the solution. These carbon fluxes lead to the accumulation of CO₂ in groundwater, with concentrations one to two orders of magnitude higher than in the atmosphere (Schlesinger & Melack 1981; Johnson et al. 2008; Macpherson et al. 2008; Macpherson 2009; Monger et al. 2015; Macpherson & Sullivan 2019). When groundwater discharges into surface water, gas efflux is driven by the concentration gradient between the water and atmosphere.

Given the extensive distribution of carbonate terrains across the world (Weary & Doctor 2014; Chen et al. 2017), the sensitivity of carbonate weathering to changes in temperature and water residence time (Sullivan et al. 2019), and the complexity of groundwater–surface water interactions, it is necessary to better quantify CO₂ efflux in carbonate-hosted headwaters. In addition, it is critical to relate the CO₂ flux to groundwater discharge to understand its contribution to global greenhouse gases. We present the results of a ∼1-year study of a headwater stream at the Konza Prairie Long-Term Ecological Research Site (LTER) and Biological Station (Konza), where merokarst (thin limestones alternating with shales) presents both the opportunity and challenge to characterize point sources of groundwater discharge (herein referred to as point sources) and resulting CO₂ efflux. We contrast the CO₂ fluxes from the stream water passing over outcrops or subcrops of limestone and shale, quantify the stream CO₂ emission within the lower half of the watershed (∼1.1 km stream length), and test how efflux is related to δ¹³C-CO₂.

Konza has a temperate, mid-continental climate (Hayden 1998), with high variability both seasonally and annually (Nippert & Knapp 2007). From 1983 to 2016 (Nippert 2017), the average daily air temperature at Konza was 12.8 °C, and the mean total annual precipitation was 839 mm. From July 2015 to August 2016, the period of this study, the average air temperature was 14.9 °C and precipitation was 974 mm (116% of the annual mean precipitation; Supplemental Table S1); 70% of precipitation fell during the growing season (March to October). On an annual basis, both 2015 and 2016 exceeded both the long-term average precipitation and air temperature.

Konza lies along the southeastern margin of the North American Prairie (Hayden 1998). Konza flora is dominated by C₄ grasses and forbs (C₃), while woody vegetation (C₃) is found in patches in riparian zones and on hillslopes (Veach et al. 2014). C₃ and C₄ plants differ in their method of carbon fixation during photosynthesis (Wang et al. 2012), which affects the δ¹³C of CO₂ produced during the breakdown of the plants.

Konza is underlain by thin soils and merokarst. The loess-based soils, primarily silty clays and silty clay loams (NRCS 2006), are thickest (∼1 m) at the base of slopes and patchy on the plateaus (Ranson et al. 1998). The merokarst is thin, interlayered beds of Lower Permian age limestones (mostly 1–2 m) and thicker shales (4–6 m; Supplemental Figure S1); discontinuous Quaternary alluvium occupies valleys. Bedrock strikes northeast-southwest and dips ∼2° northwest (Twiss 1988). Limestone outcrops form flat uplands and benches along the hillsides; slopes form over shales. Streams dissect the landscape into relatively steep-sided valleys (Macpherson 1996). The limestones and alluvium act as aquifers and shales as aquitards; many of the limestones are hydraulically connected to the streams (Macpherson et al. 2008; Hatley et al. 2023). Hydraulic conductivities of the limestones, controlled by secondary porosity, range over five orders of magnitude (10⁻⁸ to 10⁻³ m s⁻¹; Pomes 1995; Sullivan et al. 2020).

Geologic mapping, measurement sites, and sampling sites spanned 1,125 m of the South Fork of Kings Creek and were located with a Garmin Etrex Legend GPS. Stream segments were discretized, based on geology, into nine reaches (Figure 1(b)). Those underlain by shale and limestone account for 41 and 59% of the total stream length, respectively. Stream reaches underlain by shale are narrower than reaches underlain by limestone; for calculation purposes, widths of stream reaches underlain by shale were set at 1.0 m and those by limestones at 1.5 m. For the sake of brevity, units will be identified with their simplified geologic name: Crouse, Upper Eiss, Lower Eiss, Stearns, Morrill, and Florena.

To quantify the flux and fractionation of CO₂ in this system, measurements of stream discharge were made, water samples were collected and analyzed for basic water chemistry and for C isotopes, and CO₂ flux was measured in the stream. As this is an intermittent, flashy, stream, sampling was not possible when there was no water in the stream and also when the flows were so high as to be unsafe for sampling.

Stream discharge was measured using data from the triangular-throated weir, and point discharge measurements were taken at nine WP locations on two different days using a ∼30 cm long, 0.25 cm inside-diameter pitot tube. Streamflow measurements follow protocols outlined by the Environmental Protection Agency (Meals & Dressing 2008) and the United States Geological Survey (USGS 2016).

Throughout most of the year, temperature differences exist between groundwater (Macpherson 2020) and surface water (Brookfield et al. 2017; Nippert 2017; Nippert & Knapp 2007) at Konza. The average groundwater temperature was ∼16 °C for the study period, with the average surface water temperature of ∼21 °C. This temperature differential permitted the location of point-source groundwater discharge locations (point sources) by a FLIR^® T600 Thermal Imaging Infrared Camera (FLIR Camera). Seeps and springs flowing from fractures along the streambank and streambed were observed, and they were underlain by the Morrill, Eiss, and Crouse limestones. These point sources were used to evaluate degassing lengths.

Groundwater and stream water samples were collected over the study period. Groundwater was sampled from five observation wells: 3-5 Mor, 3-5-1 Mor, 3-6 Mor, 4-6 Eis1, and 4-6 Eis2 (Figure 1(a)). The results of chemical analysis (inorganic species and pH; Norwood 2020) were used to calculate aqueous pCO₂ by developing with PHREEQC Interactive 3.1.7–9213 (Charlton et al. 1997; Supplemental Table S2). Charge balances were also calculated using PHREEQC.

CO₂ efflux from the stream was measured with a suspended chamber (Crawford et al. 2014; Rawitch et al. 2019). Two chamber designs were lab tested (Norwood 2020) and used in the field. The first is a 3D-printed ABS plastic rectangular prism with rounded, triangular prisms on the front and back for streamlining; it has a small footprint (9.40 × 10⁻³m²) and was sealed with acetone vapor. The second is made with 5-mm thick plexiglass, is heavier, has the shape of a rectangular prism, and has a larger footprint (1.55 × 10⁻²m²).

For each trial, the bottom of the chamber was submerged ∼2 cm below the water surface, sealing the chamber to the water surface and preventing the influx of atmospheric CO₂. The chamber was connected to a pump and a Li-Cor LI-820 CO₂ Infrared Gas Analyzer using Tygon^® tubing; the circulating air (1 L/min) passed through a Drierite^® filter to remove moisture. CO₂ concentrations from the chamber were logged at 1-s intervals; pressure changes in the chamber were minimized by returning analyzed gas to the chamber. Chamber trials lasted at least 5 minutes, ensuring that the CO₂ emission rate became linear. Between trials, the chamber atmosphere was allowed to equilibrate with the atmosphere. Atmospheric CO₂ (ppm), temperature (°C), and relative humidity (RH, %) were measured using an AZ-77535 CO₂/temperature/RH meter.

Gas was sampled directly from the suspended chamber and stored in 0.5 L Tedlar^® Gas-Sampling Bags for C isotope determination. These samples were analyzed on a Picarro^® G2201-I Analyzer for Isotopic CO₂/CH₄ to determine δ¹³C-CO₂. Two reference standards for CO₂ gas (δ¹³C of −40.78‰ and −10.42‰, V-PDB) were used along with CO₂ gas standards (500 ppm and 1,000 ppm) for calibration. Samples were injected directly into the Picarro analyzer from the gas-sampling bags. Between analyses, nitrogen gas flushed the Picarro of CO₂.

The Konza LTER Program records stream discharge at the weir every 5 min; data are transformed into daily averages (Dodds 2021). The discharge measurements at the weir are assumed to be a reliable indicator of upstream flow in N04d; however, we also recorded discharge at upstream reaches on two dates when there was no flow at the weir. Over the study period, discharge was recorded at the weir 68% of the time. The mean discharge was 138 m³/day, and the median discharge was 402 m³/day (Supplemental Figure S2). The mean discharge and sum of daily discharge from 1 January 2016 to 2 August 2016 were almost twice as high as in 2015 (Supplemental Table S1).

Groundwater and surface water samples were collected over the 379 days between 21 July 2015 and 2 August 2016 (Supplemental Table S3). Charge balances are less than 5% for all but three analyses. Water chemistry is typical of limestone terrains, dominated by calcium and bicarbonate with slightly alkaline pH and moderate total dissolved solids. Analytical methods and results are similar to other studies done at the site (e.g., Macpherson 1996; Macpherson et al. 2008; Macpherson & Sullivan 2019). With the exception of seasonal variations in temperature and nitrate, very little spatial or temporal variability was observed at this site (Table 1) and is not further discussed.

Carbon flux measurements from the South Fork ranged from 2.2 to 214 g CO₂ m⁻² day⁻¹, with a mean CO₂ flux of 20.9 g CO₂ m⁻² day⁻¹ ± 41.4 g CO₂ m⁻² day⁻¹ (1 S.D.; Supplemental Table S4). WP-3A-SRC and WP-17 are high-flux outliers. For nearly all locations, flux measurements were higher for stream reaches underlain by limestone than those underlain by shale. The CO₂ fluxes measured in the three limestone units also varied, even between the more-permeable (Upper Eiss) and less-permeable (Lower Eiss) portions of the same limestone member (Table 2 and Supplemental Table S4). During the study period, the variability within a geologic unit, measured two to four times, was higher in high permeability units (Upper Eiss; Crouse [presumed high because of frequency of springs in this unit; Barry 2018]) than low permeability units (Lower Eiss, Stearns, Morrill). Among the limestones, the CO₂ fluxes and the coefficients of variation of the CO₂ fluxes decreased with unit thickness (Table 2). Point sources were not found in stream reaches underlain by shale.

Table 3 presents the mean CO₂ flux for each geologically constrained stream reach (stream reaches that correspond to only one geologic unit) during the 379-day study period. The average flux for all units was of the same order of magnitude with or without outliers: 14.5 ± 14.8 g CO₂ m⁻² day⁻¹ with and 10.2 ± 10.6 g CO₂ m⁻² day⁻¹ without outliers. Therefore, during the approximately 1-year study period, the South Fork emitted between 4.2 and 7.0 metric tons of CO₂ from the stream to the atmosphere.

The spatially and temporally variable CO₂ fluxes and carbon isotopes in this study reflect the geomorphologic and hydrologic heterogeneity of the South Fork, variable annual weather patterns, and vegetation. Point sources were only identified in reaches underlain by limestone and at geologic contacts between limestones and shales.

The spatial variation of CO₂ flux along the South Fork is attributed to point sources of groundwater discharge. In addition, the highest daily CO₂ emissions corresponded to the thickest (Crouse) and most permeable (Upper Eiss) limestones (Figure 3; Supplemental Table S4). Where measurements were made downstream of the point source, most locations showed a large decline in CO₂ efflux (Figure 3; Supplemental Table S4), reinforcing the complexity of characterizing the contribution of degassed CO₂. In some stream reaches, multiple point sources occur over short distances, further increasing the complexity of the degassing patterns (e.g., WP-1, Lower Eiss). These short distances are in contrast to the previous work that evaluated degassing over stream lengths of 100 s–1,000 s of meters (e.g., Mohammadi et al. 2020). In the Morrill reach, the discharge and degassing patterns may also be complicated by complex groundwater flow patterns. This includes groundwater moving upstream with respect to surface water flow, which then discharges below the weir, suggesting losing stream behavior (Sullivan et al. 2020, their supplemental information; Barry 2018); and small-scale gaining and losing stream behavior within the outcrop in the stream (Norwood 2020). The extensive outcrop of the Crouse in the upstream portion (Figure 1(b) and Table 3) has multiple point sources with different CO₂ fluxes (Supplemental Table S4): the highest CO₂ flux corresponds to a spring located just west of the stream and the second highest to the Crouse Limestone-Easly Creek Shale boundary, suggesting the importance of less-permeable units forcing groundwater discharge.

As a simplified analysis of this data, we consider the rate of changes in fluxes downstream of point sources in the same manner as one would consider the half-life of radioactive decay. We calculate these ‘decay rates’ using measurements provided at 1 and 2 m downstream of the point sources measured in the Eiss Limestone and assume a linear decay of the natural logarithms of the measurements, as the nature of the decay rate is not known. However, the measured data fit the linear functions with coefficients of determination of 0.96, 0.84, and 0.86, for August 2015, March 2016, and April 2016 measurements, respectively. We do not consider the Crouse Limestone due to the previously mentioned complexities in this thicker limestone unit. Using the values for the Eiss Limestone in Figure 4(b), the distance to reduce the CO₂ flux by half for August 2015, March 2016, and April 2016 were 0.70 m, 1.57 m, and 1.54 m, respectively. It appears that this rate is dependent upon the initial flux rate at the point-source location, where a higher initial flux results in a shorter distance to reduce that flux by half. This analysis is clearly limited by the few data points and measurement locations, but we feel it provides an interesting direction for future research.

Repeat measurements at the same locations demonstrate the variability in CO₂ flux at this site (Figure 3 and Table 2). Others have documented stream CO₂ content that may be strongly influenced by diurnal, seasonal, and annual factors (e.g., Hotchkiss et al. 2015; Marx et al. 2018). We propose that flux variability results from variations in hydraulic connectivity, as well as weather-driven recharge and groundwater flow directions, which is consistent with the recent work by Duvert et al. (2018). Table 2 shows the variability in measurements for all locations where at least two measurements were made. For measurements in limestones with lower hydraulic conductivity (Morrill, Lower Eiss), coefficients of variation of CO₂ flux are lower (1–53%) than those for limestones with higher hydraulic conductivity (Upper Eiss, Crouse; 52–155%).

The k₆₀₀ values ranged over two orders of magnitude, from ∼0.0002 to ∼0.02 m d⁻¹ (Supplemental Table S4). Excluding the highest k₆₀₀, which was measured at the Crouse spring (WP-17), the range lowers to one order of magnitude (∼0.0002 to ∼0.008), within the range of other k₆₀₀’s measured in other headwater streams with low stream velocities (e.g., Rawitch et al. 2021).

During the 2-m trials, most gas samples showed an inverse relationship between δ¹³C-CO₂ and CO₂ flux (Figures 3 and 4). Over the 2-m distance downstream of point sources, isotopic signatures are higher by ∼5–7‰ (Figure 4). Many have proposed that carbon isotope exchange with atmospheric CO₂ (δ¹³C ∼ −8‰) is partly responsible for increasing downstream δ¹³C values measured in streams and rivers (Taylor & Fox 1996; Yang et al. 1996; Atekwana & Krishnamurthy 1998; Karim & Veizer 2000; Hélie et al. 2002; Mayorga et al. 2005). However, isotopic equilibrium between dissolved CO₂ and atmospheric CO₂ can only be attained after the equilibrium of CO₂ concentration between atmosphere and stream/river water (Doctor et al. 2008). Thus, if a drive for CO₂ flux from water to the atmosphere exists, as observed in this study, fractionation of δ¹³C-CO₂ will be primarily affected by gas efflux, rather than the carbon isotope exchange that accompanies carbon mixing. The observed depletion in ¹²C during degassing therefore supports the hypothesis proposed by Doctor et al. (2008) and Venkiteswaran et al. (2014) that fractionation in the isotopic composition of DIC or CO₂ is driven by the process of gas transfer from stream to atmosphere, rather than mixing of CO₂ between the stream and atmosphere.

The largest number of measurements in a single geologic unit, point sources in the Upper Eiss Limestone (Figure 8), demonstrates variable weather and climate. Moisture, as reflected by precipitation and streamflow, is a primary control, as it affects the flux direction across the groundwater–surface water interface, residence times (reactions times), and plant functioning. Further, the general direction of groundwater flow at Konza has been shown to be opposite of the streamflow direction in the Morrill, but the same as streamflow in the Eiss (Sullivan et al. 2020). The variation in CO₂ flux and isotopic values measured in the Eiss further illustrate the complexity of characterizing CO₂ flux in headwater streams with: (1) thin limestone aquifers, with and without point sources; (2) large interannual variations in total meteoric precipitation and in timing of meteoric precipitation; and (3) large interannual variations in stream discharge. This is consistent with the study by Duvert et al. (2018) who found that CO₂ evasion is heavily influenced by spatial heterogeneities in the surface and subsurface. Recent work also highlights the influence of temporal variability in hydrogeologic and hydrologic conditions on CO₂ evasion, particularly through changes in the interactions between surface and subsurface and their respective responses to meteorology (Duvert et al. 2018; Marx et al. 2018).

The difference between CO₂ flux and δ¹³C-CO₂ collected from the Upper Eiss in 2015 and 2016 and the Crouse in 2016 illustrates the effect of hydrologic response to meteorology and limestone heterogeneity. In drier years, e.g., 2015 (preceded by an even drier year in 2014; Supplemental Table S5), groundwater residence time is longer and the buildup of CO₂ in groundwater is evidenced as higher CO₂ in point sources discharging to the stream. Higher water tables in wetter years reduce the distance between the surface and water table, with a wetter vadose zone increasing the unsaturated hydraulic conductivity, lowering groundwater residence time (Brookfield et al. 2017; Hatley et al. 2023).

This study examines the timing and extent of CO₂ transport from shallow aquifers to the atmosphere and supports findings from similar studies that headwater streams are significant contributors to local and regional carbon cycling. At Konza, 45 measurements of CO₂ efflux were taken at 18 different locations along a headwater stream. The spatial variability of CO₂ flux along the 1.1-km stream segment reflects the underlying merokarst geology and the rapid decrease in CO₂ flux downstream of point sources of discharge. Point sources, detected by water temperature, were only observed in reaches underlain by limestones; at point sources, CO₂ flux measured 2 m downstream from the point source ranged from 3 to 40% of the point-source flux, regardless of the magnitude of the point-source CO₂ flux.

The stable isotopic composition of CO₂ (δ¹³C-CO₂) was studied as a potential tracer of groundwater influx and predictor of CO₂ efflux rate. Lower δ¹³C-CO₂ values were often accompanied by larger CO₂ fluxes, but the inverse relationship is not predictive, likely because disequilibrium between stream CO₂ and atmospheric CO₂ is the main driver of efflux. Thus, δ¹³C-CO₂ can be a reliable indicator of groundwater discharge into a stream where a high contrast in partial pressure of CO₂ exists between groundwater and stream water, but it is not a reliable indicator of the CO₂ efflux rate.

The large range of CO₂ fluxes and δ¹³C-CO₂ values observed over small spatial extents reinforces the importance of point-source measurements in headwater streams, especially in those areas where karst, or other preferential flow paths, controls groundwater flow. Hence, large-scale investigations of stream degassing based on tracer tests or mass balance equations may overlook significant CO₂ contributions. The suspended chamber proved to be a simple and effective method for collecting measurements of CO₂ flux directly from the stream surface. This study provides additional reasons to consider shallow aquifers and headwater streams when accounting for carbon sinks and sources on local, regional, and global scales.

Climate and stream discharge datasets were provided by the Climate and Hydrology Database Projects (http://climhy.lternet.edu/), a partnership between the Long-Term Ecological Research program and the US Forest Service Pacific Northwest Research Station, Corvallis, Oregon. Significant funding for these data was provided by the National Science Foundation Long-Term Ecological Research program and the USDA Forest Service. B. S. N. and G. L. M. are grateful for support from the Konza Prairie LTER program (DEB-0823341 and DEB-1440484), the Geology Associates Fund of the KU Endowment Association, and the KU Department of Geology. A portion of coauthor PLS's time was supported by NSF EAR 2024388. B. S. N. thanks Mike Rawitch, Trevor Osorno, Emily Barry, Mackenzie Creamens, and Brooks Bailey for field assistance and consultation.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.