Modelling of physical and reactive processes during biodegradation of a hydrocarbon plume under transient groundwater flow conditions

Journal of Contaminant Hydrology 59 (2002) 113 – 131 www.elsevier.com/locate/jconhyd

Modelling of physical and reactive processes during biodegradation of a hydrocarbon plume under transient groundwater flow conditions H. Prommer a,b,*, D.A. Barry a, G.B. Davis c a

Faculty of Civil Engineering and Geosciences, Department of Water Management, Delft University of Technology, Delft, The Netherlands b Centre for Applied Geoscience, University of Tu¨bingen, Germany c CSIRO Land and Water, Centre for Groundwater Studies, Private Bag No. 5, Wembley, WA 6014, Australia Received 1 May 2001; received in revised form 11 March 2002; accepted 12 March 2002

Abstract Numerical experiments of non-reactive and reactive transport were carried out to quantify the influence of a seasonally varying, transient flow field on transport and natural attenuation at a hydrocarbon-contaminated field site. Different numerical schemes for solving advective transport were compared to assess their capability to model low transversal dispersivities in transient flow fields. For the field site, it is shown that vertical plume spreading is largely inhibited, particularly if sorption is taken into account. For the reactive simulations, a biodegradation reaction module for the geochemical transport model PHT3D was developed. Results of the reactive transport simulations show that under the site-specific conditions the temporal variations in groundwater flow do, to a modest extent, affect average biodegradation rates and average total (dissolved) contaminant mass in the aquifer. The model simulations demonstrate that the seasonal variability in groundwater flow only results in significantly enhanced biodegradation rates when a differential sorption of electron donor (toluene) and electron acceptor (sulfate) is assumed. D 2002 Elsevier Science B.V. All rights reserved. Keywords: Reactive transport; BTEX; Dispersion; Sorption; Simulations; Natural attenuation; MT3DMS; MODFLOW; PHREEQC; PHT3D; Dispersion; Vertical mixing; Transverse mixing

* Corresponding author. CSIRO Land and Water, Centre for Groundwater Studies, Private Bag No. 5, Wembley, W.A. 6913, Australia. Fax: +61-8-9333-6211. E-mail address: [email protected] (H. Prommer).

0169-7722/02/$ - see front matter D 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 9 - 7 7 2 2 ( 0 2 ) 0 0 0 7 8 - 5

114

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

1. Introduction Numerical modelling is widely used to analyse and predict the risk associated with using natural attenuation as a remediation scheme. Numerous models that account for flow and three-dimensional reactive transport associated with biodegradation of oxidisable organic contaminants have become available over the past few years (e.g., Clement, 1997; Frind et al., 1999). With few exceptions, their application to field sites involve a vertically averaged two-dimensional simulation model, i.e., vertical concentration gradients over the aquifer depth are assumed negligible. However, as pointed out, for example, by Grathwohl et al. (2000), transversal vertical mixing can be an even more critical (physical) factor controlling the length of a naturally attenuating contaminant plume than the transversal horizontal mixing that is typically studied. Unfortunately, most contaminated sites are not well characterised with respect to the vertical distribution of the aqueous contaminants within an aquifer, so the identification of the dominant mixing processes (and how to model them) in the field is difficult. A detailed three-dimensional study is currently under way for one of the few locations for which more detailed information of hydrochemical parameters have been intensively recorded (e.g., Thierrin et al., 1995; Davis et al., 1999). At this hydrocarbon-contaminated site in Perth, Western Australia, toluene, ethylbenzene and xylenes are mineralised under sulfate-reducing conditions in a seasonally varying groundwater flow field. Benzene has shown to be persistent. The observed long, thin contaminant plumes within the relatively homogeneous sand aquifer indicate low transversal dispersivities, providing a significant challenge to achieve an accurate numerical description. In support of this three-dimensional study, vertical cross-sectional modelling is carried out along the (flow) path of the contaminants. This paper reports some of the initial numerical studies that have been carried out in order to (a) compare the numerical error that is introduced by the use of three different numerical schemes for solving advective transport when studying the problem with a reactive transport model based on the widely used MT3DMS simulator (Zheng and Wang, 1999) and (b) further understand the site-specific role of the transient flow field on reactive processes and, in particular, to investigate whether additional transverse dispersion resulting from the transient flow field is likely to influence natural attenuation rates significantly. Mass transport of non-reactive chemicals in transient flow fields has been studied previously (e.g., Kinzelbach and Ackerer, 1986; Goode and Konikow, 1990). Kinzelbach and Ackerer (1986) define the apparent transversal dispersivity at,app resulting from fluctuations in groundwater flow direction as: 2 2 V V l t at;app ¼ at  þ al  ; VV VV

ð1Þ

where al and at are the true longitudinal and transversal dispersivities, respectively, V is the groundwater velocity and Vl and Vt are its component in the mean flow and transverse directions, respectively, with overbars indicating time averages. Goode and Konikow (1990) show that the apparent dispersivities at,app for small changes in the angle of the

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

115

flow direction ( < 5j) lies not significantly above the actual dispersivity at, except for ratios of at/al < 0.1.

