The Big Cypress National Preserve is a flat (mean surface slope = 0.03 m km⁻¹) karst landscape in southwestern Florida, and part of the Greater Everglades ecosystem (Figures 1A–C). The combination of low relief and a humid subtropical climate with mean annual precipitation of 1.33 m (McPherson, 1974; Shoemaker et al., 2011) results in significant water storage on the landscape. Thousands of small (< 2 ha) karst depressions dot the landscape in a remarkable regularly patterned landscape of distributed water storage (Quintero and Cohen, 2019). These depressions have far thicker soils (Watts et al., 2014) with far higher OC content (Zhang et al., 2019) than the surrounding uplands because they hold water for much of the year, thereby supporting hydrophytic vegetation such as pond cypress (Taxodium ascendens) and a variety of wetland taxa (e.g., Salix caroliniana, Cladium jamaicense). These cypress “domes,” so called because the cypress trees at the margins are typically shorter than those in the center due to growth rate differences (Katherine, 1988), are crucial elements in the complex mosaic that makes BICY an important biodiversity hotspot. The mosaic of uplands within which these depressions occur are dominated by cabbage palm (Sabal palmetto) and south Florida slash pine (Pinus elliottii var. densa) and, in stark contrast to the wetlands, are rarely inundated and lack any soil development.

Water levels in the wetlands follow pronounced seasonal patterns in rainfall and evapotranspiration; high water levels typically occur in the summer and fall, with drydown over the winter and the lowest water levels in late spring (Figure 1D). Hydrologic connectivity in this landscape occurs when the wetlands fill sufficiently to exceed a critical stage threshold where water spills out and generates landscape sheetflow (McLaughlin et al., 2019); below this spill threshold, wetlands hold water but do not appear to exchange water regionally. While construction of elevated roads and canals for logging, oil drilling, and recreation over the last century has resulted in localized hydrologic modification (McPherson, 1974), the hydrologic and biological conditions in BICY are far less impacted than in the adjacent Everglades system.

We studied four cypress domes across two contrasting lithologies within BICY. Two were in the region known as Raccoon Point, near the end of 11-Mile Road (RP1: 25.9896 N, 80.9279 W and RP3: 25.9144 N, 80.9392 W) on the Pleistocene Fort Thompson Formation, a more recent carbonate lithology with a high proportion of insoluble material. The RP domes were shallower, and lacked surface water during the dry season. The other two domes were located near Turner River Road (TR2: 26.1197 N, 81.2691 W, and TR3: 26.0208 N, 81.2661 W), located on the Miocene Tamiami Formation. The TR domes were deeper, and each wetland center was inundated during the entire study period. All sampling described below was performed during daylight hours, generally between 10:00 and 15:00.

The partial pressure of CH₄ (pCH₄) and stable isotopic composition (δ¹³C-CH₄) was measured in surface waters, bubbles collected from the sediment-water interface, and two porewater depths at three locations per wetland site—the edge of the wetland near the upland-wetland boundary, midway into the wetland, and the center of the wetland (Figure 1C and Figure 2). For each of these locations, nested porewater samplers were installed to the soil-bedrock interface and a depth halfway between the bedrock and sediment-water interface (Table 1). The samplers were constructed of 5.1 cm diameter schedule 80 PVC and had 10 cm long sections with 125 µm well screening at their base. Water was pumped from each sampler through an overflow cup which contained calibrated (daily) YSI ProPlus pH, DO, specific conductivity, and temperature probes. Water was pumped until all parameters stabilized, to ensure pristine porewaters were sampled.

Sampling for measurements of pCH₄ and δ¹³C-CH₄ values included overfilling 635 ml polycarbonate bottles for about three bottle volumes ensuring they contained no bubbles. The bottles were closed with butyl stoppers equipped with two stop-cocks and short/long straws to allow creation of a headspace by removing 60 ml of water and injecting 60 ml of N₂ to fill the headspace without contact with the atmosphere (Ward et al., 2016). The bottle was shaken vigorously for 2 min and the headspace was extracted into a 60 ml syringe.

After collecting the above samples and measuring CH₄ fluxes (described below), bubbles were collected from the sediment-water interface to measure non-microbially degraded δ¹³C-CH₄ values. Bubbles were collected by manually disturbing the sediment and capturing the released gas in an inverted funnel equipped with a 3-way luer valve (Sawakuchi et al., 2016). When the funnel filled with sediment bubbles, the gas was collected in 60 ml syringes. Prior to collecting the bubbles, care was taken to not disturb sediments as porewater were collected and flux measurements were made. δ¹³C-CH₄ data from the bubbles was used in our estimations of MOX (Estimations of CH ₄ Oxidation).

Gas samples were injected from the syringe into a Picarro G2201-i Cavity Ring-Down Spectrometer (CRDS) within ∼12 h of sampling to measure pCH₄ and δ¹³C-CH₄. Porewater and bubble samples were diluted 10–1,000 times with pure N₂ prior to analysis to achieve values within the CRDS’s detection range. pCH₄ values were corrected for dilution with N₂ based on the common gas law. pCH₄ values are reported in ppt units considering the high values recorded in porewaters and bubbles.

The flux of CH₄ was measured using various opaque chambers depending on the source of the CH₄ flux (Figure 2). For surface waters we used a “floating dome” design, which was an inverted high density polyethelene planter pot with polyethelene foam glued to its rim for flotation (Sawakuchi et al., 2017). The chamber was wrapped in tin foil to avoid heating of the headspace, which was sampled via an air-tight 2-way luer valve. The surface water chambers had a volume of 9.9 L, diameter of 0.33 m where the chamber meets the water, and surface area of 0.086 m² where the chamber meets the water. Five replicate chambers were deployed simultaneously for 5–6 min, after which the headspace was sampled using a 60 ml syringe. Ambient air was also sampled in the same manner.

