Tableau Input Coupled Kinetic Equilibrium Transport (TICKET) Model

Dec 28, 2007 - TICKET is a general-purpose, multispecies reactive transport model that is based on the tableau structure in MINEQL. The model can be u...
4 downloads 25 Views 150KB Size
Download PDF
Environ. Sci. Technol. 2008, 42, 838–844

Tableau Input Coupled Kinetic Equilibrium Transport (TICKET) Model K E V I N J . F A R L E Y , * ,†,| K E V I N J . R A D E R , ‡ AND BENJAMIN E. MILLER§ Civil and Environmental Engineering, Manhattan College, Riverdale, New York 10471, Civil and Environmental Engineering, University of Delaware, Newark, Delaware 19716, and Civil and Environmental Engineering, University of Massachusetts Lowell, Lowell, MA, 01854

Received October 18, 2006. Revised manuscript received September 20, 2007. Accepted September 24, 2007.

TICKET is a general-purpose, multispecies reactive transport model that is based on the tableau structure in MINEQL. The model can be used in solving problems from simple chemical equilibrium calculations for batch systems to complex onedimensional, reactive transport computations for surface water, groundwater, and water treatment systems. To streamline model input and model formulation, specifications of equilibrium speciation (including homogeneous precipitation, solid solutions, adsorption, and ion exchange) and linear and nonlinear kinetic reactions are defined directly in the tableau. In addition, the burden of accounting for appearing and disappearing solid phases is circumvented by approximating homogeneous precipitation as a solid solution (with an insoluble seed). The solution technique for the model is based on a onestep algorithm and can be applied to both steady-state and fully implicit, time-variable problems. This approach is particularly well-suited in handling chemical speciation-transport problems with fast, nonlinear reaction kinetics and transient chemical intermediates. TICKET model simulations are presented for several test cases to verify the computational scheme. A model application, which examines the effects of sorption and overlying dissolved oxygen concentration on Fe(II) and As(III) oxidation in a sediment column, is also presented to demonstrate the utility of TICKET in examining complex chemical speciationtransport behavior.

Introduction Chemical speciation and species-dependent transport processes are particularly important in determining the fate and effects of metals and metalloids in the environment. A number of speciation-transport models have been developed and applied over the years to address specific problems related to metal contamination. These have included work on the transport of cations in groundwater (1–6), the subsurface transport of radionuclides (7, 8), the propagation of mineral fronts in groundwater flow systems (9, 10), the cycling of metals and major elements in sediments (11–13), * Corresponding author phone: (718) 862-7383; fax: (718) 8628035; e-mail: [email protected]. † Manhattan College, Riverdale. ‡ University of Delaware. § University of Massachusetts. | Current Address: HydroQual, Inc. Mahwah, NJ 07430. 838

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 42, NO. 3, 2008

and the transport of metals in streams (14). More recently, general modeling codes such as PHREEQC (15) and PHAST (16) have been developed, primarily to address time-variable response of chemicals in groundwater. Corresponding development of general-purpose, chemical-speciationtransport models for surface water quality studies has been limited (17). The goal of our work, therefore, is to construct a generalpurpose, multispecies reactive transport modeling approach that can be readily applied to both steady-state and timevariable evaluations of surface water systems. In this formulation, specific attention is given to settling, resuspension, and sediment boundary exchange for streams and estuaries, and to settling and particle mixing for vertically layered lakes and sediment following approaches outlined in Thomann and Mueller (18). Other transport processes are defined in general terms to allow direct application of the model to batch systems, flow-through reactors, packed columns, and one-dimensional groundwater systems. Chemical speciation in the model is based on the tableau structure in MINEQL (19). To provide a more streamlined approach for model input and model formulation, the tableau is expanded to include both chemical equilibrium and kinetic reactions. Steady-state and time-variable solutions are obtained by simultaneous consideration of equilibrium speciation, transport, and kinetic reactions using a one-step solution algorithm. This approach is particularly well-suited for problems with fast, nonlinear reaction kinetics (20). The result, the tableau input coupled kinetic equilibrium and transport (TICKET) model, is described below and is followed by a discussion of model verification and model application. Model Development. The basic approach for constructing a chemical equilibrium tableau and for developing a numerical solution are described in Morel and Hering (21) and Westall et al. (19), respectively. In this approach, chemical components (which according to Morel and Hering (21) represent “a set of chemical entities that permits a complete description of the stoichiometry of a system”) are written across the top of the tableau, and a full set of chemical species are written down the left side. The elements in the tableau represent the stoichiometric coefficients for the chemical equilibrium reactions and are used in defining both mass action and mole balance equations. In TICKET, separate mass action and mole balance coefficients are specified using a double entry (ai,j, bi,j) for each element in the tableau (see top section of tableau in Table 1). This provides greater flexibility in constructing tableaux for more complex problems such as solid-solutions and homogeneous precipitation reactions as described below. In the double entry format, the ai,j coefficients are used in calculating the chemical species concentration for each row in the tableau as follows: n