2. Accuracy of numerical schemes 2.1. Purely advective transport at an angle to grid orientation In order to study transport and mixing processes caused by a seasonally fluctuating flow-field, i.e., transient changes of water-table and flow velocities that cause the contaminant plume centre(s) to move vertically, errors due to numerical dispersion or oscillations need to be minimised. In (bio)reactive multi-species/multi-component models, excessive numerical dispersion might results in unrealistic artificial mixing between electron donors and electron acceptors, leading to local over-prediction of biodegradation, i.e., contaminant mass removal rates. Several of the relevant fundamental numerical issues have been discussed elsewhere (Kinzelbach, 1992; Zheng and Wang, 1999; Cirpka et al., 1999). Among them, advection-dominated transport with a plume direction that differs from the numerical grid-orientation was shown to be a critical issue. For a static grid, this presents an unavoidable problem under transient flow conditions. Thus, the first (schematic) test problems described here (Fig. 1) were set up to study the errors introduced by such ‘angular’ transport. Clearly, the chosen spatial model discretisation (here Dx = 5 m and Dz = 0.25 m, see Table 1) plays a key role in the absolute error produced. However, as three-dimensional multi-component transport simulations produce a significant computational burden, our goal in this test problem is rather to verify and compare the relative error introduced by three different widely used solution schemes for a specific, acceptable (in practice) grid resolution. The three advection solution schemes investigated were the upstream finite difference (FD) method with central weighting in space, the third-order total-variation-diminishing (TVD) scheme (Leonard, 1988) and the hybrid method of characteristics (HMOC) scheme (Zheng and Bennett, 1995), all of which are available

Fig. 1. Model setup used for test cases 1 and 2 (vertical cross-section, confined aquifer). A transversal flow qtr is superposed on the general, longitudinal flow direction. In test case 2, the magnitude of qtr is reversed in a cyclic manner to create a seasonally changing, transient flow field.

116

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

Table 1 Set of model parameters used in test cases 1 and 2 (comparison of MT3DMS advection schemes) Flow and physical transport parameters

Value, unit

Model dimensions lx  ly  lz Discretisation Dx  Dz Base flux in longitudinal direction qlong Transverse flu qtr * Hydraulic conductivity Porosity Transversal dispersivity at Longitudinal dispersivity al/Transversal dispersivity at Resulting pore water velocity in longitudinal direction

400  1  8 m 5  0.25 m 0.2 m3 day  1 m  1 0.02 m3 day  1 m  1 20 m day  1 0.2 Between 0 and 4 mm 10 1 m day  1

* Magnitude of qtr is varying.

within MT3DMS (v4.0). They are also used in many MT3DMS-based, related reactive transport modelling packages such as RT3D (Clement, 1997), SEAM3D (Waddill and Widdowson, 1998), BioRedox (Carey et al., 1999) and PHT3D (Prommer, 2002; Prommer et al., submitted for publication). As in the present study, MODFLOW (McDonald and Harbaugh, 1988) is typically used for the flow simulations underlying the transport modelling studies carried out with these models. In the first test case, a transversal flow qtr pushes a steady-state plume of a non-reactive chemical downwards (perpendicular to the main, i.e., longitudinal flow direction) before it is, further downstream, moved back upwards by a transversal flow  qtr (see Figs. 1 and 2). As can be seen in Fig. 3, the comparison of the simulated concentration profiles with the appropriate solutions obtained from simulations without any transversal flow (assumed to be free of numerical error), demonstrates that only the HMOC method is capable of producing a solution that is virtually free of numerical dispersion for any of the tested transversal dispersivities at (0.004, 0.04, 0.4 and 4 mm). In particular, the FD scheme introduces significant numerical dispersion, while the TVD scheme produces an acceptable solution. It must be noted that the postulated high accuracy of the HMOC scheme will only be warranted if the number of particles tracked is sufficiently high. A huge total number of particles will, on the other hand, require substantially more computer memory than the TVD scheme. However, the resulting memory requirements will nowadays, and even less likely in the future, rarely be a limiting factor for most practical applications.

Fig. 2. Concentration contours of the fully developed, steady-state plume for the case of purely advective transport (transversal dispersivity at = 0).

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

117

Fig. 3. Steady-state concentration profiles of fully developed plumes (after 500 days simulation time) at the effluent end of the model domain for simulations without transversal flow ( qtr = 0, solid lines) and simulations with transversal flow ( qtr p 0, dashed lines with +).

2.2. Seasonally changing flow-field in a confined aquifer Based on the above case, a second test problem was set up to investigate plume dispersion in a non-steady flow field. In this example, transient flow in a confined aquifer is induced by a cyclic reversion of the direction of the transversal fluxes qtr and  qtr (after 150, 365, 515, 730 days, etc.). Thus, the total transversal movement of the contaminant plume’s passage through the model domain is significantly increased, compared to the previous case. Fig. 4 (upper plot) shows a comparison of the resulting concentration

118

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

Fig. 4. Upper plot: Comparison of plume contours at the effluent end of the model domain after 1825 days simulation time under transient flow conditions for three different transversal dispersivities. Lower plot: Comparison of simulated vertical concentration profiles at x = 50, 200 and 400 m for three different MT3DMS advection schemes.

