Reactive Transport Modeling of Subaqueous Sediment Caps and Implications for the Long-Term Fate of Arsenic, Mercury, and Methylmercury

Abstract

A 1-D biogeochemical reactive transport model with a full set of equilibrium and kinetic biogeochemical reactions was developed to simulate the fate and transport of arsenic and mercury in subaqueous sediment caps. Model simulations (50 years) were performed for freshwater and estuarine scenarios with an anaerobic porewater and either a diffusion-only or a diffusion plus 0.1-m/year upward advective flux through the cap. A biological habitat layer in the top 0.15 m of the cap was simulated with the addition of organic carbon. For arsenic, the generation of sulfate-reducing conditions limits the formation of iron oxide phases available for adsorption. As a result, subaqueous sediment caps may be relatively ineffective for mitigating contaminant arsenic migration when influent concentrations are high and sorption capacity is insufficient. For mercury, sulfate reduction promotes the precipitation of metacinnabar (HgS) below the habitat layer, and associated fluxes across the sediment–water interface are low. As such, cap thickness is a key design parameter that can be adjusted to control the depth below the sediment–water interface at which mercury sulfide precipitates. The highest dissolved methylmercury concentrations occur in the habitat layer in estuarine environments under conditions of advecting porewater, but the highest sediment concentrations are predicted to occur in freshwater environments due to sorption on sediment organic matter. Site-specific reactive transport simulations are a powerful tool for identifying the major controls on sediment- and porewater-contaminant arsenic and mercury concentrations that result from coupling between physical conditions and biologically mediated chemical reactions.

1 Introduction

Subaqueous sand caps are a remedial alternative for managing contaminated sediments that have received considerable attention in recent years (Azcue et al. 1998; Palermo 1998; Wang et al. 2004; USEPA 2005). This technology involves the placement of clean sediment or sand over contaminated material to physically isolate chemical constituents from ecological receptors. Although relatively simple and optimum in low-energy environmental conditions (e.g., lakes, estuaries, low-velocity river reaches), the effectiveness of capping can be compromised when there are diffusive and/or advective fluxes of chemical contaminants from underlying groundwater (Liu et al. 2001). Depending on the magnitude of these fluxes, successful mitigation may require (1) addition of mineral and organic substrates for adsorption (Moo-Young et al. 2001; Ying and Axe 2005; Viana et al. 2008), (2) favorable geochemical conditions that promote precipitation (or co-precipitation) of metals as crystalline and/or amorphous compounds (Sengor et al. 2007), and/or (3) addition of reactive chemical amendments to facilitate sequestration by adsorption or dissolution/precipitation (Jacobs and Forstner 1999; Yang et al. 2007; Kumpiene et al. 2008).

To simulate the effectiveness of a sediment cap for mitigating contaminants, it is important to consider porewater geochemistry and cap mineralogy because these parameters ultimately govern the partitioning of metals between the aqueous and solid phases. Studies indicate that porewater and mineral phases may evolve over time in capped sediments due to the establishment of microbial and macrofaunal populations. Redox stratification associated with biological activity can result in reductive dissolution of iron oxides, the precipitation of metal sulfides, and/or the formation of methylmercury (Slowey and Brown 2007; Himmelheber et al. 2008; Johnson et al. 2010).

Within the context of simulating sediment cap performance, geochemical and reactive transport modeling is a useful tool because it allows laboratory results to be extrapolated to long timescales. Although a number of models have been developed to understand contaminant transport and fate in subaqueous caps (Alshawabkeh et al. 2005; Lampert and Reible 2009; Viana et al. 2008; Arega and Hayter 2008; Go et al. 2009), most have focused on describing physical processes such as consolidation, bioturbation, advection, and dispersion. Geochemical and biogeochemical processes have generally been simplified in these models and have often relied on assumed distribution coefficients, rather than examining a complete set of integrated equilibrium and kinetic chemical reactions, to represent partitioning and fate.

The objective of the present study was to use a biogeochemical reactive transport model to simulate the mineralogical evolution and long-term fate of redox-active contaminants in a subaqueous sediment cap. Arsenic and mercury were selected for analysis because they are common contaminants in sediment and possess dissimilar geochemistry. Arsenic predominantly occurs as an oxyanion aqueous species, either arsenite or arsenate, and its partitioning to sediment is affected by the presence of iron oxide and arsenic sulfide phases. An important limitation on the effectiveness of capping of arsenic-contaminated sediment is its potentially high solubility under anoxic conditions (Mucci et al. 2000; O’Day et al. 2004). Mercury, in contrast, is a siderophile element that most commonly occurs in groundwaters and porewaters as dissolved inorganic or organic sulfide complexes. Although mercury can be sequestered in relatively insoluble sulfide phases [HgS(s) as cinnabar or metacinnabar], mercury can also be stabilized in solution by complexation with natural organic ligands and sulfide (Skyllberg 2008; Slowey and Brown 2007; Slowey 2010). Mercury may become methylated by sulfate- (and/or iron-) reducing bacteria that populate sediment caps (Ullrich et al. 2001), which is a key step in the production of bioavailable methyl mercury species. This work is preceded by similar studies that have simulated geochemical processes near the sediment–water interface (Dueri et al. 2003; Canavan et al. 2007; Jung et al. 2009; Couture et al. 2010) but did not specifically evaluate contaminant transport in sediment caps associated with contaminant remediation. A comprehensive review and compilation of thermodynamic and kinetic data were performed to construct the database for the reactive transport model presented here. Simulations of the 1-D transport of arsenic- and mercury-contaminated porewater into a clean (quartz) sand cap undergoing biogeochemical reduction near the sediment–water interface were used to (1) assess the relative effectiveness of an unamended cap for immobilizing metals and thus protect ecological receptors in the overlying aqueous environment and (2) identify key factors that should be considered during cap design.

