Subsurface Redox Interactions Regulate Ebullitive Methane Flux in Heterogeneous Mississippi River Deltaic Wetland

As interfaces connecting terrestrial and ocean ecosystems, coastal wetlands develop temporally and spatially complex redox conditions, which drive uncertainties in greenhouse gas emission as well as the total carbon budget of the coastal ecosystem. To evaluate the role of complex redox reactions in methane emission from coastal wetlands, a coupled reactive‐transport model was configured to represent subsurface biogeochemical cycles of carbon, nitrogen, and sulfur, along with production and transport of multiple gas species through diffusion and ebullition. This model study was conducted at multiple sites along a salinity gradient in the Barataria Basin at the Mississippi River Deltaic Plain. Over a freshwater to saline gradient, simulated total flux of methane was primarily controlled by its subsurface production and consumption, which were determined by redox reactions directly (e.g., methanogenesis, methanotrophy) and indirectly (e.g., competition with sulfate reduction) under aerobic and/or anaerobic conditions. At fine spatiotemporal scales, surface methane fluxes were also strongly dependent on transport processes, with episodic ebullitive fluxes leading to higher spatial and temporal variability compared to the gradient‐driven diffusion flux. Ebullitive methane fluxes were determined by methane fraction in total ebullitive gas and the frequency of ebullitive events, both of which varied with subsurface methane concentrations and other gas species. Although ebullition thresholds are constrained by local physical factors, this study indicates that redox interactions not only determine gas composition in ebullitive fluxes but can also regulate ebullition frequency through gas production.

• Subsurface methane productivity declines with increasing salinity gradient • Ebullitive methane flux is episodic and highly variable in coastal wetland system • Subsurface redox interactions regulate the frequency and gas composition in ebullitive flux

