Simulating transport and removal of xylene during remediation of a sandy aquifer

JOURNAL OF

Contaminant Hydrology ELSEVIER

Journal of ContaminantHydrology 19 (1995) 205-236

Simulating transport and removal of xylene during remediation of a sandy aquifer Wolfgang Sch~ifer a,*, Ren6 Therrien b a Institut fiir Umweltphysik, Universidit Heidelberg, D-69120 Heidelberg, Germany b Ddpartement de G~ologie et Gdnie G~ologique, Universit~ Laval, Qudbec, Qud. G1K 7P4, Canada

Received 27 September 1994; accepted after revision7 March 1995

Abstract Xylene, originating from a spill, is present both as a nonaqueous-phase liquid (NAPL) at residual saturation near the water table, and as a dissolved groundwater component contaminating a sandy aquifer beneath an abandoned refinery. Three remediation wells are in operation on the site to prevent further xylene migration in the groundwater. Field observations indicate that microbially-mediated xylene degradation and oxygen and nitrate reduction occur in the aquifer. To realistically simulate dissolved xylene migration at this site, a three-dimensional numerical flow and transport model incorporating biochemical multispecies interactions and xylene dissolution from the NAPL has been developed. In the calibration process the variable contact area between the NAPL and groundwater and the vertical transverse dispersivity were identified as crucial parameters controlling the fate of xylene. The simultaneous modeling of a whole set of related reactive species made it also possible to quantify the observed biodegradation. Results indicate that it contributes in the same order of magnitude to total xylene removal than does extraction by the wells. The calibrated model will be used to assist in the design of an in situ bioremediation scheme, where biodegradation in the aquifer is enhanced by injection of an electron acceptor.

I. Introduction 1.1. P r e v i o u s w o r k

Groundwater contamination by petroleum products is a widespread problem in many industrialized countries. These chemicals are considered immiscible, but will nevertheless dissolve in the groundwater and can reach concentrations that are harmful for

* Corresponding author. 0169-7722/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0169-7722(95)00018-6

206

W. Schiller, R. Therrien / Journal of Contaminant Hydrology 19 (1995) 205-236

human health. As dissolved species, they will migrate due to advection and dispersion, but will also be subject to reactions in the subsurface, such as biodegradation. The understanding of the fate of these contaminants in groundwater systems therefore requires a good knowledge of the physical and biochemical processes affecting their transport. Several models that simulate the reactive transport of organic contaminants have been developed. These models usually combine a classical description of advective-dispersive transport with a source/sink term that is used to represent the reaction between the organic compound and other dissolved species or aquifer material. Reactions can be represented by simplified zero- or first-order source/sink terms, in which case the reactive transport of a single species, dissolved organic carbon, is simulated. This approach might be suitable for some cases, but a more rigorous representation of reactive transport is often required. One conceptualization which has proven to be successful consists in describing the biodegradation of organic contaminants by microbially-mediated redox reactions. Models have been presented where the simultaneous transport of dissolved organic carbon and dissolved oxygen, acting as an electron donor and acceptor, respectively, along with the related dynamics of a bacterial population are simulated (e.g., Borden and Bedient, 1986; MacQuarrie et al., 1990). Additional reactions have also been considered, including denitrification (Widdowson et al., 1988; Kinzelbach et al., 1991) and fermentation (Kindred and Celia, 1989). There are, however, only a few reported studies where numerical models, with varying degrees of complexity, have been applied to real cases of groundwater contamination by hydrocarbons. In an early effort (GTC, 1982), a two-dimensional non-reactive transport model was used in the vertical plane to simulate groundwater contamination by hydrocarbons downstream of a special waste disposal site in Gloucester, Ontario, Canada. The authors qualitatively compared measured to simulated concentration distributions and proposed a hydraulic remediation design based on their simulations. Diersch and Kaden (1984) used a two-dimensional areal model to simulate the transport of organic contaminants in the groundwater below a manufacturing facility in the former German Democratic Republic. With the help of the model they compared several alternative remediation strategies. Kinzelbach (1985) used a two-dimensional flow model and a three-dimensional random walk transport model, including first-order decay and linear retardation, to simulate the fate of trichloroethane in an alluvial aquifer in southern Germany. He compared measured to simulated breakthrough curves and concentration distributions, and predicted the effects of a proposed hydraulic remediation measure. The studies listed up to here did not consider reactive transport or made use of the simplified source/sink approach mentioned previously. Departing slightly from this simplified approach, Borden et al. (1986) studied the reactive transport of hydrocarbons in the groundwater below a creosote site in Texas, U.S.A. They used a two-dimensional horizontal transport model linked to an improved reaction model that combined first-order decay and a very simple form of a multispecies model. They matched the model results to measured concentrations by varying the aquifer dispersivities. This represented the first effort where a calibrated model was used to study in situ bioremediation alternatives. An application of a more sophisticated model simulating microbially-mediated organic carbon degradation at the same site was

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

207

presented by Borden et al. (1989). The authors compared measured and simulated concentrations for a push-pull test to investigate the potential of in situ aerobic biodegradation. Chiang et al. (1989) simulated the transport of benzene, toluene and xylene and the aerobic degradation of these compounds with a two-dimensional flow and transport model coupled with the reaction model of Borden et al. (1986). They simulated measured concentration distributions for a contaminated aquifer in Michigan, U.S.A., and identified the initial distribution of the organic compounds and the longitudinal dispersivity to be the crucial model parameters. Another type of reactive transport model application is presented by MacQuarrie and Sudicky (1990), who simulated benzene transport in a vertical cross-section of the Borden aquifer in Ontario, Canada, with a two-dimensional flow and transport model incorporating a kinetic model of aerobic degradation. Their innovation consisted in considering the effect of aquifer heterogeneity on benzene degradation under natural flow conditions through the use of a heterogeneous distribution of the hydraulic conductivity. A more complex degradation model was used by Semprini and McCarty (1991) and Semprini and McCarty (1992) to investigate in situ biorestoration beneath the Moffet Naval Air Station in California, U.S.A. The presence of a uniform flow field allowed them to use a one-dimensional advective-dispersive transport model. Their coupled

e

e C ) R44.6 •

(~ R44.8

R 41.2

OC) R 41.4

R 41.3 •

o-xylene/ =% / • ~ / / ~ } I

•

observation well

0

remediation well

100 m I

J

Fig. t. Plan view of a part of the abandoned refinery with the position of observation wells (full circles) and remediation wells (open circles) relative to the location of the former xylene plants.

208

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

reactive model considered microbial growth, subsequent methane and oxygen consumption, cometabolic degradation of chlorinated aliphatic compounds and rate-limited sorption. As part of the monitoring of a biotreatment experiment in the field, they performed an extended comparison of measured and simulated concentrations which was aimed at the optimization of the biorestoration measure. Kinzelbach et al. (1991) used a two-dimensional flow and reactive transport model to simulate in situ bioremediation of groundwater contaminated by benzene and toluene in the Upper Rhine Valley in Germany. Their biochemical model considered aerobic and denitrifying microbial growth and subsequent substrate consumption. The modeling results indicated that the exchange process between the pore water and biophase was the limiting factor for the overall contaminant removal. In another study, Brusseau (1992) used a one-dimensional flow and transport model coupled with a model for rate-limited sorption and a two-domain approach to successfully simulate measured organic contaminant behaviour for selected field-scale experiments described in the literature. Bekins et al. (1993) simulated the methanogenic degradation of phenolic compounds in an aquifer below a wood treatment plant in Florida, U.S.A. They used independently measured microbial kinetic parameters to reproduce measured concentrations along a one-dimensional profile. The value of the constant biomass and the microbial decay constant, which were not determined in the laboratory, turned out to be important fitting parameters in this case. The review presented here indicates that, in general, there is an increasing complexity in the type of (bio)chemical models used to simulate the subsurface degradation processes. This evolution is in part due to the fact that simplified approaches are not always suitable to study reactive transport, and also a result of the steadily increasing computational power of the computers on which these models can be run.

1.2. General situation and objectives The site we examined during this study is an abandoned refinery located in the lower Rhine area in Germany. Several processing units were operated at the refinery and have caused the contamination of the surficial sand aquifer by a variety of petroleum products. Because of the spatial distribution of these former processing units, there are discrete contamination zones on the site, with each zone characterized by the presence of one specific compound or a group of organic compounds. This paper focuses on the area of the former xylene processing plant, where the subsurface has been contaminated mainly by the three dimethylbenzene isomers (o-, m- and p-xylene) and by ethylbenzene, with m- and p-xylene being by far the most abundant contaminants. For the sake of clarity, all four chemicals will be hereafter referred to as xylene.