contours (C = 0.05C0 and C = 0.95C0) after 1825 days simulation time for one particular advection scheme (HMOC) and three different transversal dispersivities (at = 0, 0.04 and 0.4 mm). In addition, the lower plot of Fig. 4 shows a comparison of the simulated vertical concentration profiles at three different locations (x = 50, 200 and 400 m) for the latter case (HMOC scheme, at = 0.4 mm) in comparison with the results of the two other tested advection schemes (TVD and FD). The vertical concentration profiles at x = 50 m are not impacted by the transversal flow and the results for all three advection schemes agree with each other. However, the results for the profiles at the locations where the plume was

Fig. 5. Second spatial moment (z-direction) along the spreading plumes for different vertical transversal dispersivities (at = 0 mm: solid lines, at = 0.04 mm: dashed lines, at = 0.4 mm: dash-dotted lines) and for different numerical solution schemes for advective transport.

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

119

exposed to transversal flow and transient flow conditions (x = 200 m and x = 400 m) show, like in test case 1, that the use of the FD scheme introduces considerable numerical dispersion. A moment analysis was used to quantify and compare the additional plume spreading caused by numerical dispersion. In the case of a model grid with a constant layer thickness, the solute plume’s first spatial moment (centre of mass) in the z-direction –z1 is: X zij cij i¼1;nlay  z 1;j ¼ X ð j ¼ 1; ncol Þ; ð2Þ cij i¼1;nlay

where zij is the spatial coordinate (z-direction) of the grid cell centre located in the ith layer and the jth column (x-direction) of a spatially discretised domain with ncol columns and nlay layers and cij is the simulated concentration within this grid cell. The computation of the second spatial moment –z 2, serving as a measure of vertical transversal plume spreading, is done by: X ðzij   z 1;j Þ2 cij i¼1;nlay  X ð j ¼ 1; ncol Þ ð3Þ z 2;j ¼ cij i¼1;nlay

In a steady state flow-field with purely parallel (groundwater) flow, the (Fickian) dispersion model, as used in this study, must lead to a linear increase of –z 2. As the simulations using the HMOC scheme show (see Fig. 5), such a linear increase can indeed be observed in both cases where transversal dispersion (at = 0.04 mm and at = 0.4 mm) was included, indicating that the transversal flow does not significantly affect plume spreading. The spatial moment remains constant for the non-dispersive case (at = 0), confirming that the HMOC scheme can accurately describe purely advective transport in the transient groundwater flow-field. The FD solution is unsuitable for correctly simulating advective transport with the chosen discretisation for the low transversal dispersivities investigated here (Table 1).

3. Reactive single-species transport (sorption) Based on the hydrogeological properties of the field site (Davis et al., 1999), a more realistic, but still schematic vertical cross-sectional flow and transport model that captures the main hydrological characteristics of the field site was developed. In contrast to the schematic model shown in Fig. 1, the constant head boundary condition at the effluent end of the model domain was changed such that the aquifer was modelled as being largely unconfined (as it is in reality). The transversal flow qtr indicated in Fig. 1 was replaced by a seasonally acting recharge to the aquifer, which is applied temporally. The model base was assumed to be a no-flow boundary. Recharge and inflow were chosen to represent schematically the typical Perth climate of a 5-month rainy period (recharge = 1.5 mm

120

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

day  1) and a 7-month dry period (no recharge). The model area was again discretised into grid elements of 0.25 m thickness (vertical extension) but with a variable discretisation in flow direction (Dx between 2.5 and 10 m). The simulated annual variation in piezometric heads (between 18.34 and 20.16 m relative to Australian Height Datum—AHD) agrees well with the measured annual variability (Davis et al., 1999). For simplicity, the piezometric head at the effluent end was kept constant although small head variations were observed in the field. Recharge and boundary inflows were averaged over one annual cycle for the corresponding steady state flow model. The HMOC scheme was used for all model runs. All (physical) parameters that define the flow-field of the modelling problem are listed in Table 2. The computed flow field is now, on average, characterised by steadily increasing flow velocities between the upstream model boundary and the effluent boundary. In the transient simulations, water table fluctuations of up to 1.66 m occur at the influent boundary. Fig. 6 shows simulated concentration (depth) profiles for the steady state and transient solutions for both, the non-reactive case and the case where sorption was included. Note that the concentration profiles resulting from the steady state flow field are not influenced by sorption. Fig. 7 shows (for one annual cycle) the computed first and second spatial moments and the normalised, vertically integrated mass along the main flow direction. In both figures, it can be seen that, as a consequence of the increasing flow velocity, the vertically integrated mass successively decreases between influent and effluent end. The inclusion of sorption reactions (linear, equilibrium) has a ‘buffering’ effect for the temporal variability of the vertically integrated dissolved mass, first spatial moment (vertical position of plume centre) and second spatial moment (plume spreading). The plots also show that for the transversal dispersivity used in these simulations (0.4 mm), the accelerating flow has a converging effect for the plume width. This convergence results from the decreasing saturated thickness of the aquifer and its effect eliminates the plume spreading caused by hydrodynamic dispersion. As can be seen in Fig. 7, the effect leads

Table 2 Set of model parameters used in the reactive transport simulations Flow and physical transport parameters

Value, unit

Model dimensions lx  ly  lz Discretisation Dx  Dz Base flux in longitudinal direction qlong in rain season Base flux in longitudinal direction qlong in dry season Average base flux in longitudinal direction qlong (steady state model) Groundwater recharge in rain season Groundwater recharge in dry season Average groundwater recharge (steady state model) Hydraulic conductivity Porosity Longitudinal dispersivity al Transversal dispersivity at Resulting pore water velocity in longitudinal direction

