A test problem for advectivedispersivereactive transport was developed by TebesStevens and Valocchi (1997) and TebesStevens and others (1998). Although based on relatively simple speciation chemistry, the solution to the problem demonstrates several interacting chemical processes that are common to many environmental problems: bacterially mediated degradation of an organic substrate; bacterial cell growth and decay; metal sorption; and aqueous speciation, including metalligand complexation. In this example, the test problem is solved with PHREEQC, which produces results almost identical to those of TebesStevens and Valocchi (1997) and TebesStevens and others (1998).
The test problem models the transport processes when a pulse of water containing NTA (nitrylotriacetate) and cobalt is injected into a column. The problem includes advection and dispersion in the column, aqueous equilibrium reactions, and kinetic reactions for NTA degradation, growth of biomass, and cobalt sorption.
The dimensions and hydraulic properties of the column are given in table 38.
TebesStevens and Valocchi (1997) defined an aqueous model to be used for this test problem that includes the identity of the aqueous species and log K s of the species; activity coefficients were assumed to be 1.0. The database file in table 39 was constructed on the basis of their aqueous model. For the PHREEQC simulation, NTA was defined as a new “element” in the SOLUTION_MASTER_SPECIES data block named “Nta”. From this point on “NTA” will be referred to as “Nta” for consistency with the PHREEQC notation. The gram formula weight of Nta in SOLUTION_MASTER_SPECIES is immaterial if input units are moles in the SOLUTION data block, and is simply set to 1. The aqueous complexes of Nta are defined in the SOLUTION_SPECIES data block. Note that the activity coefficients of all aqueous species are defined with a large value for the a parameter (1×10^{ 7} ) in the gamma identifier, which forces the activity coefficients to be very nearly 1.0.
The background concentrations in the column are listed in table 40. The column contains no Nta or cobalt initially, but has a biomass of 1.36×10^{ 4} g/L. A flux boundary condition is applied at the inlet of the column, and for the first 20 h (hours), a solution with Nta and cobalt enters the column; the concentrations in the pulse also are given in table 40. After 20 h, the background solution is introduced at the inlet until the experiment ends after 75 h. Na and Cl were not in the original problem definition but were added for charge balancing sorption reactions for PHREEQC (see Sorption Reactions).
Nta is assumed to degrade in the presence of biomass and oxygen by the reaction:
HNta^{2} + 1.62O_{2} + 1.272H_{2} O + 2.424H^{+} = 0.576C_{5} H_{7} O_{2} N + 3.12H_{2} CO_{3} + 0.424NH_{4} ^{+} .
PHREEQC requires kinetic reactants to be defined solely by the moles of each element that enter or leave the solution because of the reaction. Furthermore, the reactants should be charge balanced (no net charge should enter or leave the solution). The Nta reaction converts 1 mol HNta^{ 2} (C_{ 6} H_{ 7} O_{ 6} N) to 0.576 mol C_{ 5} H_{ 7} O_{ 2} N, where the latter is chemically inert, and its concentration can be ignored. The difference in elemental mass contained in these two reactants provides the stoichiometry of the elements C, H, O, and N in the reaction. This stoichiometry is equal to the sum of the elements on the righthand side of the equation, excluding C_{ 5} H_{ 7} O_{ 2} N, minus the sum of the elements on the lefthand side of the equation. The corresponding change in aqueous element concentrations per mole of HNta^{ 2} reaction is given in table 41 (positive coefficients indicate an increase in aqueous concentration, and negative coefficients indicate a decrease in aqueous concentration).
 

 

 
1.0 

3.12 

1.968 

4.848 

0.424 
The following multiplicative Monod rate expression is used to describe the rate of Nta degradation:
where is the rate of HNta^{ 2} degradation (mol L^{ 1} h^{ 1} , mole per liter per hour), is the maximum specific rate of substrate utilization (mol/g cells/h), is the biomass (g L1h1, gram per liter per hour), is the halfsaturation constant for the substrate Nta (mol/L), is the halfsaturation constant for the electron acceptor O_{ 2} (mol/L), and c indicates concentration (mol/L). The rate of biomass production is dependent on the rate of substrate utilization and a firstorder decay rate for the biomass:
where is the rate of cell growth (g L1h1), Y is the microbial yield coefficient (g cells/mol Nta), and b is the firstorder biomass decay coefficient (h^{ 1} ). The parameter values for these equations are listed in table 42.
where i is either Co^{ 2+} or CoNta^{ } (mol/L), s _{ i} is the sorbed concentration (mol/g sediment), is the mass transfer coefficient (h^{ 1} ), and is the distribution coefficient (L/g, liter per gram). The values of the coefficients are given in table 43. The values of K _{ d} were defined to give retardation coefficients of 20 and 3 for Co^{ 2+} and CoNta^{ } , respectively. Because the sorption reactions are defined to be kinetic, the initial moles of these reactants and the rates of reaction are defined with KINETICS and RATES data blocks; no surface definitions (SURFACE, SURFACE_MASTER_SPECIES, or SURFACE_SPECIES) are needed. Furthermore, all kinetic reactants are immobile, so that the sorbed species are not transported.
1 h^{ 1} 
5.07e3 L/g 