Fig. 2. Interpolated concentration distributions for the xylene production area shown in Fig. 1: (a) dissolved organic carbon (DOC) coming from xylene (contours: 1-20 mg L -1 , interval: 1 mg L - I ) ; (b) dissolved oxygen (DO) (contours: 1-8 mg L -1, interval: 0.5 mg L - 1); (c) nitrate (contours: 1-40 mg L - l , interval: 2 mg L - l ) ; and (d) alkalinity (contours: 0.8-7.6 mmol L-1, interval: 0.4 mmol L - 1). The measurements were made in August 1990.

W. Schlifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

(a)

209

(b) DOC [rag/I] 0

~3

0

! o

Ioo rn

®remediation well

L.u

~

na

"~ remedmtlon well

(d)

Cc)

9 I I I

loq m ~

®remediation well

l\J I I I~LIA

I/[I~UIA/II/I I I I LJ

? loo m ® remediation well 2~Ll_L[ J l l d j ~ l l L / ~ I L.t t L ! ~.L~LLI~LI

210

W. Schiller, R. Therrien / Journal of Contaminant Hydrology 19 (1995) 205-236

In May 1988, after contamination was discovered, three partially-penetrating pumping wells were installed to prevent further spreading of xylene in the groundwater. A total of 23 observation wells were also installed to monitor the extent of the contamination in the area surrounding the former plant (Fig. 1). The measured concentration distributions of some species in those wells revealed an interesting pattern (HPC, 1991). The upper part of the sand aquifer is characterized by enhanced xylene concentrations but very low dissolved oxygen (DO) and nitrate concentrations, whereas DO and nitrate concentrations are high outside the contaminated area (Fig. 2a-c). Furthermore, higher calcium concentrations and alkalinity values are measured inside and downstream of the contaminated area (Fig. 2d). These observations strongly suggest that aerobic and denitrifying bacteria are active and are using xylene as their carbon and energy source, thereby contributing to the aquifer remediation. In the spring of 1992, a large amount of xylene was detected in the form of a pure nonaqueous-phase liquid (NAPL) present at residual saturation. This pure phase xylene is located in the upper portion of the aquifer, near the water table and upstream of the pumping wells. The maximum concentrations measured were as high as 72 g of xylene per kg of soil. Based on the field measurements, the total amount of residual pure-phase xylene is estimated to be between 20 and 30 t 1. The presence of this NAPL suggests that pure-phase xylene migrated downward in the unsaturated zone after an accidental spill at the plant, and, being lighter than water, formed a body floating near the water table. The NAPL eventually reached residual saturation and became immobile. This NAPL constitutes the source for xylene contamination because of dissolution of xylene in the surrounding groundwater. The location of the NAPL measured in the soil is consistent with observations of the vertical distribution of xylene at the site, which reveal that dissolved xylene concentrations are higher near the water table and decrease with depth. The purpose of our study was to develop and use a numerical model to identify the mechanisms governing xylene transport at the former refinery site and to assess their relative importance with respect to xylene removal. The numerical model simulates the three-dimensional transport of several reactive species in the groundwater along with their biochemical interactions. We will first briefly describe the numerical model and then present the information available on the aquifer properties and input transport parameters. We will then describe our model calibration, which consisted in reproducing the concentrations measured at the remediation wells, and discuss the influence of the most sensitive parameters affecting xylene transport. We will finally comment on the importance of biodegradation for xylene removal for this specific case. The model application presented here combines for the first time a fully three-dimensional flow and transport model with a complex model incorporating microbially-mediated contaminant degradation processes. Because this field example is representative of many contamination cases by petroleum products that are lighter than water, the modeling strategy and the results we present should also be applicable to similar cases.

i 1 t = 1 metric tonne=

10 3

kg.

W. Schiller, R. Therrien/Journal of Contaminant Hydrology 19 (1995) 205-236

211

2. Modeling approach 2.1. The numerical model

The field observations for the present case dictated our choice of a conceptual model. The evidence of biologically-mediated reactions required that a multiple species approach be used along with the specific consideration of relevant reactions. Also, the irregular distribution of the contaminant in both the horizontal and the vertical dimensions and the presence of partially-penetrating wells strongly suggested the use of a three-dimensional model. The partially-penetrating wells are drawing contaminated water from the upper saturated part of the aquifer and clean water from the bottom part. This mixing effect cannot be reproduced with a more simple vertically-averaged two-dimensional model. Finally, the presence of the NAPL had to be explicitly considered in the model because it represents the source of xylene. We therefore needed to develop a model that could simulate three-dimensional groundwater flow and multispecies transport, along with microbially-mediated nitrate and oxygen reduction, as well as xylene dissolution from the NAPL present at residual saturation. The three-dimensional groundwater flow and solute transport model is a simplified version of the one developed by Therrien and Sudicky (1995). We only consider fully-saturated flow for this work. The model is based on a standard Galerkin finite-element method for the discretization of the flow and the solute transport equations. Details of the method are not given here but can be found in standard texts (e.g., Huyakorn and Pinder, 1983). An option to use a finite-difference representation for the governing equations is also implemented in the model (Panday et al., 1993). Upstream-weighting of velocities and a fully implicit discretization in time are used for the transport solution for the mobile reactive species. The resulting matrix equations are solved with an efficient preconditioned ORTHOMIN iterative technique. Verification problems that test the flow and non-reactive transport model against analytical solutions are presented in Therrien and Sudicky (1995). The biochemical activity model is an expanded version of the one described in Kinzelbach et al. (1991) and, in more detail, in SchMer (1992). The biochemical activities considered are: (1) aerobic and denitrifying microbial growth; (2) subsequent consumption of xylene (electron donor and carbon source for cell growth), DO and nitrate (alternate electron acceptors); (3) effects of microbial carbon dioxide production and proton consumption on the calcite equilibrium in the aquifer. Xylene is allowed to dissolve into the pore water from the NAPL, which is fixed to the aquifer material. Dissolved xylene can also adsorb to or desorb from the aquifer material directly. Exchange processes are also considered between the mobile pore water and an immobile biophase, where the microbially mediated redox reactions take place. The exchange processes are modeled with linear exchange terms, and the exchange rates are governed by individual exchange coefficients. Biochemical degradation activity is simulated for a total of nine reactive species. Oxygen, nitrate and organic carbon species are present in the mobile pore water and the

212

w. Schiller, R. Therrien/ Journal of Contaminant Hydrology 19 (1995) 205-236

immobile biophase. The model considers organic carbon adsorbed onto the aquifer material and organic carbon present in the NAPL. Growth or decay of the bacterial population are also explicitly simulated. The effect of organic carbon oxidation and nitrate reduction on the equilibrium of calcite in the aquifer is modeled in a subroutine of chemical equilibration. The processes considered are the dissociation of carbonic acid, the auto-protolysis of water and the dissolution of calcite. An infinite reservoir of calcite is assumed to exist and reactions are taking place in a closed system where carbon dioxide is not allowed to degas. The equations describing transport and reaction of all species in the model are presented in the Appendix. Coupling of transport and biochemistry is done via an iterative three-step procedure. In a first step advective-dispersive transport is calculated implicitly for all the mobile species, including those involved in the calcite equilibrium (Eq. A-2 in the Appendix) while concentration changes due to the reactions associated with microbial activity and calcite equilibrium are handled as explicit source terms (Sm). In the second step, the exchange processes and microbial growth and subsequent electron acceptor and donor consumption are computed implicitly (Eqs. A-6-A-15 in the Appendix), while the concentration changes from advective-dispersive transport are represented by explicit source terms with values taken from the first-step calculations. In the third step the calcite equilibrium is updated with respect to the concentration changes from physical transport (step 1) and microbial carbon dioxide production and proton consumption (step 2). The operator splitting technique described above presents some important advantages compared to a one-step scheme, where transport and reaction equations are solved simultaneously, but it also creates additional numerical errors (e.g., Valocchi and Malmstead, 1992). However, these additional errors are minimized in the model by iterating between the threesteps. During the course of one iteration, the explicit reaction term in the transport step is updated with the results from the reaction steps of the preceding iteration (Herzer and Kinzelbach, 1989). The performance of a two-step procedure for a system of transport and biochemistry is described in Sch~ifer and Kinzelbach (1990) and Sch~ifer (1992). 2.2. Aquifer properties and input parameters

The contaminated aquifer consists of calcite-rich sand and gravel. Its saturated thickness, determined from two boreholes drilled to the bottom of the aquifer, is between 14 and 16 m. The related depth to the water table varies between 6 and 8 m. The transmissivity of the aquifer was determined by several pumping tests (HPC, 1992b). Different methods of evaluation for four locations close to the xylene spill yielded transmissivities between 1400 and 3500 m 2 day-1 with an average value of 2600 m 2 day -1. No quantitative investigation has been performed to study the spatial variability of the hydraulic conductivity of the aquifer. However, the aquifer material from the aforementioned two boreholes was qualitatively assessed to consist of rather homogeneous coarse sand. Although not verified by measurements the aquifer was

W. Schiller, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

213

assumed to be moderately anisotropic in the vertical direction. The recharge rate was estimated from climatic conditions and soil type. Piezometric heads at the site are extensively monitored through a large number of observation wells. Data at these wells have been collected weekly since May 1988 and have enabled a reliable determination of the mean regional flow direction and hydraulic head gradient. The three remediation wells and the observation wells are partially penetrating the aquifer. They are screened from the water table down, with lengths varying between 2 and 6 m. It was thought at the time of construction that locating the screens at these elevations would provide a more efficient capture of the xylene plume, which is restricted to the upper part of the aquifer. Fewer measurements were available to determine input parameters for the transport of the dissolved solutes. The effective porosity, the dispersivities and the parameters of the microbial kinetics had to be estimated, since no independent tests were conducted to evaluate these parameters. A broad range of values for the specific yield factor, the maximum growth rate of the bacteria and the Monod constants are reported in the literature either for natural (e.g., Button, 1985) or contaminated systems (e.g., Goldsmith and Balderson, 1988; G~illi and McCarty, 1989; O'Reilly and Crawford, 1989). We decided to use a low value for the growth rate because of the low temperature prevailing in the groundwater and a low Monod constant because of the mainly oligotrophic nature of the groundwater system (Ishida et al., 1982). The microbial decay rate, the utilizable portion of dead bacteria and the yield factor were chosen arbitrarily, while the coefficients of oxygen and nitrate consumption were calculated from the yield factor via stoichiometric relationships. However, we found from our numerical simulations that the transport conditions were at near steady-state around the plume for the time span under consideration. Thus although microbial degradation could not be neglected, there was almost no net growth or decay of microorganisms in the model and therefore the parameters describing microbial kinetics had only a minor effect on the simulated xylene concentrations in the wells. The maximum solubility and the octanol-water partitioning coefficient for xylene were taken from Montgomery and Welkom (1990). The fractional organic carbon content of the aquifer was assumed to be low, because the aquifer material consists mainly of organic-poor coarse sand. The exchange coefficients for the transfer of xylene and the electron acceptors are, from their nature, fitting parameters. As default, the values for the three exchange parameters o~, /3 and Y (compare the Notation Section A-5 in the Appendix) were set to infinity, i.e. instantaneous exchange was assumed. During calibration we found no necessity to change these values except for the exchange between NAPL and pore water. Here it turned out that using a slight diffusion limitation (i.e. a decreased value of T) provided a better fit. The meaning of these results will be discussed in Section 4.2. The penetration depth of NAPL-xylene into the groundwater and the transverse dispersivity in the vertical direction, which had to be fitted, turned out to be the most important input parameters. They will be discussed in more detail in the section on the model calibration. The selected values for the hydraulic and biochemical model parameters are given in Tables 1 and 2, respectively.

214

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

Table 1 Values for selected hydraulic model parameters Parameter

Value

Horizontal hydraulic conductivity, k h Vertical hydraulic conductivity, k z Regional head gradient, ah/Oy Effective porosity, n e Longitudinal dispersivity, a I Horizontal transverse dispersivity, ath Vertical transverse dispersivity, atv Groundwater recharge, q Average pumping rates at wells R41.2, R41.3, R41.4, R44.6 and R44.8 Discretization x-direction, y-direction, z-direction,

172 m day-1 69 m day- l 0.07% 25% 5.0 m 1.0 m 0.01 m 200 mm yr- l 202, 167, 185, 234 and 256 m 3 day- l, respectively

in: /ix Ay Az

10-100 m 10-100 m 0.25-5 m

Extent of model in: x-direction y-direction z-direction

650 m 630 m 15m

Table 2 Values for selected biochemical model parameters Parameter Maximum aerobic and denitrifying growth rate, fil.aer and /Zde., respectively Constant microbial decay rate, /Xdec Monod constants for: oxygen, K o nitrate, K N organic carbon (from xylene), Koc

Value 0.5 day- 1 0.05 day i

0.2gm 0.5gm 2.0gm

3 3 3

Utilizable portion of dead bacteria, f Yield factor for organic carbon assimilation, Y~er, Yde, Coefficient of oxygen consumption, R o Coefficient of nitrate consumption, R N

90% 0.1 g bio-carbon/g organic carbon 31.5 g o x y g e n / g bio-carbon 48.6 g nitrate/g bio-carbon

Exchange coefficients: pore water-biophase, a pore water-aquifer material, /3 pore water-NAPL, 3'

(instantaneous exchange) oo (instantaneous exchange) 0.5 day- 1

Maximal solubility of organic carbon from xylene in pore water, OCso l Octanol-water partitioning coefficient for xylene, Kow Fractional organic carbon content of the aquifer, foc Adsorption coefficient for xylene, k d

170gm

3

1,400 0.07% 0.6 cm 3 g - 1

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

215

2.3. Discretization and boundary conditions Fig. 3 shows a horizontal view of the model grid together with the location of the NAPL as measured in the soil, the position of the wells and simulated piezometric heads. It should be noted that we considered the three remediation wells constructed especially for xylene removal (R41.2, R41.3, R41.4) and also two wells (R44.6, R44.8) from a nearby site that influence the flow field in the model area (compare Fig. 1). The inner part of the model area is surrounded by increasingly larger model cells, such that the influence of the boundaries on the flow field is minimized. The grid spacing in the inner area is 10 m in the x- and y-directions. We vertically discretized the saturated portion of the aquifer with six layers, with layer thicknesses varying between 0.2 m at the top of the domain and 5.0 m at the bottom of the aquifer. The grid is oriented such that it coincides with a constant potential on the upper and lower boundary (constant-head boundaries) and a streamline on the left- and right-hand sides (impervious boundaries). It is assumed that the constant-head boundaries are uniform in the vertical direction. There is a constant recharge at the top of the domain and the bottom of the aquifer is considered impermeable. The concentrations of the mobile species are fixed at the inlet (upstream) boundary, with values interpolated from concentration measurements in five observation wells located close to this boundary (HPC, 1991). At the outlet (downstream) boundary, the solutes are allowed to freely advect out of the domain. The recharging water has a constant and uniform DO concentration of 10.0 mg L - I and is free of nitrate and

630 m

0

650 m

Fig. 3. Horizontal view of the top layer of the model grid, The black area denotes the location of the NAPL. Also shown are the location of the wells and the simulated piezometric heads. The discretization in the central part is 10 m x 10 m. The discretization and boundary conditions are the same for all six layers.

216

w. Schiifer, R. Therrien/Journal of ContaminantHydrology 19 (1995) 205-236

xylene. It is further assumed that there is no mass flux across the lateral boundaries and the bottom of the aquifer. The mass of xylene present in the NAPL is obtained from field measurements (HPC, 1992a). 2.4. Initial conditions for transport

Data on the distribution of xylene and the other species have only been collected since the beginning of the remediation process in May 1988. Consequently, no data are available concerning the xylene spill and the plume evolution prior to remediation. Since our objective was to reproduce the available data by simulating xylene migration after the start of remediation, we needed to specify initial concentrations for all species. We had the choice between two scenarios to start our simulation. The first option consisted in specifying initial concentrations for xylene and the other species based on the field measurements made at the start of remediation. The second scenario involved performing a history modeling of the contamination by starting our simulation at some earlier time, which corresponds to the time when the spill occurred, and letting the plume form until the start of remediation in May 1988. This preliminary simulation produces a set of initial concentrations at the time remediation starts. With respect to the first scenario, there is an initial time span required in multispecies systems during which the reactive species have to adapt to each other until consistency according to the boundary conditions is reached. The problem in the field case presented here consisted in specifying dissolved xylene concentrations at the beginning of the well operation that correspond to the xylene flux from NAPL-xylene dissolution. If the initial xylene concentrations do not fit the concentrations that result from the dissolution process, a numerical xylene wave moves through the model area and disturbs the interpretation of modeling results. We found that the adaptation time for the xylene concentration with respect to dissolution was on the order of one-fifth of the total desired simulation time of 4.5 yr. We judged that this period was too long to be acceptable. Therefore, to avoid the interferences resulting from adaptation, we opted instead for the second scenario. To obtain the initial concentrations at the start of remediation using history modeling, we needed to determine the time period for which the NAPL-xylene was in contact with groundwater prior to the installation of the pumping wells. In other words, we needed to estimate the time at which the xylene spill occurred. There was no information available about the spill but we could still determine this time span by comparing initial xylene concentrations at wells R41.2 and R41.3. Well R41.2 is directly downstream of well R41.3, which is located at the downstream edge of the NAPL zone (see Fig. 3). This implies that if the NAPL had been in contact with groundwater for a very long time prior to the start of remediation at the wells, both wells would have shown nearly the same initial concentrations. We assumed here that dilution is negligible for this short distance between the wells and the proximity of the source. In reality, the initial concentration measured at well R41.2 was eight times lower than that at well R41.3. We therefore iteratively decreased the pre-pumping contact time until the modeled initial concentrations in the two wells matched the measurements. We estimated by using this procedure that the NAPL must have been in contact with the groundwater for a period of

W. Schiifer, R. Therrien/Journal of Contaminant Hydrology 19 (1995) 205-236

217

3 yr before the start of remediation. We therefore modeled the plume formation after the spill for a period of 3 yr, under regional flow conditions, to obtain the initial concentrations at the start of remediation. This allowed us to avoid spurious effects from numerical adaptation that occur when we directly start the simulation at the time the wells are put into operation. We used a time step equal to 10.0 days for transport for the 3-yr pre-pumping period and 2.0 days for the actual remediation simulation, which starts in May 1988 when the wells are put in operation. Remediation was simulated for 1642 days until November 1992. In November 1992, well R41.3 had to be shut down because of clogging and some new wells were installed. Future simulations will consider the period after November 1992.

3. Model calibration

3.1. Uniqueness of calibration results The first step in our model calibration consisted in simulating the flow system in the aquifer. We had a large amount of information on the spatial distribution of piezometric heads, well characteristics, pumping rates and hydraulic conductivity. We could therefore suitably reproduce the observed piezometric heads by using these known input parameters. Table 3 shows a comparison between measured and simulated piezometric heads for January 1992.

Table 3 Comparison of measured and simulated piezometric heads in the observation and pumping wells of the contaminated area shown in Fig. 3 for January 1992 Well No.

Measured head (m)

Obseruation wells: R5.0 20.05 R 5.1 20.(18 R5.2 20.20 R5.3 20.23 R5.4 20.21 R5.5 20.13 R 5.6 20.09 R 5. 7 20.05 R5.8 20.26 R5.9 20.22 R6.0 19.98 R6.1 20.62 R6.2 20.14 R6.3 20.18 R6.4 20.00

Simulated head (m) 20.2 20.09 20.20 20.22 20.21 20.11 20.09 20.03 20.25 20.22 19.95 19.95 20.12 20.14 20.03

Well No.

Measured head (m)

Obseruation wells (cont.): R6.5 19.8l R 6. 7 19.72 R6.8 19.99 R6.9 19.83 R20.2 19.98 R41.1 19.99 R44.0 19.97

Simulated head (m) 19.92 19.93 19.92 19.96 19.98 19.96 19.91

Pumping wells: R41.2 R41.2 R41.4 R44.6 R44.8

19.21 18.61 18.62 19.39 19.39

19.94 20.10 20.00 19.90 19.91

The units are meters above sea level, the simulated heads are taken from the top layer. The model grid is too coarse to resolve the specific effects that lead to the large drawdown measured in the pumping wells.

W. Schiller, R. Therrien / Journal of Contaminant Hydrology 19 (1995) 205-236

218

30

O0

(a)

WELL

R412

2 7 O0 o data 24 O0 21

~

--mod

O0

< ~18

O0

u15

O0

72 00 x

9 O0 6 O0

o

3 O0 O0 O0

34000

6B000 duration

102000

of

138000

remediation

770000

[d]

(b)

30 O0

R41 3

WELL 27. O0 0 data

24 O0 21

--model

O0

< o~18 O0 E u15

O0 , 0

0

12 00 x

9 O0

oo-:Vo-O-~'~_%

°

6 O0

oooo

o

o

o o °oO

o

o

O0

i

O0

o

o i

O0

i

680

duratlon 30 O0

0

o

o

i

340

Oo 0 0 0

ooo o

o o

3 O0

~7

o

(c)

WELL

L

O0

of

i

1020

_ _ i

k

O0

remediation

1360

O0

1700

O0

[d]

R414

O0 o data

24 O0

21

--model

O0

< ~18

O0

o f 5 O0 ~ 12 O0 x

9 O0

o

6 O0 o

3 O0 o

o

~T ~

O0 O0

340

O0

~ - - - L ~ ~LoooO o,oo~ o O ~ o ~ a ~ = ~ 6 8 0 O0 1020 O0 1360 O0 duration of remediation [d]

, I 7 0 0 O0

~,~

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

2190

piezometrzc

21 7 0 21 so

heads

zn

219

AS

0

DATA ~

MODEL

~ 2 ~ .so 24

10

I

o 2 0 90 20

70

20

50

20

30

20

10

19

90

I

1

o

t

I

c

I

-

O0

340

I

O0

680 duratior

'

O0 of

1 0 2 0 O0 remedzation

1360 [d]

O0

I Z00

O0

S

1988 1989 ] 1990 I 1991 1992 Fig. 5. Water table elevation measured in an observation well in the vicinity of the former xylene plant from May 1988 until November 1992. The dashed line indicates the piezometric heads used for the stepwise stationary approximation in the flow model.

Once the flow field was satisfactorily calibrated, we needed to simulate xylene migration in the aquifer. The central criterion used to assess the quality of the transport model results consisted in comparing simulated to measured concentrations in the three pumping wells, where concentrations have been recorded at weekly to monthly intervals since the beginning of the remediation in May 1988. Additionally, we checked modeled distributions of xylene, DO, nitrate and alkalinity against those measured in the pumping and observation wells during an extensive site characterization in August 1990. We cannot exclude that there exists one or more alternative sets of parameters to the ones we found during model calibration that could also suitably reproduce the measured xylene concentrations in the remediation wells. We agree with what Konikow and Bredehoefl (1992) stated for groundwater flow modeling, i.e. that "validation, per se, is a futile objective".This stems from the nature of a model being only a representation of certain aspects of the system being analyzed and not an identical copy of this system. But, as the authors further state, the process of iterative model calibration leads to a better understanding of the system under consideration. The reliability of the calibration results generally increases with an increasing database. In our case of reactive transport with a large number of parameters, the

Fig. 4. Comparison of measured and simulated xylene concentrations from beginning of remediation in May 1988 (day 0) until November 1992 (day 1642) for: (a) well R41.2; (b) well R41.3; and (c) well R41.4. A constant penetration depth of NAPL-xylene into the saturated part of the aquifer was assumed.

220

W. Schiller, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

ambiguity of the model results could be restricted by the fact that we had at our disposal detailed information on the flow in the aquifer, on the location of the NAPL source region and on xylene concentrations measured simultaneously in several wells for a period of > 4 yr. In the following we discuss the influence on the fate of xylene of the two most sensitive parameters determined from the calibration: the zone of contact between the NAPL and the groundwater and the vertical transverse dispersivity.

3.2. NAPL-xylene In the model we assumed that the only source of xylene in groundwater was dissolution from the NAPL-xylene. Thus the degree of mixing in a pumping well governs the xylene concentration in the pumped water. This mixing depends mainly on the three-dimensional capture zones of the individual wells and the penetration depth of the NAPL into the saturated zone of the aquifer (NAPL-depth). The simulation of the capture zones is based on well-known data like pumping rates, filter lengths of the wells and hydraulic conductivity of the aquifer. Although the horizontal extent of the zone where the NAPL is located was measured with a good resolution, very little information was available concerning the NAPL-depth. It was therefore subject to a fitting procedure. In a first attempt we assumed that the NAPL-depth was constant everywhere during the whole simulation period. Fig. 4 shows results for a NAPL-depth of 0.17 m. Significant deviations between model results and measurements arise, especially for the central well R41.3 which shows relatively high concentrations of xylene. For this example the model results were acceptable for the initial period, but became unsatisfactory for later times. Using a smaller value for the constant NAPL-depth led to a shift in the time span for which the model results were acceptable. We found that it was impossible to match the data over the whole simulation period with a single value for the NAPL-depth. An explanation for this can be obtained by examining Fig. 5, where the water table elevation is shown for the simulation period. Besides seasonal variations, one can notice a long-term drop of the water table from ~ 21.9 m in May 1988 to 19.9 m in early 1992. Because pure-phase xylene is less dense than water, it cannot penetrate deeply into the saturated zone of the aquifer. The appreciable water table fluctuations shown in Fig. 5 will therefore strongly influence the contact area between the flowing groundwater and the residual NAPL-xylene. To account for this effect, we decided to use 23 sub-periods (dashed line in Fig. 5) to simulate the remediation instead of using a constant average steady-state flow field for the whole simulation. We assumed steady-state conditions for each sub-period, which made it possible to incorporate the variable water table and therefore the variable

Fig. 6. Comparison of measured and simulated xylene concentrations from beginning of remediation in May 1988 (day 0) until November 1992 (day 1642) for: (a) well R41.2; (b) well R41.3; and (c) well R41.4. A variable penetration depth of NAPL-xylene into the saturated part of the aquifer according to the water table elevation shown in Fig. 5 was assumed. The concentration time series are shown for the final run of the calibration process.

W. Schii:er, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236 30 00

(a)

WELL R41

2

27 00 e

24 00 21

data

--model

O0

T. 1800 O0

u15

c 1200 ×

9 00 6 00

3 00 00

~

~

-

30 00

~

;

;

340 00

00

~ ~oo;~o~:e~o

680 00 duration of

'(b)

1020 00 cemediation

......

:--

__t:__

1360 00

,

1700 00

[d]

R413

WELL

27 00 o data 24 00

--model

21 00 c~18 00 E o 1500

~12 x

o

O0

9 00 6 00 o

oo

3 00 i

00

o i

340 00

00

30 00

o

o o

(c)

i

680 00 duration of

L

i

1020 00 remediation

i

k

i

i

1360 00

1700 00

136000

170000

[d]

WELL R 4 1 4

27 00 o 24 00 21

data

--model

00

q ~18

00

u15

O0

i

~;12 O0 × 9 O0 6 O0 o

3 O0 O0 O0

o

340 00

68000 duration of

102000 remediation

[d]

221

222

w. Schiifer, R. Therrien/Journal of Contaminant Hydrology 19 (1995) 205-236

NAPL-depth. As can be seen in Fig. 6, the consideration of the water table fluctuations allowed us to clearly improve the model results. 3.3. Dispersivities

Due to a lack of data, we did not explicitly consider a heterogeneous hydraulic conductivity distribution in the aquifer, but used instead the macrodispersion approach (e.g., Gelhar et al., 1979; Dagan, 1984), where the effects of aquifer heterogeneities on solute transport are represented by longitudinal and transverse macrodispersivities. The values of the macrodispersivities are scale-dependent and related to the degree of aquifer heterogeneity. In our case, where a quantitative information on hydraulic conductivity distribution and thus on the degree of aquifer heterogeneity was not available, we had to fit the dispersivity values during the calibration process. To allow for different values for horizontal and vertical transverse dispersivities, we use the approximation given by Burnett and Frind (1987) for the dispersion tensor. We performed a series of simulations to investigate the sensitivity of the modeled breakthrough curves in the wells to the values of dispersivity. We found that the sensitivity was greatest for the vertical transverse dispersivity and lowest for the longitudinal dispersivity. The relative insensitivity for the latter parameter could be expected since the continuous well operation creates a quasi steady-state situation for longitudinal solute transport to the wells. We therefore assumed a value of 5.0 m for the longitudinal dispersivity. This value appears reasonable for the dimensions of the domain considered and the rather homogeneous sandy aquifer material. It also fulfils the grid-P6clet criterion in the central area of the model area such that physical dispersion dominates the numerical one. The increased sensitivity for the horizontal transverse dispersivity can be explained by its influence on the horizontal extent of the contaminant plume, which affects the concentration of xylene especially in the two lateral wells R41.2 and R41.4. The calibration process yielded a value of 1.0 m for the horizontal transverse dispersivity. This relatively large value comprises the effects of the non-resolved fluctuations of the direction of regional flow. The high sensitivity of xylene concentrations in the wells to the vertical transverse dispersivity, atv, is a direct consequence of the large vertical concentration gradient that exists for xylene, with a near-saturation concentration of ~ 160 mg L -1 inside the NAPL-area and zero concentration in the uncontaminated water some decimeters below. Therefore, vertical transverse dispersion plays a key role for the mixing of contaminated and uncontaminated water. Enhanced mixing in the vertical direction increases the portion of contaminated water that enters the wells and causes an increase in xylene concentration in all three wells.

Fig. 7. Simulated concentration distributions for a horizontal plane ~ 2 m below the water table for: (a) dissolved organic carbon (DOC) coming from xylene (contours: 1-10 mg L -1, interval: 1 mg L-1); (b) microorganisms (contours: 1.3.105-1.3.106 cells/g dry soil, interval: 1.3.105 cells/g dry soil); (c) dissolved oxygen (DO) (contours: 1-8 mg L -1, interval: 1 mg L- 1); and (d) nitrate (contours: 2-10 mg L -1, interval: 2 mg L- 1). The concentrations are shown for the central part of the model aquifer which is more finely discretized (cf. fig. 3).

223

W. Sch~'fer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

(a) DOC [mg/1]

(b) microorganisms [cells/g dry soil]

~l.3xlO

1.3xlO5

(c)

(d)

6

224

W. Schiifer, R. Therrien / Journal of Contaminant Hydrology 19 (1995) 205-236

During the model calibration, we found that we needed a value of atv equal to 0.01 m to adequately reproduce the observed mixing behaviour in the wells. We could not get comparably satisfying results with a decreased value for at~ and, to compensate, an increased value for the NAPL-depth. This large value for atv is somewhat surprising since results from field studies (e.g., Sudicky, 1986) usually suggest a smaller vertical transverse dispersion, possibly being on the order of molecular diffusion. One explanation for this larger value for vertical dispersivity could be that it is an effective macroscopic parameter that represents unresolved vertical velocity components. However, this is unlikely for the case presented here as the vertical fluxes caused by the partially-penetrating wells and groundwater recharge are considered explicitly in the three-dimensional flow model. We think that a plausible explanation for the increased dispersivity fitted during the calibration results from the way we had to represent the NAPL body in our model. While the horizontal extent of the NAPL-xylene spill was known from total xylene measurements made on a 10 m × 10 m grid in the soil just above the water table (HPC, 1992a), virtually nothing was known about the vertical distribution of the residual xylene. It was not measured and could only be determined by the fitting procedure described in the preceding section. But the fitting procedure provides only a rough estimate and allows no insight into the detailed vertical xylene distribution. Consequently, we had to represent the body of NAPL-xylene by a simple compact slice with sharp upper and lower boundaries following the layering of the model. For the real vertical NAPL distribution, however, we would expect a much more irregular structure. The transition between the area completely contaminated with NAPL-xylene and the uncontaminated deeper part of the aquifer is gradual in reality while it is abrupt in the model. The increased vertical dispersivity we needed during calibration can then be thought as being an effective simulation parameter that reproduces the smooth NAPL structure transition occurring in the real aquifer.

4. Microbial activity 4.1. Effects on reactive species distribution Aerobic and denitrifying degradation of xylene have been shown to occur in the laboratory (e.g., Zeyer et al., 1986). The measured distribution pattern of the reactive species involved in xylene degradation strongly emphasizes that this type of degradation also takes place for the present field case. Direct measurement of biodegradation in the field is hard to achieve. Local measurements of rates with in situ devices (e.g., Gillham et al., 1990) are useful but are somewhat limited because of the large spatial heterogeneity that prevails in natural aquifers. The measurement of the contaminant concentration alone is also usually not sufficient to monitor degradation, because of the complex interactions between microbial activity and physical transport.

(a) R41.3

R41.2 DOC ling/l]

0

L

630 m

g:

(b) 14.5 m

R41.3 < ~ 1 . 3 x 1 0 5 ~

R41.2 1"31106~--

microorganisms [cells/g dry soil]

I I,,a

630 m Fig. 8. Simulated concentration distributions for a vertical slice passing through wells R41.2 and R41.3 for: (a) dissolved organic carbon (DOC) coming from xylene (contours: 1-10 mg L -1, interval: 1 mg L 1); and (b) microorganisms (contours: 1.3. 105-1.3 - 106 cells/g dry soil, interval: 1.3- 105 cells/g dry soil). A vertical exaggeration equal to 10 is used.

226

W. Schiifer, R. Therrien/Journal of Contaminant Hydrology 19 (1995) 205-236

In our numerical simulation we quantify degradation by modeling the organic contaminant together with the potential electron acceptors DO and nitrate, as well as the microbial effects on alkalinity and calcium concentrations. Figs. 7 and 8 show modeled concentration distributions for xylene, microorganisms, DO and nitrate ~ 2.5 yr after beginning of the well operation. Fig. 7 presents results for a horizontal plane located 2 m below the water table, which is very close to the lower boundary of the contaminant plume. Fig. 8 shows concentrations for a vertical cross-section intersecting wells R41.3 and R41.2 and parallel to the y-axis; this vertical section is oriented approximately along the flow direction. We can recognize an area of low xylene concentrations inside the plume around well R41.3 (Figs. 7a and 8a). These low concentrations are the result of vertical mixing of upper contaminated and deeper uncontaminated water in the vicinity of the partiallypenetrating well. The distribution of microorganisms shows an area of higher activity at the upstream boundary of the plume (Figs. 7b and 8b). There, the mixing of oxygen- and nitrate-rich water with the contaminated water creates favourable conditions for the growth of the aerobic and denitrifying microorganisms. Even more favourable growth conditions prevail around the wells where mixing is enhanced by the well operation. Except for the area around well R41.3, there is no aerobic or denitrifying microbial activity occurring inside the plume. This is caused by the fact that the oxygen and nitrate, which are transported by advection and dispersion, are completely consumed in the outer area (Fig. 7c and d). Under the simulated conditions the limiting factor for the microbial degradation activity is the rather low supply of the electron acceptors oxygen and nitrate. 4.2. Contribution to total xylene removal

Although designed as a purely hydraulic measure, xylene removal at the wells is actually complemented by in situ biodegradation. Results from modeling show that biodegradation decreases the xylene concentrations by about a factor of 2 in the two outer wells R41.2 and R41.4. The calibrated model also enabled us to evaluate the contribution of microbial degradation to total xylene removal. This is graphically illustrated in Fig. 9, where xylene removal by extraction and biodegradation is indicated at selected times prior to and after the start of remediation. Day zero corresponds to the start of the well operation in May 1988. Prior to that, regional flow conditions prevailed and there was of course no extraction by the wells; microbial degradation was therefore the only mechanism of xylene removal. With a very low removal rate of 0.25 kg day- 1 under regional conditions, > 100 yr would be required for the removal of the estimated 20-30 t of xylene present in the aquifer. In the initial phase of remediation, extraction dominated the effects of microbial activity, but the latter was also enhanced by the well operation. This was due to the increased influx of electron acceptors into the contaminated area induced by the discharge at the wells. Along the course of the remediation the total extraction rate at the wells continuously decreased. This was caused by the decreasing xylene concentration in the pumped water and also by the decreased discharge capacity of the central well due to

xylene removal rates 3 "~

2,5

--"

2

@

N

[ ?

I.

?: 7

•

microbial degradation

1,5 [ : extraction

~

1 g~

;~

0,5

o -100

196

504

868

1230

1642

duration of remediation [d] Fig. 9. Contributions of microbial activity and extraction to the total xylene removal rate under regional flow conditions (day - 100) and for selected times during the remediation.

bo ~n I

228

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

filter clogging. However, the decrease of the microbial activity was less pronounced such that in the end it became again the main factor for xylene removal. The controlling factor for the overall aerobic and denitrifying microbial degradation rate for the period modeled is the electron acceptor (oxygen and nitrate) flux into the contaminated area. This implies that degradation is limited by the hydraulics of the system and not by (bio)chemical kinetic processes (e.g., the time scale of net microbial growth is much smaller than that of total xylene removal). Also the time scale of matter exchange between different phases of the aquifer must be smaller than the residence time of oxygen and nitrate inside the contaminated area. These finding are supported by the field data which show: (a) a relatively static behaviour of the contaminant plume in the horizontal directions; and (b) a complete reduction of oxygen and nitrate inside the contaminated zone (see Fig. 2). Therefore the main potential of the three-dimensional flow and transport model applied for this situation consisted in providing spatially and temporally resolved oxygen and nitrate fluxes into the contaminated part of the aquifer. Based on these fluxes the microbial degradation activity could be rather easily calculated in the biochemical sub-model. In a later stage of remediation, for example, when the NAPL reservoir becomes more and more depleted, in situ conditions might change and processes such as rate-limited mass exchange between pore water and biophase or rate-limited desorption might then control the overall degradation rate. MacOuarrie and Sudicky (1990) state that the macrodispersion approach developed for non-reactive solutes, which was also used in the simulations presented here (compare Section 3.3), might overestimate organic carbon degradation. In our opinion this is mainly a problem for long-term reactive transport under natural flow conditions. In our case solute transport is dominated by advective fluxes induced by the pumping wells, such that possible shortcomings of the macrodispersion concept in conjunction with reactive transport should not cause noticeable effects on our simulated xylene degradation. The specific role of the vertical transverse dispersivity has been discussed in Section 3.3. In general, an acceptable quantitative agreement between measurements and model results could be achieved for xylene, DO and nitrate concentrations in the pumping and observation wells. However, only a qualitative agreement was found for alkalinity and for calcium. Their absolute concentration values were too small by a factor of 2 in the model. One reason for this could be that we underestimated the oxygen and nitrate fluxes into the contaminated area and therefore underestimated aerobic and denitrifying microbial degradation. But this is unlikely since the water fluxes are determined by the pumping rates of the three wells, which are well known, and the oxygen and nitrate concentrations are sufficiently resolved by the observation wells. Another factor that could contribute to an increase in microbial carbon dioxide production and therefore in alkalinity and calcium concentration is that additional types of electron acceptors are used by microorganisms. This hypothesis is supported by field measurements showing that dissolved iron and manganese concentrations were higher and sulfate concentrations were lower inside the contaminated zone than outside. These findings suggest that microbial activity goes beyond aerobic and denitrifying growth and that microbial activity under iron-, manganese- and sulfate-reducing conditions is also involved in xylene degradation. This also means that the relative contribution of total microbial

w. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

229

activity to xylene removal might be underestimated in the present model. The effects of additional electron acceptors will be investigated in the future with an extended biochemical model. Nevertheless, the presence of appreciable microbial activity in the contaminated aquifer suggests to exploit the microbial degradation potential of the system by implementing active measures from outside. Because our results indicate that the lack of electron acceptors influx into the plume area limits microbial activity, an in situ bioremediation program with addition of oxygen and/or nitrate into the subsurface is planned in the future. This appears to be a very promising way of speeding up aquifer remediation. The calibrated numerical model will be a useful tool in designing this remediation scheme.

5. Conclusions The application of a numerical model yielded insight in the mechanisms that govern xylene transport and removal in the groundwater below an abandoned refinery. Three distinct features could be identified to play key roles for the migration of xylene in this field case, which constitutes a typical example of contamination by petroleum products lighter than water: (1) Xylene in the form of NAPL at residual saturation exists in the transition area between the unsaturated and the saturated zones upstream of the wells. A large amount of xylene is entrapped in the NAPL. Its dissolution into the flowing pore water constitutes a long lasting source of groundwater contamination. (2) The interaction of five partially penetrating wells and the non-uniform vertical distribution of xylene requires that xylene migration be analyzed in all three dimensions and, finally, (3) The apparent microbial degradation activity influences xylene concentrations in the pore water. Therefore, the numerical model needed to study this field example had to consider three-dimensional flow and transport of xylene and other mobile reactive species as well as their mutual biochemical interactions. It also had to include the dissolution of xylene from the NAPL into the mobile pore water. The calibration process showed that the penetration depth of NAPL xylene in the flowing groundwater (NAPL-depth) constitutes a crucial model parameter. Changes in the NAPL-depth caused by water table variations had a prominent effect on xylene concentrations in the wells. Furthermore, we found that the vertical transverse dispersivity also greatly influences xylene migration. It controls the vertical mixing of contaminated and uncontaminated water. The unusually large vertical transverse dispersivity value obtained from calibration was attributed to its role as an effective parameter that accounts for the unresolved vertical irregularity of the NAPL body. This suggests that the characterization of the residual NAPL for similar cases would greatly help in designing remediation measures and estimating clean-up times.

230

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

The mechanistic multispecies model provided a realistic description of the xylene concentrations in the pumping wells over the modeled time period. The simultaneous consideration of the transport of the relevant species in the three-dimensional flow field and their biochemical reactions further allowed to demonstrate that the contribution of microbial activity to total xylene removal is as important as extraction by the wells. The calibrated numerical model represents a very useful tool to assist in the design of an intended in situ bioremediation, where microbial degradation potential will actively be enhanced by injection of electron acceptors.

Acknowledgements We gratefully acknowledge the collaboration of Dr. R. Samtleben from the Deutsche BP AG, Hamburg, and Dr. A. Greving, C. Heinecker and A. Voss from Harress Pickel Consult GmbH, Kassel, who provided us with all data and some additional measurements necessary to apply and calibrate the numerical model. We also acknowledge the comments and advice of Professor W. Kinzelbach during the various stages of this work. Funding for this research was provided by Koerber Foundation, Hamburg and by the Natural Sciences and Engineering Council of Canada.

Appendix The governing equations describing three-dimensional saturated groundwater flow and the transport of the reactive species in the biogeochemical model are presented. A description of all variables can be found at the end of the Appendix.

A-1. Saturated groundwaterflow Saturated groundwater flow is described according to the following equation: Oh

V(K. Vh) + q =Ss-~-t

(A-l)

A-2. Advective-dispersive transport The transport of dissolved mobile species in the pore water is governed by the following equation: ~7(0'

0Cm

~Tfrn ) -- ~ 7 ( v . Crn) -.~-am = - -at

(A-2)

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

231

The exchange between the mobile and immobile regions of the system is represented by the source term S m in Eq. A-2. Depending on the species, this source term S m will have the following form: Source term for the transport of oxygen: Sore = - o / ( 0 m - Oim )

(m-3)

Source term for the transport of nitrate: SNm = --c~(N m -Nim )

(A-4)

Source term for the transport of organic carbon: Sock, = - o ¢ ( O C

m - OCim ) - ~ ( O C m - O C m a x ) -- ' ~ ( O C m - O C s o , )

(A-5)

The source terms for the species from the calcite system are given implicitly by the new equilibrium. For comparison purposes, the transport of a non-reactive tracer is also modeled. The source term in Eq. A-2 becomes equal to zero for the tracer. A-3. Biochemical reactions

The transport of the immobile species is described by the following ordinary differential equation:

OCim at

(A-6)

= Sire

All reactions involving the immobile species are represented by the source term Sim in Eq. A-6. Depending on the species, this source term will have the following form: Reactions involving oxygen in the biophase: SO,m = '~(Om -- Oim) - R o

~

a~r

(a-7)

Reactions involving nitrate in the biophase: SNim

= c r ( N m - N i m ) - R N -0t - den

(A-8)

Reactions involving organic carbon in the biophase:

1 Soc'm

=

°~(Ofm

- Ofim)

0X

- Y'~er ~ ' -

1 aer

0X

Yden ~

+f

-~t

den

(A-9) dec

Sorption and desorption of organic carbon onto the aquifer material: SOCma t = ] ~ ( O C m - O C m a x )

where (A-10)

Ofma x = kdOC m

k a = 0.63fockow

(Karickhoffetal., 1979)

232

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

Dissolution of organic carbon from in the NAPL:

SOCNAeL= '~(OC m - OCsol)

(A-11)

Growth and decay of the bacterial population:

Sx= [~]aer-~- [ ~ t ]den-- [~]dec

(A-12)

where

[~X] [OX ]

OCim Oim aer = /'ZaerF(Oim) Koc + OCim Ko q- Oim x

(A-13)

OC im Nim den ~ /'~den[1-- g(Oim)] Koc + Oeim KN + Sim g

(A-14)

The term F(Oim) in Eqs. A-13 and A-14 represents a switching function between aerobic and denitrifying conditions and is defined according to Kinzelbach et al. (1991): F(Oim ) = 0.5 + ~.-1 arctan[(Oi m _ Ot)fs, ]

(1-16)

F(Oim) approaches a value of one if the oxygen concentration in the biophase Oim is well above the threshold concentration O t and a value of zero if Oim becomes small. The slope fst determines wether the transition between aerobic and denitrifying conditions is smooth or abrupt.

A-4. Calcite equilibrium For the calculation of the calcite equilibrium in the aquifer we consider the following six dissolved species: O H - , H +, H2C03,

HCO 3 , C032- , Ca 2+

which undergo the following four reversible equilibrium reactions:

H20 ~ H + + OH

H2CO 3 ~ HCO 3 + HCO 3 ~ C O ~ -

H+

+H +

CaCO 3 ~ Ca 2+ + CO 2 The law of mass actions provides four equations for the four reversible reactions. To come to a system of six equations to solve for the six unknown species we still need two equations representing conservation of mass conditions. A procedure to determine the

W. Schlifer, R. Therrien /Journal of Contarainant Hydrology 19 (1995) 205-236

233

conservation equations in a very general manner is described in Morel (1983). We used the following expressions: TOT H = H + - O H - - HCO 3 - 2CO32- + 2Ca 2÷ TOTH2co3 = H2CO 3 + HCO 3 + CO32- - Ca 2+ For this formulation of the conservation equations TOT n represents the total charge of the system and TOTH2co~ the total dissolved inorganic carbon species and calcium. The resulting system of six algebraic equations is solved with an Newton-Raphson method. The source terms for the six species are then calculated as the concentration difference before and after equilibration. The calcite equilibrium calculations are described in more details in SchMer (1992). A-5. Notation

Symbols used for the equations Symbol

Description

aer, dec, den

indexes for aerobic, decay and denitrifying, respectively concentration of an immobile species [M L 3] concentration of a mobile species [M L-3] dispersion coefficient tensor [L: T -1] utilizable portion of dead bacteria [-] switching function between aerobic [-] and denitrifying conditions slope of the switching between aerobic [-] and denitrifying conditions fluid source/sink [L3 L-2 T 1] hydraulic head [L] hydraulic conductivity tensor [L T 1] octanol-water partitioning [-] coefficient for xylene fractional organic carbon content [-] of the aquifer distribution coefficient for xylene [L3 M -1] half-velocity concentrations for [M L 3] oxygen, nitrate and organic carbon, respectively concentration of immobile oxygen [M L 3] threshold concentration of immobile oxygen [M L-3] concentration of mobile oxygen [M L-3] concentration of immobile nitrate [M L-3] concentration of mobile nitrate [M L 3] concentration of immobile organic carbon [M L -3 ]

Cim Cm D

f F

L, q h K kow

foc ko K o, K N, K o c

Oim Ot Om Nim Nm

OCi m

Dimension

234

W. Schiller, R. Therrien / Journal of Contaminant Hydrology 19 (1995) 205-236

Notation (continued) OC m OCmat

OCmax OCNAPL

OC~l Ro RN aim

Sm Ss t

X Y u

/3 Y /a.

concentration of mobile organic carbon concentration of organic carbon adsorbed onto aquifer material concentration of dissolved organic carbon in equilibrium with OCmat concentration of organic carbon in the NAPL solubility of organic carbon in the pore water coefficient of oxygen consumption during aerobic organic carbon dissimilation coefficient of nitrate consumption during denitrifying organic carbon dissimilation source term for an immobile species source term for a mobile species specific storage coefficient time concentration of microorganisms yield factor for organic carbon assimilation average pore-water velocity exchange coefficient between pore water and biophase exchange coefficient between pore water and aquifer material exchange coefficient between pore water and NAPL maximum growth rate of bacteria

[M L -3 ]

[M L 3] [M L -3 ] [M L 3] [ M L -3 ]

[-] [-] [M L -3 T-1] [M L -3 T-1] [L- 1] [T] [M L - 3 ]

[-] [L T- 1] [T 1] [T -1 ]

[T- 1 ] [T- 1 ]

References Bekins, B.A., Godsy, E.M. and Goerlitz, D.F., 1993. Modeling steady-state methanogenic degradation of phenols in groundwater. J. Contam. Hydrol., 14: 279-294. Borden, R.C. and Bedient, P.B., 1986. Transport of dissolved hydrocarbons influenced by oxygen-limited biodegradation, 1. Theoretical development. Water Resour. Res., 22: 1973-1982. Borden, R.C., Bedient, P.B., Lee, M.D., Ward, C.H. and Wilson, J.T., 1986. Transport of dissolved hydrocarbons influenced by oxygen-limited biodegradation, 2. Field application. Water Resour. Res., 22: 1983-1990. Borden, R.C., Lee, M.D., Thomas, J.M., Bedient, P.B. and Ward, C.H., 1989. In situ measurement and numerical simulation of oxygen limited biotransformation. Ground Water Monit. Rev. Winter, pp. 83-91. Brusseau, M.L., 1992. Transport of rate-limited sorbing solutes in heterogeneous porous media: Application of a one-dimensional multifactor nonideality model to field data. Water Resour. Res., 28: 2485-2497. Burnett, R.D. and Frind, E.O., 1987. Simulation of contaminant transport in three dimensions, 2. Dimensionality effects. Water Resour. Res., 23: 695-705. Button, D.K., 1985. Kinetics of nutrient-limited transport and microbial growth. Microbiol. Rev., 49: 270-297.

W. Schiifer, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

235

Chiang, C.Y., Salanitro, J.P., Chai, E.Y., Colthart, J.D. and Klein, C.L., 1989. Aerobic biodegradation of benzene, toluene, and xylene in a sandy aquifer - - data analysis and computer modeling. Ground Water, 27: 823-834. Dagan, G., 1984. Solute transport in heterogeneous porous formations. J. Fluid. Mech., 145: 151-177. Diersch, H.-J. and Kaden, S., 1984. Contaminant plume migration in an aquifer: Finite element modeling for the analysis of remediation strategies: A case study. Collab. Pap. from Int. Inst. Appl. Systems Anal., Laxenburg. G~illi, R. and McCarty, P.L, 1989. Kinetics of biotransformation of 1,1,1-trichloroethane by Clostridium sp. strain TCAIIB. Appl. Environ. Microbiol., 55: 845-851. Gelhar, L.W., Gutjahr, A.L. and Naff, R.L., 1979. Stochastic analysis of macrodispersion in aquifers. Water Resour. Res., 19: 161-180. Gillham, R.W., Starr, R.C. and Miller, D.J., 1990. A device for in situ determination of geochemical transport parameters, 2., Biochemical reactions. Ground Water, 28: 858-862. Goldsmith, Jr., C.D. and Balderson, R.K., 1988. Biodegradation and growth kinetics of enrichment isolates on benzene, toluene and xylene. Water Sci. Technol., 20: 505-507. GTC (Geologic Testing Consultants Ltd.), 1982. Gloucester special waste disposal site: Hydrostratigraphic interpretation and remedial measures assessment using mathematical modeling techniques. Geol. Test. Consult. Ltd., Ottawa, Ont., Rep. to Dep. Supply Services, Ottawa, Ont. Herzer, J. and Kinzelbach, W., 1989. Coupling of transport and chemical processes in numerical transport models. Geoderma, 44: 115-127. HPC (Harress Pickel Consult G.m.b.H.), 1991. Ver~inderugen der hydrochemischen Milieubedingungen dutch Zugabe einer KW-Quelle in einen Lockergesteinsaquifer. Harress Pickel Consult G.m.b.H., Kassel, Res. Rep. for Dtsch. BP A.G., Abt. VS 4241, Hamburg. HPC (Harress Pickel Consult G.m.b.H.), 1992a. Verdichtende Untergrunderkundung und SanierungsmaBnahmen im Bereich des ehemaligen Aromatenfeldes auf dem Gel~inde des Terminal Ruhr der Deutsche BP AG, Hiinxe. Harress Pickel Consult G.m.b.H., Kassel, Rep. for Dtsch. BP A.G., Abt. VS 4241, Hamburg. HPC (Harress Pickel Consult G.m.b.H.), 1992b. Bestimmung der Gesteinsdurchl~issigkeiten und Ermittlung der Entnahmeraten flir die Sanierungsbrunnen auf dem Aromatenfeld des Terminal Ruhr der Deutsche BP AG in Hiinxe. Harress Pickel Consult G.m.b.H., Kassel, Res. Rep. for Dtsch. BP A.G., Abt. VS 4241, Hamburg. Huyakorn, P.S. and Pinder, G.F., 1983. Computational Methods in Subsurface Flow. Academic Press. New York, NY. Ishida, Y., Imai, I., Miyagaki, T. and Kadota, H., 1982. Growth and uptake kinetics of a facultatively oligotrophic bacterium at low nutrient concentrations. Microb. Ecol., 8: 23-32. Karickhoff, S.W., Brown, D.S. and Scott, T.A., 1979. Sorption of hydrophobic pollutants on natural sediments. Water Resour. Res., 13: 241-248. Kindred, J.S. and Celia, M.A., 1989. Contaminant transport and biodegradation, 2. Conceptual model and test simulation. Water Resour. Res., 25: 1149-1161). Kinzelbach, W.K.H., 1985. Modeling of the transport of chlorinated solvents in groundwater: A case study. Water Sci. Technol., 17: 13-21. Kinzelbach, W., Sch~ifer, W. and Herzer, J., 1991. Numerical modeling of natural and enhanced denitrification processes. Water Resour. Res., 27:1123-1135. Konikow, L.F. and Bredehoeft, J.D., 1992. Groundwater models cannot be validated. Adv. Water Res., 15: 75-83. MacQuarrie, K.T.B. and Sudicky, E.A., 1990. Simulation of biodegradable organic contaminants in groundwater, 2. Plume behavior in uniform and random flow fields. Water Resour. Res., 26: 223-239. MacQuarrie, K.T.B., Sudicky, E.A. and Frind, E.O., 1990. Simulation of biodegradable organic contaminants in groundwater, 1. Numerical formulation in principal directions. Water Resour. Res., 26: 207-222. Montgomery, J.H. and Welkom, L.M., 1990. Groundwater Chemicals Desk Reference. Lewis, Chelsea, MI. Morel, F.M.M., 1983. Principles of Aquatic Chemistry. Wiley, New York, NY. O'ReiUy, K.T. and Crawford, R.L., 1989. Kinetics of p-cresol degradation by an immobilized Pseudomonas sp. Appl. Environ. Microbiol., 55: 866-870. Panday, S., Huyakorn, P.S., Therrien, R. and Nichols, R.L., 1993. Improved three-dimensional finite-element techniques for field simulations of variably-saturated flow and transport. J. Contain. Hydrol., 12: 3-33.

236

W. Schiller, R. Therrien /Journal of Contaminant Hydrology 19 (1995) 205-236

Schiller, W., 1992. Numerische Modellierung mikrobiell beeinflusster Stofftransportvorg~tngeim Grundwasser. Oldenbourg, Miinchen. Sch~ifer, W. and Kinzelbach, W., 1990. Modelling of pollutant transport in groundwater including biochemical transformations. In: G. Gambolati, A. Rinaldo, C.A. Brebbia, W.G. Gray and G.F. Pinder (Editors), Computational Methods in Subsurface Hydrology. Computational Mechanics Publications, Southampton, pp. 405-413. Semprini, L. and McCarty, P.L., 1991. Comparison between model simulations and field results for in-situ biorestoration of chlorinated aliphatics, Part 1. Biostimulation of methanotropic bacteria. Ground Water, 29: 365-374. Semprini, L. and McCarty, P.L., 1992. Comparison between model simulations and field results for in-situ biorestoration of chlorinated aliphatics, Part 2. Cometabolic transformations. Ground Water, 30: 37-44. Sudicky, E.A., 1986. A natural gradient experiment on solute transport in a sand aquifer; Spatial variance of hydraulic conductivity and its role in the dispersion process. Water Resour. Res., 22: 2069-2082. Therrien, R. and Sudicky, E.A., 1995. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. J. Contain. Hydrol. (submitted). Valocchi, A.J., Malmstead, M., 1992. Accuracy of the operator splitting for advection-dispersion-reaction problems. Water Resour. Res. 28: 1471-1476. Widdowson, M.A., Molz, F.J. and Beneficial, L.D., 1988. A numerical transport model for oxygen- and nitrate-based respiration linked to substrate and nutrient availability in porous media. Water Resour. Res., 24: 1553-1565. Zeyer, J., Kuhn, E.P. and Schwarzenbach, R.P., 1986. Rapid microbial mineralization of toluene and 1,3-dimethylbenzene in the absence of molecular oxygen. Appl. Environ. Microbiol., 52: 944-947.