When surface waters were not present, e.g., at the edge of the wetland during certain sampling periods, we utilized triplicate soil flux chambers constructed from white PVC pipe with a sealed endcap on one end and air-tight 2-way luer valve. Chambers were gently placed on the soil surface for 4–5 min and samples were collected in 60 ml syringes. The soil chambers had a volume of 2.0 L, diameter of 0.10 m where the chamber meets the soil, and a surface area of 0.008 m² where the chamber meets the soil.

Another chamber, made out of the same PVC materials used for the soil chambers, was designed to fit over cypress knees protruding from the water surface (Figure 2). At inundated sampling stations the chamber’s opening was placed just below the water surface so as to not touch the sediment-water interface for 5 min prior to sampling with a 60 ml syringe. One chamber was deployed in triplicate over one knee per sampling location when present; no cypress trees or knees were present at TR3. In order to calculate how much of the chamber’s volume was submerged under water and filled with the knee above the water-air interface, we measured how deep the chamber was submerged underwater, the knee height above the water surface, and the circumference at the bottom and top of the knee assuming the geometry of a tapered cylinder.

For each chamber design, a Licor-820 with external diaphragm pump was interfaced to one chamber to monitor CO₂ fluxes (not reported here) and to ensure that the chamber’s headspace was initially in equilibrium with the atmosphere and did not reach equilibrium with the surface waters/soils/knees during the deployment (i.e. gas fluxes remained linear). Water and air temperatures were measured using a YSI Pro Plus multiparameter sonde. Samples collected from the flux chambers were analyzed by direct injection on the CRDS as previously described. The rate of change ( d C H 4 d t ; atm L⁻¹ min⁻¹) was calculated as follows: d C H 4 d t = p C H 4 f   −   p C H 4 i t   (1) where p C H 4 f is the partial pressure of CH₄ inside the flux chamber at the end of the deployment (atm L⁻¹), p C H 4 i is the partial pressure of CH₄ in ambient air collected with a syringe (atm L⁻¹), and t is the length of time between placing the chamber on the water, soil, or knee and sampling the chamber’s headspace (minutes). CH₄ fluxes (F_CH4; μmol m⁻² min⁻¹) were then calculated based on the ideal gas law as follows: F CH 4 =   d C H 4 d t   x   V R  x  T  x  A     (2) where, V is the volume of the chamber in L units, R is the universal gas constant (0.082057 L atm °K⁻¹ mol⁻¹), T is temperature in °K, and A is the surface area of the chamber opening in m² units. We then converted F_CH4 values to mmol m² d⁻¹ for reporting. The above flux calculations were applied to our three chamber types to represent: (Eq. 1) the total methane flux from surface waters (F_T), which include ebullitive and diffusive fluxes, (Eq. 2) the total methane flux from soils (F_S), and (Eq. 3) the total methane flux from cypress knees (F_K) to the atmosphere.

We then calculated the relative contributions of diffusive vs. ebullitive fluxes (F_D and F_E, respectively) in the case of surface waters using gas transfer velocities similar to other river and lake studies (Bastviken et al., 2004, Bastviken et al., 2010; Sawakuchi et al., 2014). To do this, we first converted surface water and atmospheric pCH₄ from μatm L⁻¹ units to μmol L⁻¹ concentrations ([CH₄]_SW and [CH₄]_ATM, respectively) based on the temperature dependent Henry’s Law constant (K_H) (Wiesenburg and Guinasso, 1979): K H =   e [ − 68.8862   +   ( 101.4956   x   100 T )   +   ( 28.7314   x ⁡ ln ( T 100 ) ) ] (3) where T is temperature in °C units. Gas transfer velocity (k _CH4) was then calculated in m d⁻¹ units for each deployed chamber as: k CH 4 = F T [ CH 4 ] SW   −   [ CH 4 ] ATM (4) k _CH4 was then converted to k ₆₀₀ based on the following equation (Wanninkhof, 1992; Wanninkhof, 2014): k 600 = ( 600 Sc CH 4 ) 0.5   x   k CH 4 (5) where Sc_CH4 is the temperature-dependent Schmidt number for methane calculated as follows: Sc CH 4 = 1897.8 − 114.28   x   T + 3.2902   x   T 2 − 0.039061   x   T 3 (6) where T is temperature in °C. We then calculated the ratio of k ₆₀₀ in each individual chamber to the minimum k ₆₀₀ observed per deployment of five chambers (k ₆₀₀/k _600-min) to determine which chambers captured bubbles from ebullition (Bastviken et al., 2004; Bastviken et al., 2010; Sawakuchi et al., 2014). For example, a high k ₆₀₀/k _600-min value indicates a large contribution of ebullitive vs. diffusive methane fluxes. The frequency distribution of k ₆₀₀/k _600-min for all chamber deployments indicates two distinct groups of ratios, whereby a ratio of 1.75 is the point of distinction (Figure 3). We chose the threshold value of 1.75 based on change point analysis (described in Data Analysis). Thus, we assumed that chambers with a k ₆₀₀/k _600-min greater than 1.75 experienced ebullition. We first calculated the average rate of CH₄ diffusion for chambers with a k ₆₀₀/k _600-min below 1.75. We then calculated ebullition as the difference between total and diffusive fluxes. For comparison across the study domain, we calculated the average ebullition rate for all five simultaneously deployed chambers, using a value of zero for chambers with a k ₆₀₀/k _600-min below 1.75.