400  1  8 m 2.5 – 10  0.25 m 0.25 m3 day  1 m  1 0.125 m3 day  1 m  1 0.176 m3 day  1 m  1 1.5  10  3 m  1 day  1 0 6.2  10  4 m  1 day  1 22 m day  1 0.2 4 mm 0, 0.4 mm Variable, on average approximately 1 m day  1

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

121

Fig. 6. Simulated concentration for advective, dispersive, non-reactive transport (upper row) and, in comparison, advective, dispersive transport with linear equilibrium sorption (bottom row, retardation factor R = 1.88). Thick lines represent the solution for steady state flow, thin lines represent solutions from the transient simulations.

Fig. 7. First moment, vertically integrated mass and second moment for (a) advective, non-reactive transport; (b) advective, dispersive, non-reactive transport and (c, d) advective, dispersive transport with linear equilibrium sorption (retardation factors R = 1.45 and 1.88, respectively). Thick lines represent the solution for steady state flow, thin lines represent solutions from the transient simulations.

122

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

Table 3 Parameters used for the simulation of kinetically controlled dissolution and biodegradation of toluene Parameter

Value

KSulf KToluene KN Kinhib,bio vmax, vdecay Ctol,solu

2.0  10  3 mol l  1 1.0  10  3 mol l  1 1.0  10  5 mol l  1 1.0  10  2 mol l  1 17.2 day  1 0.1 d  1 1.0  10  2 day  1 3.5  10  4 mol l  1

to a reduction of the second spatial moments at the model effluent boundary, compared to the influent boundary, where the mass is undispersed (Table 3).

4. Reactive multi-component transport Finally, the influence of the transient flow field on reactive processes was studied. Biodegradation of oxidisable contaminants, in our case benzene, ethylbenzene, toluene and xylenes (BTEX), is known to occur mainly at contaminant plume fringes where both organic substrates and electron acceptors (here sulfate) are available. Thus, it might be expected that the cyclic transversal movement of the plume, caused by the seasonal water table fluctuations, would increase transversal dispersion and mixing of the substrates and electron acceptor. This would then possibly imply increased biodegradation rates, changing the shape and reduce the total mass of the contaminant plume(s). In particular, this might be expected when electron donor(s) and electron acceptor(s) have different sorption characteristics, which leads to a difference in the relative mobility of these main reactants and chromatographic mixing might occur. In order to quantify this effect, the reactive multicomponent transport model PHT3D (Prommer et al., submitted for publication) was used to compare the characteristics of contaminant plumes in the transient flow-field (similar to the one described in the previous section) with the corresponding cases of steady state plumes. 4.1. Reactive transport model 4.1.1. Model description The PHT3D model couples the three-dimensional transport simulator MT3DMS (Zheng and Wang, 1999) with the geochemical model PHREEQC-2 (Parkhurst and Appelo, 1999). As PHREEQC-2 is capable of solving simultaneously both equilibrium and kinetically controlled reactions, the coupled model can address a wide range of mixed equilibrium and/or kinetic reactive transport problems. The reactive transport equation solved by the coupled model for the (mobile) aqueous species/components (in indical notation) is:   BC B BC B ¼ Dij ðvi CÞ þ rreac ; ð4Þ  Bt Bxi Bxj Bxi

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

123

where vi is the pore velocity in direction xi, Dij is the hydrodynamic dispersion coefficient tensor, rreac is a source/sink rate due to chemical reaction and C is the total aqueous component concentration (e.g., Westall et al., 1976; Yeh and Tripathi, 1989; Engesgaard and Kipp, 1992). The (local) redox-state, pe, is modelled by transporting chemicals/ components in different redox states separately. The pH is modelled from the (local) charge balance. The coupling of the two (sub)models is achieved through a modified (Prommer et al., submitted for publication) sequential operator-splitting technique (Herzer and Kinzelbach, 1989; Yeh and Tripathi, 1989; Steefel and MacQuarrie, 1996; Barry et al., 1996a,b, 1997). The advection and dispersion terms within Eq. (4) are computed for all mobile chemicals by the transport simulator MT3DMS. The reaction term rreac in Eq. (4) is computed by PHREEQC-2 for all, i.e., mobile and immobile chemicals. The coupled model has been successfully benchmarked against a series of reactive multi-species and multi-component transport problems from the literature (Prommer, 2002; Prommer et al., submitted for publication). 4.1.2. Biodegradation of organic compounds A site-specific PHT3D reaction module for toluene degradation under sulfate-reducing conditions as the key reactive process was developed for this study. As recently discussed by Brun and Engesgaard (2002), biodegradation models can be classified, in a general way, into so-called one-step or two-step process models. The latter class of models, into which the present model can be allocated, is characterised by the assumption that the electron flow during biodegradation of an oxidisable organic compound, e.g., of toluene under sulfate-reducing conditions might be described by two (sub)steps. The overall reaction for this two-step reaction is: þ  C7 H8 þ 4:5SO2 4 þ 3H2 O þ 2H ! 7HCO3 þ 4:5H2 S:

ð5Þ