Supporting Information:
Supporting Information may be found in the online version of this article.
Methane emissions are driven by a combination of processes including production, consumption and transport, all of which can be regulated by multiple biotic and abiotic factors.Especially in coastal river deltaic wetlands, where freshwater and seawater exchange frequently, methane fluxes are highly variable across gradients of salinity and are becoming increasingly difficult to predict as they are subjected to press and pulse disturbances such as sea level rise, riverine flooding, and storm events.In addition to biogeochemical processes, methane transport within soil and across interfaces is also subject to change with local biological (e.g., vegetation, bioturbation) and physical conditions including hydrology, meteorology, and temperature which also influences methane production and consumption (e.g., Ardón et al., 2018;Bao et al., 2021;Kellner et al., 2006;Tokida et al., 2007).Among different gas transport processes, gas ebullition, unlike plant-mediated transport which depends on vegetation presence and physiology and growth factors, can occur in aquatic, mudflat and other unvegetated environments (S.Chen et al., 2021;Fechner-Levy & Hemond, 1996;Gao et al., 2013;Olsen et al., 2019).Meanwhile, compared with diffusion transport, gas ebullition is an important gas transport pathway that could account for the high spatiotemporal variability in field measurements, but it is less understood and quantified in the field due to its episodic behavior (Baird et al., 2004;X. Chen et al., 2017;Tokida et al., 2007;Villa et al., 2021).Methane ebullitive flux can be related to saturated concentration in subsurface as well as hydrostatic and/or meteorological dynamics (Maeck et al., 2014;Scandella et al., 2016).Methane ebullitive flux could have large temporal variations in coastal wetlands due to the influence in tidal flushing and wave action, leading to large discrepancies between observed and simulated fluxes.
Estimation of methane flux in coastal wetlands is usually based on empirical studies using static chambers, lab incubations, or eddy covariance flux measurements, but methods can yield highly variable measurements due to differences in the processes captured and the spatiotemporal coverage (Krauss et al., 2016;Reid et al., 2013;Shahan et al., 2022).For example, sediment core incubations can be used to measure diffusive flux (and potentially ebullition based on incubation time) but will not account for plant-mediated transport included in eddy covariance measurements.In addition, collocated measurements of redox state and other key processes may be necessary to fully understand the biogeochemical context of surface methane flux (Ardón et al., 2018;Capooci & Vargas, 2022;Chanton et al., 1989).
Considering the complicated driving factors in methane emission, process-based modeling of methane fluxes is a valuable tool for estimating fluxes across broader spatial and temporal scales than can be achieved using existing measurements, as well as for investigating underlying process interactions and potential responses of fluxes to projected changes in environmental conditions.However, most previous process-based biogeochemical modeling studies on methane flux either apply oxygen availability as the main inhibitor of methane production or use fixed production and oxidation rates for methane excluding other redox species and gas products (Kellner et al., 2006;Peltola et al., 2018; W. J. Riley et al., 2011;Y. Zhang et al., 2002).While such approaches can work well in freshwater wetlands, coastal wetlands are subject to a complex set of redox conditions due to frequent tidal water flooding and solutes exchange with land and coastal water (Doroski et al., 2019;Xiao et al., 2017).Especially, in riverine coastal wetlands, methane production and consumption are considered to be impacted by other redox species which are abundant in freshwater and seawater like nitrate, iron, or sulfate (Capooci & Vargas, 2022;Thiel et al., 2019;Ubben & Hanson, 1980).Existing methane flux models also use simplified representations of ebullition fluxes driven either by fixed concentration thresholds or by hydrostatic pressure dynamics assuming other gas species are in equilibrium with the atmosphere (e.g., Peltola et al., 2018;W. J. Riley et al., 2011).However, the complex redox conditions of coastal wetlands, along with dynamic hydrological conditions, could further modulate gas composition in ebullitive flux and regulate the episodic frequency of ebullition (Kellner et al., 2006;Maeck et al., 2014;Scandella et al., 2016).A better understanding of the correlations between complex redox reactions and ebullition fluxes in both field and model studies may help explain and constrain episodic features of ebullition fluxes which are difficult to predict.
We configured a reactive transport model that represents key biogeochemical reactions in coastal wetland sediment (Figure 1) for three purposes: to better represent complex redox reactions in coastal wetlands, to understand whether subsurface redox interactions regulate methane ebullition fluxes, and to better constrain the scale and variability of episodic methane fluxes over coastal wetland systems.The model incorporates multiple redox interactions such as sulfate reduction and acetate oxidation, that influence the production and consumption of multiple gas species including methane, nitrous oxide, and carbon dioxide in the subsurface sediment.Transport processes in the model include bioturbation, diffusion and pressure-driven ebullition for multiple gas species from subsurface.
We conducted a case study using the new model along a salinity gradient in the Barataria Basin at the Mississippi River Deltaic Plain (MRDP).As one of the fastest deteriorating sediment and freshwater receiving basins in the MRDP, the Barataria Basin is subdivided into three regions from north to south: (a) the upper basin which is mainly controlled by an engineered freshwater river diversion (Davis Pond freshwater diversion), (b) the middle basin which is impacted by both riverine water and marine water, and (c) the lower basin which is connected to the northern Gulf of Mexico and is mainly impacted by tidal and wave activities.Although few process-based biogeochemical modeling effort focuses on the methane emission from subsurface deltaic wetland, in the past few decades, intensive field effort on methane flux has been conducted including in situ chamber measurements, laboratory incubations, and real-time instrumental flux measurements (Alford et al., 1997;DeLaune et al., 1983;G. O. Holm et al., 2016;Krauss et al., 2016;Leventhal & Guntenspergen, 2004;D. Wang et al., 2021;Yu et al., 2008).Two eddy covariance towers were mounted in the west of the brackish middle bay LA1 (G.Holm et al., 2011Holm et al., -2012) ) and the freshwater upper bay LA2 (G.Holm et al., 2011Holm et al., -2013)), which were maintained by the United State Geological Survey.The combination of high temporal frequency flux data and spatially widespread historical measurements in this region provides a valuable data set for evaluating model simulations of key redox reactions and surface fluxes.In this case study, five sites were selected to represent a salinity gradient that includes a freshwater marsh with an eddy covariance tower (LA2, active), an intermediate marsh, a brackish marsh with an eddy covariance tower (LA1, retired) at bay west, a brackish marsh at bay east, and a salt marsh (Figure S1 in Supporting Information S1).The objectives of this modeling study are (a) to connect surface methane flux with subsurface methane productivity through including new complex redox reactions instead of using empirical relationships (Table 1), (b) to understand how subsurface redox interactions in coastal wetland contribute to methane ebullition flux, and (c) to constrain the variability of episodic methane fluxes from dynamic coastal wetland systems.

Complex Redox Reaction Network
To better resolve the redox complexity and multiple gas species dynamics in coastal wetland, a newly coupled reactive transport model framework was developed to simulate subsurface redox dynamics in coastal wetland.A new complex reactive network of key redox reactions were defined in the chemical reaction solver-PFLOTRAN  1 and 2. The salinity gradient spans from salt marsh, brackish marsh to intermediate marsh and freshwater marsh (Figure S1 in Supporting Information S1). ) reduction, iron (Fe) reduction, nitrification, denitrification, methanogenesis, and methanotrophy along with dissimilatory nitrate reduction to ammonium (DNRA) and abiotic mineral formation such as mono-sulfide-iron precipitation (Tables 1 and 2).Soil organic matter decomposition in PFLOTRAN follows previous implementations in land models with modifications over decomposition rates and C:N ratios (B.N. Sulman et al., 2022;Tang et al., 2016).The DNRA includes a fermentative pathway and an autotrophic pathway with hydrogen sulfide (H 2 S) or methane (CH 4 ) as the electron donors.Both aceto-clastic and hydrogenotrophic pathways of methanogenesis are included as separate reaction pathways.Aerobic methanotrophy as well as anaerobic CH 4 oxidation by other substrates (e.g., NO 3 − , SO 4 2-, Fe 3+ ) are also included as CH 4 consumption pathways, which together reflect the major aerobic and anaerobic methane redox reactions in coastal environment.Sediment organic matter (SOM) and minerals including mono-sulfide-iron, ferric hydroxide (Fe(OH) 3 ) and iron  2004) Note. "/" indicates information is not available.Temperature impacts on microbial-driven reactions (Equation S1 in Supporting Information S1) are simulated with the Q 10 (2.0) and the optimum temperature is set to 20 o C (Ma et al., 2017).Due to the lack of site-specific and high temporal resolution measurements of sediment temperature for all sites, sediment temperature was assumed to be the same as the water temperature.Although microbial reactions respond to temperature differently, detailed microbial activities measurements are limited in the study region.Thus, the same Q 10 parameter is applied to all microbial reactions in the model, although CH 4 production has a wide range of Q 10 and may be more sensitive to temperature than other reactions (Inglett et al., 2012).