Microbial oxidation of CH₄ in either oxic shallow sediment layers or surface waters results in an increase in δ¹³C-CH₄ values relative to its origin (Bastviken et al., 2002). Thus, we calculated the extent to which MOX limits diffusive surface water methane emissions by comparing surface water δ¹³C-CH₄ to both bubble and porewater δ¹³C-CH₄. We applied two different models to assess the largest range of uncertainty. We used a model that represents an open system at steady state (Eq. 7) (Happell et al., 1994) and a second Rayleigh model for closed systems (Eq. 8) (Sawakuchi et al., 2016) to calculate the fraction of CH₄ oxidized in shallow sediments and surface waters (f): f = δ SW −   δ b ( a − 1 )   x   1000 (7) ln ( 1 − f ) = ln ( δ b + 1000 ) −   ln ( δ SW + 1000 ) a − 1 (8) where δ_SW is the δ¹³C-CH₄ of surface waters, δ_b is the δ¹³C-CH₄ of bubbles and/or porewater, and a is the isotopic fractionation factor. Each calculation was made using two fractionation factor values—1.025 and 1.033—representative of the range of literature values available for similar systems (Tyler et al., 1997; Zhang et al., 2013). We also made the calculations using δ¹³C-CH₄ from bubbles and both porewater depths considering it was not always logistically possible to collect bubbles in shallow waters. We present MOX estimations as the average ± 1 SE of all calculations (i.e., two equations, two a values, and one or more belowground δ¹³C-CH₄ values).

We estimated total annual CH₄ emissions from each wetland in order to evaluate 1) the contribution of each flux pathway to wetland-scale CH₄ emissions and 2) the sustained global warming potential impact of CH₄ emissions compared to OC burial in the wetland sediments. The hydrology and topographic setting of the studied wetlands are extensively described by McLaughlin et al. (2019). In short, we measured water levels continuously for the entire study period using pressure transducers (Solinst Gold Levelogger) deployed in shallow groundwater wells situated in the deepest part of each wetland. Barometric pressure correction was performed using data collected from barometric pressure transducers (Solinst Barologger) deployed in a dry well adjacent to each water level recorder. LIDAR digital elevation models were provided by the National Center for Airborne Laser Mapping, processing of the LIDAR digital elevation models are described in detail in Quintero and Cohen (2019). Elevations are reported with an estimated mean accuracy of 5 cm across each domain, and with a spatial resolution of 5 m. The surface area of the center, intermediate, and edge sections of each wetland were determined based on elevation thresholds that aligned with shifts in vegetation coverage and inundation extent. Edge delineations capture the entirety of the wetland depression within the delineated extent while making sure to exclude any upland area. Intermediate delineations capture the bald/pond cypress extent of each wetland while making sure to exclude dwarf cypress and pine at higher elevations. Center delineations capture the center communities in our domes, these include open water areas, which typically differ from the dominant cypress community.

Surface water CH₄ fluxes were multiplied by the percentage of the year that each section was inundated and soil CH₄ fluxes were multiplied by the percentage of the year that was dry. To scale the knee surface area we used an available literature value for the density of knees (0.315 knees m⁻²) in a similar cypress swamp setting since knee density was not measured in this study (Pulliam, 1992). Pulliam (1992) determined knee density across two 1,250 m transects in the frequently inundated lower floodplain of the Ogeechee River with 1 m² plots established every 8m (i.e., 156 plots per transect).We then multiplied this knee density by the average surface area of knees that we made flux measurements on (0.57 m² knee⁻¹), and finally multiplied by the landscape surface area of each compartment. TR3 had no cypress trees or knees, whereas the other wetlands had knees in qualitatively similar abundance from the center to edge. Due to logistical constraints (i.e., time constraints), knee CH₄ flux measurements were not able to be made in each wetland section, thus, in those cases the average flux rate for a given wetland was applied to each of its sections.

Carbon burial rates were estimated for the center of each wetland based on previously reported 210Pb sediment accumulation rates, which were not measured for the intermediate and edge wetland sections (Zhang et al., 2019). 210Pb dating was performed down to 42 cm at the RP3 wetland and we applied this same rate to the other three wetlands. OC abundance in sediments (%OC) and bulk density were measured every 1 cm (Zhang et al., 2019). We used the average %OC and bulk density down to 42 cm to match the 210Pb time frame for carbon burial estimates. Carbon burial rates were determined as follows: OC Burial = % OC   x   Bulk Density   x   S ediment Accumulation Rate   x   Surface Area (9)

The ratio of methane emissions to OC burial was calculated to compare their relative rates. This ratio was then multiplied by the sustained global warming potential of CH₄ for 20, 100, and 500 yr time frames (Neubauer and Megonigal, 2015) to determine the net balance in atmospheric radiative forcing from CH₄ emissions and CO₂ uptake by vegetation and subsequent burial as OC in sediments.

All statistical analyses were performed in the statistical computing language R using R Studio version 1.3.1093 (RStudio Team, 2020). Calculations of the k ₆₀₀/k _600-min threshold for ebullition to occur described in CH ₄ fluxes From Surface Water, Soils, and Cypress Knees were performed using the “changepoint” R package. The changes in mean (cpt.mean) function was used with the At Most One Change (AMOC) method and the Modified Bayesian information criterion (MBIC) penalty. The change point of k ₆₀₀/k _600-min was determined to be 1.75 with a confidence value of 1.00.

Field data was determined to not have a normal distribution based on visual inspection of density and QQ plots. Thus, we chose the non-parametric unpaired Wilcoxon test to evaluate the statistical significance of differences in measured parameters (e.g., pCH₄, δ¹³C-CH₄, CH₄ fluxes, and MOX) between wetland sites, sampling stations (i.e., center, intermediate, and edge), measurement types (e.g., soil vs. surface water fluxes), and sampling season. For example, we used the Wilcoxon test to determine if a given flux pathway was significantly higher at the site with the highest average flux rate compared to each of the other three sites individually. Thus, reported p values represent comparisons of one group to another (e.g., TR2 vs. RP1, TR2 vs. RP3, etc.). Significant differences were considered to fall within a 95% confidence interval (i.e., p < 0.05). In cases where all comparisons yielded insignificant differences (e.g., no wetland sites were statistically distinct from one another) we report the range of computed p values as (maximum p value < p > minimum p value) for the sake of brevity.

