Integrating Tide‐Driven Wetland Soil Redox and Biogeochemical Interactions Into a Land Surface Model

Redox processes, aqueous and solid‐phase chemistry, and pH dynamics are key drivers of subsurface biogeochemical cycling and methanogenesis in terrestrial and wetland ecosystems but are typically not included in terrestrial carbon cycle models. These omissions may introduce errors when simulating systems where redox interactions and pH fluctuations are important, such as wetlands where saturation of soils can produce anoxic conditions and coastal systems where sulfate inputs from seawater can influence biogeochemistry. Integrating cycling of redox‐sensitive elements could therefore allow models to better represent key elements of carbon cycling and greenhouse gas production. We describe a model framework that couples the Energy Exascale Earth System Model (E3SM) Land Model (ELM) with PFLOTRAN biogeochemistry, allowing geochemical processes and redox interactions to be integrated with land surface model simulations. We implemented a reaction network including aerobic decomposition, fermentation, sulfate reduction, sulfide oxidation, methanogenesis, and methanotrophy as well as pH dynamics along with iron oxide and iron sulfide mineral precipitation and dissolution. We simulated biogeochemical cycling in tidal wetlands subject to either saltwater or freshwater inputs driven by tidal hydrological dynamics. In simulations with saltwater tidal inputs, sulfate reduction led to accumulation of sulfide, higher dissolved inorganic carbon concentrations, lower dissolved organic carbon concentrations, and lower methane emissions than simulations with freshwater tidal inputs. Model simulations compared well with measured porewater concentrations and surface gas emissions from coastal wetlands in the Northeastern United States. These results demonstrate how simulating geochemical reaction networks can improve land surface model simulations of subsurface biogeochemistry and carbon cycling.


Introduction
Coastal wetlands can sequester carbon at exceptionally high rates (>200 g C m 2 year 1 ), with global total C sequestration rates (up to >100 Tg C year 1 ) rivaling those of temperate, tropical, or boreal forests (Chmura et al., 2003;McLeod et al., 2011).Vegetated areas, including salt marsh, mangrove, and sea grass, have been estimated to account for half of global ocean organic carbon burial (Duarte et al., 2005).Salt marshes are also vulnerable to loss, with an estimated equivalent global emissions of 16 Tg CO 2 equivalent year 1 due to loss of salt marsh area between 2000 and 2019 (Campbell et al., 2022).Wetland emissions of greenhouse gases such as methane are highly sensitive to salinity, and particularly to sulfur cycling driven by seawater influence in coastal systems (Poffenbarger et al., 2011).While wetlands account for 40% of global natural methane emissions (Saunois et al., 2020), coastal wetlands are thought to be a relatively minor contributor (<2% of wetland emissions).However, significant methane production has been observed in some salt marshes (Capooci et al., 2024) and mangroves (Rosentreter et al., 2018).Changes in sulfur dynamics and seawater influence can contribute to peat collapse and rapid carbon loss in coastal wetland systems subject to changing sea levels (Chambers et al., 2019).Tidal wetlands represent a key challenge for existing carbon cycle modeling frameworks, due to their outsized role in the carbon cycle and the complex combination of hydrology, redox dynamics, and interactions of different chemical cycles that drive subsurface biogeochemistry in these systems (Ward et al., 2020).
Chemical interactions including pH dynamics, redox cycling, oxygen consumption, and mineral interactions are recognized as key drivers of soil carbon cycling in both oxic (Hall et al., 2018;Li et al., 2021;Sollins et al., 1996) and anoxic (Kögel-Knabner et al., 2010;Lipson et al., 2010;Sutton-Grier et al., 2011) environments.Redox interactions are particularly important in determining greenhouse gas emissions in inundated soils subject to redox fluctuations (Ginn et al., 2017;Sulman et al., 2022).However, land surface models (LSMs) that are used to simulate and project carbon and nutrient cycling as part of Earth system model (ESM) simulations typically use simplified representations of organic matter cycling that include only carbon, macronutrients (N and P), water, and energy cycling (Todd-Brown et al., 2013).When representing decomposition under saturated conditions, LSMs typically do not explicitly represent oxygen consumption or redox interactions, but instead use empirical functions of moisture to simulate organic matter decomposition rates as well as CO 2 and CH 4 production (Wania et al., 2013).These omissions could drive uncertainties and predictive errors when simulating biogeochemical responses to changing hydrological conditions or projecting carbon cycling across different soil types.
In saturated soils, the omission of redox cycling and oxygen concentrations could lead to bias in simulations of organic matter degradation as well as greenhouse gas production.Existing LSMs typically treat soil saturation as a proxy for redox state, assuming that saturated conditions translate directly to oxygen depletion (Wania et al., 2013).In reality, both organic matter decomposition and methane (CH 4 ) production are sensitive to the presence and depletion of terminal electron acceptors (TEAs) including oxygen, iron, sulfate, nitrate, and manganese (Estop-Aragonés et al., 2013;Herndon et al., 2015;Poffenbarger et al., 2011).While some models do include a temporal delay in methane (CH 4 ) production as a proxy for the depletion of TEAs (Riley et al., 2011), such proxy approaches may not be sufficient to represent variations in TEA patterns across variations in soil mineral content or in situations where flows of dissolved oxygen or plant-mediated oxygen transport are important.In frequently flooded coastal or riparian systems, such approaches may not adequately represent the addition and mixing of TEAs, and may fail to accurately predict methane fluxes in coastal systems where increasing sulfate availability suppresses methane production even as water levels rise (Kirwan et al., 2023).Redox conditions and porewater chemical concentrations can also affect plant growth.For example, plant tolerance to salinity and inundation varies widely (LaFond-Hudson & Sulman, 2023) and sulfides produced via sulfate reduction under anoxic conditions can be toxic to plants (Koch et al., 1990;Lamers et al., 2013).Thus, representing dynamics of redox-active chemical species may allow LSMs to simulate wetland carbon cycling processes and greenhouse gas emissions more accurately.
Incorporating representation of chemical interactions directly into LSM codes has been challenging due to the complexity of introducing processes specific to individual chemicals into already-complex model structures.Specialized reactive transport simulators do exist that can simulate complex biogeochemical reaction networks (Frei et al., 2012;Hammond et al., 2014;Perzan et al., 2021;Steefel et al., 2015;J. Tang et al., 2022), and simulators such as PFLOTRAN (Hammond et al., 2014) include flexible configuration systems allowing alternative reaction network structures to be represented without work-intensive changes to model code (Hammond, 2022).Previous work to couple reactive transport simulators to existing LSMs has demonstrated the feasibility of offloading biogeochemical calculations from fixed representations in LSM code to more flexible reaction network simulators, but these implementations have not previously moved beyond demonstrating that existing LSM soil C and macronutrient representations can be reproduced in the coupled codes (G.Tang et al., 2016;J. Tang et al., 2022).
Here, we coupled the Energy Exascale Earth System Model (E3SM) Land Model (ELM; Burrows et al., 2020) to the reaction network simulator PFLOTRAN (Hammond et al., 2014) via the application programming interface (API) Alquimia (Andre et al., 2013) to enable simulations of flexibly defined reaction networks and robust representation of oxygen and TEA concentrations, mineral precipitation and dissolution, and chemical interactions with organic matter cycling within an LSM.As a demonstration of the model framework, we simulate the effect of tidal cycling on subsurface oxygen and salinity concentrations as well as sulfur cycling in tidal wetland soils, and we compare simulated production and surface emissions of carbon dioxide and methane across gradients of salinity.