Transport Processes
Vertical transport of chemical species within the sediment column was simulated using a customized transport code coupled to the reaction solving capabilities of PFLOTRAN.Transport of gas species between sediment and atmosphere and/or overlaying water occurs in the model through multiple pathways (Figure 2).Because this model study is primarily configured to represent the sediment biogeochemical dynamics in coastal wetlands, only transport processes driven by diffusion, advection and ebullition are coupled with the reaction processes in PFLOTRAN.Vegetation plant functions are beyond the investigation scope, so plant-mediated gas transport, which varies among plant species and requires detailed information on plant functional types, root density and leaf area, etc., is not included in this model study.Since plant-mediated transport can be a critical pathway for O 2 , CH 4 , and other gas species to exchange between sediment and atmosphere (Harden & Chaton, 1994;Y. Li et al., 2021;Murray et al., 2019;Villa et al., 2020), the omission of this pathway means that the simulations likely underestimate the surface CH 4 gas flux but may also underestimate CH 4 consumption driven by transport of oxygen into the subsurface.
Diffusion in the model happens across interfaces and within the sediment column along concentration gradients of the dissolved ions between water/air and sediment porewater, which is governed by the Fick's second Law (Equation S2 in Supporting Information S1).The actual diffusion coefficient in the model is determined by the apparent ion diffusion coefficient of each chemical species and the sediment porosity at each depth (Table .2, Equation S3 in Supporting Information S1).Additionally, bioturbation is considered as a diffusion enhancing factor for all dissolved species and changes with water inundation depth (Equation S4 in Supporting Information S1, (Soetaert et al., 1996)).Moreover, a maximum depth exists for bioturbation activity in the sediment below which it is negligible.An estimated vertical sediment accumulation rate is included and is subjected to change with inundation depth (Equation S5 in Supporting Information S1).Diffusion and advection transport processes within sediment layers are simulated with a finite volume solution scheme.Wet and dry boundaries are considered at the surface sediment layer of the simulated sediment profile.To represent the aerobic condition during sediment exposure, gas diffusion coefficients are two orders of magnitude higher under the dry surface boundary than those under the wet surface boundary.
Ebullition is considered as the other transport pathway for the dissolved gas species in sediment porewater.Multiple factors can influence ebullition including hydrostatic pressure, temperature variation, gas production and consumption rate, and sediment properties (Kellner et al., 2006;Liu et al., 2016;Maeck et al., 2014;Scandella et al., 2016;Tokida et al., 2007).In this model, multiple gas species can be produced or consumed in the simulated sediment column.Ebullition happens when the sum of the partial pressure of each dissolved gas species exceeds the hydrostatic pressure (Equations S6-S13 in Supporting Information S1, (Sander, 2015)).Thus, ebullition flux of each gas species is determined by the fraction of its partial pressure relative to the sum of all gas species.It is assumed that the ascending gas reaches to the air interface immediately after ebullition occurs (Peltola et al., 2018; W. J. Riley et al., 2011).

