Simulated Biogeochemical Effects of Seawater Restoration on Diked Salt Marshes, Cape Cod National Seashore, Massachusetts, USA

1. Introduction

Historical blockage of tidal exchange within salt marshes through diking reduces or eliminates seawater interaction and resulted in the freshening of surface water and marsh porewater. This diking of coastal wetlands has resulted in water-quality problems such as hypoxia and acidified leachate [1,2,3], the loss of shellfish and fisheries habitat and biodiversity [4,5,6,7], loss of buffer from sea-level rise [8], and stagnated freshwater ponds that host mosquitos and midges [1]. Diking of coastal wetlands has been practiced extensively for hundreds of years in western Europe and the North Atlantic [5]. In southern New England, pre-colonial salt marshes have been ditched, drained, and diked for livestock grazing, hay production, and later insect control to the point where about 50 to 70 percent were destroyed [9,10]. In the contiguous U.S., it has been estimated that 5,000-6,000 km² of wetland or potential wetland space in the intertidal zone has impaired tidal exchange, behind dikes, levees, or transportation infrastructure [11,12,13]. These concerns with diking have led to widespread efforts in salt marsh restoration, including the removal of tidal restrictions. In addition to the reversal of water-quality problems listed above, restoration of disconnected saline tidal flows acts to substantially reduce emission of CO₂ and CH₄ [11,12,13,14]. However, restoration of diked coastal marshes from fresh to saltwater can cause substantial eutrophication from the release of freshwater sediment-derived nutrients and also the release of acidity and iron (Fe) from Fe sulfide minerals [4,15]. Therefore, understanding these diked marsh systems can help inform coordination of restoration efforts to avoid potential problems.

The response of marsh biogeochemistry to tidal restoration depends on the flood-tide depth and duration and their effects on sediment, porewater, and substrate composition [4]. Diked marshes are generally either seasonally flooded or perennially drained depending on the freshwater drainage capacity of culverts and creek systems [5]. Diked marshes with sediments that are seasonally flooded (referred to as “diked, seasonally flooded” or DF) accumulate peat of higher organic content than salt marsh peat [4,16,17] because the inorganic sediment supply normally imported by flood tides is limited [18] and the decomposition in freshwater peats, dominated by methanogens, is slow [19]. Furthermore, porewater nutrient concentrations are low [5] probably as a result of both anoxic conditions, which are conducive to denitrification, and sorption of ammonium (NH₄⁺) and phosphate (PO₄^3-) to surfaces. Seawater restoration in these diked wetlands should cause nutrient release because the increase in SO₄^2- concentrations and corresponding SO₄^2- reducing conditions will enhance organic matter degradation, which increases pH and production of bicarbonate, and subsequently leads to desorption of NH₃ [4]. Diked marshes with sediments that were drained (referred to as “diked, drained” or DD) resulted in air penetration into previously anoxic sediments, subsequent oxidation of organic matter and Fe sulfide minerals, and release of acidity, dissolved Fe, Al, and nutrients including NH₄⁺ and PO₄^3-. DF and DD sediments represent the end members in these diked marsh systems, and understanding the way they respond to seawater restoration is critical in assessing the resulting biogeochemical changes.

Acidic soils are common in estuaries, and their inundation by tidal restoration or sea-level rise can lead to significant amounts of desorption and ion exchange [4,20,21,22]. The major distinction between sorption (adsorption and absorption) and exchange is that geochemical model equations for sorption use the concentration of one chemical only, and ignore the effects of other solutes, whereas exchange involves the replacement of one chemical for another at the solid surface [23,24]. The behavior of PO₄^3- in aqueous systems, for example, can be affected by both ion exchange and sorption processes [25,26]. Geochemical models such as PHREEQC can simulate both sorption/desorption and exchange; sorption/desorption can be modeled as surface complexation reactions (such as the Dzombak and Morel [27] model) and exchange can be modeled using several conventions such as Gaines-Thomas [28]. Sorption (or desorption) is generally used herein to describe both sorption/desorption and ion-exchange reactions.

Based on studies of salt marsh restorations and microcosm experiments [4,5], immediate re-entry of seawater into diked and seasonally flooded and drained marshes results in significant die off of freshwater biomass and subsequent oxidation of organic matter coupled to the reduction of oxygen (Equation (1)), ferric oxyhydroxides (Equation (2)), and SO₄^2- (Equation (3)) and results in anoxic conditions and increased concentrations of alkalinity, Fe(II), and S^2- in sediments:

(1) CH₂O_(s) + O_2(g) → CO_2(aq) + H₂O_(l)

(2) CH₂O_(s) + 4Fe(OH)_3(s) + 7H⁺_(aq) → 4Fe²⁺_(aq) + HCO₃^-_(aq) + 10H₂O_(l)

(3) CH₃COOH_(s) + SO₄^2-_(aq) + 2H⁺_(aq) → H₂S_(aq) + 2CO_2(aq) + 2H₂O_(l).

Ferric hydroxides can also be reduced by dissolved S^2- [29]. These biogeochemical reactions in saltmarsh sediments are analogous to those described for marine sediments during diagenesis [30] and follow the well-known sequence of microbially mediated reactions that proceed from most energetically favorable to less favorable reactions. Organic matter is first oxidized by oxygen, then by nitrate, Fe(III), SO₄^2-, and finally carbon (C) through methanogenesis by the hydrogenotrophic pathway (Equation (4)),

(4) HCO₃^-_(aq) + H⁺_(aq) + 4H₂ → CH_4(aq) + 3H₂O_(l)

or by the acetoclastic pathway [31]. Fe(II) or Fe(III) can react with the S^2- to form Fe monosulfides (FeS), such as mackinawite (Equation (5)):

(5) 2FeOOH_(s) + 3HS^-_(aq) → 2FeS_mack + S⁰_(aq) + 3OH^-_(aq) + H₂O_(l)

and eventually Fe disulfide (FeS₂) such as pyrite, which can take place rapidly with H₂S as a reactant [32,33,34]. Large amounts of NH₄⁺ and PO₄^3- could be released by desorption and cation exchange as pH rises, and from reductive dissolution of Fe(III) hydroxide minerals, causing eutrophication; P is also released from sorption sites.

The aeration of peat marsh sediments can cause oxidation of pyrite (Equation (6)):

(6) FeS_2(s) + 7/2O_2(g) + H₂O_(l) → Fe²⁺_(aq) + 2SO₄^2-_(aq) + 2H⁺_(aq)

that generates highly acidic (pH <4) porewaters rich in dissolved Fe(II) and SO₄^2- and can lead to fish mortality [1,2,35]. The Fe(II) initially released may be oxidized to Fe(III), which is a more effective oxidant than oxygen [36,37,38], and can be solubilized by organic ligands [39]. The Fe²⁺ oxidizes and may precipitate Fe(III) hydroxide and Fe(III) hydroxysulfate minerals, such as goethite (Equation (7)) or jarosite minerals (Equations (8) and (9)):

(7) Fe(OH)²⁺_(aq) + H₂O_(l) ↔ FeOOH_(goethite) + 4H⁺_(aq)

(8) 3Fe³⁺_(aq) + M⁺_(aq) + 2HSO₄^-_(aq) + 6H₂O_(l) → MFe₃(SO₄)₂(OH)_6(jarosite) + 8H⁺_(aq)

(9) 3Fe³⁺_(aq) + NH₄⁺_(aq) + 2HSO₄^-_(aq) + 6H₂O_(l) → NH₄Fe₃(SO4)₂(OH)_{6(jarosite-amm)} + 8H⁺_(aq).

Jarosite is commonly found in acid- SO₄^2- soils as a byproduct of pyrite oxidation in DD sediments [5] and is well-documented to form at low pH [40]; M⁺ can include monovalent cations such as Na⁺ or NH₄⁺, as in ammonium jarosite or “jarosite-amm” (Equation (9)). Goethite is the predominant secondary phase on pyrite surfaces under alkaline conditions [41]. Cations such as Na⁺ can compete with Fe²⁺ for sorption sites on goethite and other Fe oxide minerals [42] such as Na⁺ for FeX₂ (Equation (10)):

(10) FeX_2(exch) + 2Na⁺_(aq) ↔ 2NaX_(exch) + Fe²⁺_(aq)

and can result in increased dissolved Fe²⁺ concentrations [43,44]. The biogeochemical cycling of Fe in estuaries has a strong seasonal component because of the integral role of organic ligands, which are generated and decomposed during the growing season, in the Fe cycle [39].

Fe sulfide mineral (FeS or FeS₂) oxidation (Equation (6)), Fe(II) oxidation, and jarosite precipitation (Equations (7)–(9)) generate acidity that can be harmful to aquatic life. The United States Environmental Protection Agency (USEPA) water-quality criteria for pH in freshwater suggest a range of 6.5 to 9; pH below this range physiologically stresses many species and can result in decreased reproduction, decreased growth, disease, or death [4,45]. NH₃ is toxic to saltwater aquatic life, with chronic and acute values of 0.23 mg/L (0.016 mmol/L) and 1.7 mg/L (0.12 mmol/L), respectively [46]. NH₄⁺ is the dominant species at pH <9.2 and is likely mobilized by desorption during salination [4,47]. Al is also mobilized during salination of diked marshes and is a particular concern for aquatic life, particularly fish [48].