We performed linear regressions to determine correlations between measured parameters (i.e., pCH₄ vs. δ¹³C-CH₄, fluxes vs. temperature, and fluxes vs. wetland stage) and report R² values.

All data is reported as the average ± standard error of the mean. This standard error was propagated through our estimations of total wetland methane emissions. The distribution and variability of field data is shown in box plots (Figures 4–6).

Methane measurements were made during two dry seasons (May 2015 and 2016) and one wet season (October 2015). Wetland stage was lower immediately before and after the May 2015 sampling, whereas wetland stage was relatively high during the May 2016 dry season compared to 2015 and 2017 (Figure 1D). The partial pressure of CH₄ in surface waters varied from 0.148 to 4.53 ppt with an average of 0.805 ± 0.128 ppt throughout the study period (Figures 4A–C). The δ¹³C-CH₄ of surface waters ranged from -69.2 to -42.0‰ with an average of -50.3 ± 0.8‰ throughout the study period (Figures 4D–F). Surface water pCH₄ was highest on average at RP1 (1.31 ± 0.128 ppt) compared to RP3 (1.00 ± 0.284 ppt; p = 0.60), TR3 (0.654 ± 0.239 ppt; p = 0.06), and TR2 (0.477 ± 0.068 ppt; p = 0.05) (Figure 4A). The most enriched average δ¹³C-CH₄ values (−46.7 ± 0.7‰) were observed TR2, which also had the lowest pCH₄ and this difference in δ¹³C-CH₄ was significant (p = 0.00 when comparing TR2 to each wetland individually) (Figure 4D). Average δ¹³C-CH₄ values were more similar to each other at RP1 (−50.2 ± 0.5‰), RP3 (−53.7 ± 1.8‰), and TR3 (−53.0 ± 2.6‰) with no statistically significant differences (0.79 < p > 0.16 when comparing each wetland to one another). There was also no significant difference between the center, midway point, or edge of the wetlands for surface water pCH₄ (0.79 < p > 0.16) and δ¹³C-CH₄ (0.74 < p > 0.58) (Figures 4B,E). Finally, there were no statistically significant differences in surface water pCH₄ (0.66 < p > 0.28) or δ¹³C-CH₄ (0.45 < p > 0.26) when comparing each sampling season (Figures 4C,F).

Porewater CH₄ concentrations ranged from 0.748 to 75.2 ppt with an average of 16.1 ± 1.85 ppt. The δ¹³C-CH₄ of porewaters ranged from −69.1‰ to −43.3‰ with an average of −57.0 ± 0.6‰ throughout the study period. Average pCH₄ in pore waters was 20 times higher than surface waters (p = 0.00) and porewater δ¹³C-CH₄ was more 6.7‰ more depleted than surface waters (p = 0.00) (Figure 4). Porewater collected at mid-depth had an average pCH₄ that was 1.1 times higher than at the bedrock interface, but this difference was not significant (p = 0.35). Likewise, average porewater δ¹³C-CH₄ at mid-depth (-57.7 ± 0.7‰) was slightly more depleted than at the bedrock interface (-56.5 ± 0.9‰; p = 0.28), but this difference was not significant.

Seasonally, the highest average porewater pCH₄ occurred during the October 2015 high wetland stage period (17.9 ± 2.32 ppt) compared to the low wetland stage May 2015 period (14.4 ± 7.35 ppt; p = 0.03) and moderate wetland stage May 2016 period (13.2 ± 2.63 ppt; p = 0.07) (Figure 4C). δ¹³C-CH₄ was the most depleted in May 2016 (−59.1 ± 1.3‰) compared to May 2015 (−55.7 ± 1.0‰; p = 0.09) and October 2015 (−56.3 ± 0.8‰; p = 0.07) (Figure 4F). Porewater pCH₄ at the bedrock interface was highest at TR3 (28.0 ± 3.75 ppt) compared to TR2 (16.0 ± 2.82 ppt; p = 0.02), RP3 (14.3 ± 3.65 ppt; p = 0.02), and RP1 (16.0 ± 3.65 ppt; p = 0.00) (Figure 4A). Following the same trend, δ¹³C-CH₄ at the bedrock interface was 9.0–10.5‰ more depleted at TR3 compared to the other wetlands (p = 0.00) (Figure 4D). Porewater pCH₄ at mid-depth was even more elevated at TR3 (46.7 ± 7.92 ppt) compared to RP3 (11.1 ± 2.18 ppt; p = 0.00), TR2 (8.30 ± 1.40 ppt; p = 0.00), and RP1 (3.29 ± 0.917 ppt; p = 0.00) (Figure 4A). The difference in δ¹³C-CH₄ between sites was not as large for mid-depth porewaters; TR3 was 2.7–5.2‰ more depleted at TR3 compared to the other wetlands (0.23 < p > 0.00). In general, pCH₄ was not correlated with δ¹³C-CH₄ for either surface or porewaters (R² = 0.12 and 0.03, respectively), except for surface waters at TR3, where pCH₄ was negatively correlated with δ¹³C-CH₄ (R² = 0.89).

The pCH₄ of bubbles collected from the sediment-surface water interface ranged from 55.0 to 1,000 ppt (i.e., pure CH₄) for all sites with an average value of 376 ± 71.1 ppt. The average δ¹³C-CH₄ of bubbles was −61.3 ± 0.68‰, which was 4.3% more depleted than porewaters (p = 0.00). Because bubbles were not always able to be collected (e.g., when surface water level was too low for the sampling funnel), we did not attempt statistical comparisons between wetland sites, seasons, or sampling stations.