2 Methods

2.1 Model Description

1-D reactive transport simulations were performed using the USGS-supported geochemical software PHREEQC (Parkhurst and Appelo 1999). Chemical processes included in the model were heterogeneous and homogeneous equilibrium speciation reactions, kinetic-based reactions describing biodegradation of organic carbon, reduction–oxidation (redox) transformations, mercury methylation, and demethylation. Physical processes included porewater diffusion, with and without porewater advective flow.

The model simulated a 1-meter-thick quartz sand cap emplaced on top of a layer of contaminated sediment in either an estuarine (salt water) or freshwater setting. The model domain was discretized into 100 grid cells, and constant concentration boundary conditions were applied at the ends (Fig. 1). The lower boundary consisted of an anaerobic porewater with elevated arsenic (10^−3.9 M, or 10 mg/L) and mercury (10^−7.3 M, or 10 µg/L), and the upper boundary and initial porewater composition consisted of either salt water or freshwater with no contamination (Table 1). The initial conditions within the cap consisted of entrained water with a chemical composition identical to oxidized surface water (either salt or fresh). The cap was also assumed to include a maximum of 1 % sediment organic matter (SOM) at the sediment–water interface, with the SOM concentration decreasing exponentially with depth in the top 0.15 m of the cap to simulate a habitat layer (Canavan et al. 2006). Porewater was allowed to diffuse and/or advect through the cap for a period of 50 years (Table 2). In some model simulations, the effect of background iron oxide concentrations on contaminant transport and fate was investigated by assuming an initial coating of goethite on the quartz sand in the cap (0.01 mol goethite/kg_sediment) (Table 2). Model simulations were also performed using influent porewater with a low concentration of arsenic (10^−5.2 M, or 0.5 mg/L) and mercury (10^−8.6 M, or 0.5 µg/L) to examine model response to input concentration (Table 3 summarizes the five different estuarine and freshwater model scenarios).

Fig. 1

PHREEQC model domain showing grid discretization, upper and lower boundary conditions, and sediment cap composition

Table 1 Chemical composition of aqueous solutions used in cap model

Table 2 Physical and chemical description of quartz sand cap

Table 3 Description of model scenarios

2.2 Thermodynamic Database

Equilibrium speciation calculations in PHREEQC utilized the default LLNL thermodynamic database (Delany and Lundeen 1990) with the following modifications. (1) Addition of an internally consistent thermodynamic database for arsenic aqueous species and minerals (Nordstrom and Archer 2003; Vlassopoulos et al. 2010), (2) inclusion of recently compiled thermodynamic values for inorganic (Powell et al. 2005) and organic mercury aqueous species (Skyllberg 2008), and (3) inclusion of equilibrium constants for FeS(s) and FeS(aq) (Rickard and Luther 2007).

2.2.1 Arsenic

Equilibrium constants for arsenic used in this study are reported in “Appendix 1.” Values for arsenic hydroxide and oxyanion species were taken from Nordstrom and Archer (2003), after adjusting the reported values for consistency with the LLNL database (Delany and Lundeen 1990). The sulfide speciation scheme and corresponding equilibrium constants of Vlassopoulos et al. (2010) were used for dissolved arsenic sulfide complexes. Adsorption of dissolved arsenate and arsenite to iron oxide phases was modeled using surface complexation constants and default surface site concentrations reported in the LLNL database for ferrihydrite from Dzombak and Morel (1990). Because model results indicated that goethite and magnetite were thermodynamically stable rather than ferrihydrite, the concentration of surface sites was reduced by a factor of 10 to account for differences in reactive surface area of the crystalline iron oxides (Dixit and Hering 2003). No adsorption to iron sulfide minerals was included in the model because prior studies indicated that precipitation of a realgar-type mineral phase is the dominant mode of arsenic attenuation in the presence of mackinawite (Gallegos et al. 2007) and because of a lack of surface complexation constants for sorption on sulfides consistent with the database. Realgar and orpiment were allowed to precipitate at thermodynamic equilibrium using the solubility constants reported in Vlassopoulos et al. (2010).

2.2.2 Mercury

Equilibrium constants for reactions with mercury are reported in “Appendix 2.” Prior studies have discussed the uncertainty in the value of the equilibrium constant for HgS(aq) (written as the species HOHgSH⁰ in “Appendix 2”). As discussed in Skyllberg (2008), experimental conditions designed to measure the formation constant for this species likely included colloidal mercury, which results in an overestimation of the equilibrium constant for dissolved mercury sulfide. Although the equilibrium formation constant for this species has also been estimated theoretically (Dyrssen and Wedborg 1986), this value has been shown to underestimate methylmercury concentrations in the presence of HgS(s) (Drott et al. 2007). Model results of this study predicted that metacinnabar was a stable phase. Therefore, the experimentally derived value of Skyllberg (2008) was used for HOHgSH⁰ in order to obtain a more accurate representation of methylmercury concentrations.

Sorption of mercury and methylmercury on SOM was included in the model by modifying equilibrium reactions for the formation of dissolved mercury-organic complexes (“Appendix 2”). This modification was accomplished by converting the aqueous complexation reactions reported in Skyllberg (2008) to exchange constants that conform to the Gaines-Thomas convention (Appelo and Postma 2005). This conversion depends on the concentration of total thiol (–SH) surface sites on SOM and, to a lesser extent, the protonated fraction of these sites. The reported values in “Appendix 2” are conditional stability constants assuming a concentration of thiol sites associated with SOM of 1.9 × 10⁻³ mol/L water. The abundance of thiol-type functional sites associated with organic matter was set to 0.047 meq/g organic matter (meq/g_OM) for both dissolved organic matter (DOM) and solid-phase SOM (Skyllberg 2008).