Sediment Column Set-Up
Data from sediment cores collected by the Coastwide Reference Monitoring System (CRMS, https://cims.coastal.la.gov/monitoring-data/) were used to configure sediment columns in the model.Each simulated column has a total depth of 31 cm and was divided into six layers (Figure 2), although a longer core length can also be easily set up in this model frame.A core length of 31 cm was configured based on field core length of CRMS soil cores.Sediment properties used in the model include porosity and bulk density (Figure S2 in Supporting Information S1).
The simulated sites were treated as saturated all the time due to the lack of high temporal resolution field saturation data.However, a saturated status can be easily reached at the study area due to the local humid climate, frequent precipitation and flooding conditions.The surface deposition rate of organic matter to the sediment column used by this model is assumed to be constant over time with a fixed C:N ratio (Table .3),although the deposition rate and C:N ratio could vary with time or space in reality (S.Li et al., 2021;Upreti et al., 2021).The deposited sediment organic matter is further buried downward within the sediment column.Considering the shallow depth of 31 cm of the sediment column, sediment compaction and groundwater influences are neglected in this model setting.

Physical Forcing
Hourly hydrological monitoring data from the CRMS program was used to drive the reactive transport model.Flooding (wet) or exposure (dry) boundary conditions were determined by water level relative to the marsh elevation.Water temperature and water level were used to inform the sediment temperature and hydrostatic pressure at each site in order to derive the Henry's Law constant for each dissolved gas species and saturated concentration in sediment porewater at each depth.The ebullition process in this study assumes that sediment temperature, which can affect the Henry's Law constant, is the same as that of the overlying flood water, and also assumes that gas Note."/" indicates information is not available.

Table 3
Sediment Column Set-Up Parameters Used by the Model concentration in the flooding water is in equilibrium with the atmosphere, which means that the escaping gasses are directly added to the air interface regardless of flooding or exposure conditions.
To calculate concentration boundary conditions for exchange of dissolved ions through the water-sediment interface, some dissolved ion concentrations in overlying water were converted to molarity concentrations from surface water salinity according to fraction relative to chloride (Table .4,Equations 1 and 2 in Supporting Information S1, (Culkin, 1965)).Concentrations of other dissolved species which are not directly associated with salinity are provided with values from local field measurements and are assumed to be constant in the surface water.When the simulated sites were exposed to air, diffusion coefficients of dissolved gasses through the air-sediments interface were increased by two orders of magnitude to represent an aerobic condition.
[Cl − ] = Sal 0.00180665 (1) Cl − is the chloride mass concentration in the unit of mg⋅L −1 , Sal is the water salinity in unit of ppt, [Ion] is the ion mass concentration in the unit of mg⋅L −1 , and f ion is the mass concentration of each ion per unit chloride mass in water.

Model Evaluation
After spun up for 22 years, an 11-year model simulation (2009-2019) was analyzed for five selected sites along a salinity gradient from freshwater to salt marsh.Because the model was primarily configured to simulate subsurface redox reactions, model validation was based on available field measurements of porewater salinity from CRMS and porewater dissolved CH 4 and CO 2 at selected sites (He et al., 2022).Constrained by data availability, porewater salinity and porewater dissolved CH 4 and CO 2 measured at different dates and time (see Figures 3 and 4 for detailed dates) were compared with the model results.Yearly mean porewater salinity was calculated for each site, and compared with field measurements of year 2006 (Figure 3).Monthly measurements of dissolved CH 4 and CO 2 concentrations at adjacent sites were also compared with model simulations at one freshwater (salinity: 0-0.5 ppt), one intermediate (salinity: 0.5-5 ppt), two brackish (salinity: 5-18 ppt) and one salt marsh (salinity: 18-30 ppt) (intermediate marsh: Figures 4 and 5; brackish marsh at bay east: Figures S4 and S5 in Supporting Information S1; salt marsh: Figures S6 and S7 in Supporting Information S1).Marsh type and salinity range are derived from CRMS vegetation and salinity zonation (Visser et al., 2002).
Field eddy covariance fluxes implicitly include fluxes influenced by a more complete suite of processes and sources including fluxes from wetland sediment, tidal water, plant-mediated transport, and wind or atmospheric pressure/ temperature effects unlike the simulated fluxes in this model which only considers gas produced in subsurface redox reactions and includes only diffusion and ebullition as their major transport pathways from subsurface to surface interfaces.Thus, we do not expect model simulations to reproduce measured eddy covariance flux patterns at high temporal resolution.Because the model uses a highly simplified estimate of plant carbon inputs to the system, we also do not expect the model to reproduce the observed magnitude of total carbon fluxes.Moreover, field recordings are usually interrupted by field conditions which adds more challenges to high temporal model-field data comparison.Therefore, we focus our evaluation of simulated surface fluxes on monthly mean fluxes and normalized seasonal patterns.To account for gaps in measured flux data, only days with at least 12 hr of valid measurements and only months with at least 14 days of valid measurements were included in the model evaluation.Because of the high number of gaps in the measured fluxes from the brackish marsh site (LA1, Figure S3 in Supporting Note."/" indicates information is not available.

Table 4 Molarity Concentration of Dissolved Species in Flooding Water
Information S1), the monthly and seasonal comparisons focused on the freshwater site (LA2) fluxes, and the brackish site fluxes were used for qualitative comparison of mean fluxes across contrasting salinity environments.

Porewater Chemistry
Overall, our model simulated porewater salinity profiles are comparable to the measured profiles, especially at depths below the top 5 cm (Figure 3).However, the model overestimated salinity at one of the brackish sites.Simulated porewater salinity varied across years, with interannual variability declining with depth and increasing with average site salinity (spatial salinity gradient).However, these interannual variations are out of the scope of field measurements because only 1 year of measurements was available for comparison with the model.
Simulated porewater CH 4 and CO 2 concentrations and their vertical gradients also compared well to the field measurements (Figures 4 and 5).Moreover, the simulated porewater CH 4 concentrations also captured the variation among seasons shown in the field instantaneous measurements (concentration increased from winter to fall, 2019, Figure 4).Porewater CO 2 concentrations did not have a clear seasonal pattern in either simulations or measurements, although there was variability among different sampling dates (Figure 5).
Overall, the vertical and temporal differences in dissolved species at different sites clearly reflect the major redox reactions/relations configured in the redox network which is informed by previous studies (Ardón et al., 2018;Doroski et al., 2019;King & Wiebe, 1980).The simulated porewater chemistry presents high heterogeneity with depth, time, and across the selected salinity gradient (Figure 6).Different sites share similarities in the vertical distributions of some dissolved chemical species.For example, the simulated oxygen and sulfate concentrations in porewater decrease from top to bottom at each site, while CH 4 concentration is much lower at the surface and peaks at intermediate depths.These vertical distributions follow the redox reaction network configured in the model, where oxygen and sulfate reduce CH 4 concentration both by suppressing CH 4 production (through competition for acetate) and by serving as terminal electron acceptors for CH 4 consumption.Porewater sulfate concentration greatly increases with site salinity which results in much lower CH 4 concentrations in the salt and brackish marshes than in the freshwater and intermediate marshes.Thus, the performance of the configured redox network and the concentration differences in porewater CH 4 among sites are consistent with the common understanding that saline wetland environments produce less CH 4 .

Surface Methane Fluxes
Overall, the seasonal pattern of measured CH 4 flux is captured by the model, with fluxes highest in summer and lowest in winter (Figure 7; Figure S3 in Supporting Information S1).The simulated CH 4 flux is higher at the freshwater marsh than the brackish marsh (Figures 8a and 8b), matching the qualitative contrast seen in measured fluxes (G.O. Holm et al., 2016).However, the model substantially underestimated the magnitude of CH 4 flux (by up to an order of magnitude) when it was compared to field measurements, which is within expectation.The ratio of modeled to measured fluxes is lower during the warm season and higher during the cold season.

Subsurface Methane Productivity
In general, surface CH 4 flux increases with increasing water level and temperature at all sites (Figures 8a and 8b).Slightly positive relationships exist between monthly mean water level and surface CH 4 flux at all sites, although different sites have different hydrological conditions due to their elevation and distance to tidal channels and/or river channels (Figure 8a, Figure S8 in Supporting Information S1).Strong positive relationships exist between surface CH 4 flux and temperature for all sites (Figure 8b).In general, high porewater CH 4 concentration leads to high surface flux along the salinity gradient, but large variations exist in surface flux under high porewater CH 4 concentrations (Figures 8a and 8b).The positive correlation between CH 4 porewater concentration and surface flux highlights that surface CH 4 flux is largely driven by subsurface CH 4 productivity.However, patterns of porewater concentration were also visibly different than patterns of surface CH 4 flux, indicating some decoupling between concentration and flux driven by vertical transport and CH 4 oxidation processes.
Similar to the surface CH 4 flux, porewater CH 4 concentration positively increases with rising water level at all sites (Figure 8c).Since the impacts of temperature on the redox reactions are based on Q 10 (2.0), strong positive relationships exist between temperature and porewater CH 4 (Figure 8d).However, under similar inundation and temperature conditions, porewater CH 4 concentrations were highly variable along the salinity gradient (Figures 8c  and 8d).Although inundation patterns are similar among sites, the fresher sites have much higher porewater CH 4 than the saltier sites.Moreover, the increasing temperature has less effect on porewater CH 4 in the saltier marshes than that in the freshwater marshes (Figure 8d), indicating that there are other constraints on CH 4 production.
Porewater O 2 concentration at the sediment surface is sensitive to inundation-exposure conditions driven by water level fluctuations.Because the diffusion rate of oxygen through saturated sediment layers and surface water is configured about 4 orders of magnitude lower than diffusion through air and unsaturated sediments in the model (Nuttle & Hemond, 1988), sediment porewater O 2 declines rapidly during flooding and gets replenished when the sediment surface is exposed.Since porewater O 2 concentration declines quickly with depth in the sediment column (Figures 6a-6d), only O 2 in the surface sediment is discussed here.In general, a negative relationship is expected between porewater oxygen and CH 4 , because aerobic conditions tend to inhibit CH 4 production and anaerobic conditions usually favor CH 4 production.This negative relationship is clearly reflected in the vertical profiles of O 2 and CH 4 concentration in porewater (Figure 6), but it becomes less obvious if depth-averaged porewater O 2 and CH 4 from all simulated sites are considered (Figure 8e).Instead, CH 4 concentrations are very scattered from low to high oxygen conditions among different sites.One reason is that depth-averaged CH 4 concentration mostly reflects the high CH 4 concentration in deep sediment (anaerobic condition) whose variations are beyond the control of surface oxygen concentration.Thus, for a specific environment (i.e., freshwater wetland), the temporal variations of oxygen in surface sediment may not properly predict porewater CH 4 dynamic in subsurface.Second, over a large salinity gradient, the relationship between porewater CH 4 and oxygen can also be greatly altered by other chemical species such as sulfate from flooded water in coastal system.This is supported by the strong negative relationship between porewater sulfate and CH 4 (Figure 8f).Increasing porewater sulfate leads to low CH 4 concentration in porewater (Figure 8f).This further indicates that surface oxygen dynamics could have limited control over the variations in deep CH 4 production and consumption, and surface oxygen dynamics may not be an effective temporal and spatial predictor for CH 4 dynamic in deep sediment.Instead, as salinity increases in coastal wetland system, sulfate reduction may be a dominant factor that controls the temporal and spatial variations of CH 4 productivity in deep sediment within a coastal wetland system (Figures 8e and 8f).

Episodic Ebullitive Fluxes and Methane Fraction in Total Ebullitive Fluxes
In addition to variable gas production in the subsurface, gas escaping through ebullition drives large fluctuations in surface fluxes but is poorly constrained in both field and model studies.In our model simulations, ebullition accounted for less than 20% of total CH 4 flux at monthly time scale, but ebullition flux was usually distributed within high flux ranges and thus drove episodic events of very high CH 4 flux rates (Figures 9a and 9b, Table S1 in Supporting Information S1).Most importantly, the fractions of ebullitive contribution to surface CH 4 fluxes are more positively related to the increasing porewater CH 4 in freshwater marsh, but are more variable and scattered in the brackish and salt marshes (Figure 9b).Averaged over monthly scale, ebullition event frequency is more variable and less related to porewater CH 4 concentration (<75 μM) at brackish and salt marshes than at freshwater marsh (Figure 9c).At freshwater marsh, ebullition event frequency increases with porewater CH 4 concentration (>75 μM).Unlike ebullition frequency, CH 4 fraction in the total ebullitive fluxes increases more steadily with increasing porewater CH 4 until it reaches to an upper bound 0.18 (Figure 9d, Figures S9 and S10 in Supporting Information S1).

Subsurface Methane Production in Heterogeneous Wetland Systems
A major challenge in estimating and projecting ecosystem CH 4 gas budget is their high spatiotemporal heterogeneity within and across ecosystems.Surface fluxes depend on the physical and biological processes that alter the production, consumption, or transport pathways for gases in the subsurface.Here, we used a novel model framework to explore the impacts of physical and biogeochemical conditions on the spatiotemporal heterogeneity of surface fluxes and subsurface concentrations of CH 4 in deltaic wetland.This study focused on fluctuations in water level and temperature as drivers of hourly to seasonal changes.In the study area, hydrological conditions are controlled by river discharge, wave activities, and tidal currents, which determine the seasonal inundation-exposure patterns as well as the positive seasonal relationship between water level and temperature at each site.Since the study area is influenced by both freshwater and saltwater, wetlands are distributed along a salinity gradient.Thus, our study investigated the interaction of physical and biogeochemical factors by simulating variations in CH 4 production patterns across selected sites representing a range of salinities within the Mississippi Delta region.
Our model simulations substantially underestimated total CH 4 emissions relative to field measurements.The model bias and its seasonal variation in simulated surface flux were likely due to multiple reasons.Eddy covariance measurements provide high temporal fluxes at large spatial scale, which are the results of a complete suite of transport processes across the ecosystem (Krauss et al., 2016).However, our model is configured to represent small spatial scale processes.Moreover, plant-mediated gas transport also contributes to field scale fluxes but is not included in this model.The plant-mediated gas transport mechanism is neither well understood nor well quantified in field studies and most numerical models, although a couple of approaches (e.g., concentration gradient driven vs. partial pressure driven) have been applied in some wetland models but still have not significantly improved model performance (W.J. Riley et al., 2011;Y. Zhang et al., 2002).Although this complex redox network includes the major subsurface redox reactions, there are other important reactions which may influence the ratio of a specific gas in the bubbles but are not included in this network (i.e., H 2 ).Additionally, this model configuration does not include pH impacts on the redox reactions, but pH dynamics can be critical processes or controlling factors for subsurface soil and coastal water environments and it is suggested to evaluate pH dynamics in future modeling effort (B.N. Sulman et al., 2022).Finally, our model framework used a highly simplified representation of plant carbon inputs to the subsurface, which were simulated as a constant input rate of sediment organic matter in the top sediment layer.Labile organic matter inputs from fresh plant litter and root exudates can be readily converted to CH 4 (Girkin et al., 2018;Hoyos-Santillan et al., 2016), so underestimates of plant-derived labile C inputs in the model could drive underestimates of total CH 4 production.This study therefore represents a first step toward incorporating more complex redox chemistry into carbon cycle models, and improved site-specific mechanisms and carbon input estimates will likely be necessary for carbon cycle estimates and quantitative comparisons with field measurements.Future mechanistic understanding on roots activities, microbial temperature and salinity sensitivity may greatly improve similar model framework.Besides, high temporal (daily or monthly or seasonal) porewater chemistry, site-specific sediment surface flux (i.e., autonomous chamber flux), sediment temperature and moisture measurements can help resolve model disparities in CH 4 budget estimates and improve our knowledge on the interactions between subsurface and surface processes in carbon cycle.Since this model primarily focuses on representing redox interactions, a detailed sensitivity analysis of parameters is not included in the current analysis.However, improved parameter constraints and sensitivity information can be helpful when applying this model framework to quantitative predictions of carbon cycling.
In coastal wetlands, both water level fluctuations and temperature could influence subsurface CH 4 production and consumption either through alternating anoxic/oxic conditions or regulating microbial activities (Bao et al., 2021;Boon et al., 1997;Hondula et al., 2021;Roslev & King, 1996;Segers, 1998).However, our simulations suggest that surface oxygen dynamics have limited control over CH 4 variations in deep sediment, and surface oxygen dynamics may not be an effective temporal and spatial predictor for CH 4 dynamic in deep sediment.In saline systems, flooding also introduces sulfate to the system.High sulfate concentration may regulate CH 4 production by allowing sulfate reducers to compete with methanogens for organic substrates, and may also increase subsurface CH 4 consumption by serving as a terminal electron acceptor for CH 4 oxidation (Koebsch et al., 2019;La et al., 2022).Additionally, frequent flooding increases the connectivity between interior wetland and coastal water, and thus lateral export of dissolved gases (i.e., CH 4 ) can be high in some coastal system (Bogard et al., 2020;Santos et al., 2021;Tamborski et al., 2021).Thus, to better constrain the CH 4 budgets in coastal wetland systems, it is critical for field and model effort to focus on the redox interactions in deep sediment including different anaerobic reactions as well as possible oxygen inputs by roots or (submarine) groundwater.This may also apply to other greenhouse gases budgets in similar complex redox systems (Figures S9 and S10 in Supporting Information S1).

Subsurface Redox Interactions Contribute to Methane Ebullition
Gas ebullition is usually considered as driven primarily by pressure and temperature (Baird et al., 2004;Tokida et al., 2007;Villa et al., 2021).However, in coastal wetlands, subsurface redox interactions may become a first-order driving factor in the ebullition process by producing multiple gas species and creating steep subsurface gas gradients.The redox interaction network included in this model allows it to simulate the consumption and/ or production of multiple gas species across a range of redox conditions, including O 2 , N 2 , H 2 , CO 2 , CH 4 , N 2 O, and H 2 S (Tables 1 and 2).Thus, the model framework applied here allows a detailed investigation of controls on ebullition flux of CH 4 and other greenhouse gases (Figures S9 and S10 in Supporting Information S1), and allows further evaluation of how strongly subsurface redox interactions regulate ebullition frequency and ebullitive gas composition.In this study, although multiple gas species are included, we only focus our analysis on CH 4 .
Our simulations suggest that ebullitive CH 4 fluxes are more variable and less related to porewater CH 4 concentration at brackish and salt marshes than at the freshwater marsh.This contrast across the salinity gradient highlights the spatiotemporal heterogeneity of ebullition contribution to surface flux for specific gas species and highlights that the contribution of ebullition to surface gas flux is not driven only by dissolved gas concentration gradients.Ebullition flux for specific gas species is determined by event frequency and gas composition.In coastal wetland systems, physical factors such as hydrostatic pressure usually constrain the physical thresholds and the event frequency, while complex redox reactions not only determine gas composition but also affect event frequency through gas production.Our simulations suggest for specific gas species such as CH 4 , its gas composition ratio may have an upper bound regardless of the porewater concentration gradient.Hence, a parameterized porewater concentration threshold for ebullition flux in models may result in disparities between model and field measurements.Moreover, subsurface redox interactions could lead to high ebullition frequency under a high gas production environment.Thus, in a coastal wetland system that spans over a large salinity gradient, ebullition events may be observed at similar frequencies across the salinity gradient, but higher ratios of CH 4 in ebullitive gas may be expected in freshwater marshes than saltier marshes and could thereby drive higher CH 4 emissions in freshwater marshes.The coastal wetland system plays key roles in carbon budgets within the land-to-ocean aquatic continuum (Regnier et al., 2022).However, whether the coastal wetland system is a net sink of buried carbon or a net source of escaping carbon and how its role in the carbon budget could shift with changing climate are still unknown.Coastal wetland systems are impacted by more complex and frequent physical and biogeochemical conditions compared with land or ocean systems, which adds challenges in resolving their physical and biogeochemical processes in models.Considering the complex physical and biogeochemical conditions characteristic of coastal wetlands, this model study investigated the connection between subsurface redox interactions and CH 4 emissions from these coastal systems.Our model results suggest that in dynamic coastal wetland systems, the first-order factors that drive greenhouse gas release from the subsurface especially through ebullition frequently shift from physical forcings to redox conditions, which could increase the difficulty of measuring and interpreting greenhouse gas emissions in field studies and raise challenges for accurate projections of coastal wetland carbon cycling in Earth system models.Ebullition flux is highly episodic and variable over time, which highlights the importance of including subsurface redox interactions and their effects on ebullition in future carbon budget studies for coastal wetland systems.Field measurements of ebullitive CH 4 is very important for future model improvement.This study further suggests that caution should be taken when temporally limited field data is used to estimate greenhouse gas budgets in coastal wetland systems because episodic ebullition fluxes could introduce large variations in flux measurements.Thus, it is critical to consider in-situ physical conditions (e.g., tide or wave activities, weather condition, distance to channels) as well as subsurface redox environments (e.g., salinity gradient or wet-dry condition) in data preparation and interpretation.

Conclusions
Through coupling a reactive biogeochemical model-PFLOTRAN with a customized transport module via Alquimia interface, this model study investigates wetland CH 4 production, consumption, and emission along a salinity gradient in the submerging Mississippi River delta.A complex redox network is built to include biogeochemical cycles of carbon, nitrogen, and sulfur, and production of multiple gas species including CH 4 , N 2 , N 2 O, and CO 2 .For CH 4 specifically, the total surface flux responds to subsurface CH 4 productivity, transport, and oxidation.Subsurface CH 4 productivity is controlled by the redox interactions directly (e.g., methanogenesis, methanotrophy) and/or indirectly (e.g., competition for substrate with sulfate reduction) under aerobic and/or anaerobic conditions, which modifies the relationships between CH 4 flux and oxygen as well as temperature.Although the ebullition threshold is constrained by physical factors, this study also suggests that complex redox reactions determine gas composition in ebullition flux and can also affect ebullition frequency through redox-dependent rates of gas production.Thus, complex redox conditions are important considerations in mechanistic studies on greenhouse gas budget in the field/laboratory.Additionally, proper representation of ebullition process in Earth system models is important in estimating greenhouse gas budget, especially for ecosystems that are subject to complex redox interactions, such as coastal wetland systems.

Figure 1 .
Figure 1.Major redox zones and salinity gradient considered in this study.Specific redox pathways, reaction rates and parameters for each chemical species used by the model are listed in Tables1 and 2. The salinity gradient spans from salt marsh, brackish marsh to intermediate marsh and freshwater marsh (FigureS1in Supporting Information S1).
hydroxide (Fe(OH) 2 ) are considered as solid species in the model.Mineral precipitation processes are critical for an iron-rich and sulfate rich coastal wetland environment which are part of the iron and sulfate cycles and could directly or indirectly impact methane dynamics via microbial activities (i.e., competition between sulfate reducers and methanogen) or iron oxidation (availability of iron by including mineral precipitation).Gas species are represented as dissolved gas phase in sediment porewater, including O 2 , N 2 , H 2 , CO 2 , CH 4 , H 2 S, and N 2 O.Other chemical species are treated as dissolved ions in both flooding water and sediment porewater in the model.

Figure 2 .
Figure 2. Diagram of transport processes across the sediment water interface (SWI) and within single sediment column.

Figure 7 .
Figure 7.Comparison between model simulation and Ameri Flux observation during 2011-2013 at the freshwater marsh CRMS3166 (eddy covariance tower: US-LA2).(a) Ratio between model simulation and field data of monthly mean methane flux.(b) Normalized monthly mean flux from model simulation (orange) and field data (light blue).

Figure 8 .
Figure 8. Relationships between modeled monthly mean surface CH 4 flux, depth-averaged porewater CH 4 concentration and relevant driving factors along the salinity gradient (freshwater: ; intermediate: ; brackish (bay west): ; brackish (bay east): ; salt: ).(a) Surface CH 4 flux versus water level; (b) CH 4 flux versus temperature; (c) porewater CH 4 versus water level; (d) porewater CH 4 versus temperature; (e) depth-averaged porewater CH 4 versus porewater oxygen concentration; (f) porewater CH 4 versus porewater sulfate.Note: zero water level (gray line) means marsh surface, and color bars from top to bottom represent depth-averaged porewater CH 4 , surface CH 4 flux and water level, respectively, and all porewater solute concentrations are depth-averaged if it is not otherwise noted.

Figure 9 .
Figure 9. Spatiotemporal heterogeneity in modeled gas ebullition.(a) Hourly mean CH 4 flux versus hourly mean CH 4 ebullition flux and gray line represent the 1:1 line (note: hourly mean data is adapted to show the large fluxes contributed by CH 4 ebullition flux; different symbol/color represents different sites along the salinity gradient), and (a') hourly mean total CH 4 flux versus hourly mean ebullition flux at the freshwater site from 05/15/2012 to 05/31/2012; (b) monthly mean porewater CH 4 concentration versus fractional contribution of ebullition to total CH 4 flux; (c) monthly mean porewater CH 4 concentration versus ebullition events per hour; (d) monthly mean porewater CH 4 concentration versus CH 4 ratio in total gas released in the ebullition events.Note: all data are monthly mean data if not otherwise noted, and each symbol represents one salinity range.

Table 1
Complex Redox Reactions, Chemical Species and Reaction Rates Used by the PFLOTRAN Model (Note: Reaction Rates (mol • L −1 • s −1 ) Used in the Model Were Converted From Rates in References Which Are Usually in the Unit of mol • m −2 • s −1 )

Table 2
Parameters of Chemical Species Used by the PFLOTRAN Model