The total flux rate of CH₄ from surface waters to the atmosphere (F_T) for all floating chamber deployments ranged from 0.17 to 137 mmol m⁻² d⁻¹ with an average of 11.3 ± 1.45 mmol m⁻² d⁻¹. The highest F_T rates were observed in May 2016 (16.2 ± 3.22 mmol m⁻² d⁻¹) compared to May 2015 (11.4 ± 3.22 mmol m⁻² d⁻¹; p = 0.38) and October 2015 (8.19 ± 1.39 mmol m⁻² d⁻¹; p = 0.14) (Figure 5C). The TR2 wetland had the highest F_T rates (17.3 ± 3.86 mmol m⁻² d⁻¹) compared to RP1 (12.0 ± 2.91 mmol m⁻² d⁻¹; p = 0.77), TR3 (10.7 ± 2.48 mmol m⁻² d⁻¹; p = 0.65), and RP3 (3.95 ± 0.64 mmol m⁻² d⁻¹; p = 0.05) (Figure 5A).

The gas transfer coefficient (k ₆₀₀) for surface water fluxes ranged from 0.002 to 18.0 m d⁻¹ (average = 0.711 ± 0.14 m d⁻¹). The k ₆₀₀/k _600-min ratio for each deployment of five chambers reached as high as 148 in the case of high ebullition, with an average of 4.82 ± 0.91 (Figure 3). Based on k ₆₀₀/k _600-min ratio, we separated out F_T into surface water fluxes derived from diffusion and ebullition (F_D and F_E, respectively) for each deployment. Average CH₄ ebullition rates for all four wetlands combined (7.79 ± 1.37 mmol m⁻² d⁻¹) were 2.2 times higher than diffusion rates (3.50 ± 0.22 mmol m⁻² d⁻¹; p = 0.00); F_T accounted for 69% of the total surface water CH₄ flux rate on average. Similar to F_T, the ebullitive flux F_E was highest in May 2016 (12.5 ± 3.01 mmol m⁻² d⁻¹) compared to May 2015 (6.87 ± 3.15 mmol m⁻² d⁻¹; p = 0.38) and October 2015 (4.92 ± 1.33 mmol m⁻² d⁻¹; p = 0.03) (Figure 5C). In contrast, the diffusive surface water flux F_D was highest during the May 2016 sampling (4.57 ± 0.44 mmol m⁻² d⁻¹) compared to May 2015 (3.72 ± 0.41 mmol m⁻² d⁻¹; p = 0.03) and October 2016 (3.26 ± 0.28 mmol m⁻² d⁻¹; p = 0.03) (Figure 5C).

Across sites, F_D was highest at RP1 (4.72 ± 0.63 mmol m⁻² d⁻¹) compared to TR3 (3.68 ± 0.43 mmol m⁻² d⁻¹; p = 0.14), TR2 (3.21 ± 0.27 mmol m⁻² d⁻¹; p = 0.69), and RP3 (2.31 ± 0.32 mmol m⁻² d⁻¹; p = 0.00) (Figure 5A). In contrast, F_E was highest at TR2 (14.0 ± 3.75 mmol m⁻² d⁻¹) compared to RP1 (12.0 ± 2.91 mmol m⁻² d⁻¹; p = 0.14), TR3 (7.04 ± 2.48 mmol m⁻² d⁻¹; p = 0.34), and RP3 (3.95 ± 0.64 mmol m⁻² d⁻¹; p = 0.08). F_T at the wetland edge stations were 1.6 and 2.6 times higher than the center (p = 0.39) and intermediate stations (p = 0.00), respectively (Figure 5B). Likewise, F_D at the wetland edge stations were 1.5 and 3.9 times higher than the center (p = 0.05) and intermediate stations (p = 0.00), respectively. FE was also 1.6 and 2.3 times higher at the edge stations compared to center (p = 0.32) and intermediate stations (p = 0.73), respectively, but this difference was less statistically significant likely due to higher variability (Figure 5B).

The average flux rate of CH₄ from soils (F_S) that were not inundated during the time of sampling (18.4 ± 5.14 mmol m⁻² d⁻¹) was 1.6 times higher than the total flux from surface waters, but this difference was not statistically significant (p = 0.77). However, there was a more significant difference between F_S and F_D and F_E when considered separately (p = 0.06 and 0.00, respectively). F_S was substantially higher in May 2016 (34.8 ± 8.31 mmol m⁻² d⁻¹) compared to both May 2015 (3.56 ± 1.01 mmol m⁻² d⁻¹; p = 0.02) when water levels were substantially lower and October 2015 when only non-inundated soils were present at RP1 (2.85 ± 1.69 mmol m⁻² d⁻¹; p = 0.29) (Figure 5C). F_S was substantially higher at TR3 (44.8 ± 13.9 mmol m⁻² d⁻¹) compared to RP1 (18.7 ± 7.01 mmol m⁻² d⁻¹; p = 0.10), RP3 (4.41 ± 1.39 mmol m⁻² d⁻¹; p = 0.03), and TR2 (0.718 ± 0.035 mmol m⁻² d⁻¹; p = 0.06) (Figure 5A).