Mercury polysulfide complexes were not included in the model because (1) elemental sulfur was not predicted to form in the cap and (2) formation of these complexes is very weak compared with HgS(aq) in the presence of metacinnabar (Benoit et al. 1999; Jay et al. 2000; Drott et al. 2007). Mercury adsorption to iron sulfide minerals was excluded because the predominant mode of mercury attenuation in the presence of mackinawite is the precipitation of a HgS mineral (in this case, metacinnabar) rather than sorption (Skyllberg and Drott 2010).

2.3 Biogeochemical Reactions and Rates

Biodegradation of two organic carbon fractions was included in the model: (1) DOM, originating in both surface water and porewater, and (2) SOM, originating from surface water and mixed within the top 0.15 m of the sediment cap by bioturbating organisms (Boudreau 1998; Canavan et al. 2006). Concentrations are given in Tables 1 and 2, respectively.

The Monod-type description of biogeochemical reactions and rates in the model, including both primary reduction and secondary oxidation reactions (Table 4), was based on prior studies (Van Cappellen and Wang 1996; Hunter et al. 1998; Canavan et al. 2006). The overall degradation rate of organic matter (OM), including either SOM or DOM, is the sum of individual reaction rates [R _i ] of successive terminal electron acceptors (EA _i ):

$$ \begin{gathered} R_{i} = k_{\text{OM}} \left( {1 - \sum\limits_{0}^{i - 1} {f_{i} } } \right)\frac{{\left[ {{\text{EA}}_{i} } \right]}}{x}\quad {\text{and}} \hfill \\ \left\{ {\begin{array}{*{20}c} {x = 1} & {{\text{for}}\quad \left[ {{\text{EA}}_{i} } \right] > K_{{{\text{EA}},i}} } \\ {x = K_{\text{EA}} } & {{\text{for}}\quad \left[ {{\text{EA}}_{i} } \right] < K_{{{\text{EA}},i}} } \\ \end{array} } \right\} \hfill \\ \end{gathered} $$

(1)

where k _OM is the rate constant for OM degradation, f _i is the fraction of electrons consumed by the ith primary reduction half-reaction, [EA _i ] is the concentration of the ith terminal electron acceptor species [EA], and K _EA,i is the limiting concentration for the respective electron acceptor. As shown in Eq. 1, when the concentration of [EA] is below its limiting value K _EA, the corresponding primary reduction rate is reduced. Rate constants (k _i ) and limiting concentrations (K _EA) used in the model are reported in Table 5. The biodegradation rate constants for SOM and DOM were set to 0.002 and 0.0001 year⁻¹, respectively, based on literature review (Canfield et al. 1993; Hulthe et al. 1998; Fossing et al. 2004; Arzayus and Canuel 2005; Wallmann et al. 2006; Komada et al. 2004). Although the biogeochemical rate constants used in this study were not explicitly calibrated, all model parameters were selected to be within the range of values observed at field sites in order to be generally applicable.

Table 4 Kinetic reactions used in model simulations

Table 5 Parameter values for kinetic reactions (shown in Table 4)

Arsenate reduction and secondary arsenite oxidation reactions were added to the suite of reactions compiled from prior studies (Table 4). Arsenate reduction was assumed to precede reduction of ferrous to ferric iron in the sequence of terminal EA based on a larger (more negative) change in the Gibbs free energy of the half-reaction for arsenate reduction compared with goethite reduction. Values for the limiting concentrations and secondary reaction rates for iron sulfide oxidation by arsenate were estimated (Table 5) from rates reported for arsenate reduction by sulfide (Rochette et al. 2000).

The rate of methylmercury formation was calculated using the following equation developed by Gilmour et al. (2008):

$$ {\text{rate}}_{{{\text{CH}}_{ 3} {\text{Hg}}^{ + } }} = \lambda \times \left[ {\text{cells}} \right] \times \left[ {{\text{C}}_{1} } \right] \times \left( {\left[ {\text{OHHgSH}} \right] + \left[ {{\text{Hg}}\left( {\text{SH}} \right)_{2} } \right]} \right) $$

(2)

where the parameter λ represents the fraction of sulfate-reducing bacteria that methylate mercury at the same rate as Desulfobulbus propionicus (a value of 0.2 was used both in Gilmour et al. (2008) and here); [cells] is the number of cells of sulfate-reducing bacteria calculated from the model-determined sulfate reduction rate (M_SO4 year⁻¹, Table 4) and a reported cell density corresponding to this rate (1.63 × 10¹² cells year M ⁻¹_SO4 ) (King et al. 2000); the methylation rate parameter [C ₁] was assumed to be 8.54 × 10⁻⁸ M_MeHg L cell⁻¹ M_ΣHgSneutral ⁻¹ year⁻¹ (Gilmour et al. 2008); and the total concentration of dissolved neutral mercury sulfide complexes ([OHHgSH] + [Hg(SH)₂]) (M_HgS L⁻¹) was computed by PHREEQC. The equivalent expression re-written as Eq. 20 in Table 4 by combining terms in Eq. 2 is:

$$ {\text{rate}}_{{{\text{CH}}_{ 3} {\text{Hg}}^{ + } }} = k_{20} \times R_{5} \times \left( {\left[ {\text{HOHgSH}} \right] + \left[ {{\text{Hg}}\left( {\text{SH}} \right)_{2} } \right]} \right) $$

(3)