The Herring River estuary in Wellfleet, Massachusetts (MA) (Figure 1) is an example of such a diked marsh system, and the third and largest salt marsh restoration project (445 hectares) is currently underway. The ocean tides were blocked from the Herring River following construction of the Chequesset Neck Road Dike in 1909, causing a transition from salt marsh to freshwater wetlands. Partial restriction on inflowing tides has resulted in lower mean sea level landward of the dike so that surficial sediments in wetland areas adjacent to the Herring River main stem have been dewatered (Figure 2). In 2017-18, the typical tide range on the upstream side of the Herring River tide-control structure was about 0.67 m, whereas on the downstream side (Wellfleet Harbor), it was about 3.1 m [49]. Restoration of the Herring River marsh involved demolition of the existing dike and removal of tidal restrictions in a phased approach by outfitting the new bridge with sluice gates that will allow water and salinity levels to rise gradually over the coming decade, enabling the ecosystem to grow back more naturally. As of March 2025, tidal flow has begun to be transitioned from the old dike culverts to the tidal gates of the new Chequessett Neck Road Bridge and Water Control Structure, but the restricted tidal flow remains the same and the restoration process will not begin until construction is complete in another year [50]. Considerable concern over the marsh restoration of the Herring River exists because (1) a stressed river herring population that includes alewife (Alosa pseudoharengus) and blueback herring (Alosa aestivalis) [51] migrates through the river into upstream kettle ponds; (2) acidic waters with low dissolved oxygen (DO) and high Fe and nutrient concentrations, common to diked marsh creeks, are released during seawater interaction and intense rain events [3,52]; and (3) these waters are received by Wellfleet Harbor, which contains productive shellfish aquaculture that is sensitive to water-quality degradation.

Reflooding with seawater during the Herring River restoration would likely result in rapid infiltration into unsaturated DD sediments, as seawater infiltrates pore spaces and replaces soil gas, but relatively slowly in the DF sediment, where freshwater will eventually be replaced by saltwater. Out of concern that this marsh restoration program could severely affect the Wellfleet Harbor Area of Critical Environmental Concern [53], sediment cores were collected from areas near and in the area of anticipated restoration and used in marsh sediment microcosm salination experiments to assess the biogeochemical effects of seawater intrusion, as described in Portnoy and Giblin [4],Portnoy and Giblin [5].

In this study, we used biogeochemical modeling to describe seawater reintroduction to diked marsh sediments based on microcosm experiments using sediment cores from Herring River marsh and the North Sunken Meadow salt marsh in Eastham, MA (Figure 1 and Figure 2) [4,5]. The reactive-transport model simulations help relate the biogeochemical changes observed in the microcosms after salination [4] to chemical reactions and reaction rates so that the timing and sequence of these reactions (including desorption, and a redox shift to SO₄^2- reduction) can be simulated over longer time periods for the resalination of the Herring River, and potentially other diked marsh settings.

The USGS computer program PHREEQC (“pH-redox-equilibrium model in the C programming language”; Parkhurst and Appelo, 2013) was used to simulate one-dimensional reactive transport in diked salt marsh sediments. The models account for spatial and temporal changes in seawater flowing through freshwater sediments, transport of associated chemical constituents, and a set of chemical reactions including desorption of constituents, organic matter decomposition, and dissolution and precipitation of iron-bearing minerals. The online supplementary information (SI) for this paper (include link here) provides further details on the PHREEQC [28] model input and output files, along with aquatic life water-quality criteria for simulated pH and solute concentrations. The PHREEQC model input files, and instructions on running the model, are also available as a USGS data release [54].

Figure 1. Location of the Herring River study area, Wellfleet, Massachusetts, and location of the North Sunken Meadow core sample collection site. Modified from [49].

Figure 1. Location of the Herring River study area, Wellfleet, Massachusetts, and location of the North Sunken Meadow core sample collection site. Modified from [49].

Figure 2. Herring River watershed basin boundary with the major soil series as determined by the Natural Resources Conservation Service [55]. The Freetown and Swansea soil series (green) are primarily diked, flooded (DF) soils that occur primarily on upgradient arms of the marsh system, where groundwater levels are presumably higher in elevation. The Maybid soil series (red) are primarily diked, drained sediments with acidic pH. Modified from [56].).

Figure 2. Herring River watershed basin boundary with the major soil series as determined by the Natural Resources Conservation Service [55]. The Freetown and Swansea soil series (green) are primarily diked, flooded (DF) soils that occur primarily on upgradient arms of the marsh system, where groundwater levels are presumably higher in elevation. The Maybid soil series (red) are primarily diked, drained sediments with acidic pH. Modified from [56].).

2. Methods

Geochemical models were developed and calibrated with PHREEQC [28] using data from [4] microcosm experiments (Appendix A.1), then “basecase” models were used to simulate a single flow velocity over longer time periods (Table 1) and with varying Fe(III) source minerals. PHREEQC is ideal for biogeochemical simulations such as tidal marsh restoration because it implements several types of aqueous models that have capabilities for calculations including speciation and saturation index, and 1D transport with reactions including aqueous, mineral, gas, solid-solution, surface-complexation (used to model sorption and desorption) and ion-exchange equilibria, and kinetic reaction [28]. It also has a built-in BASIC interpreter that allows for complex reaction rate laws and output post-processing options [28]. PHREEQC was used with the phreeqc.dat thermodynamic database; input files, concentration-time plots, and selected output are provided in Supplementary Information (SI) Tables S1–S10 and in Brown [54].

2.1. Geochemical Models of Microcosm Experiments

This section discusses how previous microcosm experiment data [4] were evaluated and used to develop and calibrate one-dimensional models using PHREEQC [28] to simulate the resalination of DF and DD sediments at the Herring River after dike removal. The Portnoy and Giblin [4] microcosm experiments are discussed in Appendix A and summarized in Table A1. The DF and DD marsh areas at Herring River are shown in Figure 2. The DF marsh areas that have remained saturated with water consist mostly of Freetown and Swansea series soils, which are deep, poorly drained peat soils formed in decomposed organic material in freshwater environments [4,5,57]. Air exposure of the DD sediments has resulted in extensive oxidation of pyrite and acid SO₄^2- soil conditions that generally correspond to the areas mapped as Maybid soil series (Figure 2) [57]. Data from [4] were evaluated to determine the simulated rates of organic matter oxidation coupled to the reductive dissolution of Fe(III), and SO₄^2- reduction, and also curve matching was used to calibrate pH and solute concentrations of NH₄⁺ and PO₄^3-.

Two microcosm experiment models were created, one for salination of each sediment type (DF and DD) used in the experiments. The models assume that fresh porewater initially fills the sediment cores (Figure 3) and is replaced by an influx of seawater for 20 months, to be consistent with Portnoy and Giblin [4] experiments. In simulations of DF and DD sediments, seawater passed through the columns after about 8 hrs and 5 hrs, respectively. The models had similar initial solute concentrations of freshwater, but the pH, Fe(II), and ammonium concentrations differed. Seawater concentrations (Table 1) were the same for both microcosm models.

2.1.1. Sediment Cores

The microcosm core dimensions were incorporated into a model grid that uses 45 1-cm cells for the 45-cm core length (Figure 3), with a core volume of 7.95 liters (L). During the Portnoy and Giblin [4] experiments, the DF and DD sediment cores were initially flushed with 10 L of water to exchange the freshwater for saltwater over a 12-hour period. Subsequent flow rates from saltwater additions to make up for sampling water and evaporation were 1 to 2 orders of magnitude slower. The same amounts of saltwater were added in both the DD and DF cores; because of different porosities in DD and DF, flow rates varied between the cores (Table 1).

For the two PHREEQC models (Tables S1a, S1b, S2a and S2b), the fresh porewater initially in the columns is defined as solutions 1-45 and equilibrated with the DF and DD sediments, respectively. Establishing the initial chemistry involved assigning surface exchange sites, solid phase reactant concentrations, kinetics, and solution concentrations for each cell. The infilling seawater solution (“Solution 0”) was then introduced to the column, and transport was modeled by “shifting” solution 0 to cell 1, the solution in cell 1 to cell 2, and so on until solution 45 is shifted to cell 45, which is equilibrated with exchange 45, surface 45, equilibrium phases 45, and kinetics 45. The keywords, including “exchange”, “surface”, “equilibrium phases”, and “kinetics”, define the input for PHREEQC. “Exchange” defines the exchange assemblage composition, “surface” defines the composition of surfaces, “equilibrium phases” defines a singular phase that reacts reversibly, and “kinetics” specifies kinetic reactions and parameters for reactive-transport calculations. The number of cell shifts and the time per shift used in the simulation are shown in Table 1. At every time step, the speciation of the solution and its reaction with sediments was computed.

2.1.2. Organic Matter Reaction Rates

In the model simulations, seawater (Solution 0) is passed through the 45 sediment core cells, bringing O₂ and SO₄^2- electron acceptors into contact with the Fe(III)- and organic matter-bearing sediments, which oxidizes organic matter, Fe(II), and H₂S, and displaces the fresh porewater. Reactions defined in the conceptual model required the calculation and (or) matching of the solute curves from the Portnoy and Giblin [4] to the simulated curves based on (1) the amount of Fe that is reactive for reductive dissolution, (2) fitting the rates of organic carbon reaction coupled to either DO reduction, Fe(III) reductive dissolution, SO₄^2- reduction, or methanogenesis (Table 3 and Table 4), and (3) surface reactions between ions and sorption sites (Table 4). Solid-phase chemistry for each model included organic matter, FeS, jarosite, Al(OH)_3(a), clays with cation-exchange capacity, and hydrous-ferric-oxide surfaces. Organic matter, which consists of roots and litter from modern terrestrial vegetation, is most abundant in surface soils [4]. Both jarosite (KFe₃(SO4)₂(OH)₆) and goethite (FeOOH) commonly occur in root channels above the continually waterlogged peat containing intact salt marsh rhizomes [58]. Model simulations included decomposition of organic matter assuming a C:N:P molar ratio of 106:16:1 (Table 1) [59,60]. Coastal salt marshes typically contain a higher C ratio because of contributions from vascular plants, but simulations using higher percentages of C and a higher N:P ratio [61,62] than the Redfield ratio did not improve the model fit.