Ch,i ) Kh,i

∏ (X

ai,j h,j)

(1)

j

The notation in eq 1 is written in a slightly more general format for use in multicell, speciation-transport problems. For this purpose, h, i, and j are used as the cell index, species index, and component index, respectively, and n is the number of components in a cell. Accordingly, Ch,i is the concentration of the ith species in cell h; Kh,i is the conditional mass action chemical equilibrium constant for the ith species in cell h; and Xh,j is the activity of the jth component in cell h. For this calculation, Kh,i is determined from the intrinsic mass action chemical equilibrium coefficient and the ap10.1021/es0625071 CCC: $40.75

 2008 American Chemical Society

Published on Web 12/28/2007

TABLE 1. TICKET Input Tableau for Chemical Equilibrium and Kinetic Reactionsa Xh,1 species Ch,1 Ch,2 • • • kinetic reactions Rh,1 Rh,2 • • • a

Xh,2

•••

a1,1;b1,1 a2,1;b2,1 • • •

a1,2;b1,2 a2,2;b2,2

•••

rxn rxn a 1,1 ;b1,1 rxn rxn a 2,1 ;b2,1 • • •

rxn rxn a 1,2 ;b1,2 rxn rxn a 2,2 ;b2,2

•••

K K1 K2 • • • k k1 k2 • • •

Mq

Me

Mw

Mk′

Mq1 Mq2 • • •

Me1 Me2

Mw1 Mw2

Mk1 Mk2

h is used in multicell problems as the cell index.

propriate activity coefficients (e.g., as calculated from the Debye–Hückel or Davies equations). Information for kinetic reactions is added to the bottom section of the tableau (Table 1) by also specifying two entries in each element of the tableau. In this approach the first rxn, corresponds to the exponents in the kinetic rate entry,ai,j law, and the second entry, brxn i,j , corresponds to stoichiometric coefficients in the kinetic reaction. By convention, negative rxn denote the loss of reactants and positive values values of bi,j denote reaction products. The overall rate for each reaction can then be expressed as follows: n

Rh,i ) kh,i

∏

balance equation for the total component concentration (TOTX) for the jth component in cell h is given as follows: Vh

m

Qh-1,h +

rxn

(2)

j

where Rh,i is the overall reaction rate for the ith reaction in cell h; and kh,i is the reaction rate coefficient. In this format, separate forward and backward reactions are specified for reversible reactions that are kinetically controlled; and Michaelis–Menten reactions can be expressed as the combination of two reactions as described in Morel and Hering (21). Information on chemical equilibrium constants and kinetic rate coefficients are included in the tableau (Table 1). In addition, species mobility factors that define the relative effect of flow, diffusion/dispersion, settling/burial, and interfacial mass transfer processes on each chemical species are specified as Mqi, Mei, Mwi, Mki’, respectively. Values for Mqi, Mei, Mwi, and Mki’ are typically set equal to one for species that are transported by flow, diffusion/dispersion, settling/burial and/or interfacial mass transfer processes, and set equal to zero for immobile species. Corresponding information for transport in multicell problems is given in Table 2 for the simple case of constant transport through a one-dimensional system. This information includes cell volumes (Vh) along with parameters to describe transport in the x-direction between connecting cells h and h + 1: the volumetric flow rate (Qh,h+1), the crosssectional area (Ah,h+1), the diffusion/dispersion coefficient (Eh,h+1), the mixing length (∆Lh,h+1), and the settling/burial velocity (wh,h+1). In addition, transfer of chemical in the y or z directions (e.g., to simulate lateral inflows, flow diversions, settling, gas exchange or other boundary exchange processes in a stream) can also be added to Table 2 by specifying lateral inflow (Qbound,h), flow diversion (Qh,bound), interfacial area ′ ), (Ah,bound), interfacial mass transfer rate coefficient (kh,bound and settling/burial velocity (wh,bound) between cell h and the cell boundary. Transport information (Table 2) is combined with information on equilibrium and kinetic reactions given in the tableau (Table 1) through dynamic mole (or mass) balance equations. For a one-dimensional system, the dynamic mole

m

∑ Mq b i

Ch-1,i -Qh,h+1

i,j

i

+ (Xh,j)a i,j

d(TOTXh,j) ) dt

Eh-1,h Ah-1,h ∆Lh-1,h Eh,h+1 Ah,h+1 ∆Lh,h+1

[ [∑

∑ Mq b i

i,j

Ch,i

m

m

∑ Me b i

i,j

Ch-1,i -

i m

∑ Me b

Ch,i

∑ Me b

Ch,i

i

i,j

i m

Mei bi,j Ch+1,i -

i

i

i,j

i

m

+ wh-1,h Ah-1,h

∑ Mw b i

i,j

m

∑ Mw b

Ch-1,i - wh,h+1 Ah,h+1

i

i

i,j

bound Ch,i - Qh,bound

[∑

∑ Mq b

bound Mki′ bi,j C h,i -

i m

- wh,bound Ah,bound

∑ Mw b i

i

i

i,j

Ch,i

Ch,i

i m

m

′ + kh,bound Ah,bound

i,j

m

∑ Mq b i

i

i

m

+ Qbound,h

] ]

i

i,j

∑ Mk b ′ i

i r

Ch,i + Vh

∑b

rxn i,j

i,j

]