where the rate constant [k ₂₀] is 4.2 × 10³ M M ⁻¹_SO4 M ⁻¹_HgS year⁻¹ (M_HgS represents the total concentration of the two neutral mercury sulfide complexes) (Table 5), R ₅ is the rate of sulfate reduction (Table 4), and [HOHgSH] and [Hg(SH)₂] are the concentrations of the two neutral mercury sulfide aqueous complexes included in the model (“Appendix 2”). Demethylation was modeled by using a first-order rate constant with respect to the total MeHg concentration (dimethylmercury plus monomethylmercury species).

Performance of the model in determining mercury methylation was evaluated by simulating the controlled microcosm experiments of Johnson et al. (2010). Concentrations of methylmercury were calculated by the model based on the chemical description of the sediment and porewater from Johnson et al. (2010), the estimated organic carbon biodegradation rate from that study, and methylation/demethylation implemented in the PHREEQC model with the parameters given in Table 5. The model predicted a total MeHg concentration in the sediment after a 120-day simulation period of 230 pg/g, which is equivalent to 0.07 % of the total mercury concentration. This concentration compares favorably with reported concentrations of 215–232 pg/g (equivalent to 0.06–0.08 % of total mercury present as methylmercury) at the deepest measurement point in the sediment microcosms (0.03 m) (Johnson et al. 2010).

3 Results

Model results were assessed for the following: (1) the evolution of solid-phase contaminant sequestration, either by mineral precipitation or surface adsorption, in the subaqueous sediment cap (reported as total moles of solid/L of water to facilitate direct comparison with the aqueous concentrations), (2) the flux of contaminants to surface water over time, and (3) geochemical environments where capping would be most effective based on model outcomes. The latter criterion was assessed by comparing simulated aqueous and solid-phase concentrations with water and sediment quality screening values that might be employed during remedial design. Table 3 summarizes the model scenarios.

3.1 Arsenic

The simulated evolution of cap mineralogy and porewater geochemistry for the case of estuarine boundary conditions and diffusion-only porewater flux (scenario E1) are shown as a function of time and depth in Fig. 2a. The first 30 years of simulation results are illustrated, after which steady-state conditions are established. After the first decade, pH (~6.5–6.7) and dissolved oxygen (DO ~ 10⁻⁸ M) are constant within the sediment cap, with the exception of a thin oxidized zone (~0.02 m) at the sediment–water interface. Dissolved arsenic is highest at the base of the cap from porewater influx and decreases toward the sediment–water interface. The concentration of As(V) decreases with time in the lower half of the cap and dissolved arsenic is dominated by As(III) at steady state. Goethite precipitates at the oxidized sediment–water interface as a result of dissolved Fe(II) oxidation by diffusion of oxygen from surface water. Goethite also precipitates at the base of the cap, where more oxidized porewater (due to surface water initially entrained in the cap) mixes with influent, anaerobic porewater containing Fe(II) (Table 1). Goethite precipitation near the base of the cap (0.9–1.0 m depth) decreases upward as Fe(II) introduced by porewater is depleted. A zone of reduced iron sulfide minerals (mackinawite and pyrite) forms between ~0.02 and 0.3 m depth from the reduction of sulfate supplied by diffusion of seawater. These minerals replace goethite with time as the primary iron solid phases in the upper part of the cap. However, arsenic sulfide minerals do not form under these conditions. For the diffusion-only scenario, the total dissolved arsenic concentration in equilibrium with goethite in the top 0.1 m at the end of the model simulation (50 years) is 10^−5.5 M (0.26 mg/L) (Table 6).

Fig. 2

Concentrations of Fe and As mineral and dissolved species (in log moles/L of water) as a function of depth and time in the sediment cap for the case of diffusion only, log [As]_tot = −3.9 and log [Hg]_tot = −7.3 in influent porewater, and the following environments. a Estuarine (scenario E1), and b freshwater (scenario F1). Rows 1 and 3 (left–right): pH; dissolved oxygen; total dissolved As(III); total dissolved As(V). Rows 2 and 4 (left–right): goethite; mackinawite and pyrite; As adsorbed to goethite; realgar and/or orpiment

Table 6 Depth-averaged results (from 0 to 0.1 m) for cap simulations after 50 years

In the freshwater scenario (F1) (Fig. 2b), steady-state concentrations are not established in most of the cap until after ~20 years. Similar to the estuarine case, reduced conditions are established throughout the cap, but the steady-state pH is lower (~6) in the freshwater than in the estuarine simulation. Dissolved arsenic is dominated by As(III) at steady state, but As(V) is elevated from the base of the cap to ~0.15 m below the surface in the habitat layer. Arsenic sequestration by solid phases in the freshwater scenario is distributed between adsorption by goethite and a zone of realgar (AsS) precipitation at 0.1–0.15 m depth. Greater quantities of both minerals are precipitated in the upper part of the cap compared with the estuarine scenario. For goethite, this outcome is directly related to lower total sulfur concentration, and therefore lower sulfide levels, in freshwater. Lower total sulfur affects goethite precipitation in two ways: (1) no pyrite or mackinawite forms in freshwater, which increases the concentration of dissolved Fe(II) available for oxidation by surface water, and (2) less sulfide is available to reduce Fe(III) (and thereby dissolve goethite). As a result, an interval of goethite forms between ~0.15 and 0.3 m below the surface in addition to zones above the base of the cap and at the sediment–water interface (Fig. 2b). For realgar, the effect of lower sulfide in the freshwater scenario is to decrease the concentration of As-sulfide complexes and thus increase the concentration of dissolved As(OH)₃°, exceeding reaglar solubility. Realgar precipitation, however, does not greatly decrease dissolved arsenic concentrations in the cap because of its higher equilibrium solubility compared with other sulfide minerals. At 50 years, the total dissolved arsenic concentration at the top of the cap is predicted to be the same (10^−5.5 M) as for the estuarine scenario (Table 6), although the distribution of arsenic species is different throughout the simulation period.