The average flux of CH₄ from cypress knees (F_K; 42.0 ± 6.33 mmol m⁻² d⁻¹) was 3.7 and 2.3 times higher than F_T (p = 0.00) and F_S (p = 0.00), respectively. F_K was 2.1 times higher in May 2016 (58.2 ± 10.8 mmol m⁻² d⁻¹) compared to October 2016 (27.6 ± 5.48 mmol m⁻² d⁻¹; p = 0.01), and was not sampled in May 2015 (Figure 5C). F_K was highest at TR2 (65.1 ± 8.46 mmol m⁻² d⁻¹) compared to RP1 (33.5 ± 7.31 mmol m⁻² d⁻¹; p = 0.01) and RP3 (31.0 ± 11.2 mmol m⁻² d⁻¹; p = 0.00) (Figure 5A). Cypress knees are not present at TR3. The higher cypress knee flux rates at TR2 may be linked to elevated porewater pCH₄ at that site (Figure 4A). The highest average F_K rates were observed at intermediate sampling stations (54.7 ± 7.76 mmol m⁻² d⁻¹) compared to edge stations (45.9 ± 13.1 mmol m⁻² d⁻¹; p = 0.18) and center stations (12.1 ± 4.67 mmol m⁻² d⁻¹; p = 0.00) (Figure 5B), but there was not a significant difference between porewater concentrations by station (0.99 < p > 0.36) (Figure 4B). CH₄ fluxes were not correlated to either temperature or wetland stage for F_T, F_D, F_E, F_S, or F_K (0.17 < R² > 0.00).

On average, we estimate that 26 ± 3% of the CH₄ diffused from soil porewaters into surface waters is oxidized to CO₂ before evasion to atmosphere, resulting in an observed surface water diffusive CH₄ flux (i.e., F_D) that is 26% lower than the rate of CH₄ input to surface waters from sediments. The extent of MOX was higher in the TR wetlands than the RP wetlands (Figure 6A); TR2 had the highest extent of MOX on average (40 ± 5%), which was 1.4 times higher than TR3 (28 ± 7%; p = 0.46), 2.0 times higher than RP1 (20 ± 5%; p = 0.01), and 2.7 times higher than RP3 (15 ± 7%; p = 0.04). There was no significant variability between sampling period for the wetlands when all wetlands were considered together (1.00 < p > 0.37) (Figure 6C). MOX was highest in the center sampling stations (35 ± 6%) compared to intermediate stations (25 ± 5%; p = 0.17) and edge stations (18 ± 7%; p = 0.11) (Figure 6B), but these differences were not significant.

The low elevation center of each dome comprised 3–12% of the surface area of the studied wetlands. The intermediate and edge (i.e., mid to high elevation) sections comprised of 25–47% and 49–88% of each wetland, respectively (Table 2). The center, intermediate, and edge components were non-inundated (i.e., no surface water) 9–17%, 15–27%, and 16–46% of the study period, respectively. The total flux of CH₄ from each wetland ranged from 239 ± 70 kg C yr⁻¹ for the small RP3 wetland to 1741 ± 220 kg C yr⁻¹ for the larger RP1 wetland. Diffusion and ebullition of methane from surface waters to the atmosphere contributed to 14 ± 2% (range = 8–17%) and 25 ± 6% (range = 12–34%) of the total spatially-scaled flux estimates from all four wetlands, respectively (Table 2). Methane emissions from non-inundated soils contributed to 34 ± 10% of the total spatially-scaled flux estimates from all four wetlands, but had a much larger range (1–71%) for each individual wetland compared to surface water fluxes (Table 2). Cypress knees contributed to 26 ± 5% of the total spatially-scaled flux estimates from all four wetlands. Knees were not present at TR3. When only considering the other three wetlands, knees accounted for 39 ± 5% of the total scaled methane fluxes with a range of 24–75% for each individual wetland. These totals do not include fluxes from the stems and leaves of the cypress trees, which is likely to increase the total.

OC burial rates over the last century were estimated for the dome centers based on a sediment accretion rate of 0.45 cm yr⁻¹ (Zhang et al., 2019). %OC averaged across the top 42 cm ranged from 13.5% to 19.1% for each wetland center and bulk density ranged from 1.20 to 1.37 g cm⁻³ (Table 3). Carbon burial rates ranged from 160 kg C yr⁻¹ for the small RP3 wetland to 1,263 kg C yr⁻¹ for TR3. The molar ratio of CH₄ emissions to OC burial ranged from 0.03 to 0.14, or, in other words, spatially-scaled OC burial rates were 7–36 times greater than total spatially-scaled CH₄ emissions from the wetland centers. When ecosystem CH₄ emissions are sustained (as opposed to a one-time pulse), CH₄ has a sustained global warming potential of 96, 45, and 11 times that of CO₂ on 20 yr, 100 yr, and 500 yr timescales (Neubauer and Megonigal, 2015). Over the 20 and 100 yr time frames, the wetland centers of all four wetland sites were a net source of warming potential to the atmosphere; i.e., the methane they emitted to the atmosphere had a larger greenhouse gas impact than the CO₂ removed from the atmosphere by vegetation that gets buried in sediments as OC (Table 3). The ratio of methane warming potential to carbon burial ranged from 2.7–13.3 and 1.3–6.3 for 20 yr and 100 yr time frames, respectively (Table 3). Over the 500 yr time frame the wetland centers at three wetland sites were net sinks of warming potential (i.e., OC burial outweighed the impact of CH₄ emissions) with the ratio of methane warming potential to carbon burial ranging from 0.3–0.9. TR2 was a net source of warming potential even over 500 yr time scales with a warming potential ration of 1.5 (Table 3).