Ch,i

Rh,i + Sh,j (3)

i

where m is the total number of chemical species in each cell; r is the total number of reactions in each cell. For terms in eq 3, line 2 represents the transport of the total chemical component between connecting cells by flow (where flow is assumed to be in the direction of increasing cell number); lines 3 and 4 represent the diffusive/ dispersive transport of the total chemical component between connecting cells; line 5 represents the settling of the total chemical component from one cell to the next (e.g., to describe settling in a vertically layered lake or reservoir); line 6 represents additional input of the total chemical component from lateral inflows and the loss of total chemical component by flow diversions; line 7 represents the mass transfer of the total chemical component across a boundary interface (e.g., to describe gas exchange); line 8 represents the loss of the total chemical component by settling to a boundary (e.g., to describe settling from stream waters to the underlying sediment), the kinetic reactions, and direct sources of the total chemical component into the cell h, respectively. If needed, other transport processes can readily be added to the formulation. Steady-State Solution. The steady-state solution for TICKET can be obtained by setting the time derivative, d(TOTXh,j)/ dt, in eq 3 equal to zero. The remaining terms are rearranged following the approach of Thomann (22) and Thomann and Mueller (18) to obtain the following: VOL. 42, NO. 3, 2008 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

839

TABLE 2. TICKET Input Parameters for Constant, One-Dimensional Transporta cell 0 1 2 • • • a

volume (h)

flow (h,h+1)

area (h,h+1)

disp. coef. (h,h+1)

mixing length (h,h+1)

sett./ burial (h,h+1)

boundary inflow (bound,h)

boundary outflow (h,bound)

boundary area (h,bound)

boundary transfer (h,bound)

boundary sett./burial (h,bound)

V1 V2 • • •

Q0,1 Q1,2 Q2,3 • • •

A0,1 A1,2 A2,3 • • •

E0,1 E1,2 E2,3 • • •

∆L0,1 ∆L1,2 ∆L2,3 • • •

w0,1 w1,2 w2,3 • • •

Qbound,1 Qbound,2 • • •

Q1,bound Q2,bound • • •

A1,bound A2,bound • • •

k′1,bound k′2,bound • • •

w1,bound w2,bound • • •

Time-variable and multi-dimensional transport can be specified but they require a more complicated input format.

m

inputsh,j ) -

∑

m

m

upstr f h,i bi,j Ch-1,i +

i

∑

cell f h,i bi,j Ch,i

inputsh,j ) -

∑ i

i

-

∑f

r

down h,i

bi,j Ch+1,i - Vh

i

∑

rxn bi,j

+ Rh,i

(4)

m

∑ [Q

′ bound Mqi + kh,bound Mki′ ] C h,i

bound,h

i

upstr f h,i )

cell f h,i )

down f h,i )

[

Eh-1,h Ah-1,h Mei ∆Lh-1,h +wh-1,h Ah-1,h Mwi

Qh-1,h Mqi +

[

Qh,h+1 Mqi + Qh,bound Mqi +

+

[

]

Eh-1,h Ah-1,h Mei ∆Lh-1,h

Eh,h+1 Ah,h+1 Mei + wh,h+1 Ah,h+1 Mwi ∆Lh,h+1

′ +wh,bound Ah,bound Mwi + kh,bound Mki′

]

cell h,i

i

]

Eh,h+1 Ah,h+1 Mei ∆Lh,h+1

-

jj

down h,i

jj

r

∑b

rxn i,j

i

n

m

TOTXh,j )

∑b

i,j

Ch,i

(6)

i

Final solution of the steady-state, chemical speciationtransport equations is obtained by substituting the chemical equilibrium mass action equations (eq 1) and the overall reaction rate laws (eq 2) into the dynamic mole balance equation (eq 4). The resulting equations are as follows: 840

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 42, NO. 3, 2008

rxn

(Xh,jj)ai,jj

kh,i

(7)

jj

These are then solved simultaneously for each component in each cell, Xh,j, using a Newton–Raphson iteration technique. Steady-state species concentrations, Ch,i, are then obtained directly from the chemical equilibrium mass action equations (eq 1). The method is identical to the solution in MINEQL and is outlined in detail in Westall et al. (19). Additional information that is needed to solve steady-state problems includes specification of direct sources (Sh,j) and boundary bound concentrations (C h,i ) for mobile and/or reacting components; total component concentrations (TOTXh,j) for immobile, nonreacting components; activities of fixed components (Xh,j); and initial activity guesses for “non-fixed” components (Xh,j)o. Time-Variability Solution. A time-variable solution for eq 3 can be determined in a similar fashion. Using a first-order, implicit solution scheme, eq 3 can be rearranged and expressed as follows:

+

m

Vh ∆t

∑b

i,j

t-∆t Ch,i )

Vh ∆t

∑b

i,j

t Ch,i -

i m

m

i

∑f

+

∑f

cell h,i

upstr h,i

t bi,j Ch-1,i

i m

m

For problems where immobile, nonreactive components are present (e.g., immobile sorbent phases or electrostatic corrections for surface complexes), the dynamic mole balance equation (eq 4) is replaced by the static mole balance equation that is used in chemical speciation programs like MINEQL:

(Xh+1,jj)ai,jj

bi,j Kh,i

i

inputsh,j + (5)

n

∑f

- Vh

(Xh,jj)ai,jj

bi,j Kh,i

m

In this equation, inputsh,j represents all external inputs of the jth component into cell h from direct sources and other upstr is associated with the transport flux boundary inputs; f h,i cell th of the i species from the upstream cell into cell h; f h,i is th species out of cell associated with the transport flux of the i down h; and f h,i is associated with the transport flux of the ith species from the downstream cell into cell h as expressed in eq 5:

inputsh,j ) Sh,j +

n

∑f

i

(Xh-1,jj)ai,jj

jj

m

m

( ∑ ) ( ∑ ) ( ∑ ) ( ∑ ) n

upstr f h,i bi,j Kh,i

t bi,j Ch,i -

i

∑f

down h,i

t bi,j Ch+1,i

i

r

-Vh

∑b

rxn i,j

t Rh,i

(8)

i

where Ct-∆t denotes known concentrations at the beginning of a time step, and Ct, Rt denote unknown concentrations and overall reaction rates at the end of a time step. The chemical equilibrium mass action equations (eq 1) and the overall reaction rate laws (eq 2) are then substituted into the mass transport equation (eq 8). The resulting set of equations are solved simultaneously at each time step using the Newton– Raphson solution. In addition to information cited above for steady-state applications, specification of the initial total component concentrations for mobile and/or reacting components (TOTXh,j)o is also required for time-variable problems.

TABLE 3. TICKET Input Tableau for a Solid Solution of ML(s) and M’L(s)a M species M M′ L ML(s) {ML(s)} M’L(s) {M’L(s)} a

M′

L

TS

1,1 1,1 1,1 1,0 1,1 1,0

1,1 1,1 1,0 1,1 1,0

1,0 0,1 1,0 0,1

K 0 0 0 1/Ksp1 1/Ksp1 1/Ksp2 1/Ksp2

Mq

Me

Mw

Mk′

1 1 1 0 0 0 0

1 1 1 0 0 0 0

0 0 0 0 0 0 0

0 0 0 0 0 0 0

TS is the total moles of solid in the solid solution.

Special Case: Solid Solutions and Homogeneous Precipitation. For tableau-based solutions, the concentration of a species can be determined directly from the system components in the tableau using eq 1. Solid solutions, however, present a special problem. To demonstrate this point, consider two metals, M and M′, precipitating with ligand L to form a solid solution of ML(s) and M’L(s). For this case, the molar concentration of ML(s), which would be given as {M} × {L} × (1/Ksp) × ([ML(s)] + [M’L(s)]) can not be represented directly in the tableau structure. Solid solutions however can be handled by modifying the tableau following the approach of Farley et al. (23). As shown Table 3, a new component, TS, is introduced in the tableau to represent the total molar concentration of solid in the solid solution (i.e., TS ) [ML(s)] + [M’L(s)]). In addition, both activities and molar concentrations of each solid are considered in the list of species. Mass action coefficients, ai,j, are then specified for the activity of each solid based on a rearrangement of the solubility product relationship {ML(s)} ) {M} {L}

1 Ksp

(9)

Mass action coefficients, ai,j, for the molar concentration of each solid are determined by specifying the activity of the solid as its mole fraction (i.e.,{ML(s)} ) [ML(s)] ⁄ TS ) which after rearranging the equation and substituting the solubility product relationship (eq 9) yields [ML(s)] ) {ML(s)}

TS ) {M} {L}

1 TS Ksp

(10)

For molar concentrations of the solids, ML(s) and M’L(s), the mole balance coefficients, bi,j, are set equal to the mass action coefficients, ai,j, in all but the “TS” column to properly account for the amount of the solids in the metal and ligand mole balance equations. The mole balance coefficients, bi,j, for the activities of the solids, {ML(s)} and {M’L(s)}, are set equal to “1” in the “TS” column. Lastly, “Total TS” which now represents the sum of activities of solids in the solid solution is, according to solid solution theory, set equal to “1”. This basic approach for solid–solution precipitation can be extended to homogeneous precipitation by specifying the second solid as an imaginary, highly insoluble seed. The error introduced in this approximation is on the order of the concentration of the seed and can be controlled by setting the seed concentration at an extremely low value. By treating homogeneous precipitates in this way, we avoid the extra burden of accounting for appearing and disappearing solid phases that can potentially occur during TICKET model simulations. Solution Algorithm. In many chemical equilibriumtransport models, time-variable solutions are based on a twostep or split operator scheme (7, 17, 20). In this approach,