Table 2. Rates of organic matter decay, and mineral dissolution and precipitation reactions used in the resalination models of diked, flooded (DF) sediments, and diked, drained (DD) sediments. Note that oxygen reduction and denitrification are not significant to oxidation of carbon; PHREEQC simulations show that dissolved oxygen and any existing nitrate are reduced very quickly after sediments are filled with seawater.

Reaction number	Reaction	Reaction type	Equation	Rates (1/s)-determined by calibration for DF sediments	Rates (1/s)-determined by calibration for DD sediments
1	Oxidation of organic matter by reduction of Fe(III) [ammonium jarosite is source of the Fe(III) used in the model]	Redox, kinetically controlled	DOM1 + 4x Fe3+4xOH- = xCO2(g) +4x Fe2+ + 3/xH2O + y NH3 + z H3PO4	Rate = 3.5e-9	Rate = 7e-9
2	Oxidation of organic matter by reduction of SO42-	Redox, kinetically controlled	DOM1 + (x/2)SO4-2 + (y-2z)CO2 + (y-2z)H2O --------> (x/2)H2S + (x+y-2z)HCO3- + yNH4+ + zHPO42-	IF(tot(“Fe(+3)”) <= 2e-6, THEN rate = 2.5 e-10	IF(tot(“Fe(+3)”) <= 2e-6, THEN rate = 1.76e-9
3	Oxidation of organic matter by methanogenesis	Redox, kinetically controlled	DOM1 + (y-2z)H2O -------> x/2CH4 + (x-2y+4z/2)CO2 + (y-2z)HCO3- + yNH4+ + zHPO42-	IF (tot(“Fe(+3)”) <= 2e-6, AND (tot(“SO4-2)”) <= 2e-6, THEN rate = 3 e-11	IF (tot(“Fe(+3)”) <= 2e-6, AND (tot(“SO4-2)”) <= 2e-6, THEN rate = 1.76e-10
4	FeS precipitation	Equilibrium	FeS = Fe+2 + S-2	NA	NA
5	Al(OH)3 (amorphous) precipitation and dissolution	Equilibrium	Al(OH)3 + 3H+ = Al+3 + 3 H2O	NA	NA

^(a)DOM = (CH₂O)x(NH₃)y(H₃PO₄)z, where x, y, z represents the C:N:P ratios [65].

Table 2. Rates of organic matter decay, and mineral dissolution and precipitation reactions used in the resalination models of diked, flooded (DF) sediments, and diked, drained (DD) sediments. Note that oxygen reduction and denitrification are not significant to oxidation of carbon; PHREEQC simulations show that dissolved oxygen and any existing nitrate are reduced very quickly after sediments are filled with seawater.

Reaction number	Reaction	Reaction type	Equation	Rates (1/s)-determined by calibration for DF sediments	Rates (1/s)-determined by calibration for DD sediments
1	Oxidation of organic matter by reduction of Fe(III) [ammonium jarosite is source of the Fe(III) used in the model]	Redox, kinetically controlled	DOM1 + 4x Fe3+4xOH- = xCO2(g) +4x Fe2+ + 3/xH2O + y NH3 + z H3PO4	Rate = 3.5e-9	Rate = 7e-9
2	Oxidation of organic matter by reduction of SO42-	Redox, kinetically controlled	DOM1 + (x/2)SO4-2 + (y-2z)CO2 + (y-2z)H2O --------> (x/2)H2S + (x+y-2z)HCO3- + yNH4+ + zHPO42-	IF(tot(“Fe(+3)”) <= 2e-6, THEN rate = 2.5 e-10	IF(tot(“Fe(+3)”) <= 2e-6, THEN rate = 1.76e-9
3	Oxidation of organic matter by methanogenesis	Redox, kinetically controlled	DOM1 + (y-2z)H2O -------> x/2CH4 + (x-2y+4z/2)CO2 + (y-2z)HCO3- + yNH4+ + zHPO42-	IF (tot(“Fe(+3)”) <= 2e-6, AND (tot(“SO4-2)”) <= 2e-6, THEN rate = 3 e-11	IF (tot(“Fe(+3)”) <= 2e-6, AND (tot(“SO4-2)”) <= 2e-6, THEN rate = 1.76e-10
4	FeS precipitation	Equilibrium	FeS = Fe+2 + S-2	NA	NA
5	Al(OH)3 (amorphous) precipitation and dissolution	Equilibrium	Al(OH)3 + 3H+ = Al+3 + 3 H2O	NA	NA

^(a)DOM = (CH₂O)x(NH₃)y(H₃PO₄)z, where x, y, z represents the C:N:P ratios [65].

Table 3. Solution chemistry of diked, flooded (DF) and diked, drained (DD) sediments and artificial seawater used in PHREEQC simulations of Portnoy and Giblin (1997a) microcosm experiments. Freshwater chemistry from Portnoy and Giblin except as noted, using freshwater and column-bound concentrations from the same types of cores and at the same depth. Concentrations in mmol/L.

Constituent or property	DF solution 1-45	DD solution 1-45	Seawater composition, solution 0
Temperature °C	25	25	25
pH	6.7	4	8.5b
pe	--	--	8.45
Na	0.026a	0.026a	468d
Ca	0.004a	0.004a	10.2d
Mg	0.0015a	0.0015a	53.2d
K	0.001a	--	10.2d
Cl	0.040	0.715b	545d
S(6), sulfate	--	--	28.2d
Alkalinity as HCO3-	4b	0.1b	2.3d
N(-3), ammonium	0.001	0.075b	-
S2-	0.1b	0b	--
P	0.001c	0.002b	--
Fe(II)	0.0001b	0.1b	--
Fe(III)	0.0001b	0.001	--
Si	--	--	0.07d
O(0), diss. oxygen	0.01	0.01	0.75
Al	0.01b	0.3b	--

^a Generally based on Portnoy and Giblin [4],Portnoy and Giblin [5]. ^c fit parameter-- for phosphorus DF and aluminum DF. ^d Nordstrom, et al. [63].

Table 3. Solution chemistry of diked, flooded (DF) and diked, drained (DD) sediments and artificial seawater used in PHREEQC simulations of Portnoy and Giblin (1997a) microcosm experiments. Freshwater chemistry from Portnoy and Giblin except as noted, using freshwater and column-bound concentrations from the same types of cores and at the same depth. Concentrations in mmol/L.

Constituent or property	DF solution 1-45	DD solution 1-45	Seawater composition, solution 0
Temperature °C	25	25	25
pH	6.7	4	8.5b
pe	--	--	8.45
Na	0.026a	0.026a	468d
Ca	0.004a	0.004a	10.2d
Mg	0.0015a	0.0015a	53.2d
K	0.001a	--	10.2d
Cl	0.040	0.715b	545d
S(6), sulfate	--	--	28.2d
Alkalinity as HCO3-	4b	0.1b	2.3d
N(-3), ammonium	0.001	0.075b	-
S2-	0.1b	0b	--
P	0.001c	0.002b	--
Fe(II)	0.0001b	0.1b	--
Fe(III)	0.0001b	0.001	--
Si	--	--	0.07d
O(0), diss. oxygen	0.01	0.01	0.75
Al	0.01b	0.3b	--

^a Generally based on Portnoy and Giblin [4],Portnoy and Giblin [5]. ^c fit parameter-- for phosphorus DF and aluminum DF. ^d Nordstrom, et al. [63].

Table 4. Ion-exchange and surface complexation sites for diked, flooded (DF) and diked, drained (DD) sediments used in PHREEQC simulations of Portnoy and Giblin [4] microcosm experiments.

Parameter	DF biogeochemistry	DD biogeochemistry
Exchange concentrations (mol/L)	0.50	0.85
Surfaces sites (mol/L)	0.027	0.2
Equilibrium phases-- ammonium jarosite concentration (mol/L)	0.0001 throughout (0.005 in basecase)	Varies with profile (Figure 4)
Rate of reaction	DF rates (Table 2)	DD rates (Table 2)
Solutions	DF freshwater; artificial seawater (Table 3)	DD freshwater; artificial seawater (Table 3)

Table 4. Ion-exchange and surface complexation sites for diked, flooded (DF) and diked, drained (DD) sediments used in PHREEQC simulations of Portnoy and Giblin [4] microcosm experiments.

Parameter	DF biogeochemistry	DD biogeochemistry
Exchange concentrations (mol/L)	0.50	0.85
Surfaces sites (mol/L)	0.027	0.2
Equilibrium phases-- ammonium jarosite concentration (mol/L)	0.0001 throughout (0.005 in basecase)	Varies with profile (Figure 4)
Rate of reaction	DF rates (Table 2)	DD rates (Table 2)
Solutions	DF freshwater; artificial seawater (Table 3)	DD freshwater; artificial seawater (Table 3)

Figure 4. Total iron concentration in DF (diked, flooded), DD (diked, drained), and NA (natural) marsh cores [5], with reactive Fe based on a fit parameter. The water table fluctuated seasonally at DF and DD sites; the DF site was waterlogged to the surface from December through March. At the DD site, groundwater was always below 15 cm and regularly dropped below 40 cm during the summer [5]. The NA marsh was constantly flooded during the winter, and waterlogged to the surface for about half of the spring-neap cycle during the growing season [5]. Microcosm experiments were not performed on the NA marsh cores.

Figure 4. Total iron concentration in DF (diked, flooded), DD (diked, drained), and NA (natural) marsh cores [5], with reactive Fe based on a fit parameter. The water table fluctuated seasonally at DF and DD sites; the DF site was waterlogged to the surface from December through March. At the DD site, groundwater was always below 15 cm and regularly dropped below 40 cm during the summer [5]. The NA marsh was constantly flooded during the winter, and waterlogged to the surface for about half of the spring-neap cycle during the growing season [5]. Microcosm experiments were not performed on the NA marsh cores.

The organic matter decay rates were varied iteratively to match the [4] microcosm experiment solute curves for alkalinity, Fe, and SO₄^2- until the occupied sites simulated matched as closely as possible. Rates of organic degradation were estimated by fitting the simulated peak of Fe(II) resulting from Fe(III) dissolution with the Portnoy and Giblin [4] observed peak and assuming that the Fe(II) peak is reached when the supply of reactive Fe(III) mineral (jarosite-amm) is exhausted. As soon as the peak dissolved Fe(II) is reached, the rate decreases because continued SOM degradation switches to SO₄^-2 reduction, which is thermodynamically less favorable than Fe(III) reduction; the S^2- generated reacts with Fe(II) and precipitates as FeS(ppt). The rate of SO₄^-2 reduction was adjusted so that the decrease in dissolved Fe(II) generally matched the measured decrease observed in the Portnoy and Giblin [4] experiments. The rate of methanogenesis could not be similarly fitted because SO₄^-2 was not exhausted in the microcosms, so the rate of methanogenesis was set to a factor of about 10 less than the rate of degradation by SO₄^-2 reduction.

2.1.3. Solution Chemistry, Ion Exchange, and Surface Complexation (Sorption)

Solutions used in the models include porewater solutions, which differ between the DF and DD sediment core models, and seawater solution, which was identical for both models (Table 3). The porewater solution chemistry used for DF and DD sediment simulations is based on values provided in Portnoy and Giblin (1997a) that represent constituent concentrations or properties measured in DF and DD cores; however, values were modified to include fit parameters for P and Al. The seawater solution was based on artificial seawater composition used in Portnoy and Giblin [4] experiments and mean seawater composition [63] (Table 3). PHREEQC version 3 [28] was used to calculate redox conditions based on the DO concentration that was assigned to the solution and (or) that varied during the simulations. PHREEQC assumes a redox condition, which is expressed in terms of pe (negative log of the electron activity) of 4, unless redox couples such as “H₂O/O₂” are specified. The pe estimates for the “H₂O/O₂” redox pair were computed using both the standard equilibrium expression (2H₂O = O₂ + 4H⁺ + 4e⁻; log K −86.08) and the relation of Sato [64] (log K −45.54), which indicates a lower equilibrium potential at a given pH than the standard log K and is generally consistent with that for the reaction H₂O₂ = O₂ + 2H⁺ + 2e^-.

Interactions between solutes and solid surfaces were modeled following Parkhurst and Appelo [28] (Example 14) using ion exchange and surface complexation, which can be collectively referred to as sorption or desorption. The cation exchange equilibrium reaction that is common to seawater intrusion of freshwater sediments is shown by Equation (11):

(11) Na⁺_(aq) + X^-_(exch) ↔ NaX_(exch),

where X- is a cation exchange site (e.g., filled by H+, K+, or 1/2Ca2+) common in clays and Na⁺ is an abundant cation in seawater.

Cation exchange parameters were determined from data on cation extraction from freshwater DD cores [4]. Calculating the number of exchange sites for cations was an iterative process that involved adjusting the number of sites based on the difference between measured and simulated concentrations on the surface in equilibrium with fresh porewater concentrations (Table 4), using literature values of cation exchange equilibrium constants from the PHREEQC database. All cations that may compete at exchange sites were included, and the process began with 1 mole of sites per liter pore water. At each iteration, the simulation-equilibrated amount of Fe(II) and NH₄⁺ on the solids was compared to measured amounts leached from the solids by KCl extractions [4], and the number of sites adjusted for the next iteration was decreased if simulated bound concentrations were greater than measured or increased if less than measured. This iteration process continued until the occupied sites simulated matched as closely as possible the measured amount of exchangeable Fe(II) and NH₄⁺.

Phosphorus sorption was modeled as a surface-complexation process ([26,27,28] that proceeds by the reaction:

(12) SiteOH + PO₄^3-_(aq) + 2H⁺_(aq) ↔ SiteHPO₄^- + Fe²⁺_(aq) + H₂O_(aq)

The default database for PHREEQC contains thermodynamic data for hydrous ferric oxide (Hfo) on a particle surface; the values are derived from Dzombak and Morel [27], who defined a strong binding site (Hfo_s) and a weak binding site (Hfo_w) and used 0.2 M per mol Fe for weak sites and 0.005 M for strong sites, a surface area of 5.33x10⁴ m²/M Fe, and a gram-formula weight of 89 g Hfo per mol Fe. The surface complexation model of Parkhurst, et al. [26] was used to simulate interactions for surface complexation of H⁺ and PO₄^3-, two constituents that react at surface sites and are of interest for eutrophication and water quality. Following Parkhurst, et al. [26], the number of sites for H⁺ sorption for DD sediments was varied until the simulated pH matched the measured pH at the end of the [4] experiment; the number of surface sites for DF were adjusted for the difference in porosity.

2.2. Basecase Geochemical Models

Two basecase models also were developed using the same reaction rates but a constant porewater velocity was used for each of the DF and DD sediment cores, unlike in the microcosm models (Table 1). The PHREEQC simulations provide estimates of speciation of nutrients, metals, and other constituents over the same period as the Portnoy and Giblin [4] experiments and for the basecase models over 12-yr and 100-yr periods (Tables S2–S4 and S10a,b). Whereas the microcosm models required several velocities to simulate the Portnoy and Giblin [4] experiments, constant velocities similar to those used in the microcosm experiments were used for the basecase model (1.21x10^-8 m/s for DF sediments and 1.97 x10^-8 m/s for DD sediments). These velocities are on the low end of the range observed in tidal wetlands by Harvey, et al. [66] of 1.2 x10^-5 to 2 x10^-7 m/s, after converting from hydraulic conductivity using porosities of 90 and 55 percent for DF and DD cores, respectively (Table 1) [4], and vertical hydraulic gradients of 1. Sensitivity analyses were performed by running scenarios under varying types and strength of Fe(III) source minerals.

3. Results and Discussion

3.1. Geochemical Modeling of Microcosm Experiments

The models for resalination of both DF and DD marsh sediments assume that fresh porewater initially fills the sediment cores (Figure 3) and is replaced by an influx of seawater for about 20 months. As part of the model calibration, curve fitting was used to develop both the DF and DD models, as described in Section 2.1.

3.1.1. Calibration of Microcosm Models

The organic matter decay rates were varied iteratively to match the [4] microcosm experiment solute curves for alkalinity, Fe, and SO₄^2- (Table 2). A rate of organic matter decay of 3.5 x 10^-9 mol L^-1s^-1 was used under Fe(III)-reducing conditions but was estimated to slow to a 0.25 x 10^-9 mol L^-1s^-1 if Fe(III) was not present (assumed SO₄^-2 reducing conditions) (Table 2). The organic matter decay rates for Fe(III) reduction, SO₄^2- reduction, and rates of methanogenesis used for the DD cores were greater (7 x 10^-9, 1.76x10^-9, and 1.76 x10^-10 mol L^-1s^-1, respectively) than those used for DF cores and result from less reducing (methanogenic) conditions in DD sediments. Furthermore, because of the redox gradient in DD sediments beneath the subaerial exposed surface, 9 layers of different equilibrium phases were necessary in the DD model to best fit the Fe(II) solute curve of Portnoy and Giblin [4] (Table S1b). Organic matter (OM) decay rates from Table 2 were normalized to the dry sediment weight for DF and DD sediment column samples (0.7952 and 3.181 grams as carbon, respectively) for comparison with the organic carbon (C) mineralization rates of other studies. Sediment-weighted organic C mineralization rates were 3.8x10^-9 mol kg^-1s^-1 for DD sediments and 8.9x10^-9 mol kg^-1s^-1 for DD sediments. These rates are similar to those reported by Neubauer, et al. [67] who measured rates of anaerobic C metabolism for brackish marshes in Chesapeake Bay, Maryland (USA) of 8.1x10^-9 to 4.9x10^-8 mol kg^-1s^-1 for the months of June and July. The slightly lower rates estimated in this study could result in part from the fact that these are fitted rate estimates for the Portnoy and Giblin [4] entire microcosm column experiments, compared to those from [67], which were determined during the growing season.

As discussed in Section 2.1.3, calculating the number of exchange sites for cations was also an iterative process. Simulated values of 0.5 and 0.85 moles of exchange sites per liter of pore water were used for DF and DD sediment cores, respectively (Table 4) as shown in the PHREEQC input files (Tables S1a and S1b). Anion sorption is important for PO₄^3-, which could be important to eutrophication in Wellfleet Harbor. The number of sites and the log k for H⁺ were adjusted until the simulated pH best fit the Portnoy and Giblin [4] data (Table 4; Tables S1a,b).

Sorption (or desorption) of PO₄^3- was modeled using the surface complexation reaction (Equation (12)) as described in Section 2.1.3. The number of sites for H⁺ sorption in DD sediments was varied until simulated and measured pH matched at the end of the experiment at 0.2 mol/L; the log k for H⁺ was adjusted to 7.5 until the beginning pH best fit the Portnoy and Giblin [4] data.

Surface sorption sites were adjusted for DF according to difference in porosity (Equation (13)),

(13) AnionSites_DF = 0.2[0.55/(1-0.55)][(1-0.9)/0.9]= 0.027.

The immediate dissolved Fe increase observed in the microcosm DF and DD cores after seawater flushing indicates that Fe(II) was mobilized by desorption and reductive dissolution of Fe(III) minerals [4]. Exchangeable Fe²⁺ (FeX₂) would have been largely removed by exchange with Na⁺ or desorption during salination so that reductive dissolution must be largely responsible for the post-flushing Fe build up.

Fe(II) concentrations measured in marsh cores were roughly constant with depth in DF and natural (NA) cores but increasing with depth in DD cores [5] (Figure 4). For the DF simulation, a reactive Fe (jarosite-amm) value of 0.1 mmol/L was used as a fit parameter such that the simulated peak in dissolved Fe(II) provided a reasonable match with the measured value and was used for a jarosite-amm equilibrium phase in the PHREEQC input file (Table S1a). In the DD sediment simulation, total Fe and reactive Fe (upper x-axis) distribution was represented by a depth profile that increased towards the air/water interface because of ferrous iron migration and precipitation (Figure 4) [5]. The ratio between total Fe and reactive Fe was established as a fit parameter in the DD calibration but was distributed in proportion to the profile of total Fe, and a range of 0.015 to 0.038 mol/L for a jarosite-amm equilibrium phase was used in the PHREEQC input file (Table S1b).

Na jarosite was used as the Fe(III) oxide phase in initial simulation trials, but was changed to jarosite-amm (eq 9) because of an otherwise unaccounted for increase in NH₄⁺ concentration observed by Portnoy and Giblin [4] that generally correlated with Fe dissolution. As with Fe(II), much NH₄⁺ would be exchanged from sediment during the initial salination, followed by an increase attributed to reductive dissolution. No reasonable ratio of N:P such as the Redfield ratio (16:1) for organic matter in coastal wetlands [60,62] would account for the buildup of NH₄⁺. When jarosite-amm was substituted for Na jarosite, the simulated NH₄⁺ increase matched the measured NH₄⁺ increase. Additional support for the NH₄⁺ form of jarosite is the concentration of NH₄⁺ (0.080 mmol/L) in the freshwater DD sediments as compared to Na at 0.058 mmol/L (average in precipitation, Truro, MA) and the greater chemical stability of jarosite-amm as compared to Na jarosite [68].

The simulated DF and DD core output chemistry plotted over time (Figure 5 and Figure 6) generally corresponds with those measured during the Portnoy and Giblin [4] microcosm experiments, as discussed below.

3.1.2. Biogeochemical Model of Diked, Flooded (DF) Sediments

Because of the accumulation of organic matter and acidity over decades of diking, seawater flooding could potentially release a large amount of Fe(II), as well as nutrients such as of NH₃ (including NH₄⁺) and PO₄^3- [4,26] by desorption and (or) mineral dissolution. Portnoy and Giblin [4] attributed a 40- to 60-fold NH₄⁺-N increase in DD porewater to desorption during resalination experiments. Aerated DD soils were found to retain high N and P concentrations, with NH₄⁺ adsorbed on silts and clays and PO₄^3- likely associated with abundant Fe and Al oxides [2].

PHREEQC simulations of resalination show that metals and nutrients are mobilized from DF sediments by desorption. Nutrients and dissolved inorganic carbon (DIC) are also mobilized by the decay of organic matter. Whereas healthy tidal marshes are characterized by SO₄^2- reducing redox conditions [8], DF sediments have low amounts of Fe(III) and SO₄^2-, and so the predominant redox process is methanogenesis. The geochemical model of DF sediments, therefore, includes organic matter, equilibrium phases including FeS, a small amount (0.1 mmol/L) of Fe(III) source mineral (i.e., jarosite-amm), and Al(OH)_3(a), and exchangeable Fe(II), NH₄⁺, and PO₄^3-(SI Table S1a). The response of pH and solute concentrations is shown in Figure 5, and the PHREEQC output file and more detailed (shorter time interval) plots are provided in SI Table S3.

The PHREEQC model calibrated to the Portnoy and Giblin [4] experiment data shows that seawater passes through the DF sediment column after about 8 hrs. The DO initially present in seawater is immediately reduced after reacting with DF sediments that contain abundant organic matter, Fe(II), and S^2-, and Fe(II) increases sharply in the first 9 hrs. The increasing front of dissolved solids concentration is accompanied by a brief decrease in pH as hydrogen ions (H⁺), Fe(II), and other cations are desorbed from sediments, and a small amount of jarosite-amm is dissolved (Figure 5a–d). The seawater is oxic, as evidenced by the pe (Sato) peak in column effluent from 8 to 14 days, before decreasing again as dissolved Fe(II) begins to increase. The desorption and exchange of Na⁺ for exchangeable cations including FeX₂ (Equation (8)) and HX results in an initial Fe(II) concentration peak at 0.4 days (Figure 5c). Dissolved Fe(II) in porewater increases again and peaks at 3.4 mmol/L at 5 days (Figure 5c) in response to the reductive dissolution of jarosite-amm (Equation (2); Figure 5c,d). Once jarosite-amm is depleted after 10 days, SO₄^2- reduction begins and S^2- is generated after 14 days (Figures 5d, 5f; SI Table S3). The pH and alkalinity increase with the introduction of alkaline, SO₄^2--rich seawater and the generation of alkalinity from SOM oxidation coupled to Fe(III) mineral reduction and SO₄^2- reduction (Equations (2) and (3); Figure 5d,e). As the concentration of SO₄^2- is reduced, S^2- increases, reacts with dissolved Fe(II), and begins to form FeS at 33 days (Figure 5d–f). FeS mineral precipitation peaks after 118 days, then decreases with the concentration of dissolved Fe(II) and begins to dissolve after 149 days (Figure 5d,e).

Both alkalinity and NH₄⁺ increases steadily throughout the DF sediment column simulation, along with a steady decrease in SO₄^2- consistent with SO₄^2- reduction coupled to SOM oxidation (Figure 5e,g). Much of the total NH₃ increase (includes both the NH₄⁺ + NH₃ forms) observed in the DF column experiments seems to result from the weathering of jarosite-amm and degradation of SOM rather than strictly desorption as surmised by Portnoy and Giblin [4]. NH₄⁺ is the dominant and least toxic form of total NH₃ in the moderately acidic DF porewater throughout most of the simulation. NH₃ is the most toxic form to aquatic life and, at elevated pH, dominates over the ionized form, NH₄⁺. Although P is commonly considered the “limiting nutrient” in freshwater ecosystems, it is generally not limiting in estuarine environments unless N inputs are sufficient, in which case PO₄^3- inputs can cause eutrophication and harmful algal growth [69]. Phosphorus concentrations (nearly all in the form of PO₄^3-) were less than 0.002 mmol/L as P at the beginning of the simulation and steadily increased after about 80 days, reaching 0.003 mmol/L after 582 days (Figure 5g).

Simulated Al concentrations increases sharply as a result of desorption, to 15 mmol/L during the onset of salination, and as pH decreases at 9 hrs, some Al(OH)_3(a) precipitates (1.1 mmol/L) (Figure 5a,h). By 12 hrs, Al levels off at 0.62 mmol/L and gradually declines below the original concentration by about 180 days (SI Table S3).

3.1.3. Biogeochemical Model of Diked, Drained (DD) Sediments

Whereas DF sediments have low amounts of Fe(III) and SO₄^2-, DD sediments have abundant Fe(III) minerals from the oxidation of FeS minerals, and so Fe(III)-reducing conditions prevail below the upper sediments that are subaerially exposed. The geochemical model of DD sediments includes organic matter, the equilibrium phases jarosite-amm, FeS, and Al(OH)_3(a), and exchangeable Fe(II), NH₄⁺, and PO₄^3-(SI Table S1b). The presence of Fe(III) and SO₄^2- in these sediments (products of FeS oxidation) provides electron acceptors for Fe(III) reduction and SO₄^2- reduction coupled to the oxidation of organic matter. The responses of pH and solute concentrations to salination of DD sediments are shown in Figure 6, and the PHREEQC output file and more detailed (shorter time interval) plots are provided in SI Table S4.

The simulated pH is much lower at the beginning of the simulation than the DF sediment because of the acidity produced by Fe sulfide oxidation and desorption of cations including H⁺ and Fe(II) (Figure 6a–c). Jarosite-amm was included as an equilibrium phase that varied with depth (SI Table S1b) to best match the Fe(II) concentration observed by Portnoy and Giblin [4]. After about 77 days of salination, the pH rises above 5 and alkalinity continues to increase, both from the high alkalinity in seawater and from the oxidation of SOM coupled to SO₄^2- and Fe(III) reduction (Figure 6a–e; SI Table S4). At the beginning of the simulation, DO decreases to low concentrations after reacting with DD sediments that contain organic matter, Fe(II), and S^2- (SI Table S4.) The greater organic matter decay rates used for DD sediments result from the energetically more favorable reductive dissolution of jarosite-amm compared to the Fe(III)-poor DF sediments (Table 2; Figures 5c and 6c); this redox shift is supported by the higher alkalinity concentration (Figure 6e). The jarosite-amm is not depleted until 126 days (as compared to 12 days until depletion for the DF core), at which point the pe(Sato) drops to -12.6, SO₄^2- reduction becomes the dominant electron-accepting process, and S^2- begins reacting with dissolved Fe(II) to form FeS minerals (Figure 6c,d). The SO₄^2- concentrations are higher in DD sediments than in DF sediments due to the SO₄^2- remaining from FeS oxidation, and S^2- concentrations are lower (Figure 6g), due to the greater amount of FeS mineral being formed (Figure 6d). The higher alkalinity concentrations in DD porewater are consistent with organic matter mineralization rates from the microcosm model fittings, which are greater for the DD sediments compared to the more reducing DF sediments, as discussed in Section 2.1.2 (Table 2).

Total NH₃ concentrations are much higher in DD sediments than in DF sediments, which demonstrates the greater decay rates of SOM, the greater abundance of jarosite-amm, and the acidic soils that allow for more sorbed NH₄⁺ (Figure 6e,g). Simulated NH₃ concentrations in porewater are similar to those observed in the Portnoy and Giblin [4] column experiments, ranging as high as 7.1 mmol/L (99 mg/L as N). PO₄^3--P concentrations are initially much higher (up to 0.12 mmol/L) in the DD sediment cores than in the DF cores (0.0019 mmol/L), probably due to the large sorbed fraction on the Fe oxides (jarosite-amm) in the DD cores that is released with reductive dissolution of jarosite-amm; PO₄^3--P concentrations decreased to 0.006 mmol/L once the Fe oxide is depleted. Simulated inorganic PO₄^3--P concentrations are initially higher than those observed in the Portnoy and Giblin [4] column experiments but decrease to background levels, contrary to measured inorganic PO₄^3--P concentrations in the experiment that increased to 0.15 mmol/L. This discrepancy may result from a Fe(III) phosphate mineral, such as vivianite, that is present in DF sediments but was not accounted for as an equilibrium phase in the model input. In acid SO₄^2- soils, Krairapanond, et al. [70] found that PO₄^3- sorption was mostly associated with Fe oxides and secondarily with Al oxides. In summary, smaller amounts of NH₄⁺, PO₄^3-, and Fe(II) were mobilized from the DF sediments as compared to the DD, because (1) less NH₄⁺ is sorbed to these sediments, and (2) dissolved Fe(II) mobilized by seawater is precipitated as Fe sulfide minerals.

Al concentrations increase significantly, from 0.09 mmol/L to briefly peak at 200 mmol/L after 7 hrs of salination (Table S4), then gradually declines below the original concentration after 79 days (Figures 6h and 8b). A large amount of dissolved Al (244 mmol/L) is exchanged off the acidic sediments during the first 88 days and some of this Al (1.4 mmol/L) is precipitated as Al(OH)_3(a); Portnoy and Giblin [4] observed a 6-fold increase in porewater Al after salination that was attributed to cation exchange in the acidic Herring River marsh sediments.

3.2. Basecase Simulations of Seawater Flooding in Diked Marsh Sediments

The biogeochemical responses of DF and DD porewaters to salination in the 12-yr run basecase models (Figure 7; Tables S6 and S7) are generally similar but extende over a longer period compared to solute responses of the microcosm experiment simulations, which began with a high rate of flushing in the first 12 hrs (Table 1). Because the initially high flow rate was not used in basecase models, these simulations can be used as a more conservative estimate of porewater movement through low permeability marsh sediments. Besides using jarosite-amm, additional Fe(III) source minerals were used in a sensitivity analysis to assess the effects of their varying reactivities and solubilities on simulated porewater biogeochemistry (Tables S2a–j). In addition to jarosite-amm (log k = -19.02), other Fe(III) oxide minerals include schwertmannite (another Fe(III) oxyhydroxy-sulfate mineral, log k = 2.25), Fe(OH)_2(a) (log k = 4.89) and goethite (log k = -1.0) [Thermoddem: 71] that were included in separate simulations (Section 3.1.3). In addition, jarosite-amm was simulated at a higher initial concentration (0.11 mol/L), an order of magnitude greater than was used in the basecase scenarios. Finally, extended basecase simulations that were run for over 100 yrs provide more long-term understanding of conditions and processes that affect constituents of concern.

The porewater pH in the basecase models during the first two years was acidic in both cores, dipping briefly to 4.8 at DF, and remaining below pH 5 at DD for over 2 yrs. The pH of DD porewater was much lower at the beginning of the simulation than the DF sediment because of the acidity produced by Fe sulfide oxidation and cation exchange that mobilizes cations including H⁺ and Fe(II). The low Fe(II) concentration at DF (Figure 7c) is consistent with the relatively small amount of jarosite-amm dissolution (Figure 7d,f), whereas the more abundant Fe(III) source (jarosite-amm) at DD resulted in accompanying increases of Fe(II), NH₃, Na, and PO₄^3- concentrations (Figure 7a–h). Alkalinity concentrations in the basecase simulations increased up to 0.01 mol/L in DF porewater and up to 0.05 mol/L in DD porewater after 5 yrs (Figure 7e,f) and demonstrate the greater mineralization of organic C from Fe(III) reduction, SO₄^2- reduction, and methanogenesis.

Concentrations of Al species increase with the drop in pH during both microcosm and basecase simulations (Figure 8 and Figure 9; Tables S3, S4, S6 and S7). Al concentrations are dominated by Al⁺³, AlSO₄⁺, and AlOH⁺² at both DF and DD sites. Al concentrations peaks much earlier in the microcosm simulations, because of the relatively high rate of flushing used in the beginning of the experiments. For the DF microcosm simulation, Al peaks only briefly at 0.36 days, coincident with the 3.6-hour dip in pH from 6.7 to 4.6, before decreasing to background levels by 0.5 days (Figure 8a). Because of the sustained acidic conditions at DD, porewater concentrations of Al species are much higher (0.2 mol/L) and extend over a longer period (80 days) compared to the DF sites (Figure 7bnew; Tables S3 and S4). For the basecase simulations, concentrations of Al⁺³, AlSO₄⁺, and AlOH⁺² also dominate the Al species. Al⁺³ has the highest concentration and AlOH⁺², the most toxic species, has the lowest concentration. Al in porewater at DF peaks between 319 days to 425 days and reaches 0.01 mol/L (Figure 9a; Table S6). The DD microcosm simulation has the same relative concentrations of Al species at DD, peaking at 0.19 mol/L with a coincident drop in pH to 4.2 at 315 days. Al concentrations doesn’t decline until the rise of pH after 2 yrs (Figure 9b; Table S7).

Extended (100-yr) basecase simulations show that S^2- concentrations in the more reduced DF sediments reach a stable concentration of 0.0034 mol/L after a year of resalination; after 1.6 yrs, S^2- concentrations at DD increased as Fe(III) was depleted and FeS minerals stop forming and begin to dissolve (Figures 7c,d), but slowly and still have concentrations of 0.04 mol/L at the end of the 100-yr simulations (Table S10b). Alternatively, S^2- in DD sediments are higher (0.02 mol/L) but decrease to near zero by 19 yrs (Table S10a,b). Simulated FeS in DF marsh sediments is much lower in abundance (by 15x) than that of DD, oxidizes more rapidly, and is depleted by 33 yrs (Table S10a). NH₃ concentrations in DF porewater increase at about 9 months, reach 1 mol/L after 4 yrs, and remain steady for 100 yrs. NH₃ concentrations in DD porewater reach 0.078 mol/L at just under 5 yrs and declines with the gradual rise in pH and the depletion of jarosite-amm to a negligible concentration after 24 yrs (Table S10b).

The 12-yr basecase scenario was then modified in a sensitivity analysis to assess the effects of varying types and amounts of Fe(III) source minerals on pH and solute concentrations. The dissolution of Fe(III) hydroxide minerals varied; the more soluble Fe(OH)_2(a), schwertmannite, and goethite were generally depleted earlier than jarosite-amm (Figure 10a–d; SI Tables S8a,b,c,d and S9a,b,c,d). All of the Fe(III) sources in the simulations were reduced by organic C oxidation and not by S^2- because they were depleted by the time S^2- concentrations increased after seawater was introduced and SO₄^2- reduction became dominant. Fe(III) minerals in DF sediments are much less abundant than in DD sediments and are generally depleted within 6 months, with the exception of the simulation including high (10x) jarosite-amm, which remained until about 7.7 yrs. Porewater concentrations of Fe(II) in DF sediments are also low–0.036 mol/L for 10x jarosite-amm, compared to 0.22 mol/L in DD sediments (Figure 10c,d). FeS minerals precipitation increases with the availability of Fe(II) from Fe(III) mineral dissolution (Figure 10c,d) and S^2- (reduced from SO₄^2- in seawater). Although jarosite-amm has a lower solubility than the other Fe(III) sources, simulations indicate that Fe(II) concentrations peaks higher and lasts longer (Figure 10c,d). Once the Fe(III) minerals are depleted, S^2- concentrations increase and remain high as SO₄^2- reduction becomes the dominant electron acceptor for organic C mineralization. The amount of FeS mineral in DF sediments remains steady or slightly declining over the 12-yr simulation, higher in jarosite-amm than in simulations with other Fe(III) minerals (Figure 10e,f; SI Table S8). The high (x10) concentration of jarosite-amm used in simulations also results in higher pH (greater than 8) compared to the default concentration and reflects the consumption of acidity during both Fe(III) mineral dissolution and FeS precipitation (Tables S8d and 9d).

3.3. Biogeochemical Implications of Seawater Restoration

The microcosm experiment models and the basecase models can be used to better understand the biogeochemical implications of seawater restoration on the Herring River and Wellfleet Harbor by incorporating speciation for constituents of concern (Al), varying the Fe(III) hydroxide source and strength, varying the flow rate, and varying the time interval. These models could also be applied to other freshwater marshes along coastal areas that are vulnerable to changing climatic conditions and sea-level rise. Aquatic organisms are adversely affected by low pH, suboxic or anoxic conditions, and high concentrations of dissolved Fe, S^2-, and Al. Fluxes of solutes such as NH₄⁺ and PO₄^3- to coastal water could cause eutrophication and hypoxia.

The sediment column salination simulations (and the column experiments) are a 1-D simplification of the actual field conditions that contain temporal and spatial variability associated with the seasons, tidal effects, precipitation, interaction with groundwater, and the heterogeneities in biogeochemistry, sediment texture, organic matter content, bioturbation, and flow rates. For example, organic matter decay rates and organo-metallic complexation are likely to be higher during warm seasons [39,72]; Al toxicity increases with temperature; the acidity and concentrations of nutrients, Fe, S^2-, and Al vary with organic matter content [39,73]; and nutrient fluxes will vary with tide amplitude because of varying dilution by seawater. These results can be used to determine and help isolate biogeochemical processes, estimate rates, and determine potential concerns to aquatic life. In particular, the models can help understand the biogeochemical processes, the interaction of constituents of concern with pH, redox, and complexation, and the general timing and sequence of the contaminant releases and reactions. Future work could include coupling the PHREEC models with two- or three-dimensional groundwater flow and solute transport models to simulate progressive mixing of tidally varying saltwater into the freshwater marsh sediments.

Several constituents derived from the DF and DD marsh sediments are a concern to both freshwater and saltwater aquatic health in the Herring River and Wellfleet Harbor (Table S5). Since 2021, biogeochemical conditions similar to those expected during marsh restoration have been observed episodically at Herring River monitoring stations at Duck Harbor due to overwash that carries saltwater over the dune crest into the Herring River watershed (Figure 2; Timothy Smith, National Park Service written communication, 2023). These overwash events, which are observed almost monthly during spring tides, have caused a significant die off of vegetation and have resulted in low pH and elevated concentrations of Fe(II) and NH₄⁺ in water samples and/or continuous monitors in the Herring River (Sophia Fox, National Park Service written communication, 2023).

In the microcosm experiment simulations for DF sediments, pH was as low as 4.7 after about 8 hrs, which is well under the freshwater-quality criteria range (6.5-9), until waters became saline at 8.7 hrs. At 8.7 hrs and later, the pH was well under the acceptable saltwater-quality criteria range (6.5-8.5) and remained so until about 400 days; pH varying outside those ranges approach the lethal limits of some species [74] (Table S5). The basecase models showed delayed changes in solutes, and the pH drop began over 300 days into the simulation.

Dissolved Fe in microcosm simulations were below the freshwater health criterion of 1,000 g/L for both DF and DD cores. Although there are no saltwater health criteria for Fe [45], Fe can be detrimental to the health of fish, invertebrates, and periphyton and can cause oxidative injury to various organs and physical damage to the gills [75,76]. Fe concentrations were well over 5 mmol/L (300 µg/L) in the DD core after 10 hrs until near the end of the simulation (Table S5).

Table S5 shows both the EPA’s 2018 and the 1988 Al guidance criteria for freshwaters [74,77]. The early peak in Al concentrations for the DF core is above EPA’s 2018 guidance, but below the Al guidance range that considers pH and hardness [45]. Al concentrations in saline water exceed 24 μg/L from 8.5 hrs until 58 day of the DF microcosm experiment simulation. For DD sediments, the pH and temperature affect the proportion of the neutral NH₃ species and therefore are used to determine the freshwater aquatic criteria for NH₃. There are two exceedances of the total NH₃ saltwater aquatic criteria for the DD microcosm simulation, 19 mg/L (1.4 mmol/L) and 16 mg/L (1.1 mmol/L) near the end of the 589 day simulation after pH had risen and NH₃ was more dominant (Table S5) [46,78]. Freshwater criteria for NH₃ are pH, temperature, and life-stage dependent, whereas those criteria for saltwater are only pH and temperature dependent.

Al is an example of a constituent of concern that can vary significantly in porewater with biogeochemical conditions during seawater restoration, as illustrated in differences between the microcosm experiment models and the basecase models (Figs 5 and 6). High Al concentrations are most common in freshwater because of the greater solubility of Al under the often acidic (<7 pH) conditions. Al toxicity is reduced by increases in cations including Ca²⁺ and Mg²⁺ that compete with Al³⁺ for uptake by aquatic organisms [79,80] and can be measured as total hardness [45]. Al toxicity also decreases when Al complexes with dissolved organic carbon (DOC), thereby reducing the bioavailability to aquatic organisms [79,81,82]. Al concentration, speciation, and toxicity in estuaries varies with pH and with temperature, DOC concentrations, and base cations, including Ca and Mg. Al minerals are moderately soluble below pH 5, with toxicity increasing up to a pH of about 6.0 (at 25° C), at which point amorphous Al(OH)_3(a) or gibbsite [83,84] dissolve and aqueous Al peaks. Al species were included in the microcosm and basecase salination models of the Herring River sediment porewater (Figure 8 and Figure 9; Tables S3, S4, S6 and S7). The toxic Al species Al(OH)₂⁺ dominates between pH 5 and 6 [83]. Concentrations of Al(OH)₂⁺ exceeded the freshwater acute Al value (varies with temperature, pH, DOC, and hardness) in DD sediments until saline conditions ensues after 10 days, and thereafter exceeds the chronic Al value (24 µg/L; [85]) for more than 500 days (SI Table S5). With the basecase models, Al concentrations increases at both sites with the drop in pH, peaking after about one year, and consisted mostly of Al⁺³, AlSO₄⁺, and AlOH⁺² (Figure 9). As mentioned previously, Al toxicity decreases when Al complexes with DOC, which varies seasonally, although DOC complexation was not simulated here. In addition, cold water will have a greater proportion of toxic species at higher pH values than warm water [86]. As a result of these factors, seasonality could have important effects on Al toxicity in the Herring River, and salination during colder months could be problematic to aquatic life.

Another potentially important factor not studied here but directly relevant to post restoration biogeochemical conditions is the compaction of sediments resulting from accelerated aeration, organic decomposition, and subsidence that occurs following seawater restoration [4,5]. Portnoy and Giblin [4] observed 6-8 cm of subsidence in the DF microcosm cores. The diking has also blocked the flood-tide sediment source to the marsh and led to further subsidence [5]. Because of this subsidence and exacerbated by sea-level rise, tidal restoration can lead to a greater extent of seawater inundation onto marsh sediments compared to natural conditions [87,88]. The prolonged seawater flooding resulting from this subsidence could therefore lead to a greater extent of desorption and redox chemistry changes, and larger pulses of acidity, metals, and nutrients that are harmful to aquatic ecosystem, as described in this section.

The persistence of FeS (and presumably also FeS₂) in extended simulations, particularly in DD sediments, reflects the tidally flooded conditions common to normal salt marsh soils in which Fe, S, and acidity are generally immobilized. Burton, et al. [20] found that tidal reflooding of wetland acid SO₄^2- soils resulted in enrichment of Fe²⁺ and SO₄^2- and pyrite reformation within the subtidal zone and did not result in rapid restoration. However, extended simulations of Herring River sediments show that pH increases and concentrations of most solutes in DF and DD sediments converge over time (Table S10a,b). Simulated concentrations of S^2- in DD porewater remain elevated for several years (Figure 7) before decreasing to near zero after almost 20 yrs (Table S10a,b); S^2- toxicity is less of a concern in DD sediments because the high concentrations of Fe(II) mobilized by seawater precipitate sulfides. S^2- toxicity can cause reduced plant uptake [89], and various fish species are sensitive to unionized H₂S [74].

These biogeochemical modeling simulations illustrate the complexities and improve understanding of how diked or otherwise freshwater wetlands and their globally important carbon reservoir respond biogeochemically to dike removal and sea-level rise. Organic C mineralization and its effect on biogeochemical cycling of Fe, S, CH₄, N, and P is intimately connected to microbe-soil-plant interactions in sediments [90]. Both the DF and DD sediments have NH₄⁺ sorbed on silts and clays and PO₄^3- likely associated with abundant Fe and Al oxides (Dent 1986). As shown in the model simulations, salination has the potential to mobilize a large pulse of NH₄⁺, PO₄^3-, H⁺, and Fe(II) (as already observed in part at the Duck Harbor overwash area on the Herring River; Figure 2), as well as potentially toxic levels of Al. The large amount of S^2- stored in these diked sediments, even the drained aerobic sediments [58], can result in a sharp increase in acidity but will eventually moderate by seawater buffering. Therefore, gradual restoration and careful monitoring of tidal restorations in diked wetlands could help mitigate the immediate negative effects of biogeochemical changes associated with restoration activities. With the resalination of diked coastal wetlands such as the Herring River, SO₄^2- reducers will likely overtake methanogens in organic matter decomposition, which increases thermodynamic energy yields and subsequently increases the turnover of N and P.

4. Summary and Conclusions

Biogeochemical models were developed and calibrated to greenhouse microcosm experiments that were conducted to assess the biogeochemical effects of seawater restoration on the historically diked Herring River marsh for diked, flooded (DF) sediments and diked, drained (DD) sediments Portnoy and Giblin [4]. Models for 1-D reactive transport of solutes in these marsh sediment porewaters were developed using PHREEQC, and matching the simulated solute curves with those results from Portnoy and Giblin [4] microcosm experiments. Specifically, rates of organic C mineralization, and timing and magnitude of Fe(II) and S^2- concentration peaks and pH were established as fit parameters. The sediment porosities, amounts of reactive Fe(III), number of sorption sites for PO₄^3-, and the cation exchange parameters for NH₄⁺ and Fe(II) were determined from data from freshwater marsh sediment cores [4]. Conditions for the DF and DD cores were simulated by reacting solute concentrations from atmospheric deposition data with solids including FeS, jarosite-amm, and Al(OH), and then reacting these cores with several volumes of seawater over a nearly 600-day period. Basecase models were then modified from the microcosm experiment models to simulate salination over longer time periods using single flow velocities and varying Fe source mineral types and strength to better understand biogeochemical conditions. Desorption of Fe²⁺, S^2-, H⁺, NH₃, and Al was significant and simulated using both ion exchange and surface complexation reactions.

Based on the geochemical modeling of coastal marsh salination presented herein, aquatic organisms would initially (ranging from days to up to 2 yrs, depending on sediment permeability and porewater velocity) be adversely affected by acidic pH (primarily in DD sediments) and elevated concentrations of Fe²⁺, S^2-, NH₃, and Al released from porewaters [4]. Furthermore, fluxes of solutes such as NH₃ and PO₄^3- could cause eutrophication and hypoxia in coastal waters. The DF and DD sediments evolved from primarily methanogenic and Fe(III)-reducing (oxic where subaerially exposed) redox conditions, respectively, to a system dominated by the SO₄^2-- S^2- pair. The pH rose to neutral-alkaline levels with infiltration of seawater.

The reintroduction of seawater to marsh sediments in the microcosm simulations showed flushing of acidic porewater with high concentrations of Fe(II) (up to 3.4 mmol/L in DF sediments). Once DF sediments became saline at 8.7 hrs, the pH was below the acceptable saltwater-quality criteria range (6.5-8.5) until about 400 days, approaching the lethal limits of some aquatic species. The simulated pH in DD porewater was as low as 4.2, resulting from Fe sulfide mineral oxidation and desorption that mobilized H⁺ and Fe(II), and well under the freshwater-quality criteria range (6.5-9) until waters became saline. Simulated redox conditions in DF porewater shifted from methanogenic to SO₄^2- reducing after 10 days; concentrations of S^2- increased at 14 days and reacted with Fe(II) to form FeS beginning at 33 days, lowering the dissolved Fe(II) level. In DD sediments, the predominant redox process quickly became reductive dissolution of Fe(III), and Fe(II) concentrations increased to 40 mmol/L after 130 days as a result of both Fe(II) desorption and jarosite dissolution. NH₄-N and PO₄^3--P concentrations in DF porewater increased steadily throughout the microcosm simulation, ranging up to levels (0.52 mmol/L and 0.0033 mmol/L, respectively) that could cause eutrophication and harmful algal growth [69,77; Table S5]; NH₃-N exceeded the chronic saltwater aquatic criteria and could harm aquatic health [77]. NH₃-N concentrations in DD porewater increased to 7 mmol/L and exceeded the total NH₃ saltwater aquatic criteria near the end of the simulation. NH₃ was much higher in DD porewater than in DF sediments because of the greater decay rates of SOM and the acidic soils that allowed for more NH₄⁺ sorption and represent a greater concern for aquatic life. High concentrations of Al (16 mmol/L) peaked early from desorption, then decreased to 5.5 mmol/L by 10 hrs and 0.012 mmol/L by the end of the microcosm experiment simulation. After conditions became saline, Al in DF porewater exceeded the 24 μg/L saltwater-quality guideline [91] from 8.5 hrs until 58 days of the simulation. Al concentrations in DD porewater exceeded the 24 µg/L saltwater aquatic criteria [91] for DD sediments once the porewater became saline (10 hrs) and thereafter exceeded the chronic saltwater Al value (24 µg/L) for more than 500 days. Although Al and most metals would be much less soluble once seawater infiltrates marsh sediments and pH increases, the toxic species Al(OH)₂⁺ dominated the total Al concentrations under the initially acidic conditions and exceeded the USEPA freshwater acute Al value in DD sediments (Table S5). All of these effects pose a potential threat to aquatic organisms. Although the exact outcomes of resalination are difficult to predict, these simulations suggest that there could be potential short- and long-term harm associated with porewaters of both DF and DD sediments of the Herring River. Because of their high peat content, DF sediments are also affected by significant subsidence following organic matter mineralization and as a result could be more prone to prolonged seawater flooding, desorption, and a larger pulse of acidity, metals, and nutrients.

The basecase simulations (run over 12 yrs) of DF and DD porewaters showed similar biogeochemical responses to salination as the microcosm simulations, but the lower flow rates represent a conservative flow estimate that might be observed in low permeability marsh sediments. These basecase reactive transport models, which are based on calibrated real-world diked marshes, allow simulations of varying flow and biogeochemical conditions, and speciation of contaminants such as Al. Understanding the types, concentrations, and distribution of Fe(III) source minerals in coastal wetland sediments is critical to aquatic life and coastal fisheries because they can release large amounts of exchangeable Fe(II) as well as NH₃ and PO₄^3- [4,21,92], but also provide sources for FeS formation that sequester S^2-. Jarosite-amm, a likely Fe(III) source in Herring River sediments, has a lower solubility than other possible Fe(III) sources [93] and simulations indicate that Fe(II) concentrations peaked later and longer than for Fe(OH)_2(a), goethite, or schwertmannite. Higher simulated amount of jarosite-amm resulted in higher pH compared to those in the basecase models and showed a greater consumption of acidity during both Fe(III) mineral dissolution and FeS precipitation [41]. All of the Fe(III) sources in the simulations were reduced by organic C oxidation and not by S^2- because they were depleted by the time seawater was introduced, SO₄^2- reduction became dominant, and S^2- was sufficiently available. Alkalinity concentrations in the high jarosite-amm simulations increased up to 0.01 mol/L in DF porewater and up to 0.05 mol/L in DD porewater after 5 yrs and demonstrate the greater mineralization of organic C from Fe(III) reduction, SO₄^2- reduction, and methanogenesis.

Extended basecase simulations were run for over 100 yrs to provide long-term understanding of conditions and processes that affect constituents of concern. For example, S^2- concentrations in the more reduced DF sediments reach a stable concentration of 0.0034 mol/L after a year of resalination, whereas S^2- in DD sediments are higher (0.02 mol/L) but decrease to near zero by 19 yrs (Table S10a,b). Simulated FeS in DF marsh sediments was much lower in abundance (by 15x) than that of DD, oxidized more rapidly, and was depleted by 33 yrs (Table S10a). NH₃ concentrations in DF porewater increased rapidly at about 9 months, reached 1 mol/L after 4 yrs, and remained steady for 100 yrs. NH₃ concentrations in DD porewater reached 0.078 mol/L at just under 5 yrs and declined to a negligible concentration after 24 yrs with the gradual rise in pH and the depletion of jarosite-amm.

This study was undertaken to help inform National Park Service staff on the level and duration of expected risk and the best strategies associated with dike opening. Although there is significant uncertainty in 1-D geochemical simulations, this study utilizes site-specific data that are common to data collected in other tide-restricted marsh systems worldwide that contain low pH and harmful or toxic levels of acidity and constituents including S^2-, NH₄⁺, Al, and Fe(II) [8]. In addition, these findings are consistent with recently observed spring tide overwash events at the Herring River that have led to locally acidic and Fe-fouled conditions (Sophia Fox, National Park Service written communication, 2023). The availability of data in similar systems and the consistency of results to overwash events demonstrate that geochemical simulations could potentially be used in other areas.

Restoration of tidal flushing to the Herring River would provide twice-daily infusions of seawater from Cape Cod Bay and Wellfleet Harbor, which would increase the pH and concentrations of alkalinity and DO; without regular flooding with aerobic seawater, the main stem of the diked Herring River can be hypoxic to anoxic during summer months. Further, most organisms in the Herring River would not be exposed directly to porewater discharge, and there would be some level of dilution with river and estuary water. However, tidal marsh restoration after 110+ years of diking could mobilize an initial pulse of high acidity, Fe(II), and Al, along with more gradual increases of NH₄⁺, S^2-, and PO₄^3-. Conducting resalination slowly and monitoring it carefully could help reduce the severity of the pulses of acidity, Fe(II), and Al [5]. Seasonality can affect the toxicity of constituents released from Herring River sediments that could be considered together with other factors such as herring migrations, algal growth, and shellfishing. For example, tidal restoration during colder months could minimize risk to herring, algae, and shellfish because biological activity and nutrient levels are lower than during summer months [39,94]. Alternatively, tidal restoration during colder months could be problematic for metals such as Al because of its decreased organic complexation that limits toxicity and also because the cold-water temperatures host a greater proportion of toxic Al species at higher pH values than warm water [86]. Therefore, stakeholders likely need to consider multiple factors when making decisions about restoration activities.

Future work could include coupling the PHREEC models with two- or three-dimensional groundwater flow and solute transport models to simulate progressive mixing of tidally varying saltwater into the freshwater marsh sediments. These results can be used to help understand the general timing and sequence of the contaminant releases and reactions (including complexation, desorption, and degradation of organic matter), and the potential concerns to aquatic life.