Model outcomes for arsenic with advective transport in addition to diffusion in both estuarine (scenario E2) and freshwater (scenario F2) systems are shown in Fig. 3a and b, respectively. One effect of advection is to decrease the time required for the establishment of steady-state dissolved profiles to a period of <10 years for most species. A second is to introduce greater quantities of arsenic into the cap and therefore generate higher dissolved arsenic concentrations throughout the cap and near the sediment–water interface. As shown in Table 6, dissolved arsenic concentrations at the top of the cap are ~10 times higher than in the diffusion-only scenarios. The introduction of higher arsenic concentrations into the habitat layer (where sulfide is generated from the reduction of sulfate) causes more precipitation of realgar where the upward-propagating arsenic flux meets the downward-propagating sulfide flux. Arsenic is also sequestered in higher amounts by the thin goethite layer at the sediment–water interface. The advection scenarios highlight the competition between arsenic sequestration by adsorption to goethite and precipitation of arsenic sulfide minerals.

Fig. 3

Concentrations of As mineral and dissolved species (in log moles/L of water) as a function of depth and time in the sediment cap for the case of advecting groundwater, log [As]_tot = −3.9, log [Hg]_tot = −7.3 in influent porewater, and the following environments. a Estuarine (scenario E2), and b freshwater (scenario F2). Rows 1 and 2 (left–right): total dissolved As(III); total dissolved As(V); As adsorbed to goethite; realgar and/or orpiment

In order to examine model sensitivity to iron concentration in the cap and to contaminant flux from upwelling porewater, additional scenarios were examined for the case of an initial coating of goethite on quartz in the sand cap (scenarios E3 and F3) (Knapp et al. 2002; Kent and Fox 2004) and for lower influent contaminant concentrations (scenarios E4 and F4). As shown in Table 6, the effect of goethite coatings is to reduce the mass flux of arsenic through the cap and thus the dissolved and solid concentrations at the top of the cap. Lower influent porewater arsenic concentrations result in concentrations at the top of the cap similar to those observed with more goethite coatings because attenuation is primarily controlled by arsenic adsorption to goethite in these scenarios.

3.2 Mercury

For the estuarine, diffusion-only sediment cap (scenario E1) (Fig. 4a), dissolved mercury concentrations are highest at early times and near the base of the cap as Hg complexed to DOM. Mercury speciation changes to complexation with dissolved sulfide within the first 5–10 years as reduced conditions are established and sulfate reduction produces dissolved sulfide. The primary solid-phase sequestration mechanism for mercury is precipitation as metacinnabar [HgS(s)] once its solubility is exceeded with sufficient sulfide production. Metacinnabar precipitation initially occurs in two areas of the cap: (1) near the base where a metacinnabar front propagates upward and (2) a depth between 0.4 and 0.6 m where a metacinnabar front propagates downward. Dissolved sulfide concentrations build over time in the upper third of the cap, which results in higher concentrations of dissolved Hg–sulfide complexes and shifts the depth of the precipitation front for metacinnabar downward. In the upper habitat layer (0–0.15 m), Hg is attenuated in the cap by sorption to SOM.

Fig. 4

Concentrations of S and Hg mineral and dissolved species (in log moles/L of water) as a function of depth and time in the sediment cap for the case of diffusion-only, log [As]_tot = −3.9, log [Hg]_tot = −7.3 in influent porewater, and the following environments. a Estuarine (scenario E1), and b freshwater (scenario F1). Row 1 (left–right): total dissolved sulfide (all species); Hg-DOM complexes; total dissolved Hg-sulfide complexes; dissolved methylmercury (mono- and di-methyl). Row 2 (left–right): metacinnabar (HgS); total mercury adsorbed to SOM, total MeHg adsorbed to SOM

Dissolved and solid-phase mercury profiles in freshwater (scenario F1) are generally similar to the estuarine scenario (E1) except for (1) the occurrence at a shallower depth of the metacinnabar precipitation front and (2) greater mercury complexation to SOM in the top 0.15 m of the cap (Fig. 4b). Differences in mercury speciation and partitioning can be attributed to less total sulfur, and therefore less sulfide, in freshwater versus estuarine systems. In the freshwater case, more time is required for sulfide concentrations to build to levels where metacinnabar precipitates in the cap. Thus, the concentration of Hg–sulfide complexes is lower at the top of the cap, resulting in a shallower depth of precipitation of metacinnabar and greater partitioning to SOM. The steady-state total dissolved mercury concentration in the top 0.1 m in the estuarine sediment cap is slightly less than in the freshwater cap (Table 6) because metacinnabar precipitation occurs closer to the sediment–water interface in the latter case.

The effect of porewater advection (scenarios E2 and F2) is to shift the metacinnabar precipitation front upwards in the cap, which increases dissolved mercury concentrations at shallower depths (Fig. 5a, b). For estuarine systems, steady-state dissolved concentrations in the top 0.1 m of the cap increase from 10^−11.2 M (1.2 ng/L) to 10^−10.6 M (5.0 ng/L) (Table 6). For freshwater, steady-state concentrations increase from 10^−11.0 M (1.9 ng/L) to 10^−10.6 M (5.5 ng/L). These concentrations are generally within the range reported for estuarine surface water (10^−9.8–10^−12.4 M) (Fitzgerald et al. 2007). Unlike arsenic, an initial coating of goethite in the cap causes little change in dissolved mercury concentrations because mercury adsorption to iron oxide minerals was not included in the model and is expected to be less important than sorption to OM (Feyte et al. 2010).