the chemical equilibrium calculation is typically performed at each time-step, and resulting species concentrations are then used in solving the kinetic and mass transport equations (eq 3). The advantage of split-operator schemes is that the algorithms are computationally faster for a given time step and have lower computational storage requirements (7). This is in contrast to the “one-step” solution scheme in TICKET where chemical equilibrium, transport, and kinetics are solved simultaneously (e.g., using a Newton–Raphson iteration solver). The advantage of the one-step method is that it is directly applicable to both steady-state and time-variable solutions. In addition, the one-step approach provides a more robust solution, particularly for time-variable problems with fast, nonlinear reaction kinetics (e.g., redox problems and/ or problems with short-lived reaction intermediates). For this class of problems, the one-step method allows timevariable simulations to be performed with larger (and fewer) time steps for improved overall computational performance (20). Initial development and testing of TICKET was performed in FORTRAN90 (24). Since that time, the model has been revised to include a more complete description of mass transport processes, a more straightforward solution technique, and a more refined scheme for including solid solutions and homogeneous precipitation. The current version of the model is available in Microsoft Excel with a programmed macro written in Visual Basic for Applications (VBA). A preprocessor VBA macro is also available to help generate the input tableau from chemical equilibrium and kinetic databases. The Excel/VBA version of TICKET has also been translated to Microsoft Visual Basic (VB) to take advantage of improved computational speed of the compiled VB code. Model Verification and Application. TICKET model simulations were performed for several case studies to verify the computational scheme and to highlight the general applicability and utility of the model. These included verification studies where TICKET model results were compared directly to (1) MINEQL+ (25) batch chemical equilibrium model results; (2) analytical solutions for simple one-dimensional transport; (3) field and model results for subsurface transport of major ions affected by ion exchange reactions (1), and (4) ACUCHEM (26) batch kinetic model results for fast, nonlinear reaction kinetics that have been proposed for Fe(II)-catalyzed oxidation of As(III) (27). The latter two cases are presented in the Supporting Information to demonstrate the general applicability of TICKET in handling chemical equilibrium-transport problems and fast, nonlinear kinetic problems, respectively. An additional verification study is presented below for the propagation of calcite and dolomite fronts in a packed column. This case is presented to elaborate on the use of our solid solution formulation in approximating homogeneous precipitation. Finally, an application study is presented to examine the effects of sorption and overlying dissolved oxygen concentrations on the transient response for iron and arsenic oxidation in sediment (28). Verification Study. Calcite and Dolomite Precipitation in a Packed Column. This problem was designed by Engesgaard and Kipp (10) to test MST1D, their one-dimensional geochemicaltransportmodel,againsttheprogramCHEMTRNS by Noorishad et al. (9). For this analysis, a 50 cm column of saturated porous media with a porosity of 0.32 and a bulk density of 1800 kg/m3 initially contains 21.7 µmoles of calcite per kg of porous media (i.e., 39 mmoles of calcite per m3 of bulk volume) that is evenly distributed throughout the column. A step input of 1 mM MgCl2 solution is injected into the column. The pore water velocity is given as 9.37 × 10-6 m/sec and the longitudinal dispersivity is given as 0.0067 m. VOL. 42, NO. 3, 2008 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

841

FIGURE 2. Fe-catalyzed As(III) oxidation mechanism circumneutral pH systems. Based on Hug and Leupin (27).

FIGURE 1. TICKET-calculated results (solid lines) and MST1D model results (dashed lines) for spatial distributions of (a) dolomite and calcite, (b) Mg2+ and Ca2+, and (c) pH, HCO3– , and CO32– concentrations at 21000 s. Based on the numerical study of Engesgaard and Kipp (10). The problem was modeled using 50, 1 cm cells. The time step was set at 200 s and the total simulation time was 21000 s. In this TICKET run, calcite(s) and dolomite(s) are modeled as equilibrium, homogeneous solids by specifying each solid as part of a solid solution with a highly insoluble seed solid, S_calcite(s) and S_dolomite(s), respectively. Following the approach outlined in Table 3, new components, TS_calcite and TS_dolomite, are introduced in the tableau to represent the total molar concentration of solid in the calcite and dolomite solid solutions (i.e., TS_calcite ) [calcite(s)] + [S_calcite(s)] and TS_dolomite ) [dolomite(s)] + [S_dolomite(s)], respectively). Both activities and molar concentrations of the solids and their seeds are considered in the list of species and entries in the tableau are made according to our solid solution approach (see Table 3). For the seed solids, K values are set equal to 1 (or equivalently, log K ) 0) and the total seed concentrations, TOT S_calcite and TOT S_dolomite, are set to small values (e.g., 10-20 mol/L) to control the accuracy of the approximation for homogeneous precipitation. Activity corrections for the solution species are made throughout the simulation using the Davies equation (29) at the beginning of each iteration, by first calculating ionic strength using species concentrations from the previous Newton–Raphson solution. The input tableau for the problem is given in Supporting Information Table S3. The TICKET and MST1D model results for dolomite and calcite (Figure 1a), Mg2+ and Ca2+ (Figure 1b) and pH, HCO–3, and CO32– (Figure 1c) are presented as spatial distributions 842

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 42, NO. 3, 2008