Our results show that the total flux of CH₄ from cypress dome surface waters was not significantly different from non-inundated soil CH₄ fluxes on the wetland edge or during periods of low water levels (Figure 5). Likewise, our observed methane flux rates from soils and surface waters (diffusion + ebullition) were nearly equivalent to other chamber-based measurements made in similar subtropical wetlands in Florida such as Corkscrew Swamp, FL (11.2 mmol m² d⁻¹ compared to 11.3 mmol m² d⁻¹ in this study) (Villa and Mitsch, 2014; Pereyra and Mitsch, 2018). The deep slough of the Corkscrew Swamp riverine cypress strand, which accumulates peat, is most similar to the cypress domes we studied, whereas the other four topographic settings studied by Villa and Mitsch (2014) were more similar to the unstudied upland settings surrounding our study domain and had an order of magnitude lower flux rates. Our average flux rates for surface water diffusion were also comparable to average rates reported for subtropical wetlands in a synthesis of global wetland studies, whereas our observed flux rates for total surface water emissions (diffusion + ebullition) was 2.8 times higher, soil emissions were 4.6 times higher, and knee emissions were 10.5 times higher (Turetsky et al., 2014). Only two subtropical wetlands were included in this synthesis, highlighting the paucity of comprehensive GHG data needed to constrain global wetland CH₄ budgets, which represent nearly half of all natural sources of CH₄ to the atmosphere (Saunois et al., 2020).

Cypress knees were the most rapid pathway for methane release in our study with flux rates 3.7 and 2.3 times higher than total surface water and soil flux rates, respectively (Figure 5). The highest knee flux rates were observed at the wetland site with the highest porewater pCH₄ (TR2) suggesting that subsurface methane content is an important driver of emissions through emergent root structures and perhaps also tree stems (Barba et al., 2019a), which were not studied here. Cypress knees represent a relatively small feature of the landscape in a two-dimensional plane, however, unlike the soil and water surface, gas exchange occurs from its entire three-dimensional structure. This attribute, along with the high measured flux rates, scales to represent a substantial percentage (39 ± 5%) of wetland CH₄ emissions at the three studied wetland sites with cypress trees present. This potentially important emission pathway has only been addressed by several previous studies. Observations made in the floodplain of a southeastern, USA blackwater river concluded that cypress knees contributed to less than 1% of the total methane emissions from the floodplain (Pulliam, 1992). The small contribution of cypress knees compared to floodplain surface CH₄ emissions in this case was likely due to substantially lower porewater CH₄ concentrations that resulted in average knee CH₄ flux rates that were ∼45 times lower than our measured rates. Knee CH₄ flux rates measured in the Sabine River floodplain in S.E. Texas were similar in magnitude to our measurements and it was reported that vegetation contributed to 64–96% of total ecosystem methane release (Bianchi et al., 1996). In that case the measured porewater CH₄ concentrations at depths greater than ∼10 cm were the same order of magnitude as our observations, whereas surface soils had nearly undetectable CH₄ levels, and surface waters had measurable rates of CH₄ oxidation.

Understanding the mechanisms that drive methane cycling through each pathway is a prerequisite for determining how effectively wetlands act as sources or sinks of GHGs and how they feedback and interact with changes to the Earth system. Water level is considered to be the primary factor that controls soil methane emissions in wetlands, particularly in tropical and subtropical wetlands that have wet and dry seasons (Whalen, 2005; Mitsch et al., 2010). Mitsch et al. (2010) observed that tropical wetlands with seasonal pulses in water level (e.g., monsoon cycles) had higher total annual methane emissions than permanently flooded wetlands, and that the highest seasonal emissions occurred when water levels were between 15–30 cm above the ground surface. This relationship to depth has also been observed in ombrotrophic peatlands and in both cases was attributed to enhanced MOX in shallow soils when water levels are low enough to allow soil aeration and enhanced MOX in surface waters when surface water depths are greater (Jauhiainen et al., 2005), which aligns with our observations of higher MOX in the wetland centers compared to edge (Figure 6B).

The open water environments in BICY are small and shallow, which favor CH₄ ebullition due to low hydrostatic pressure and plant mediated fluxes from the predominantly anoxic porewater environments. MOX decreased the diffusive flux from open water by 26 ± 3% in the BICY wetlands (Figure 6). However, the extent of MOX observed in our study area was considerably lower than reported in a large tropical floodplain lake in the Amazon (Barbosa et al., 2018), which found that in some cases up to 100% of the methane delivered from sediments to surface waters was oxidized. Future studies are needed to resolve the drivers of the large variability in MOX across different wetlands, and thereby better predict the conditions under which MOX significantly limits wetland methane emissions. Studies in rivers and lakes have identified some factors that influence the extent of MOX. For example, a lab incubation study showed that methane formation rates in lake sediments increased with increasing temperatures, whereas MOX potential was more closely related to CH₄ concentrations rather than temperature or sediment type (Duc et al., 2010). In rivers, NH₄ ⁺ concentrations above 20 μmol L⁻¹ and light exposure have both been shown to inhibit MOX (de Angelis and Scranton, 1993; Dumestre et al., 1999), whereas MOX generally increases with higher suspended sediment concentrations (Weaver and Dugan, 1972). The high light exposure to the shallow waters at BICY, low suspended sediment concentrations (mean = 5–7 mg L⁻¹; 1974–1994), and historic mean NH₄ ⁺ concentrations of 13–20 μmol L⁻¹ (Lietz, 2000) are all factors that may result in low MOX compared to other wetland systems.

Methane emissions from woody vegetation, which has only recently received significant attention (Barba et al., 2019a), are linked to several mechanisms relevant to our measurements. One mechanism for methane emissions from woody plants appears to be similar to herbaceous plants, whereby they take up soil-derived CH₄ by the root system (Covey and Megonigal, 2019). Stem methane emissions are not related to transpiration (Barba et al., 2019b), suggesting that diffusion through the plant’s cells is the physical mechanism for methane exchange. For example, though plant adaptations such as aerenchyma and lenticels for increased O₂ diffusion in anoxic or water-logged environments can facilitate root-mediated transport of soil-produced methane (Pangala et al., 2013), transport of soil-derived CH₄ through tree stems has also been found in plants without these structures (Maier et al., 2018). Thus, regardless of whether or not cypress knees actively provide O₂ to roots in anoxic soil layers, diffusion of methane from soils to the atmosphere through knees is likely to occur, resulting in the high flux rates we observed. Methane can also be produced inside the biomass of woody vegetation by methanogenic archaea (Zeikus and Ward, 1974; Yip et al., 2019). This phenomenon is typically attributed to trees in upland forests with high CH₄ flux rates even when soils exhibit net uptake of methane (Barba et al., 2019b), though it is certainly possible that trees in wetland environments both transport methane from soils and produce it internally depending on the redox state and microbiome of woody biomass.