Fig. 5

Concentrations of Hg mineral and dissolved species (in log moles/L of water) as a function of depth and time in the sediment cap for the case of advecting groundwater, log [As]_tot = −3.9, log [Hg]_tot = −7.3 in influent porewater, and the following environments. a Estuarine (scenario E2), and b freshwater (scenario F2). Rows 1 and 2 (left–right): total dissolved Hg; metacinnabar (HgS); total dissolved methylmercury; methylmercury adsorbed to SOM

3.3 Methylmercury

Methylmercury formation is predicted to occur within the upper 0.15 m of the cap in both estuarine and freshwater scenarios (E1 and F1) (Fig. 4a, b). Dissolved methylmercury concentrations are approximately 1–4 % of the total mercury, which is typical for sediment environments (Ullrich et al. 2001). The primary attenuation mechanisms for methylmercury are demethylation and adsorption of methylmercury to SOM. In these model scenarios, the mass of adsorbed methylmercury generally exceeds the dissolved mass (Figs. 4, 5). Dissolved methylmercury concentrations are slightly lower at the top of the freshwater cap than the estuarine cap because the concentrations of neutral mercury sulfide complexes available for methylation are lower (owing to lower total sulfur) (Table 6). In contrast, methylmercury concentrations in sediment are higher in fresh than salt water because concentrations of dissolved methylmercury-sulfide complexes are lower and more methylmercury is available for sorption to SOM.

Advective transport increases the steady-state concentration of methylmercury in both estuarine and freshwater systems (Fig. 5). The highest dissolved concentrations occur in the estuarine sediment cap with advection, and the highest sediment-bound methylmercury levels occur in freshwater. Concentrations are predicted to be lower than those typically reported for contaminated sediments (Bloom et al. 1999) and are within the range of estuarine and riverine systems (Conaway et al. 2003; Marvin-DiPasquale et al. 2000; Marvin-Dipasquale and Oremland 1998; Goulet et al. 2007; Fitzgerald et al. 2007).

3.4 Comparison with Screening Levels for Water and Sediment Quality

Model outcomes at 50 years were compared with two quality criteria: criterion continuous concentration (CCC) for dissolved contaminants in saltwater and freshwater (USEPA 2009), and sediment quality guidelines (SQG) using the effects range-low (ERL) values from (NOAA 1999). The CCC is an estimate of the highest concentration of a material in surface water to which an aquatic community can be exposed indefinitely without resulting in an unacceptable effect. The ERL guidelines are intended to represent sediment concentrations of a contaminant below which adverse effects on biota rarely occur and are derived from the 10th percentile values of compiled studies below which adverse effects occurred. Comparison of model outcomes (Table 6) showed that for all high influent arsenic concentration scenarios, dissolved arsenic concentrations at the top of the cap (depth-averaged from 0 to 0.1 m) exceeded CCC levels for both freshwater (10^−5.7 M) and saltwater (10^−6.3 M) (Fig. 6a). When influent arsenic concentration was reduced (20 times lower, scenarios E4 and F4), depth-averaged concentrations in both estuarine and freshwater sediment caps were below CCC levels (36, and 150 µg/L, respectively). In the presence of iron oxyhydroxides (as goethite) on quartz (scenario E3), concentrations approached but did not exceed the CCC in the top 0.1 m of the cap (Fig. 6a). Dissolved mercury concentrations (Table 6) were several orders-of-magnitude below the CCC for all scenarios (data not shown).

Fig. 6

Comparison of predictions from model scenarios (Table 3) to water and sediment quality screening levels at 50 years. a Dissolved arsenic concentrations compared with CCC (yellow line) for estuarine model results. Solid green diffusion, low concentration (E4), dashed blue diffusion, high concentration (E1), stippled black diffusion with goethite coating (E3), dashed red advection + diffusion, high concentration (E2). Comparison of b arsenic and c mercury sediment concentrations in freshwater and estuarine conditions with SQG (SQG = yellow line) in diffusion-only (solid green) and diffusion plus 0.1-m/year advection (dashed red) scenarios

Sediment arsenic concentrations were below ERL levels for both high and low concentration scenarios (Fig. 6b). For mercury, sediment ERL concentrations were exceeded only within the lower 0.2 m of the cap where high concentrations of HgS(s) were precipitated and were below screening levels at the top (Fig. 6c). Methylmercury concentrations in sediments (not shown) are predicted to be below regulatory screening levels.

4 Discussion

4.1 Factors Controlling Cap Effectiveness

The analysis of mineralogical and porewater changes within a sediment cap by reactive transport modeling highlights primary factors that control the migration of dissolved contaminants from the cap into surface waters. An important observation from this analysis is that different but interrelated factors are responsible for arsenic or mercury attenuation in the cap. From the standpoint of (bio)geochemical reactions, controlling processes are (1) the concentration of DOM and SOM in the system, which influences (a) the amount of iron and sulfate reduction, and (b) the extent of mercury and methylmercury complexation with DOM in solution and sorption to SOM in sediments. (2) The total amount of iron in the system from either sediment Fe(III) oxyhydroxides or groundwater influx of dissolved Fe(II). Total iron controls (a) the precipitation of Fe(III)-oxide minerals and thus arsenic sorption in more oxidized zones, and (b) the concentration of dissolved Fe(II) in more reduced zones. (3) The total amount of sulfur in the system, examined in this study by the difference between typical estuarine and freshwater systems. Total sulfur determines the maximum amount of sulfate reduction and thus, (a) the precipitation of iron, arsenic, and mercury sulfide phases and (b) the concentration of aqueous sulfide complexes. In addition, the amount of porewater advection, compared with diffusion only, will determine the flux of contaminants moving through the cap and thus the ability of (bio)geochemical processes to attenuate them.