1 h^{ 1} 
5.33e4 L/g 
When modeling with PHREEQC, kinetic reactants must be charge balanced. For sorption of Co^{ 2+} and CoNta^{ } , 1 mmol of NaCl was added to the solution definitions to have counter ions for the sorption process. The kinetic sorption reactions were then defined to remove or introduce (depending on the sign of the mole transfer) CoCl_{ 2} and NaCoNta, which are charge balanced. To convert from moles sorbed per gram of sediment ( s _{ i} ) to moles sorbed per liter of water, it is necessary to multiply by the grams of sediment per liter of water, 3.75×10^{ 3} g/L.
See Input file for example 15. shows the input file derived from the preceding problem definition. Although rates have been given in units of mol L1h1, rates in PHREEQC are always mol/s (mole per second), and all rates have been adjusted to seconds in the definition of rate expressions in the input file. The density of water is assumed to be 1 kg/L.
The 10meter column was discretized with 10 cells of 1 meter each. The first two SOLUTION data blocks (table 44) define the infilling solution and the initial solution in cells 1 through 10. The solutions are copied to solution 100 and 101, to be used later in the 20cell model.
The RATES data block defines the rate expressions for four kinetic reactions: HNta2, Biomass, Co_sorption, and CoNta_sorption. The rate expressions are initiated with start , defined with numbered Basiclanguage statements, and terminated with end . The last statement of each expression is SAVE followed by a variable name. This variable is the number of moles of reaction over the time subinterval and is calculated from an instantaneous rate (mol/s) times the length of the time subinterval (s), which is given by the variable “TIME”. Lines 30 and 20 in the first and second rate expressions and line 10 in the third and fourth rate expressions adjust parameters to units of seconds from units of hours. The function “MOL” returns the concentration of a species (mol/kgw), the function “M” returns the moles of the reactant for which the rate expression is being calculated, and “KIN” returns the moles of the specified kinetic reactant. The functions “PUT” and “GET” are used to save and retrieve a term that is common to both the HNta2 and Biomass rate expressions (see also See ReactionPath Calculations).
The KINETICS data block defines the names of the rate expressions that apply to each cell; cells 1 through 10 are defined simultaneously in this example. For each rate expression that applies to a cell, the formula of the reactant ( formula ) and the moles of the reactant initially present ( m , if needed to be different from the default of 1 mol) are defined. It is also possible to define a tolerance ( tol ), in moles, for the accuracy of the numerical integration for a rate expression. Note that the HNta2 rate expression generates a negative rate, so that elements with positive coefficients in the formula are removed from solution and negative coefficients add elements to solution. The biomass reaction adds “H 0.0”, or zero moles of hydrogen; in other words, it does not add or remove anything from solution. The assimilation of carbon and nutrients that is associated with biomass growth is ignored in this simulation. Also, the kinetics block is copied to KINETICS 101, for use in the 20cell model.
The SELECTED_OUTPUT data block punches the molalities of the aqueous species Nta3, CoNta, HNta2 and Co+2 to the file ex15.sel . To each line in the file, the USER_PUNCH data block appends the time (in hours), the sorbed concentrations converted to mol/g sediment, and the biomass. USER_GRAPH plots the concentrations as symbols without lines to facilitate the comparison with the 20cell model.
The first TRANSPORT data block defines the first 20 h of the experiment, during which Nta and cobalt are added at the column inlet. The column is defined to have 10 cells ( cells ) of length 1 m ( lengths ). The duration of the advectivedispersive transport simulation is 20 time steps ( shifts ) of 3,600 seconds ( time_step ). The direction of flow is forward ( flow_direction ). Each end of the column is defined to have a flux boundary condition ( boundary_conditions ). The dispersivity is 0.05 m ( dispersivities ) and the diffusion coefficient is set to zero ( diffusion_coefficient ). Data are written to the selectedoutput file only for cell 10 ( punch_cells ) after each shift ( punch_frequency ), and data are written to the output file only for cell 10 ( print_cells ) after each fifth shift ( print_frequency ).
At the end of the first advectivedispersive transport simulation, the initial column solution, which was stored as SOLUTION 101, is copied to SOLUTION 0, to become the influent for the next transport simulation. The second TRANSPORT data block defines the final 55 h of the experiment, during which Nta and cobalt are not present in the infilling solution. All parameters are the same as in the previous TRANSPORT data block; only the number of transport steps ( shifts ) is increased to 55.
With advectivedispersivereactive transport simulations, it is always necessary to check the numerical accuracy of the results. In general, analytical solutions will not be available for these complex simulations, so the only test of numerical accuracy is to refine the grid and time step, rerun the simulation, and compare the results. If simulations on two different grids give similar results, there is some assurance that the numerical errors are relatively small. If simulations on two different grids give significantly different results, the grid must be refined again and the process repeated. Unfortunately, doubling the grid size at least quadruples the number of solution calculations that must be made because the number of cells doubles and the time step is halved. If the cell size approaches the size of the dispersivity, it may require even more solution calculations because the number of mix steps in the dispersion calculation will increase as well.
To test grid convergence in this example, the number of cells in the column were doubled. All keyword data blocks that defined compositions for the range 1 through 10 were changed to 1 through 20. In addition, the parameters for advectivedispersive transport were adjusted to be consistent with the new number of cells. The final TRANSPORT data block in table 44 defines the 20cell model. The number of cells and number of shifts are doubled; the cell length and time step are halved. To print information for the same location as the 10cell model (the end of the column), the punch_cells and print_cells are set to cell 20. To print information at the same time in the simulation as the 10cell model, punch_frequency is set to every 2 shifts, print_frequency is set to every 10 shifts, and the time step for going from the cellmidpoint to the columnend is halved on line 10 in USER_PUNCH.
The distributions of aqueous and immobile constituents in the column at the end of 75 h are shown in figures See Dissolved concentrations and pH values at the outlet of the column for Nta and cobalt transport simulations with 10 (symbols) and 20 cells (lines). and See Concentrations of sorbed species and biomass at the outlet of the column for Nta and cobalt transport simulations with 10 (symbols) and 20 cells (lines). for the 10 and 20cell models. In the experiment, two pore volumes of water with Nta and cobalt were introduced to the column over the first 20 h and then followed by 5.5 pore volumes of background water over the next 55 h. At 10 h, HNta^{ 2} begins to appear at the column outlet along with a rise in the pH (fig. 16). If Nta and cobalt were conservative and dispersion were negligible, the graph would show square pulses that increase at 10 h and decrease at 30 h. However, the movement of the Nta and cobalt is retarded relative to conservative movement by the sorption reactions, and small concentrations arrive early because of dispersion. The peak in Nta and cobalt concentrations occurs in the CoNta^{ } complex between 30 and 40 h. The peak in Co^{ 2+} concentration is even more retarded by its sorption reaction and does not show up until near the end of the experiment.
In figure 17, solidphase concentrations are plotted against time for concentrations in the last cell of the column. The sorbed CoNta^{ } concentration peaks between 30 and 40 h and lags slightly behind the peak in the dissolved concentration of the CoNta^{ } complex. Initially, no Nta is present in the column and the biomass decreases slightly over the first 10 h because of the firstorder decay rate for the biomass. When the Nta moves through the cells, the biomass increases because the Nta becomes available as substrate for microbes. After the peak in Nta has moved through the column, biomass concentrations level off and then begin to decrease because of decay. The K _{ d} for cobalt sorption gives a greater retardation coefficient than the K _{ d} for CoNta^{ } sorption, and the sorbed concentration of Co^{ 2+} appears to be still increasing at the end of the experiment.
Both the 10cell and the 20cell models give similar results, which indicates that the numerical errors in the advectivedispersive transport simulation are relatively small; furthermore, the results are very similar to results given by TebesStevens and Valocchi (1997) and TebesStevens and others (1998). However, TebesStevens and Valocchi (1997) included another part to their test problem that increased the rate constants for the sorption reactions from 1 to 1,000 h^{ 1} . The increased rate constants generate a stiff set of partial differential equations, which are equations that describe processes that occur on very different time scales. The stiff problem, with very fast sorption reactions, proved intractable for the explicit RungeKutta algorithm, but can be solved with the cvode algorithm. With cvode, the calculation of the slow sorption column takes about four times longer than with RungeKutta and calculation of the fast sorption column is about six times longer than for slow sorption. Another way to solve the stiff problem is to introduce the fast kinetic sorption reaction as an equilibrium process, which constitutes instantaneous rates. However, even with equilibrium sorption, grid convergence was computationally much more intensive; it was necessary to use 100 cells or more to arrive at a satisfactory solution. As an estimate of relative CPU times, the 20cell model took 2.7 times the CPU time of the 10cell model. A 200cell model took approximately 600 times the CPU time of the 10cell model.