Woody vegetation plays an important role in CH₄ emissions in the case of environments with elevated porewater concentrations and high amounts of oxidation in shallow soils and/or surface waters. For example, average porewater pCH₄ was 20 times greater than surface water pCH₄ in our study (Figure 4). The roots of trees, which sample water and solutes from deeper soils than herbaceous plants, act as a direct conduit for methane to vent from deep soil porewaters with limited exposure to settings where oxidation occurs. On the other hand, O₂ may enter sediments from tree roots, oxidizing the rhizosphere and enabling MOX to occur deeper in the sediment than one might expect. This may explain several of the relatively high porewater δ¹³C-CH₄ values observed (i.e., maximum of −43.3‰), while the average porewater δ¹³C-CH₄ values (−57.0 ± 0.6‰) fall in the range of acetoclastic rather than hydrogenotrophic methanogenesis. Tree stems may contribute even more than knees to whole ecosystem methane fluxes because of their greater surface area and height, though few studies have shown how CH₄ fluxes change along the radial and vertical axes of a stem (Jeffrey et al., 2020). Tree-mediated methane emissions contribute 9–27% and 44–65% of the total ecosystem CH₄ flux in forested temperate and tropical wetlands, respectively (Pangala et al., 2017). Similar to methane emissions, the soil depths where methane oxidation can occur are linked to water level (Whalen and Reeburgh, 2000). Thus, the extent of wetland tree CH₄ emissions are likely linked to the interplay between rooting depth, water table depth, and redox zonation along the soil profile.

Are freshwater wetlands net sources or sinks of atmospheric GHGs and over what timescales? This critical question has been addressed in previous studies using the molar ratio between wetland CH₄ emissions and CO₂ uptake. These values ranged from 0.05–0.06 for two subtropical wetlands, 0.09–0.11 for two temperate wetlands, and 0.13–0.20 for three high-latitude wetlands (Whiting and Chanton, 2001). Our estimates were in the same range as subtropical wetlands studied by Whiting and Chanton (2001) (0.04–0.08) with the exception of TR2 (0.14) (Table 3), which had the highest knee flux rates and highest rates of ebullition from sediments (Figure 5A). Whiting and Chanton (2001) used global warming potential values of CH₄ based on the lifetime of a single pulse of CH₄ in the atmosphere (20 yr = 21.8, 100 yr = 7.6, 500 yr = 2.5) (Lelieveld et al., 1993) to reach the conclusion that all of their studied wetlands were net sinks of warming potential over a 500 yr time frame and the subtropical wetlands were also net sinks over a 100 yr time frame. We used updated global warming potentials based on a sustained flux of methane from an ecosystem, which are considerably higher (20 yr = 96, 100 yr = 45, 500 yr = 11; (Neubauer and Megonigal, 2015) to reach the conclusion that, over 100 yr timescales, the subtropical wetlands we studied are net sources of warming potential in contrast to prior studies.

There is not a clear consensus on what timescales are most important to consider with respect to global warming potential (Myhre et al., 2013; Kirschbaum, 2014). The 100 yr time frame represents the lifetime of CO₂ in the atmosphere—an order of magnitude longer than the residence time of CH₄—and the time it generally takes terrestrial ecosystems to reach successional steady state (Odum, 2014). The 500 years time frame is similar to the time period over which ocean circulation exerts an important influence on climate change (Lelieveld et al., 1998). For context, the BICY patterned landscape appears to have begun forming in the middle to late Holocene (Chamberlin et al., 2018) and the dominant vegetation shifted from herbaceous to woody vegetation ∼3,000 years before present as precipitation and hydroperiods increased (Zhang et al., 2019). Radiocarbon dating of long-chain fatty acids showed that sediment accretion rates in the BICY wetlands were an order of magnitude slower from 1,600 to 3,500 years before present (0.03 cm yr⁻¹ vs. 0.45 cm yr⁻¹ present-day) (Zhang et al., 2019). Although we have no way of measuring past methane emissions, it is possible that carbon burial and methane emission rates have maintained a similar relative balance throughout the various stages of the wetland system’s development, i.e., as carbon burial increases, so do CH₄ emissions. On the other hand, if woody vegetation is in fact a dominant pathway for deep soil methane to bypass oxidation, the wetlands may have shifted from a net sink to net source of greenhouse potential after the vegetation shift towards woody species around 3,000 years ago.

Results from this study build upon a growing body of work demonstrating the importance of wetlands for methane emissions and thus their contributions to climate change (Zhang et al., 2017). Model simulations suggest that natural wetland methane emissions may be a positive feedback for warming. This feedback may be amplifying anthropogenic radiative forcing by up to 5% by 2100, similar to increases in anthropogenic methane emissions over the same time scales (Gedney, 2004). Our results further demonstrate that emergent vegetation is an important flow path for releasing methane from wetland soils as it limits the amount of methane oxidized in surface waters and shallow soils. While the role of herbaceous plants on wetland methane emissions is represented in current global scale Earth system models, the role of woody vegetation is crudely parameterized (Riley et al., 2011). Adding mechanistic detail to predictive models is paramount for fully constraining the role of wetlands as a positive feedback for warming and, likewise, predicting the impact of natural and anthropogenic changes to hydrological regimes on wetland carbon cycling.