In the first step, the oxidation of an organic substrate that has diffused into a biofilm, electrons are transferred to electron carriers such as reduced nicotinamide adenine dinucleotide (NADH) (Brock et al., 1994; Rittmann and van Briesen, 1996). For example, the complete oxidation (mineralisation) of toluene coupled to the reduction of NAD + is: C7 H8 þ 18NADþ þ 21H2 O ! 18NADH þ 25Hþ þ 7HCO 3:

ð6Þ

The electrons gained in this step can now be further transferred (Brock et al., 1994) 2 either to extracellular electron acceptors (O2, NO 3 , SO4 , etc.) or, alternatively, used for the formation of new biomass. For the former case, this involves the (re)oxidation of NADH, e.g., for sulfate: þ  þ SO2 4 þ 4NADH þ 5H ! HS þ 4H2 O þ 4NAD þ free energy:

ð7Þ

The free energy gained in these reactions can then be stored as adenosine triphosphate (ATP) and, together with NADH, reinvested to generate new biomass.

124

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

For the latter case, assuming a simplified chemical composition for biomass (C5H7O2N), the reaction in which the electron flow is diverted towards the generation of new biomass is described by (Rittmann and van Briesen, 1996): þ 10NADH þ 5H2 CO3 þ NHþ 4 þ 9H þ free energy

! C5 H7 O2 N þ 10NADþ þ 13H2 O:

ð8Þ

The efficiency of the microorganisms determines the ratio at which electrons gained in the oxidation step (6) are diverted between the electron-accepting step (7) and biomass generation according to Eq. (8). For an assumed efficiency em ¼ organic carbon converted to cell carbon=organic carbon degraded ¼ 0:1; ð9Þ and combining Eqs. (6), (7) and (8) leads to: 2 þ C7 H8 þ 0:14NHþ 4 þ 4:15SO4 þ 2:58H2 O þ 1:86H

! 6:30HCO 3 þ 0:14C5 H7 O2 N þ 4:15H2 S:

ð10Þ

Note that Eq. (10) reflects the stoichiometry during microbial growth. Comparing Eqs. (5) and (10) demonstrates that less sulfate is needed under growth conditions (4.15 versus 4.5 mol of sulfate per mole of toluene degraded). To explain this discrepancy, we look at the valence and concentrations of the redox-sensitive reactants before and after the reaction during the complete mineralisation of toluene. Electron balance requires: YC7 H8 ovC7 H8 þ YSO2 ovSO2 ¼ YHCO3 ovHCO3 þ YH2 S ovH2 S ; 4 4

ð11Þ

and YH2S are the stoichiometric factors for toluene, sulfate, where YC7H8, YSO42, YHCO 3 bicarbonate and hydrogen sulfide, respectively, and ovC7H8, ovSO2 , ovHCO3 and ovH2S are 4 the appropriate operational valences (Parkhurst et al., 1980; Brun, 1997; Brun and Engesgaard, 2002) of these species. Mass balance for sulfur requires that YSO2 = YH2S 4 and mass balance for carbon requires that YHCO3 = 7, i.e., the number of carbon atoms within toluene. For the second case, ammonium and bacteria need now also be considered in the electron balance: YC7 H8 ovC7 H8 þ YSO2 ovSO2 þ YNHþ4 ovNHþ4 4 4 ¼ YHCO3 ovHCO3 þ YC5 H7 O2 N ovC5 H7 O2 N þ YH2 S ovH2 S ;

ð12Þ

where YNH4+ , YC5H7O2N, ovNH4+ and ovC5H7O2N are the stoichiometric factors for ammonium, the stoichiometric factor for bacteria, the valence of ammonium and the valence of bacteria, respectively. However, as nitrogen balance requires that YNHþ4 ¼ YC7 H5 O2 N

ð13Þ

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

125

and bacteria and ammonium both have the same valence, the respective terms are equivalent. Thus, the difference, DYSO2 , in electron acceptor (sulfate) consumption 4 between the two cases described by Eqs. (5) and (10) can be predicted from DYSO2 ¼ 4

ð7  6:3ÞovHCO3 0:7  4 ¼ 0:35: ¼ 6  ð2Þ ovSO2  ov H2 S 4

ð14Þ

However, as we will see, the electron acceptor ‘‘saving’’ from this step, will be consumed by microbial decay. The reaction describing this decay and the appropriate consumption of electron acceptors is, for sulfate:   þ þ C5 H7 O2 N þ 3H2 O þ 2:5SO2 4 ! 5HCO3 þ 2:5HS þ 1:5H þ NH4 :

ð15Þ

It can be seen that the decay of 1 mol of bacteria consumes 2.5 mol of sulfate. Thus, the decay of the 0.14 mol of bacteria that were produced during growth, as described by Eq. (6), consumes 0.14  2.5 = 0.35 mol of sulfate, which equals the difference between sulfate consumption in Eqs. (5) and (10). Accordingly, only in a system with a steady state or a quasi steady state microbial concentration, i.e., where microbial growth and decay are balanced, does the stoichiometry resulting from Eq. (5) apply exactly. In contrast, it is not necessarily a correct approximation in cases with continuous local and temporal changes of bacterial activity. Postma and Jakobsen (1996) argue that if the first step (the oxidation step) is the ratelimiting step, the second step, the electron-accepting and biomass-generating step, can be modelled as an equilibrium reaction. Given this assumption, kinetic biodegradation reactions of oxidisable organic substances can be easily incorporated within the PHREEQC-2 framework as only the organic substance and the bacterial population need to be treated as kinetic species, whereas components that serve as electron acceptors or reaction products such as CO2 are assumed to be in geochemical equilibrium. The reaction kinetics for the bacterial population, i.e., the concentration of sulfate-reducing bacteria (SRB) is based on the mass balance equation: BXgrowth BXdecay BX ¼ þ ; Bt Bt Bt