In estuarine systems, total sulfur is in excess because of diffusion of sulfate from overlying surface water and from porewater influx of sulfate (Table 1). Because dissolved sulfide is generally higher throughout the cap than in the freshwater case, iron sulfide minerals (mackinawite and pyrite) form in the upper part of the cap, and the formation of iron oxide minerals (modeled here as goethite) is less important and limited to the bottom of the cap. Therefore, arsenic sorption to iron oxides is lower overall. No precipitation of arsenic sulfide minerals (realgar or orpiment) is observed because of the higher solubility of arsenic sulfides compared with iron sulfides or metacinnabar at the pH of these systems. In the freshwater simulations with lower total sulfur, iron sulfide minerals do not form, a small amount of realgar precipitates within the top 0.15 m, and goethite is stable below the top 0.15 m. Arsenic sorption to goethite occurs within the cap and at the oxidized sediment–water interface. However, dissolved arsenic concentrations in the top 0.1 m are similar to the estuarine case because the sorption capacity of iron oxides is the limiting factor for arsenic attenuation. Thus, dissolved arsenic flux from the cap for scenarios with high porewater concentration or with advection is higher than water quality screening levels because of insufficient attenuation by sorption to iron oxides and higher solubility of arsenic sulfide minerals compared with mercury sulfide.

In real systems, metastable Fe(III)-oxide phases such as ferrihydrite, or metastable Fe(II, III) oxides such as green rust-type phases, may form instead of the thermodynamically stable phases of goethite and magnetite used in this model. These metastable iron phases generally have higher surface areas than their stable counterparts and thus may be more effective sorbents of arsenic. On timescales longer than those modeled, continued precipitation of iron oxides may create a natural sorption barrier for arsenic (Jung et al. 2009), but only if the rate of oxide precipitation and generation of reactive surface sites exceeds the rate that arsenic is introduced into the system. Because sorption is generally more effective for arsenic attenuation than sulfide precipitation, cap performance depends on arsenic concentration and flux through the cap, the concentration and type of iron oxide minerals and their sorption sites, and porewater pH within the cap.

Mercury and methylmercury concentrations are predicted to remain below water and sediment quality criteria levels in both estuarine and freshwater systems once sulfide concentrations increase after the first 5–10 years because of precipitation of metacinnabar in the lower parts of the cap. Due to the very low solubility of mercury sulfide minerals, equilibrium dissolved mercury concentrations are low. In the upper part of the cap without metacinnabar, predicted methylmercury concentrations are slightly higher in the estuarine than the freshwater setting because of higher rates of sulfate reduction. The distribution of mercury and methylmercury between aqueous and solid phases in this zone is controlled by sorption to SOM, aqueous complexation with DOM, and aqueous complexation with sulfide. Therefore, the amount and distribution of OM are the primary controlling factors for mercury and methylmercury flux from the top of the cap. The concentration of thiol sites associated with SOM is the primary factor for mercury attenuation by sorption in sediments, which competes with sulfur binding sites in DOM as a possible mode of mercury transport in the aqueous phase.

4.2 Model Uncertainties

Important sources of uncertainty in any chemical model are the values of equilibrium constants for reactions assumed to reach thermodynamic equilibrium, and in the rate constants for reactions assumed to be controlled by kinetic processes. For the majority of equilbrium reactions in the LLNL database and those that were added in this study (see “Appendices 1 and 2”), equilibrium constants are generally accepted and robust, with the following notable exceptions. As discussed above, there is a large uncertainty regarding equilibrium constants for aqueous mercury complexation (Skyllberg 2008), which affects model-predicted methylmercury levels (Drott et al. 2007). There is considerable uncertainty regarding both the identities and stability constants of arsenic and mercury mono- and poly-sulfide complexes, which may be important in high sulfide systems where elemental sulfur is stable (Rickard and Luther, 2006). Although mercury and mercury-DOM complexes may adsorb to quartz and iron oxides (Tiffreau et al. 1995; Backstrom et al. 2003), sorption to oxides is likely secondary to complexation with dissolved sulfide and DOM, and sorption to SOM, and therefore was not included. Another simplifying assumption was the use of an average diffusion coefficient for all aqueous species, which neglects the effects of porosity and tortuosity on the diffusion rates of specific ions or complexes. A greater source of uncertainty, however, probably results from the kinetic model and rate constants selected for abiotic and biotic chemical reactions. The mechanisms and rates of mercury methylation and demethylation are complex because (1) iron-reducing bacteria may contribute to methylmercury concentrations (particularly in freshwater caps; see Fitzgerald et al. 2007; Fleming et al. 2006), and (2) methylmercury production and degradation depends on temperature, pH, and other environmental variables (Ullrich et al. 2001). Uncertainties in these model parameters highlight the need for site-specific concentrations and rates in order to model remediation scenarios at contaminated sites.