for

at 21000 s. TICKET results closely match previously published numerical modeling results from MST1D (10), with the possible exception of HCO–3. Since TICKET and MST1D responses for pH and CO32– coincide, the discrepancy in the HCO–3 response is most likely due to small differences in the HCO–3 acidity constant or small differences in the ionic strength activity adjustments in the two models. (Additional batch equilibrium calculations were performed for the calcite system using TICKET and MINEQL+ (25) to confirm that TICKET was properly handling ionic strength activity adjustments for the HCO–3 calculation.) A slight offset in the calcitedolomite front is also noted and is attributed to differences in the TICKET and MST1D transport algorithms. Model Application. Fe(II)-Catalyzed As(III) Oxidation Column Study. A chemical pathway for Fe(II)-catalyzed oxidation of As(III) is described in detail by Hug and Leupin (27). As shown in Figure 2, the sequence begins with a series of one-electron transfers as Fe(II) is oxidized to a radical species, Fe(IV). The Fe(IV) is then utilized in oxidizing As(III) to the intermediate As(IV) which is subsequently oxidized to As(V). Because Fe(IV) is also consumed in a competitive reaction with Fe(II), only a small amount of As(III) was oxidized in their batch studies (27). Bisceglia et al. (28) reasoned that in sediment, where continual oxidation of Fe(II) occurs, As(III) oxidation could potentially proceed to completion. This effect was demonstrated in a sediment column experiment which showed 28% As(III) oxidation in 17 days (compared to 8% As(III) oxidation in an equivalent batch system). TICKET model simulations were also performed. For final model simulations, all kinetic reactions were taken from Hug and Leupin (27) and equilibrium reactions for the sorption of Fe(II) and As(III) to sand and for the sorption of As(V) and As(III) to precipitated Fe(III) oxide were determined through model calibration (28). As a follow-up to this work, model simulations were performed to further examine the effects of sorption and overlying dissolved oxygen concentration on the transient response for iron and arsenic in sediment. All model parameters and coefficients were taken directly from Bisceglia et al. (28) and are described in the Supporting Information. For the base case, sorption reactions were not considered and the overlying dissolved oxygen concentration was set at 250 µmol/L. Simulation results for day 10 (Figure 3a and b) show that oxygen has diffused from the overlying water into the sediment. Concomitantly, Fe(II) and As(III) were oxidized, producing a gradient in Fe(II) and As(III) concentrations and inducing diffusion of Fe(II) and As(III) to the oxic-anoxic interface. The combination of diffusive transport and oxida-

FIGURE 3. TICKET-calculated depth profiles of total Fe(II), total Fe(III), dissolved oxygen, total As(III), and total As(V) for day 10: (a, b) base case with no sorption and an overlying dissolved oxygen concentration of 250 µmol/L, (c, d) Fe(II) sorption to sand, and (e, f) overlying dissolved oxygen concentration of 125 µmol/L, and (g, h) As(III) sorption to sand. Concentrations are given in mmoles L(pore water)-1. Based on the sediment column experiment of Bisceglia et al. (28). tion processes had resulted in a build-up of an immobile Fe(III) precipitate in the oxic layer, and a more diffuse distribution of the mobile As(V) in the sediment. Overall, 46% Fe(II) and 40% As(III) were oxidized in the 10 day period. This trend continued such that by day 60, 100% Fe(II) and 95% As(III) were oxidized. The effect of sorption was then considered by first specifying Fe(II) sorption onto sand (with an equivalent partitioning of 9:1 between sorbed and dissolved phases). In this case (Figure 3c, d), sorption reduced Fe(II) concentration in the pore water and retarded Fe(II) diffusive transport by a factor of 10. The overall rate of Fe(II) oxidation was, therefore, slower (29% Fe(II) oxidation by day 10). During this time, oxygen penetrated farther into the sediment and a much smaller build-up of Fe(III) oxide was observed in the oxic layer. This occurred because the diffusive transport of Fe(II) up to the oxic–anoxic interface is now much slower than the downward diffusion of oxygen. Although Fe(II) sorption has a significant effect on Fe(II) oxidation in sediment, the overall rate of As(III) oxidation was similar to the base case with 38% As(III) oxidation by day 10. This result at first seems counterintuitive because the slower rate of Fe(II) diffusion and Fe(IV) production at the oxic–anoxic interface should slow the As(III) oxidation. Lower dissolved Fe(II) concentrations at the oxic–anoxic interface however decreased the Fe(II) competitive pathway (Figure 2) and led to more efficient oxidation of As(III). Ultimately 100% Fe(II) and 95% As(III) were oxidized by day 110. In similar fashion, the effect of oxygen transport was examined by reducing the overlying dissolved oxygen concentration by a factor of 2. In this case (Figure 3e, f), oxygen penetration was reduced, Fe(II) and As(III) oxidation occurred closer to the water–sediment interface, and a larger build-up of immobile Fe(III) oxide was observed in near surface sediment. As expected, overall rates of Fe(II) and As(III) oxidation were slower (41% Fe(II) and 34% As(III) oxidation by day 10). This occurred because oxygen limited the rate of Fe(IV) production in sediment and affected the Fe(II) and As(III) oxidation rates proportionally. Ultimately, 100% Fe(II) and 95% As(III) oxidation was achieved in 100 days. Finally, for different water chemistry (e.g., lower pH), As(III) sorption may become dominant. Effect of As(III)