ð16Þ

with the growth term in the present study being BXgrowth Ctolu Csulf CN ¼ vmax Yx Ibio X Bt Ktolu þ Ctolu Ksulf þ Csulf KN þ CN and the decay term within Eq. (16) being BXdecay ¼ vdec X ; Bt

ð17Þ

ð18Þ

where X, Ctolu, Csulf and CN are bacterial, toluene, sulfate and (total) nitrogen concentrations, vmax is an asymptotic maximum specific uptake rate, Ktolu, Ksulf and KN are the toluene, sulfate and nitrogen half-saturation constants, Yx is a stoichiometric factor, vdec is

126

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

a decay rate constant and Ibio is an inhibition term for microbial growth (Zysset et al., 1994). It is computed from Ibio ¼

hbio  X ; hbio

ð19Þ

whereby hbio is a maximum bacterial concentration. The toluene removal rate and the bacterial growth rate are proportional and the former can be computed from BXgrowth BCtolu ¼ Ytolu : Bt Bt

ð20Þ

The second, electron-accepting/cell growth step can be conveniently modelled in PHREEQC-2 as irreversible reaction (Parkhurst and Appelo, 1999; Brun and Engesgaard, 2002). For the discussed example (toluene oxidation under sulfate-reducing conditions), if the computed amount of toluene to be mineralised during a time step Dt is 1 Amol (for a particular grid-cell of a spatially discretised domain), an appropriate amount of C7H8 would be added to the aqueous composition and 0.14 Amol of C5H7O2N would be removed. During this reaction the aqueous composition and mineral concentrations will be adjusted, thereby consuming the most favourable electron acceptor present (here sulfate). The removal of 0.14 Amol of C5H7O2N from the aqueous solution represents the mass transfer of the elements C, H, O and N from the aqueous solution to the immobile ‘‘biophase’’, i.e., their incorporation into biomass. Accordingly, the bacterial mass will be increased by 0.14 Amol. In a similar way, the decay of bacteria is modelled as a release (addition) of C, H, O and N to the aqueous solution, accompanied by a decrease of the bacterial mass. 4.1.3. Initial and boundary conditions An anaerobic, equilibrated and charge-balanced groundwater composition, as sampled upstream of the contaminated area, provided the initial concentrations (uncontaminated aquifer) for the aqueous components included in the simulations (Table 4). The same water composition was assumed for the upstream (vertical) model boundary. At the field site, oxygen was occasionally observed in close proximity to the water table during times of recharge. Thus, for the modelling the recharge water composition was determined by equilibrating the initial water composition with, additionally, 3.2 mg l  1 of oxygen. As discussed above, sulfate-reducing bacteria and toluene were treated as kinetic species whereas for all other aqueous components chemical equilibrium was assumed. The complexation reactions and the appropriate reaction constants for the components included in the model simulation were taken from the standard PHREEQC-2 database. Benzene, ethylbenzene and xylene were excluded in these simulations as their influence on the electron-balance was expected to be small. 4.1.4. Contamination source Data from a multi-port borehole within the source zone indicated that residual NAPLs (nonaqueous phase liquids) over a zone 1 –1.5 m thick were located below the water table

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

127

Table 4 Initial concentrations of aqueous components and minerals of the uncontaminated aquifer used in all simulations Aqueousa component

C (mol l  1)

C (mol l  1)

pH pe O(0)b S(VI) S(  II) Fe(II) Fe(III) N(5) N(3) N(  3) C(IV) Ca Mg Na K Cl SRB

5.14 8.01 0.0 7.85  10  4 0.0 5.20  10  5 1.46  10  5 2.29  10  8 4.32  10  6 9.58  10  5 7.96  10  3 8.42  10  4 5.67  10  4 5.25  10  3 2.59  10  4 6.46  10  3 1.00  10  7

4.86 15.6 2.00  10  4 7.85  10  4 0.0 1.82  10  11 6.66  10  5 1.00  10  4 0.0 0.0 7.96  10  3 8.42  10  4 5.67  10  4 5.25  10  3 2.59  10  4 6.46  10  3 n.a.

a b

Except SRB, which is not part of the aqueous phase. Values in parentheses indicate valence state.

(Davis et al., 1999). These data showed also that the zone from which elevated BTEX concentrations dissolved generally followed the movement of the water table. However, in periods with a rising water table it appears that the lower fringe of this smearing zone

Fig. 8. First moment, vertically integrated mass and second moment of toluene for reactive multi-component transport without sorption of toluene (left side) and reactive multi component transport with sorption of toluene (R = 1.88, right side). Thin lines represent different stages/results from the transient simulations during one annual cycle, thick lines represent the solution for steady state flow.

128

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

moved slower than the water table. In the modelling study, the latter effect was neglected whereas a ‘boundary condition’ that allows kinetically controlled NAPL dissolution from a defined, constant thickness below the moving water table (1 m) was incorporated using: 