Biogeochemical Reaction Network
We implemented a network of reactions including soil organic matter (SOM) decomposition and aqueous redox chemistry in PFLOTRAN (Hammond et al., 2014), building on previous PFLOTRAN implementations of redox biogeochemistry applied to Arctic soils (Sulman et al., 2022) and coastal wetland sediments (O'Meara et al., 2024;Wang et al., 2024).SOM decomposition reactions were implemented in the PFLOTRAN Reaction Sandbox (Hammond, 2022) using the same decomposition kinetics used in ELM SOM and litter decomposition calculations (G.Tang et al., 2016).ELM litter and SOM decomposition follows pseudo-first-order kinetics with nutrient limitation according to a "converging trophic cascade" (CTC) framework (Burrows et al., 2020;P. E. Thornton et al., 2002) (Figure 1, upper right).Litter and coarse woody debris (CWD) pools decompose into SOM pools with fixed decomposition time scales (modified by temperature and moisture) and fixed C:N ratios.N mineralization or immobilization is determined by the relative C:N ratio of successive pools and the fraction of the pool C that is converted to CO 2 during a decomposition transition.Organic matter pools decompose as a solidstate process, transforming from one solid organic matter pool to the next with associated production of mineralized N and CO 2 .While ELM can simulate phosphorus (P) as well as N cycling (Yang et al., 2014), the current PFLOTRAN reaction network for decomposition omits P. Inorganic N pools (NO 3 and NH 4 + ) are also tracked, accounting for N mineralization and immobilization as well as plant root N uptake.Root N uptake rates are calculated based on plant N demand and a Michaelis-Menten function of inorganic N availability.
To incorporate dissolved oxygen consumption and aqueous-phase redox reactions into the reaction network (Table 1; Figure 1, lower portion), decomposition of litter and SOM pools was modified so that the decomposed fraction previously converted directly to CO 2 was converted instead to DOM with a fixed C:N ratio of 20.Rate constants for decomposition of litter and SOM pools to DOM under anoxic conditions were assumed to be 10% of the default values under oxic conditions, representing the decreased activity of hydrolytic enzymes under anoxic conditions (Kristensen et al., 1995).Multiple aqueous-phase chemical reactions were added representing alternative pathways of DOM decomposition, with liberated N, Fe, and sulfate content of organic matter included based on fixed stoichiometry of DOM (C:N:S:Fe = 2000:100:20:1), based on measurements of C, S, and Fe content Spartina alterniflora litter from Massachusetts sites (Breteler et al., 1981) and a global synthesis showing a median plant litter Fe concentration of 0.2 g kg 1 (Peng et al., 2023).In addition to aerobic decomposition of DOM, which consumes oxygen, anaerobic reactions including fermentation, iron reduction, sulfate reduction, and methanogenesis are included in the reaction network.Methane oxidation reactions using oxygen, iron, and sulfate are also included.Following previous applications of this framework (Sulman et al., 2022), redox reactions are implemented as multi-Monod type reactions that can include both substrate and inhibition interactions: where R is reaction rate (mol (L H 2 O) 1 s 1 ), V max is temperature-dependent maximum reaction rate (mol L 1 s 1 ), N is the set of reactant species (including substrates and terminal electron acceptors), M is the set of inhibiting species, C S N is the concentration of the Nth substrate, K S N is the half-saturation constant of the Nth substrate, C I M is the concentration of the Mth inhibiting species, and K I M is the inhibition constant of the Mth inhibiting species.Reaction stoichiometries, rates, half-saturations, and inhibition species are shown in Table 1.
Inhibition is used to prevent anaerobic reactions from occurring in oxic soil layers, and to represent the dependence of fermentation on pH and buildup of acetate concentrations.The reaction network does not include direct inhibition of redox reactions by the presence of alternative electron acceptors (e.g., inhibition of iron reduction by sulfate or inhibition of methanogenesis by Fe(III)), apart from oxygen.Rate constants and half-reaction parameters built on values used for the earlier implementation of the reaction network (Sulman et al., 2022) or used literature values where available, as specified in Table 1.Values for parameters that could not be directly constrained using literature data were estimated based on rates relative to similar reactions in the network.Because translation of rate constants from specific laboratory or field measurements to a complex system can be inexact, some calibration adjustments were applied to parameters based on the observed patterns of solute concentrations and surface fluxes used for model comparison (Section 2.7).
All aqueous reactions have a temperature-dependent reaction rate via an Arrhenius relationship: where V 0 is the maximum reaction rate at reference temperature, Ea is activation energy, R is the ideal gas constant (8.314J mol 1 ), and T is temperature in C. Ea was set to 80 kJ mol 1 (approximately a Q10 of 3.0 at 20°C) for sulfate reduction and methanogenesis, and 50 kJ mol 1 (approximately a Q10 of 2.0 at 20°C) for other reactions reflecting the higher temperature sensitivity of methanogenesis and sulfate reduction relative to aerobic respiration (Inglett et al., 2012).
PFLOTRAN solves for the mass balance of each component according to the stoichiometric relationships defined for all reactions, including kinetic (Table 1) and equilibrium reactions.pH is tracked dynamically from the appropriate proton balance of aqueous-phase biogeochemical reactions and mineral precipitation/dissolution, incorporating aqueous speciation as part of the solution, for example, CO 2 /HCO 3 and H 2 S/HS partitioning.The biogeochemical conceptual model incorporates key aqueous complexation (e.g., carbonates, sulfides, etc.) and mineral precipitation-dissolution (pyrite, Fe oxides, etc.) reactions that buffer the system with respect to pH.Fermentation has a net acidifying effect due to proton release, as do sulfide oxidation and pyrite dissolution.Fe (III) reduction causes a net increase in pH, because it is coupled to proton-absorbing dissolution of Fe oxide minerals.Note that while Fe(III) reduction is modeled as an aqueous-phase reaction, the very low solubility of Fe (III) means that Fe(III) reduction effectively occurs as reductive dissolution of Fe(OH) 3 and/or Goethite mineral species (Sulman et al., 2022).

Coupling via the Alquimia Interface
PFLOTRAN is coupled to ELM via the Alquimia interface (Andre et al., 2013), which is designed as a standardized application programming interface (API) for incorporating existing third-party biogeochemistry codes within environmental transport models.Alquimia has previously been used to connect the Advanced Terrestrial Simulator (ATS) model with PFLOTRAN for watershed-scale reactive transport simulations (Jan et al., 2021;Molins et al., 2022;Xu et al., 2022).Alquimia organizes key chemical information into mobile and immobile (sorbed) concentrations of solutes, as well as volumetric fractions of minerals.The API also includes functions for initialization, equilibration of initial and boundary conditions, and time stepping the geochemical model.Here, we implemented the Alquimia API within ELM.Alquimia initialization and initial condition equilibration subroutines were added to the ELM initialization code, and the Alquimia time stepping subroutine was added to the ELM code as described below.PFLOTRAN input and database files are read as part of the initialization process to specify the chemical species, reaction network, and reaction parameters such as rate constants, inhibition factors and thermodynamic equilibrium constants.
ELM represents key carbon and nitrogen pools including multiple litter and SOM pools as well as soil nitrate and ammonium.These pools are all represented in the PFLOTRAN reaction network used in these simulations, building on previous work to represent ELM decomposition processes in PFLOTRAN (G.Tang et al., 2016).We modified the Alquimia interface to treat solid-state SOM pools as immobile chemicals within the Alquimia data structure, allowing transparent data transfer of SOM pools from ELM to PFLOTRAN and back via the interface.This structure allows the ELM decomposition processes to be fully replaced by equivalent or modified calculations on the PFLOTRAN side, updating C and N concentrations and maintaining C and N mass balance while enabling interactions with reaction networks of arbitrary complexity as determined by the PFLOTRAN input file, provided that all ELM SOM C and N pools are included in the PFLOTRAN reaction network.
Coupling within the ELM-PFLOTRAN framework is modular, with ELM storing the state variables (e.g., concentrations) while PFLOTRAN calculates chemical transformations.However, only data that are directly relevant to ELM state (primarily organic matter and nutrient pools) are translated into ELM data structures that are visible to other model components.This allows representation of different reaction networks to have minimal effects on other parts of ELM code, and allows simulation of different reaction network configurations and complexities without any changes to ELM code specific to a particular reaction network configuration.

Vertical Gas and Solute Transport
ELM-PFLOTRAN employs operator splitting for reactive transport: ELM simulates the one-dimensional gas and solute transport within vertical columns and calls PFLOTRAN (through Alquimia) to solve the zero-dimensional biogeochemistry for each layer in the one-dimensional column.Vertical advection-diffusion is implemented using the finite volume approach of Patankar (1980).The current gas diffusion implementation does not divide soluble gases into dissolved and gas phases, but instead treats them as solutes with a higher diffusion rate in unsaturated soil layers.Diffusion coefficients are set separately for gas and non-gas solutes.Gas diffusion coefficient decreases with increasing water saturation based on (Fan et al., 2014): where D a is aqueous diffusion coefficient (m 2 s 1 ).At the beginning of the column calculation, gas concentrations in the top layer are assumed to be in equilibrium with the upper boundary layer concentrations if the top layer is unsaturated.Vertical advection of solutes is calculated by assuming that vertical flow downward from each layer is equal to the subsurface drainage flow rate of the column as calculated by ELM.Solute concentrations in downward vertical flows into the top layer are determined by the upper boundary condition.Vertical flow out of the bottom of the deepest soil layer is assumed to be zero.
Ebullition is included as a transport pathway for dissolved gases, following (Wang et al., 2024).Pressure in each layer is calculated using the weight of water in layers above, including atmospheric pressure.Partial pressure of each dissolved gas is calculated based on a temperature-dependent Henry's law relationship with a gas-specific Henry's Law constant (see Table S1 in Supporting Information S1): where P g is partial pressure of gas g, C g is concentration of gas g (mol m 3 ), H g is the Henry's Law constant for gas g (mol m 3 Pa 1 ), and H T,g is the temperature dependence of solubility for gas g (K 1 ).If partial pressure of a dissolved gas exceeds the ambient pressure, the excess concentration is removed from the layer, reducing the gas concentration in the lower layer to the saturation value.90% of the excess is moved upward one layer, thus assuming that bubbles can be re-dissolved in unsaturated upper layers.The remaining 10% is emitted to the atmosphere, representing a fraction of bubbles that move more rapidly to the surface.This process is conducted starting in the bottom layer and moving up the profile.

Time Stepping Approach
Vertical transport and chemical reactions are calculated with an operator splitting approach using Strang splitting (Strang, 1968) to reduce truncation error related to operator splitting (Carrayrou et al., 2004).A variable time stepping approach is used to account for failure of the chemical reaction simulator to converge to a valid solution when the simulated time step is too long compared to the time scale of chemical reactions, or when consumption of gases (e.g., O 2 ) is high enough that transport calculations at that time step will underestimate gas concentrations.One half time step of vertical transport is calculated first, and gas concentrations in the surface soil layer are equilibrated with the upper boundary condition.Next, chemistry is updated via Alquimia/PFLOTRAN for each soil layer, starting at the top.If any soil layer fails to converge to a valid solution, then concentrations in all layers are reset and the time step is cut in half.When the top layer is unsaturated, a reduction of greater than 25% in dissolved oxygen concentration in the top layer (which is assumed to be near equilibrium with the atmosphere) is also treated as a nonconvergence condition, because it indicates that the current time step length cannot accurately capture the rate of oxygen consumption and/or transport in the column.The column reactive transport calculations, and potential shortening of the time step, are repeated recursively until chemistry in all layers can be successfully updated.Then, the second half time step of vertical transport is calculated.The shortened time steps are repeated appropriately to ensure that the total integration matches the ELM time step (60 min in our simulations) because ELM does not natively support flexible time stepping.

Tidal Forcing
Tide-driven lateral flows into and out of the soil column built on previous work focused on boreal peatland microtopography (Shi et al., 2015) and initial implementation of tides in coastal systems that used hydrologically coupled soil columns and a sinusoidal tidal pattern (O'Meara et al., 2021).Lateral flows and tidal-driven exchange of water and solutes in this approach use a hydrological boundary condition determined by the relative height of water in a tidal channel compared with water table height in the wetland soil column.We extended the previous lateral flow implementation, which used a single lateral flow time scale, to include rapid horizontal flow when the water table or tidal water level was above the soil surface to equilibrate the surface water depth in the wetland to the tide height.
where Q surf is horizontal surface water flow into the wetland column (mm s 1 ), z tide is height of water in the tidal channel (mm, relative to wetland soil surface height), z surf is surface water height in the wetland (mm, defined as zero when water table is below the surface), and k surf is a rate constant representing the time scale of surface water transfer as a function of the difference in surface water height, set to a rapid flow to so that surface water level is close to equilibrium with the tidal forcing (7 × 10 5 s 1 ).Consistent with the previous lateral flow implementation, a slower drainage flux allows water to flow into or out of the soil column during low tide conditions when the water table and tide height are below the surface: where Q subsurf is horizontal net subsurface water flow into the wetland column, z WT is water table depth in the wetland subsurface (defined as < 0), and k subsurf is the rate constant for subsurface net flow, calculated using the mean saturated hydraulic conductivity of the column (Shi et al., 2015).In addition, ELM calculates a subsurface drainage flow rate as a function of water table depth: where Q drain is net subsurface drainage rate and f ice is an increasing function of mean column ice fraction (accounting for decreased drainage through frozen layers).
Horizontal flows in ELM (including Q surf , Q subsurf , and Q drain ) are currently calculated using a "bucket" approach that is not fully integrated with vertical flow.Vertical flows are calculated first, according to a Richards Equation approach.Next, total horizontal water outflow during the time step is removed from the column by subtracting water content from each layer one at a time, moving downward starting from the water table.Conversely, water flowing into the column is added to the layer above the water table until it reaches saturation, with the process repeated moving upwards by layer until the appropriate total amount of water has been added to the column.Because horizontal flows were not fully integrated into the ELM calculations for vertical flow within the column, the combined hydrology did not yield reasonable results for solute transport.Therefore, we represented vertical transport of solutes assuming that vertical flows balanced subsurface drainage: where Q vert is vertical flow out of the layer, z is layer depth, and z max is depth of the bottommost layer.
Lateral inflow as well as infiltration during flooded conditions are assumed to have the solute concentrations of the tidal boundary condition, which is supplied as salinity concentration in an external forcing data set.Sulfate concentration is assumed to equal 14% of the concentration of chloride (on a per mass basis).pH is calculated using a linear approximation of pH = 6.0 for fresh water and pH = 8.0 for saltwater with a salinity of 30 ppt.
Because comprehensive concentration data for all compounds in the reaction network were not available for the tidal boundary condition, salinity and sulfate are exchanged horizontally via tidal flows while other solutes (including nitrogen) are assumed to stay primarily in the soil column.Specifically, when calculating vertical and lateral transport 10% of the mass of solutes without a defined freshwater/saltwater boundary condition (i.e., excluding pH, salinity, and sulfate) was available for transport and leaching while the remaining 90% remained in the soil layer.This estimated soluble fraction approach was necessary to prevent excessive leaching of nutrients out of the subsurface.Excessive nitrogen leaking could potentially be addressed by incorporating dissolved nitrogen into the boundary condition or representing sorption of ammonium on soil surfaces in the reaction network.

ELM Simulations
Simulations used a standard ELM spinup process (Thornton & Rosenbloom, 2005) of 100 years of accelerated decomposition spinup followed by 300 years of regular spinup and 150 years of transient (historical) simulation.

Journal of Advances in Modeling Earth Systems
10.1029/2023MS004002 SULMAN ET AL.
Atmospheric forcing, including temperature and precipitation, used downscaled Global Soil Wetness Project Phase 3 (GSWP3) meteorology for the Plum Island Ecosystems site, repeated as necessary for spinup.Tidal forcing used sinusoidal tide constituents available from NOAA Tides and Currents for the Plum Island low marsh site (Station ID 8441241), with reference height corrected so that tidal height was defined relative to the marsh surface.
To test the role of salinity and associated S cycling on biogeochemistry and greenhouse gas fluxes in the model, we compared three model configurations for simulating the low marsh ecosystem.All models included the same tidal hydrology patterns.In the Fresh configuration, salinity in tide water was set to zero.In the Saline configuration, salinity in the tide water used measured concentrations from the tidal forcing data set, which ranged from 24 to 35 ppt.In the Saline + reduced GPP configuration, the same saline tide water concentrations were used, and gross primary production (GPP) was additionally reduced as a function of tidal salinity level to represent the impact of saline conditions on plant productivity: where f(s) is the salinity effect on root water uptake resistance (varying between 0 and 1), s is tidewater salinity (ppt), μ is the optimal salinity (0 ppt), and λ is the salinity tolerance (20 ppt), based on observed salinity responses of Spartina alterniflora (LaFond-Hudson & Sulman, 2023;Vasquez et al., 2006).This parameterization yielded a 30%-35% reduction in mean daily GPP when salinity was taken into account.

Comparison With Measurements
Model simulations of geochemical processes were evaluated by comparing simulated profiles of salinity, sulfide, and DOC concentrations to measurements from low marsh sites in the Plum Island Ecosystems Long Term Ecological Research (PIE LTER) monitoring program (A. Giblin et al., 2021;Morris & Sundberg, 2023).The Law's Point (LP) site is a low salt marsh dominated by Spartina alterniflora along the Rowley River, Rowley, MA ( 70.84246 E, 42.73174 N) where porewater measurements were collected monthly using porewater diffusion samplers 12 m from the creek bank from May to October 1999-2022.Shad Creek is a Spartina dominated salt marsh, also along the Rowley River ( 70.8381692 E, 42.7344368 N) where porewater measurements were collected monthly in May-October 2017 using diffusion samplers 10 m from the creek bank.The Typha site is also located along the Parker River ( 70.914799 E, 42.750869 N) and represents a brackish marsh dominated by Typha sp, where diffusion samplers were installed 4 m from the creek bank.
In addition, porewater carbonate chemistry data was collected in a creek bank dominated by S. alterniflora in the Nelson Creek catchment in July 2014 (Data Set S1).Water was pumped from depths of 10 and 20 cm using a piezometer and a peristaltic pump.Samples for dissolved inorganic carbon were placed in vials, overflowed several times, capped, stored at 4°C and run on an Apollo SciTech DIC analyzer within 24 hr from sampling.Sample pH was immediately measured in the field using a combination electrode.Measurements of conductivity were taken in the lab within a day.Simulated soil organic matter concentration profiles were compared with measured profiles from PIE LTER marsh sites (Spivak, 2020;Spivak et al., 2018).Uncertainty ranges in measured SOM profiles were calculated using the standard error of the mean over three replicate profiles.

Simulated Hydrology
The tidal-driven hydrological model, paired with a time series of hourly tide height and salinity, allowed the model to simulate patterns of surface water depth across diel tidal fluctuations as well as longer-term variations in tide height (e.g., spring/neap tide cycles; Figures 2a and 3a).The hydrological forcing represents the low marsh flux tower site, where the water table generally stays close to the surface even during low-tide conditions (Figure 2b).Simulations showed unsaturated soil conditions in the top 50 cm of the soil profile during low tides (Figures 2a and 3b), while observed water table during low tides generally stayed within 5-10 cm of the surface.However, the deeper unsaturated layers in the model remained quite wet (>85% of saturation below 20 cm), limiting oxygen concentrations below 5-10 cm depth (Figure 3d).The tidal inundation cycle was visible in soil oxygen concentrations, with oxygen in the subsurface depleted rapidly during flooded periods.
Salinity increased with depth in the top 50 cm and varied with tidal fluctuations, reflecting the influence of salinity inputs at the soil surface and transport into the subsurface, combined with the solute concentrating effect of water removal from the root zone driven by transpiration (Figures 2c and 3c).Simulated salinity matched well with the observed range for the Law's Point low marsh site (Figure 2d), while the Shad site ranged between 10 and 15 ppt and the Typha site salinity generally stayed below 10 ppt.Simulated subsurface temperatures at the site ranged from below freezing during winter to 20°C during summer (Figures 2e and 2f), with the greatest range of temperatures occurring near the surface.Simulated soil temperatures were about 5°C lower than the observed range, likely due to mismatch between the GSWP3 gridded atmospheric forcing data and local site conditions along with the model's lack of heat transfer by water flow into the model soil column.

Simulated Redox and Sulfur Cycling
The reaction network in the model connected carbon, sulfur, and oxygen cycling in the subsurface (Figure 1) and responded to seasonal and tidal cycles.Sulfate concentration in water entering the soil profile through infiltration or lateral flows was assumed to be proportional to salinity, leading to sulfate profiles that qualitatively resembled salinity profiles, both increasing with depth and peaking at 1 m (Figures 4a and 4c; Figures 2c and 2d).Sulfate reduction produced sulfide in anoxic layers, driving a sulfide concentration profile that peaked at a shallower 30 cm depth (Figures 4d and 4f).Simulated sulfide concentrations were within the range of observed values in shallow layers, with observed profiles ranging from 0.1 to 4 mM and simulated profiles in the saline configuration ranging up to 1.5 mM.However, simulated sulfide concentrations declined with depth below 30 cm and did not match the continuing increase in observed sulfide concentrations at depths down to 1.5 m.Sulfate concentrations were much lower in the freshwater configuration than in the saline configuration, driving differences in subsurface biogeochemistry (Figures 4b and 4e) which were consistent with sulfide measurements from the Typha site.
Sulfate reduction consumed DOC, lowering DOC concentrations in layers below 10 cm depth in saline relative to fresh simulations (Figure 4i).Simulated DOC concentrations in the saline simulation were within the range of measured concentrations (0.1-1 mM), although the peeper measurement technique may underestimate DOC concentrations (A.Giblin, personal communication).While measured DOC concentrations in the brackish Typha site were higher than those in the saline site, the model overestimated freshwater DOC concentration, especially in shallow layers.DOC concentrations increased in spring and declined in summer as sulfate reduction and methanogenesis rates increased.
DIC concentrations (including dissolved CO 2 and CH 4 ; Figures 4j-4l) increased in July and August as DOC was depleted.The magnitude of DIC concentration was higher in saline than in freshwater simulations.DIC concentrations in the saline simulation were higher than salt marsh measurements but within the same order of magnitude, and the simulated increase in DIC concentration from 10 to 20 cm depth was consistent with observations.Subsurface methane concentrations (Figures 4m-4o) also increased over time as DOC was depleted.Methane concentrations were low near the surface and increased in deeper layers, reflecting the predominance of methane production in more reducing subsurface layers and the consumption of methane in more oxidizing layers.Along with oxygen and Fe(III), sulfate also served as a substrate for methane oxidation, which lowered subsurface methane concentrations in the saline simulation compared to the fresh simulation.pH was lower in the freshwater Observed water levels were shifted by 9 days to account for a temporal mismatch in tide timing between model forcing and observations.simulations than in the saline simulation, and declined as DIC concentrations increased (Figures 4p-4r).The model substantially overestimated porewater pH in the saline simulation.
Peak methane efflux was 40 times higher in the fresh simulation than in the saline simulation (Figure 5a), with a seasonal cycle increasing rapidly in spring and continuing through the fall.The magnitude of methane efflux was consistent with chamber measurements from low-salinity sites in Cape Cod.The very low methane fluxes from the saline simulations were consistent with the magnitude of fluxes measured from the low marsh flux tower as well as chamber measurements from the Cape Cod sites (Figure 5b).Saline simulations in which vegetation productivity was also reduced had 43% lower surface methane emissions and were closer in magnitude to observed fluxes than saline simulations without reduced vegetation productivity.Soil CO 2 fluxes (excluding autotrophic root respiration) were slightly higher for the saline simulation than for the fresh simulation (Figure 5c).CO 2 fluxes from the reduced GPP simulation were 25%-50% lower than for the saline simulation without reductions in GPP, indicating the impact of reduced C inputs to the system.Simulated loss of DIC through lateral tidal flows were about 15%-20% the magnitude of surface CO 2 efflux, and were similar in magnitude to methane emissions in the freshwater simulation and much higher than methane emissions from the saline simulations (Figure 5d).Soil organic carbon (SOC) concentration profiles (Figure 6) declined with depth below about 10 cm, reflecting more rapid decomposition in periodically unsaturated surface layers and reduced litter inputs in deeper layers.C concentrations were slightly higher in the saline configuration compared to the fresh configuration, despite higher DIC production in the saline simulations.When salinity was combined with reduced GPP, SOC concentrations were about 50% lower.Observed SOC concentrations were generally lower than simulated SOC concentrations, but the simulation with high salinity and reduced GPP approached the magnitude of observations at depths below 50 cm.Unlike simulated SOC, observed SOC concentrations did not decline significantly with depth.
In addition to carbon and sulfur, the reaction network also included iron redox cycling, allowing simulations of iron oxide and iron sulfide mineral precipitation and dissolution (Figure 7).Fe(II) concentrations were much higher in the freshwater than in the saline simulation (Figure 7b), reflecting the rapid precipitation of Fe(II) into iron sulfide minerals in the more sulfidic saline simulation.Fe(II) concentrations in the fresh simulation peaked around 1 m depth.Fe(II) concentrations were generally consistent with measured concentrations from a freshwater mesocosm study and a saline marsh (Fettrow, Jeppi, et al., 2023;Fettrow, Vargas, & Seyfferth, 2023), although observed salt marsh Fe(II) concentrations were somewhat higher than the near-zero concentrations in the saline simulation and we were not able to compare across the full range of depths.Iron oxide concentrations were high in the more oxic surface layers and lower in deeper layers (Figure 7d).The fresh simulation showed net iron oxide loss to dissolution in anoxic layers (Figure 7c), while the saline simulation had iron oxide accumulation in shallower layers.Surface iron sulfide concentrations were higher in the saline than in the fresh simulation, but iron sulfide concentrations in the fresh simulation were higher deeper in the profile (Figure 7f).Net iron sulfide precipitation rates were highest in shallower layers in the saline simulation, reflecting higher available iron (Figure 7e).While both FeS and pyrite were included as iron sulfide minerals in the reaction network (Table 1), in model simulation results the accumulated iron sulfide minerals were entirely composed of pyrite.

Biogeochemical Insights for Modeling Coastal Wetlands
Our model framework successfully incorporated tidal-driven hydrology, redox biogeochemistry, and pH dynamics into a full-featured LSM, allowing simulations of variations in subsurface biogeochemical cycling driven by rapid hydrological fluctuations in the context of carbon cycling.Model simulations connected higher sulfate concentrations in saline wetlands to lower DOC and higher DIC concentrations along with greatly reduced methane emissions.The higher fluxes from the freshwater simulation were consistent with (Sanders-DeMott, Eagle, Kroeger, Wang, et al., 2022), who found that methane flux increased strongly along a saline to freshwater gradient in coastal wetlands.The ability to simulate suppression of methane production under salinization is key to accurately predicting coastal wetland greenhouse gas balance (Kirwan et al., 2023).Simulated SOC concentrations were slightly higher in the saline simulation than in the fresh simulation, despite the role of sulfate as a terminal electron acceptor.The difference is likely due to a slightly (2%) higher simulated GPP in the saline simulation compared to the fresh simulation, possibly driven by biogeochemical interactions with nutrient availability.When the impact of salinity on GPP was taken into account by reducing GPP, SOC concentrations were much lower.This result suggests that sulfate reduction alone may not be sufficient to explain differences in soil carbon patterns between saline and freshwater wetlands, and that plant-soil feedbacks may be necessary to explain contrasts.Plant feedbacks have been hypothesized to play a major role in peat collapse associated with salinization (Chambers et al., 2019).However, there is also evidence that seawater additions can enhance SOC mineralization (Chambers et al., 2011).While the reaction network included both FeS and pyrite as iron sulfide minerals, only pyrite accumulated in simulated sediments.This was consistent with observed rapid pyrite formation in salt marshes (A.E. Giblin, 1988) as well as theoretical modeling indicating that pyrite can form in as little as 3 hours and does not necessarily require intermediate precipitation of FeS (Rickard, 2019).

Value of Simulating Detailed Biogeochemical Interactions in LSMs
While many existing LSMs, including the E3SM Land Model, do include methane production and emission calculations (Riley et al., 2011;Wania et al., 2013), our simulations highlight the potential importance of more complex interactions in determining decomposition and greenhouse gas production.Previous studies have identified substrate limitation as a driver of seasonal patterns in methane production (Chang et al., 2020).pH dynamics can also influence methane production, both by direct impacts on microbial physiology (Wagner et al., 2017) and by changing the solubility of alternative terminal electron acceptors such as iron (Marquart et al., 2019;Sulman et al., 2022).Here, we demonstrate the capability of simulating substrate dynamics, pH changes, oxygen depletion, and their influences on methane emissions within a full-featured LSM.While we focus on methane production in this analysis, other important processes that this model framework can enable include the phytotoxic effect of sulfide in soils (Koch et al., 1990;Lamers et al., 2013), impacts of drought-driven increases in soil salinity concentration on vegetation and microbial communities, and interactions of pH dynamics with subsurface biogeochemistry.pH is widely considered a critical environmental variable, affecting carbon storage, microbiology, plant growth, and nutrient availability in environmental systems (Fierer & Jackson, 2006;Neina, 2019).Yet, dynamic pH is not included in current LSM frameworks.Thus, the ability to simulate dynamic pH in an LSM represents a significant step forward.That said, the current model implementation overestimated pH in salt marsh sediments compared to observations.Accurately simulating pH dynamics is challenging due to multiple interacting processes that affect pH, and additional model improvements such as better representation of soil pH buffering capacity will be necessary to improve the accuracy of the pH component.
The incorporation of subsurface DIC concentrations and DIC loss in runoff is an important step forward toward representing the carbon balance of coastal wetland ecosystems, where lateral export of DIC and total alkalinity can be an important component of the net carbon balance, with total inorganic alkalinity export representing a long-term carbon sink in the ocean (Reithmaier et al., 2021;Yau et al., 2022).However, the current biogeochemical parameterization has not been evaluated in detail for the accuracy of DIC speciation (i.e., what fraction of DIC is in the form of bicarbonate vs. carbonate and aqueous CO 2 ) and will need attention to other elemental cycles such as calcium to produce accurate estimates of total inorganic alkalinity production.In any case, the coupling of ELM to a detailed reaction network simulator provides the technical capability for incorporating total alkalinity production and balance into a land surface model.
While this paper focused on simulating redox dynamics, the reactive transport framework used in the model implementation builds the groundwork for a wide range of applications.PFLOTRAN includes a broad set of geochemical reaction capabilities, including microbially-mediated as well as abiotic aqueous reactions, dissolution and precipitation of different types of minerals, and sorption of solutes onto mineral surfaces (Steefel et al., 2015).The direct coupling of PFLOTRAN chemistry into ELM means that any geochemical reactions implemented in PFLOTRAN can be directly incorporated into land model simulations with minimal edits to land model code.Thus, this framework could be easily adapted to facilitate various applications including testing different SOM decomposition reaction networks, simulating dynamics of inorganic carbon storage and release from carbonate minerals, and cycling of micronutrients within the soil.PFLOTRAN's Reaction Sandbox, which allows for customized geochemical formulations to be implemented in PFLOTRAN code (Hammond, 2022), opens broad possibilities for testing geochemistry and biogeochemical interactions within a coupled ELM-PFLOTRAN system.Furthermore, the implementation of the ELM coupling using the Alquimia API opens the possibility of coupling ELM other reactive transport codes that are compatible with the API and may have different reaction simulation capabilities.
Coupling with the vegetation component of LSMs can also support new insights from subsurface biogeochemical modeling.For example, subsurface redox state is strongly dependent on vegetation productivity, oxygen transport, and organic matter deposition (Noyce et al., 2023), and biogeochemical responses of coastal wetlands to sea level rise are sensitive to plant-soil interactions (Mueller et al., 2020;O'Meara et al., 2024).Quality and quantity of litter inputs are strong controls on subsurface biogeochemical cycling that can vary with productivity and shifts in plant functional types along key gradients such as the marsh-mangrove ecotone (LaFond-Hudson & Sulman, 2023;Steinmuller et al., 2020).High rates of transpiration in coastal wetlands can draw down water table, create aerated zones, and concentrate salts (Xin et al., 2022), with halophytic plants like mangroves having particularly high potential to salinize sediments (Passioura et al., 1992;Reef & Lovelock, 2015).Direct coupling of a biogeochemical subsurface model like PFLOTRAN to a LSM capable of simulating evapotranspiration and its dependence on plant function allows exploration of aboveground-belowground interactions between plantmediated surface water balance and subsurface biogeochemical cycling.

Areas for Improvement of Model Implementation
The current model lacks a full set of boundary conditions for solutes in tidal flows, which are currently limited to salinity, sulfate, and pH.A major limitation of this approach is lack of nutrient inputs from surface water, which could lead to underestimated vegetation productivity.The lack of full solute boundary conditions including major cation and anion concentrations also makes it difficult to accurately quantify pH, DIC, and DOC dynamics of the simulated wetland.Future applications of this model framework would benefit from developing a full set of solute boundary conditions in river and tidal waters.
The one-dimensional representation of subsurface hydrology in ELM posed challenges for directly integrating reactive transport into ELM.The ELM hydrology model was designed primarily for simulating grid cell water balance and water limitation of vegetation.Lateral flows, including subsurface drainage and tide-driven lateral flows, are not fully integrated into the hydrological solver.Rather, the model calculates vertical redistribution using a Richards equation approach and afterward adds or removes water associated with lateral flows using a filling/emptying bucket approach.This causes calculated lateral and vertical flow rates to be inconsistent with the full water balance, leading to unrealistic salinities due to flow convergence and high flow velocities within the column when using ELM-simulated water flow rates directly for reactive transport calculations.In the current study, we ultimately replaced the internally calculated vertical flow rates with approximate flows that were consistent with subsurface drainage.The lack of full solute boundary conditions could also have contributed to unrealistic results when using internally calculated flow rates.Further work in this area could benefit from fully integrating lateral flows into the ELM hydrological model and simulating lateral flows using hydraulic head boundary conditions rather than height differentials.As an intermediate step, 3-dimensional simulations of hydrologic flows in coastal wetland sediments could be used to inform the parameterization of column-scale hydrological exchanges in ELM.Highlighting the importance of subsurface flows, observed PIE LTER solute concentration profiles were highly dependent on the hydrological flow regime, particularly flushing of the subsurface.For example, sulfide concentrations measured at 4 m from the Shad site creek bank were an order of magnitude lower than those measured 10 m from the creek bank due to more frequent tidal flushing (A.Giblin et al., 2021).
The current lateral flow implementation imposes hydrological flows as boundary conditions on the ELM column and does not fully integrate hydrological exchanges or solute flows with other components of the E3SM, such as the river model (MOSART) and ocean model (MPAS-Ocean) (Golaz et al., 2019).Fully integrating coastal wetland processes into a fully coupled ESM will require coupling hydrological and solute exchanges across model components so that water, carbon, and other quantities can be conserved in large-scale simulations.The prescribed boundary condition approach used here builds the groundwork for incorporating these exchanges into the ESM coupler framework.
The model underestimate of soil temperature relative to measured site conditions (Figure 2) also indicates potential issues with meteorological forcing data.We used a global gridded atmospheric forcing product (GSWP3) to drive model simulations.Such spatially coarse reanalysis products often require correction to match finer scale variations in meteorological conditions consistent with specific field sites (Berg et al., 2003).Future model-data comparisons could be made more robust by applying site-specific corrections to meteorology or using higherresolution products that assimilate weather stations (e.g., Daymet; Thornton et al., 2021).For these initial simulations, we chose to use a readily available, gap-free product.Given the strong temperature dependence of biogeochemical processes in the model, particularly sulfate reduction and methanogenesis, the underestimate of temperature likely drove underestimates of maximum biogeochemical process rates during the summer.Based on the model's temperature dependence (Equation 2), a temperature increase from 17 to 22°C would increase the rate of sulfate reduction or methanogenesis by 75%, or the rate of other reactions in the model by 42%.Note, however, that increasing rates of multiple reactions would likely not express this full change in the total rates due to substrate limitation.Some aspects of the biogeochemical reaction network omit or simplify potentially important processes.First, the reaction network does not explicitly simulate changes in microbial community structure, biomass, or metabolic capabilities.Instead, the model assumes that key parameters of biogeochemical reactions (i.e., maximum rate constants and half saturation constants; Table 1) do not change, and that actual reaction rates are driven by instantaneous solute concentrations.Previous measurements have demonstrated that wetland microbial communities can shift in response to changes in salinity, redox state, and substrate availability, potentially driving shifts in microbial community-level metabolic capabilities over space and time (Dang et al., 2019;Lambais et al., 2008).In our model framework, such shifts could be expressed either by dynamically changing reaction rate and half saturation parameters in response to integrated soil conditions, or by explicitly simulating microbial community structure by assigning microbial functional types to different reactions and coupling those reactions to dynamic microbial biomass.
In addition, the PFLOTRAN implementation assumes that redox reactions occur in the aqueous phase, which is not accurate for reductive dissolution of Fe oxides via direct microbial interactions with mineral surfaces (Kappler et al., 2021;Roden & Zachara, 1996).The very low solubility of Fe(III) means that in practice Fe(III) reduction is tightly coupled to Fe oxide mineral dissolution; however, this approach could cause modeled Fe(III) reduction to be overly sensitive to pH and its effect on Fe oxide solubility, as has been suggested in a previous application of this PFLOTRAN framework (Sulman et al., 2022).Updates to the PFLOTRAN framework to allow direct reduction of Fe oxides would help to alleviate this issue.The erroneously high porewater pH in saline simulations also lowers Fe oxide solubility, potentially biasing simulated Fe cycling in salt marsh sediments.Improvements to pH in the model should be a priority for future coastal wetland applications.
The implementation of gas transport in the soil column could also be improved.The current implementation includes moisture-dependent gas diffusion as well as a simple implementation of ebullition but does not include plant-mediated gas transport.Plant-mediated transport can be an important pathway for both methane transport out of the soil and oxygen transport into the soil, especially for aerenchymous plants (Colmer, 2003;Jeffrey et al., 2019;Noyce et al., 2023).Planned work on this model framework will include plant-mediated gas transport, with dependence on plant traits such as aerenchymous tissues and rooting depth distributions at the plant functional type level (LaFond-Hudson & Sulman, 2023).The implementation of ebullition also uses a simple approach that calculates partial pressure separately for each dissolved gas.This approach may underestimate ebullition flux when multiple dissolved gases are produced in the subsurface (e.g., methane, CO 2 , and H 2 S).An improved approach would incorporate gas mixing in bubbles, and we plan to move toward that approach in ongoing work.The current gas diffusion implementation does not explicitly divide soluble gases into dissolved and gas phases, but instead differentiates dissolved gases from non-gas solutes using diffusion coefficients.A two phase (gas and aqueous) transfer scheme that tracked the dissolved fraction of gases in each layer could lead to improved gas transport simulations.
Model parameterization is also a challenge, particularly for increasingly complex biogeochemical reaction networks.Our model parameterization does incorporate field and laboratory measurements of reaction rates and solute concentrations where possible (Table 1), but some parameters are inevitably difficult to constrain.In this initial study, we focused on demonstrating the feasibility of simulating reaction network interactions within a land surface model, and therefore did not evaluate modeled rates in detail.Applications of this framework to predictive modeling of biogeochemical cycling will benefit from additional detailed evaluation of reaction rates and concentrations in the context of porewater concentration and flux measurements.Additional parameterization of soil column hydraulic properties could also help to improve the accuracy of simulated hydrology (Figure 3a).

Conclusions
We coupled a biogeochemical reaction network solver (PFLOTRAN) to a land surface model (ELM) and implemented vertical solute and gas transport as well as tidal-driven inputs of salinity and sulfate.We applied the model to simulate biogeochemical cycling in Massachusetts tidal marshes under either saline or freshwater tidal boundary conditions.The coupled model framework allowed simulations of multiple redox reactions, pH dynamics, oxygen consumption, and methane production and oxidation to be fully integrated within a land surface model.Sulfate supplied in the saline simulation drove high levels of sulfate reduction, which reduced DOC, increased DIC, and greatly lowered subsurface methane concentrations and surface methane emissions.This new model framework builds the foundation for simulating multicomponent biogeochemical interactions in land surface models and demonstrates how directly simulating redox reactions in inundated soils can improve model simulations of organic matter decomposition and greenhouse gas production while building the groundwork for explicit geochemical representation in larger-scale land surface model and Earth system model simulations.

Figure 1 .
Figure 1.Diagram of the biogeochemical reaction network used in the simulations.Pools are shown as circles, color coded by type of pool (including coarse woody debris [CWD], soil organic matter [SOM], and dissolved organic matter [DOM] along with dissolved gases and ions).Arrows indicate transformations via the reactions shown in white ovals.Note that multiple methane oxidation pathways involving oxygen, sulfate, and iron are shown as separate reactions.Nitrogen pools and related reactions are omitted from the diagram for clarity.
estuary complex and nearby low salinity (8-15 ppt) and high salinity (20-30 ppt) salt marsh sites with unrestricted tidal flows in Cape Cod, MA.

Figure 2 .
Figure 2. Water levels, soil salinity, and soil temperature over May-November.Profile plots (panels d and f) show the time-averaged profile (solid line) and the 10th and 90th percentile values (shaded region).Symbols and error bars in panels d and f show observations with 10%-90% percentile ranges.

Figure 3 .
Figure 3. Simulated surface water depth (a), volumetric water content (b), salinity (c), and oxygen concentration (d) over an approximately 1 month time period.Observed water level is also shown as a black line in panel (a).Observed water levels were shifted by 9 days to account for a temporal mismatch in tide timing between model forcing and observations.

Figure 4 .
Figure 4. Simulated biogeochemistry over a growing season.The left column shows profiles over time for the saline simulation, and the middle column shows profiles over time for the freshwater simulation.The right column shows mean profiles (solid lines) for the saline (orange) and fresh (blue) simulation, respectively.Shaded regions show the 10th to 90th percentile range of values.Symbols in (f), (i), and (l) show measured values from the Law's Point (LP), Shad Creek (Shad), and Typha sites.

Figure 5 .
Figure 5. Simulated surface greenhouse gas fluxes over a growing season.(a) Methane fluxes.(b) Same as a, but magnifying the vertical axis so measured and simulated methane flux from the saline site are visible.(c) Soil CO 2 flux.Note that autotrophic respiration is excluded.(d) Lateral flux of DIC out of the soil column.

Figure 6 .
Figure 6.Profiles of simulated soil organic carbon concentrations for the three simulations.Line colors show the three simulations (Fresh, Saline, and Saline with reduced GPP).Shaded regions show the 10th to 90th percentile range of values over one year.Observed data from Spivak (2020).Error bars show standard error across three replicate profiles.

Figure 7 .
Figure 7. Simulated iron and iron sulfide cycling.(a) Fe(II) concentration as a function of depth and time in the freshwater simulation.(b) Mean profiles and 10-90 percentile of Fe(II) in fresh and saline simulations.(c) Mean rate of change of iron oxide mineral concentration.(d) Iron oxide mineral concentration profile.(e) Rate of change of iron sulfide mineral concentration.(f) Iron sulfide mineral concentration profile.

Table 1
Biogeochemical Reactions and Parameters Included in the Reaction Network

Journal of Advances in Modeling Earth Systems
SULMAN ET AL.

of Advances in Modeling Earth Systems
Fettrow, Jeppi, et al. (2023)e concentration measurements available at the PIE LTER field site, we used porewater concentration measurements of Fe(II) collected by Fettrow (2023a) andFettrow, Jeppi, et al. (2023)at the St. Jones National Estuarine Research Reserve in Dover, Delaware as well as DOC and Fe(II) from mesocosm incubations using 8 cm deep soil cores collected from the same site(Fettrow, 2023b; Fettrow, Vargas, & Seyfferth, 2023).
water table, we show model watSanders-DeMott, Eagle, Kroeger, Wang, et al., 2022)ses 85% of saturation.To supplement the PIE LTER flux tower measurements, which were only located in a saline low marsh system, surface methane fluxes were also compared with chamber measurements collected by (Sanders-DeMott, Eagle,Kroeger, Brooks, et al., 2022;Sanders-DeMott, Eagle, Kroeger, Wang, et al., 2022).These chamber flux measurements were collected from high salinity (8-12 ppt) and low salinity (<5 ppt) marsh sites within a diked Journal SULMAN ET AL.