sorption was, therefore, examined by specifying As(III) binding onto sand (with an equivalent partitioning of 9:1 between sorbed and dissolved phases). In this case (Figure 3g, h), sorption reduced As(III) concentration in the pore water and retarded As(III) diffusive transport by a factor of 10. This had a minor effect of the Fe(II) oxidation rate (47% Fe(II) oxidation by day 10), but slowed the rate of As(III) oxidation significantly (19% As(III) by day 10). Ultimately, only 61% As(III) oxidation was attained (with 100% Fe(II) oxidation by day 60). The Fe(II)-catalyzed As(III) oxidation column application demonstrates the utility of TICKET in examining complex chemical speciation–transport behavior. Results presented above indicate that Fe(II) sorption, overlying dissolved oxidation concentration and As(III) sorption (alone or in various combinations) will affect the vertical distribution of Fe in sediment and the rate and ultimate extent of As(III) oxidation. In other applications, TICKET is being used in assessing the effects of site-specific water chemistry and surface water transport processes on steady-state and timevariable responses of metals in lakes. This work, which is briefly described in the Supporting Information, demonstrates the importance of metal binding affinity to organic matter and organic matter cycling on the behavior of Cu and Cd in an idealized lake. Lastly, since the modeling framework (as given in eqs 1–10) is relatively straightforward, the TICKET approach can readily be used by others who wish to develop more advanced speciation-transport models (e.g., which include linked three-dimensional hydrodynamic, sediment transport and/or nutrient enrichment/organic carbon calculations).

Acknowledgments This research was funded in part by the National Institute of Environmental Health Sciences through a Superfund Basic Research Program Grant (P42ES10344), the U.S. EPA Center for Metals in the Environment, and the International Council for Mining and Metals (ICMM). Additional support for Benjamin Miller was provided through a National Science Foundation Graduate Research Fellowship. The comments and suggestions of Dominic M. Di Toro (University of Delaware) and Richard F. Carbonaro (Manhattan College) are greatly appreciated. VOL. 42, NO. 3, 2008 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

843

Supporting Information Available Verification study for subsurface transport of major cations: input tableau (Table S1), model results (Figure S1); verification study for Fe(II)–As(III) oxidation in a batch system: input tableau (Table S2), model results (Figure S2); input tableau for calcite and dolomite precipitation in a packed column (Table S3); parameters for Fe(II)-catalyzed As(III) oxidation column model application: input tableau (Table S4); model study for metal cycling in an idealized lake: model contruct (Figure S3), organic matter and sulfide cycling (Figure S4), model results for Cu and Cd (Figure S5). This material is available free of charge via the Internet at http://pubs.acs.org.