BCtolu;napl ¼ -ðCtolu;solu  Ctolu Þ; Bt

ð21Þ

where Ctolu,napl is the immobile (NAPL) concentration, Ctolu is the aqueous concentration, Ctolu,solu is the estimated multi-component solubility of toluene and - is a rate transfer constant. 4.2. Steady state versus transient simulations For both the steady state and the transient case, the simulated reactive system produces thin plumes of toluene and a relatively thin zone of sulfate depletion, as observed at the field site by Davis et al. (1999). The seasonally changing flow and, in particular, the ‘moving’ contaminant source lead to a pronounced variability of contaminant dissolution, flow path (Fig. 8, first spatial moments), plume spreading (Fig. 8, second spatial moments) and resulting biogeochemical reaction rates. Other characteristic field observation that are

Fig. 9. Contour plots for toluene, sulfate and sulfate-reducing bacteria for annual cycle under transient flow condition.

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

129

Fig. 10. Relative total mass of dissolved toluene in the model domain during the first 1500 days simulation time.

well reproduced by the transient model are (i) the dipping of the toluene plume and the sulfate-depleted zone below the watertable during times of recharge and (ii) the strong seasonal fluctuations of the vertically integrated mass. This former effect can only be reproduced when very small transversal dispersivities are used. The (simulated) microbial activity is largely confined to the vicinity of the contamination source and, further downstream, to the mixing zone where sulfate and toluene occurrence overlaps. Fig. 9 shows, for a transient case, the simulated concentration contours of toluene, sulfate and sulfate-reducing bacteria over one annual cycle. As a result of the seasonally changing hydrology, the computed total (dissolved) toluene mass present in the model domain varies considerably during an annual cycle, as can be seen in Fig. 10. The intra-annual changes in total dissolved mass are, however, decreasing with increasing sorption of toluene. From this plot it can also be seen that it is the chromatographic effect resulting from differential sorption of electron donor and acceptor which leads to enhanced biodegradation of toluene. The total mass in the transient ‘no sorption’ case is oscillating around the (constant) value for the steady case whereas in the transient ‘sorption’ cases the average total dissolved mass is up to 30% below the total dissolved mass in the corresponding steady state case.

5. Conclusions A modelling study was carried out to further understand the role of transverse mixing caused by seasonal hydrological changes and, in particular, its role in enhancing biodegradation at a field site in Perth, Western Australia. Initially, a series of numerical experiments were carried out to assess the most suitable (MT3DMS) advection scheme to be used in order to minimise numerical dispersion. The results of these experiments show that, for the (i) particularly low transverse dispersivities occurring at the study site and (ii) the practical discretisation constraints, the HMOC scheme is the only one that does not

130

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

exhibit notable numerical dispersion. The scheme is computationally more demanding than the alternative FD or TVD schemes (Zheng and Wang, 1999). However, as the total computational cost of multicomponent transport problems is dominated by the reaction term solution, this appears to be a minor disadvantage. Next, non-reactive and equilibrium sorption-only simulations were carried out for a flow field that mimics the dominating hydrological characteristics of the field site. It was found that vertical plume spreading, defined by the second spatial moment, was decreasing as a result of the accelerating flow field, i.e., the reduction of the vertical extent of the aquifer along the flow-path. The subsequent inclusion of equilibrium sorption reactions caused a reduction in vertical plume movement (temporal change of first spatial moment) and, to a smaller extent, a reduction in the vertical plume spreading. For the subsequent reactive multicomponent simulations, a site-specific reaction module that includes microbially mediated toluene degradation and relevant secondary reactions was formulated for the PHREEQC-2 based reaction simulations (and used for the computation of the reaction term in the governing transport equation). The results of these multi-component simulations indicate that the seasonal groundwater fluctuation might cause an enhancement of biodegradation. However, it is the combined effect of reactive processes—in particular the inclusion of toluene sorption—rather than an increase in the physically based (apparent) transversal dispersion for toluene, which could lead to notably enhanced biodegradation. References Barry, D.A., Bajracharya, K., Miller, C.T., 1996a. Alternative split-operator approach for solving chemical reaction/groundwater transport models. Adv. Water Resour. 19, 261 – 275. Barry, D.A., Miller, C.T., Culligan-Hensley, P.J., 1996b. Temporal discretisation errors in non-iterative splitoperator approaches to solving chemical reaction/groundwater transport models. J. Contam. Hydrol. 22, 1 – 17. Barry, D.A., Miller, C.T., Culligan, P.J., Bajracharya, K., 1997. Analysis of split operator methods for nonlinear and multispecies groundwater chemical transport models. Math. Comput. Simul. 43, 331 – 341. Brock, T.D., Madigan, M.T., Martinko, J.M., Parker, J., 1994. Biology of Microorganism. Prentice Hall, NJ. Brun, A., 1997. Reactive transport modeling of coupled inorganic organic processes in groundwater. PhD thesis, Dept. of Hydrodynamics and Water Resources, Series Papers, 65, Technical University of Denmark, Lyngby, 184 pp. Brun, A., Engesgaard, P., 2002. Modelling of transport and biogeochemical processes in pollution plumes: model formulation and development. J. Hydrol. 256 (3 – 4), 211 – 227. Carey, G.R., Van Geel, P.J., Murphy, J.R., 1999. BioRedox-MT3DMS V2.0: A Coupled Biodegradation-Redox Model for Simulating Natural and Enhanced Bioremediation of Organic Pollutants—User’s Guide ConestogaRovers and Associates, Waterloo, Ontario, Canada. Cirpka, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol. 40, 159 – 182. Clement, T.P., 1997. RT3D—A modular computer code for simulating reactive multi-species transport in 3dimensional groundwater systems. Battelle Pacific Northwest National Laboratory Research Report, PNNL-SA-28967. Davis, G.B., Barber, C., Power, T.R., Thierrin, J., Patterson, B.M., Rayner, J.L., Wu, Q., 1999. The variability and intrinsic remediation of a BTEX plume in anaerobic sulphate-rich groundwater. J. Contam. Hydrol. 36, 265 – 290. Engesgaard, P., Kipp, K.L., 1992. A geochemical transport model for redox-controlled movement of mineral fronts in groundwater flow systems: a case of nitrate removal by oxidation of pyrite. Water Resour. Res. 28, 2829 – 2843.