Physical processes not included in the model may additionally affect contaminant distribution and fate within the cap. For example, this study did not explicitly include processes such as consolidation, bioturbation, bioirrigation, and sedimentation. Consolidation can temporarily increase the porewater flux through the cap, particularly at early times, thereby requiring a thicker cap to maintain low contaminant concentrations at the cap-water interface. Also, bioturbation and bioirrigation at the cap surface can affect the concentrations of both arsenic and mercury by creating a more diffuse zone of oxygenated conditions within the bioturbated layer, which could shift the sulfate/sulfide redox boundary deeper into the sediment (Dueri et al. 2003; Benoit et al. 2006). Inclusion of bioturbation and bioirrigation would likely increase the zone of Fe(III) oxide precipitation at the sediment–water interface, and thus the amount of iron oxide available for arsenic adsorption. The impact on arsenic breakthrough over long timescales would depend on the amount of Fe(II) available for oxidation and the arsenic flux. For mercury, greater introduction of oxygen into the habitat zone may decrease methylmercury production by slowing the rate of sulfate and/or iron reduction (Benoit et al. 2006). Finally, uncertainty is introduced by the imposed boundary conditions. For example, the constant concentration boundary conditions employed in this study give rise to sharp fronts of metacinnabar and goethite precipitation near the base of the cap. These sharp fronts are model artifacts that occur because of the prescribed boundary conditions. In reality, dissolved species may also migrate to some extent into the underlying contaminated sediment, thereby shifting the zone of precipitation downward.

4.3 Cap Design Implications

The sediment cap simulations illustrate the fate of arsenic and mercury in various contaminated sediment–water environment scenarios and highlight some of the key factors to be considered for cap design. Unamended subaqueous sand caps are not expected to be effective for reducing porewater dissolved arsenic concentrations below ecological screening levels when concentrations in upwelling porewater are high (Fig. 6a). As shown in Fig. 6b, there is a potential over time for the habitat layer of the sediment cap to become contaminated in excess of SQG by either adsorption to iron oxides or precipitation of arsenic sulfides under sulfate-reducing conditions. The former is expected to be a feature of the uppermost few millimeters of the cap near the sediment–water interface where exposure to oxygenated water promotes formation of iron oxides. Arsenic sorption may be enhanced if the rate of ferrihydrite or goethite precipitation exceeds the arsenic flux into the cap, or if there is a pre-existing pool of iron oxyhydroxides in the sediment cap material (scenarios E3 and E4, respectively). Precipitation of arsenic sulfides, on the other hand, is more likely in estuarine/marine porewater environments where the downward diffusive flux of sulfate and the upward flux of arsenic provide dissolved reactants and bacterial sulfate reduction promotes the precipitation of solid-phase sulfides. The depth and extent of the zone of accumulation depends on the porewater upwelling velocity and the competitive precipitation of iron sulfide minerals. Due to the modest solubility of arsenic sulfides, they may occur closer to the top of the active sulfate-reducing zone where dissolved sulfide concentrations are greatest (Fig. 6b).

In contrast, capping with unamended sand is a potentially viable option for reducing mercury flux under all of the conditions modeled. The effectiveness is a function of (1) the rate of sulfate reduction, (2) porewater upwelling rates, and (3) cap thickness. As shown in Fig. 6c, mercury contamination is restricted to the base of the cap due to precipitation of insoluble mercury sulfide. Therefore, cap thickness appears to be a key design parameter for preventing mercury contamination of the habitat layer for a given flow regime because the upper depth limit of mercury accumulation is enhanced by porewater advection. The minimum thickness will need to be greater than the depth of the habitat layer in order to create a zone of separation between ecological receptors and mercury accumulation (Simpson et al. 2002). For the range of conditions examined in our study, the maximum thickness is expected to be <1 m, which is consistent with field studies (Moo-Young et al. 2001). The effectiveness of a sand cap for limiting mercury exposure to biota will depend on a number of site-specific physical, hydrologic, and biogeochemical factors that can be examined with reactive transport modeling as illustrated here.

5 Conclusions

Reactive transport simulations show that the effectiveness of sand caps for remediating arsenic and mercury in subaqueous sediments, as measured by solid and aqueous-phase concentrations in the habitat layer, depends on (1) the rate of biologically mediated sulfate reduction, (2) dissolved contaminant concentrations and flux, (3) porewater advection rates, and (4) cap thickness. Arsenic attenuation depends primarily on pH-dependent sorption by iron oxide minerals and secondarily on arsenic sulfide precipitation. Therefore, prediction of cap performance with time depends on arsenic concentration and flux through the cap, the concentration and type of iron oxide mineral sorption sites, and pH. The primary controlling factor for mercury attenuation within the cap is mercury sulfide precipitation. In the habitat layer, however, attenuation of mercury and methylmercury are controlled by the amount and distribution of DOM and SOM, and in particular the concentration of thiol (–SH) groups associated with OM. Unamended caps are expected to be more effective for mercury attenuation than for arsenic attenuation in estuarine settings because of the lower solubility of mercury sulfide solid phases compared with arsenic sulfide phases. Dissolved and solid-phase contaminant concentrations are higher in cases where advective transport is important because introduction of additional contaminant mass can overwhelm attenuation mechanisms, particularly for arsenic sorption to iron oxides. The depth of mercury sulfide precipitation relative to the sediment–water interface can be partially controlled by changing cap thickness. The minimum thickness to reduce exposure to ecological receptors will be the depth of the habitat layer; however, the maximum required thickness may be as little as one meter. Site-specific quantification of the potential impact of biota on organic carbon is an important consideration during cap design, as are other physical processes such as compaction, sedimentation, bioturbation, and bioirrigation that influence the depth of the redox front below the sediment–water interface. Given the complex interplay between chemical, physical, and biological processes, reactive transport modeling provides a quantitative framework that can guide site-specific cap design, as well as lend insight into coupled physical-biogeochemical processes in subaqueous sediment systems.

Appendices

Appendix 1

See Table 7.

Table 7 Thermodynamic equilibrium constants for arsenic added to default LLNL database

Appendix 2

See Table 8.

Table 8 Thermodynamic equilibrium constants for mercury added to default LLNL database