Literature Cited (1) Valocchi, A. J.; Street, R. L.; Roberts, P. V. Transport of ionexchanging solutes in groundwater: chromatographic theory and field simulation. Water Resour. Res. 1981, 17, 1517–1527. (2) Miller, C. W.; Benson, L. V. Simulation of solute transport in a chemically reactive heterogeneous system. Water Resour. Res. 1983, 19, 381–391. (3) Cederberg, G. A.; Street, R. L.; Leckie, J. O. Groundwater mass transport and equilibrium chemistry model for multicomponent systems. Water Resour. Res. 1985, 21, 1095–1104. (4) Van Ommen, H. C. ‘Mixing-cell’ concept applied to transport of non-reactive and reactive components in soils and groundwater. J. Hydrol. 1985, 78, 201–213. (5) Charbeneau, R. J. Multicomponent exchange and subsurface solute transport: characteristics, coherence, and the Riemann problem. Water Resour. Res. 1988, 24, 57–64. (6) Appelo, C.A.J. Some calculations on multicomponent transport with cation exchange in aquifers. Ground Water 1994, 32, 968– 975. (7) Yeh, G. T.; Tripathi, V. S. Critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water. Resour. Res. 1989, 25, 93– 108. (8) Yeh, G. T.; Tripathi, V. S. Model for simulating transport of reactive multispecies components: model development and demonstration. Water. Resour. Res. 1991, 27, 3075–3094. (9) Noorishad, J.; Carnahan, C. L.; Benson, L. V. Development of the Nonequilibrium Reactive Chemical Transport Code CHEMTRNS,Rep. LBL-22361;Lawrence Berkeley Laboratory, University of California: Berkeley, 1987. (10) Engesgaard, P.; Kipp, K. L. Geochemical transport model for redox-controlled movement of mineral fronts in groundwater flow systems: a case of nitrate removal by oxidation of pyrite. Water Resour. Res. 1992, 28, 2829–2843. (11) Van Cappellen, P. and Wang, Y. Metal cycling in surface sediments: modeling the interplay of transport and reaction. In Metal Contaminated Aquatic Sediments; Allen, H. E., Ed.; Ann Arbor Press: Chelsea, MI, 1995; pp 21-64. (12) Van Cappellen, P.; Wang, Y. Cycling of iron and manganese in surface sediments: a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese. Am. J. Sci. 1996, 296, 197–243.

844

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 42, NO. 3, 2008

(13) Giambalvo, E. R.; Steefel, C. I.; Fisher, A. T.; Rosenberg, N. D.; Wheat, C. G. Effect of fluid-sediment reaction on hydrothermal fluxes of major elements, eastern flank of the Juan de Fuca Ridge. Geochim. Cosmochim. Acta 2002, 66, 1739–1757. (14) Runkel, R. L.; Kimball, B. A. Evaluating remedial alternatives for an acid mine drainage stream: application of a reactive transport model. Environ. Sci. Technol. 2002, 36, 1093–1101. (15) Parkhurst, D. L.; Appelo, C.A.J. User’s guide to PHREEQC (Version 2)––A computer program for speciation, batch-reaction, onedimensional transport, and inverse geochemical calculations,Water-Resources Investigations Report 99–4259; U.S. Geological Survey: Reston, VA, 1999. (16) Parkhurst, D. L.; Kipp, K. L.; Engesgaard, P.; Charlton, S. R. PHAST––A program for simulating ground-water flow, solute transport, and multicomponent geochemical reactions,Techniques and Methods 6-A8;U.S. Geological Survey: Reston, VA, 2004. (17) Dzombak, D. A.; Ali, M. A. Hydrochemical modeling of metal fate and transport in freshwater environments. Water Pollut. Res. J. Canada 1993, 28, 7–50. (18) Thomann, R. V.; Mueller, J. A. Principles of Surface Water Quality Modeling and Control. Harper & Row, NY NY, 1987. (19) Westall, J. C., Zachary, J. L. and Morel, F.M.M. MINEQL: a computer program for the calculation of chemical equilibrium composition of aqueous systems. Ralph M. Parsons Laboratory Technical Note No. 18, MIT, Cambridge MA, 1976. (20) Saaltink, M. W.; Carrera, J.; AyoraC. On the behavior of approaches to simulate reactive transport. J. Contam. Hydrol. 2001, 48, 213–235. (21) Morel, F.M.M.; Hering, J. G. Principles and Applications of Aquatic Chemistry. John Wiley & Sons: New York, 1993. (22) Thomann, R. V. Systems Analysis and Water Quality Management; McGraw-Hill: New York, 1972. (23) Farley, K. J.; Dzombak, D. A.; Morel, F.M.M. Surface precipitation model for the sorption of cations on metal oxides. J. Colloid Interface Sci. 1985, 106, 226–242. (24) Miller, B. E. Development of a General Aquatic Multi-component Reactive Transport Computer Model, with Application to a Wetland Sediment. PhD Dissertation, Clemson University, Clemson, SC, 1997. (25) Schecher, W. D.; McAvoy, D. C. MINEQL+: A chemical equilibrium modeling system, Version 4.0 for Windows; Environmental Research Software: Hallowell ME, 1998. (26) Braun, W.; Herron, J. T.; Kahaner, D. K. Acuchem: Computer program for modeling complex chemical reaction systems. Int. J. Chem. Kinet. 1988, 20, 51–62. (27) Hug, S. J.; Leupin, O. Iron-catalyzed oxidation of arsenic(III) by oxygen and by hydrogen peroxide: pH-dependent formation of oxidants in the Fenton reaction. Environ. Sci. Technol. 2003, 37, 2734–2742. (28) Bisceglia, K. J.; Rader, K. J.; Carbonaro, R. F.; Farley, K. J.; Mahony, J. D.; Di Toro, D. M. Iron(II)-catalyzed oxidation of arsenic(III) in a sediment column. Environ. Sci. Technol. 2005, 39, 9217– 9222. (29) Stumm, W.; Morgan, J. J. Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters ,3rd ed.; John Wiley & Sons: New York, 1996.

ES0625071