H. Prommer et al. / Journal of Contaminant Hydrology 59 (2002) 113–131

131

Frind, E.O., Molson, J.M., Schirmer, M., 1999. Dissolution and mass transfer of multiple organics under field conditions: the Borden emplaced source. Water Resour. Res. 35, 683 – 694. Goode, D.J., Konikow, L.F., 1990. Apparent dispersion in transient groundwater flow. Water Resour. Res. 26, 2339 – 2351. Grathwohl, P., Klenk, I.D., Eberhardt, C., Maier, U., 2000. Steady state plumes: Mechanisms of transverse mixing in aquifers. In: Johnston, C.D. (Ed.), Contaminated Site Remediation: From Source Zones to Ecosystems. Proc. 2000 Contaminated Site Remediation Conference, Melbourne, 4 – 8 December, pp. 459 – 466. Herzer, J., Kinzelbach, W., 1989. Coupling of transport and chemical processes in numerical transport models. Geoderma 44, 115 – 127. Kinzelbach, W., 1992. Numerische Methoden zur Modellierung des Transports von Schadstoffen in Grundwasser Oldenbourg Verlag, Mu¨nchen. Kinzelbach, W., Ackerer, P., 1986. Modelisation de la propagation d’un contaminant dans un camp d’ecoulement transitoire. Hydrogeologie 2, 197 – 205. Leonard, B.P., 1988. Universal limiter for transient interpolation modeling of the advective transport equations: the ULTIMATE conservative difference scheme. NASA Technical Memorandum 100916 ICOMP-88-11. McDonald, J.M., Harbaugh A.W., 1988. A Modular 3D Finite Difference Ground-Water Flow Model. Technical report, US Geological Survey techniques of Water-Resources Investigations. Parkhurst, D.L., Appelo, C.A.J., 1999. User’s guide to PHREEQC—A computer program for speciation, reaction-path, 1D-transport, and inverse geochemical calculations. Technical Report 99-4259, US Geol. Survey Water-Resources Investigations Report. Postma, D., Jakobsen, R., 1996. Redox zonation: equilibrium constraints on the Fe(III)/SO4-reduction interface. Geochim. Cosmochim. Acta 60, 3169 – 3175. Prommer, H., 2002. PHT3D—A reactive multi-component transport model for saturated porous media. Version 1.0 User’s Manual, Technical report, Contaminated Land Assessment and Remediation Research Centre, The University of Edinburgh. Prommer, H., Barry, D.A., Zheng, C., 2002. MODFLOW/MT3DMS-based reactive multi-component transport modelling. Ground Water, in press. Rittmann, B.E., van Briesen, J.M., 1996. Microbiological processes in reactive modeling. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media. Mineral Society of America. Steefel, C.I., MacQuarrie, K.T.B., 1996. Approaches to modeling of reactive transport in porous media. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media. Min. Soc. of America, pp. 83 – 129. Thierrin, J., Davis, G.B., Barber, C., 1995. A ground-water tracer test with deuterated compounds for monitoring in situ biodegradation and retardation of aromatic hydrocarbons. Ground Water 33, 469 – 475. Waddill, D.W., Widdowson, M.A., 1998. SEAM3D: A numerical model for three-dimensional solute transport and sequential electron acceptor-based bioremediation in groundwater, Technical report, Virginia Tech., Blacksburg, Virginia. Westall, J.C., Zachary, J.L., Morel, F.M.M., 1976. MINEQL: A computer program for the calculation of chemical equilibrium composition in aqueous systems. Tech. Note 18, Water Qual. Lab., Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge. Yeh, G.T., Tripathi, V.S., 1989. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour. Res. 25, 93 – 108. Zheng, C., Bennett, G.D., 1995. Applied Contaminant Transport Modeling: Theory and Practise Wiley, New York, 440 pp. Zheng, C., Wang, P.P., 1999. MT3DMS: A modular three-dimensional multispecies model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; Documentation and User’s Guide. Contract Report SERDP-99-1, US Army Engineer Research and Development Center, Vicksburg, MS. Zysset, A., Stauffer, F., Dracos, T., 1994. Modelling of reactive groundwater transport governed by biodegradation. Water Resour. Res. 30, 2423 – 2434.