Cellular Autometa Matlab algorithm

Computers and Mathematics with Applications 56 (2008) 645–656 www.elsevier.com/locate/camwa

Mathematical modeling of bioremediation of trichloroethylene in aquifers

Benito M. Chen-Charpentier a,∗, Hristo V. Kojouharov b

a Department of Mathematics, University of Wyoming, Laramie, WY 82071-3036, United States b Department of Mathematics, University of Texas at Arlington, Arlington, TX 76019-0408, United States

Abstract

Trichloroethylene (TCE) is a very common contaminant of groundwater. It is used as an industrial solvent and is frequently poured into the soil. There exist bacteria that can degrade TCE. In contrast with most cases of bioremediation, the bacteria that degrade TCE do not use it as a carbon source. Instead the bacteria produce an enzyme to metabolize methane. This enzyme can degrade other organics including TCE. In this paper we model in situ bioremediation of TCE in an aquifer by using two species of bacteria: one that forms biobarriers to restrict the movement of TCE and the second one to reduce TCE. The model includes flow of water, transport of TCE and the nutrients, bacterial growth and degradation of TCE. Nonstandard numerical methods are used to discretize the equations. Some results are presented. c© 2008 Elsevier Ltd. All rights reserved.

Keywords: Multi-species biofilms; Trichloroethylene; TCE; Porous media; Bioremediation

1. Introduction

Trichloroethylene (TCE) is a very widely used industrial and commercial solvent. It is an unsaturated aliphatic chlorinated hydrocarbon chemical compound. Before its toxic properties were discovered it was widely used as an anesthetic and analgesic. Today it is still widely used by the dry cleaning industry and as a degreaser in many applications, including decaffeinating coffee and extracting vegetable fats. Unfortunately TCE finds its way to groundwater deposits, where it mixes with water and is a serious source of contamination. Even without having to deal with the added complication of the porous medium, the elimination of TCE from the environment is an expensive and complicated process. Recently [1] the use of biofilters was proposed to degrade this harmful compound. Many researchers have reported that cometabolism is responsible for TCE degradation. Different species of microorganisms, more specifically methanotrophs (bacteria that metabolize methane), use methane as a substrate. In order to metabolize the methane, these bacteria produce an enzyme called MMO (methane mono-oxygenase). However, this enzyme can also degrade other organic compounds, including TCE. So, by stimulating the bacteria with air and methane, they produce mono-oxygenase that will also degrade the TCE [1]. However, the bacteria can also be stimulated with phenol or toluene and also ammonia oxidizers, instead of methane, in order to produce enzymes that degrade TCE.

∗ Corresponding author. E-mail addresses: bchen@uwyo.edu (B.M. Chen-Charpentier), hristo@uta.edu (H.V. Kojouharov).

0898-1221/$ - see front matter c© 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2008.01.007

646 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

Under aerobic conditions microorganisms utilizing a primary substrate such as methane, phenol or toluene and also ammonia oxidizers produce non-specific oxygenases that degrade TCE [2–4]. Experiments carried out in anaerobic conditions also revealed the possibility of TCE degradation to ethane by means of methanol-degrading bacteria [5].

Biofilm bacteria have been successfully used for forming biobarriers to control the propagation of contaminants in groundwater. In [6] we have simulated the growth of biofilms and the formation of biobarriers in porous media using nonstandard numerical methods. Experiments and numerical models show that while subsurface biobarriers substantially control the movement of contaminants, they do not reduce it to zero as is really desirable in practice. Recently there have been some experiments [7–9] where two different types of bacteria are combined to get better results. One type is a strong biofilm-forming bacteria and the other is a bacteria that reacts with the contaminant transforming it into harmless substances. The biofilm-forming bacteria are needed to form the biobarrier so the contaminant transport is reduced and also allows the contaminant-degrading bacteria to establish themselves in the biobarrier and therefore be almost immobile and efficiently destroy the contaminant as it flows by. In many cases the contaminant-degrading bacteria uses the contaminant directly as a nutrient. It may also use the nutrient used to grow the biofilm as an additional nutrient. Several different cases of the use of two bacteria to control pollutants in groundwater were simulated numerically in [10].

As mentioned above, one of the most common contaminants in groundwater is TCE. Now we have a different situation because it has been reported that cometabolism is responsible for TCE degradation. The microbes do not use TCE as a carbon or energy source, but the amount of TCE degraded depends directly on the amount of enzymes produced, which depends directly on the amount of methane or phenol consumed. According to [8], the toluene ortho- mono-oxygenase (TOM) pathway is responsible for the degradation of TCE by Burkholderia cepacia and is located on the TOM plasmid. The authors also calculated the rate of TCE degradation per substrate consumption. There have also been some large scale demonstrations of the use of in situ bioremediation of TCE done by Los Alamos National Laboratory [11,12].

In this paper we model the water flow, the transport of nutrients and TCE as well as the growth of biofilm-forming microbes and biodegradation microbes. We only consider the limiting nutrient in the simulations and assume that all the other necessary nutrients are present in enough quantities. We will study the interaction between the two types of bacteria, the TCE and the nutrients. We are specially interested in the persistence of the biobarrier with the degrading bacteria. Under many conditions two species of bacteria using the same resources, one species will dominate and the second will eventually disappear. We want to see whether this undesirable situation will happen. The paper is organized as follows. In the next section we present the mathematical model. The model includes flow of water, transport of nutrients and TCE and growth and degradation. In Section 3 we show qualitative results of some dual- species biobarrier simulations. The purpose and value of the numerical simulations is to guide future multi-species biofilm experiments that could lead to the design of more effective bioremediation strategies. In the last section we present some conclusions and future research directions.

2. The mathematical model

In order to model multi-species biofilm interactions in porous media we consider a three-phase mixture consisting of a liquid phase, a solid rock phase and a biofilm phase. Even though the biofilm can be considered to be part of the solid phase, it is simpler to take it as a separate phase. The seven molecular species present in the porous medium are the contaminant, for example, trichloroethylene (TCE), the the biodegradation microbe unable to form a significant biofilm, which could be, for example, Burkholderia cepacia microorganism (a TCE-degrading bacteria unable to form a significant biofilm), the strong biofilm-forming bacteria, which could be, for example, Klebsiella oxytoca, the two substrates (which could be methane, phenol or toluene), and the water and rock species. The formulations used in the rest of the paper are valid for diverse species of microbes and chemical substances. In the rest of the paper, we will usually refer to these particular microbes but, of course, the formulations are valid for other microbial species (and substances).

The biofilm-forming microbes use the nutrients to grow and to form extra-polymer substances (EPS). The EPS is used to link the microbes together to form the biofilm. One modeling possibility is to consider the EPS an additional species, but this method increases the number of equations and the complexity of the problem. Since the biofilm microbes use a percentage of the available nutrients to form EPS, a second alternative is to consider the microbes and

B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656 647

the EPS as one species. The growth parameters for these microbes have to be adjusted accordingly. We will follow the second way and consider that the mass of the biofilm-forming microbes includes the EPS.

The fundamental equation for saturated transient groundwater flow of constant density, in horizontal direction, can be written in the form [13]:

Ss ∂h

∂t −

∂

∂ x

( K

∂h

∂ x

) = f. (1)

The single fluid-flow equation (1) arises from the mass balance law

Ss ∂h

∂t +

∂v

∂ x = f, (2)

when we substitute for the specific discharge vector v using Darcy’s law

v = −K ∂h

∂ x . (3)

Here h(x , t ) denotes the hydraulic head, Ss is the specific storage, K is the saturated hydraulic conductivity, and f (x , t ) represent sources or sinks. The specific discharge vector v(x , t ), called superficial or Darcy velocity, represents the speed of the water. We assume that there are no sources and sinks for the fluid, therefore f = 0 in Eq. (1). We also assume a piecewise steady-state fluid flow, due to the relatively slow changes in the porous media properties [14]. Also we are modeling very short cores with uniform biofilm distribution so we can take the velocity to be independent of x [14]. Invoking the above simplifying assumptions to Eq. (1) yields:

− ∂

∂ x

( K

∂h

∂ x

) = 0. (4)

The transport and reaction of nutrients and the contaminant, and the growth of the two microbial species are governed by a system of partial differential equations [13]. We assume that the two types of microbes are immobile, as part of the dual-species biofilm structure. Since the rock phase does not change we assume that the solid rock matrix is stationary and that the diffusion of the two microbial, nutrient and contaminant species in the solid phase is negligible. Therefore we can work only with the liquid and biofilm phases:

∂

∂t (φ

Bi o ρB ) = r B (ρB , ρK , ρN , ρM , ρT ),

∂

∂t (φ

Bi o ρK ) = r K (ρB , ρK , ρN , ρM , ρT ),

∂

∂t (φ

L ρN ) +

∂

∂ x (vρN ) −

∂

∂ x

( D

∂ρN

∂ x

) = r N (ρB , ρK , ρN , ρM , ρT ),

∂

∂t (φ

L ρM ) +

∂

∂ x (vρM ) −

∂

∂ x

( D

∂ρM

∂ x

) = r M (ρB , ρK , ρN , ρM , ρT ),

∂

∂t (φ

L ρT ) +

∂

∂ x (vρT ) −

∂

∂ x

( D

∂ρT

∂ x

) = rT (ρB , ρK , ρN , ρM , ρT ).

(5)

Here ρi , i = B, K , N , M, T , represents the intrinsic mass density of the biodegradation microbes, the strong biofilm-forming microbes, the two nutrients, and the contaminant, respectively. For a single-fluid flow, the quantity φ L = VL /(VL + VBi o) and the quantity φ

Bi o = VBi o/(VL + VBi o), where VL and VBi o represent the volumes

occupied by the liquid and by the biofilm, respectively, are the portions of the void space occupied by the biofilm and the liquid, D is the hydrodynamic dispersion coefficient, and ri represents the total rate at which species i is produced via reactions and sources.

Let X B and X K be the current biodegradation and biobarrier-forming microbial concentrations, respectively, then X̃ B = X B /ρB and X̃ K = X K /ρK are the corresponding normalized microbial concentrations. We assume that the growth and accumulation of both microbial species in the pore spaces cause changes in the porous media properties. It follows that the changes in porosity and in saturated hydraulic conductivity, for small initial biobarrier-forming

648 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

microbial concentrations [15], are given by

φ(X̃ B , X̃ K ) = φ0(1 − X̃ B − X̃ K ),

K (X̃ B , X̃ K ) = K0(1 − X̃ B − X̃ K ) nk ,

(6)

respectively, where φ0 is the clean surface porosity, K0 is the initial hydraulic conductivity and nk is an experimentally determined parameter which takes values around 3 [15]. Furthermore, we assume that the two microbial death rates are proportional to the size of the corresponding microbial populations.

For simplicity, from now on we drop the tilde from the normalized microbial concentration. We examine four different models based on the types of interactions between the microbial and nutrients species.

2.1. Dual-species biobarrier model I

First, we assume that only one species of nutrients, M , is present in the system, for example methane, phenol, or toluene, that is needed in order for the biodegradation microbes, B, to produce the TCE-degrading enzyme. Also, direct interactions in the system are assumed to occur only between the strong biofilm-forming microbial and the nutrients species, and also between the biodegradation microbial and the nutrients and contaminant species.

Incorporation the above simplifying assumptions yields the following governing system of differential equations:

∂ X B ∂t

= µ B (SM )X B G(X B + X K ) − k B X B

∂ X K ∂t

= µ K (SM )X K G(X B + X K ) − k K X K

∂ SM ∂t

+ v ∂ SM ∂ x

− ∂

∂ x

( D

∂ SM ∂ x

) = −

1 YB

µ B (SM )X B G(X B + X K ) −

1 YK

µ K (SM )X K G(X B + X K ),

∂ ST ∂t

+ v ∂ ST ∂ x

− ∂

∂ x

( D

∂ ST ∂ x

) = −

R

YB µ

B (SM )X B G(X B + X K ),

(7)

where R is the rate-constant of TCE degradation per substrate, M , consumed by the biodegradation microbes B, µ j ( j = B, K ) describes the kinetics of microbial transformations of nutrient and/or contaminant species, YB and YK are the yield rate coefficients [16], X K is the normalized concentration of biobarrier-forming microbes, X B is the normalized concentration of biodegradation microbes, SM is the concentration of the limiting nutrient, and ST is the concentration of the contaminant. The coefficients k B and k K represent the first-order microbial decay rates and account for the decay processes that diminish the active biomass and for the desorption of the microbes caused by the shear force of the surrounding liquid.

The function G represents the fraction of daughter cell of adherent bacteria that find attachment sites [17]. It is introduced to restrict the growth of microbes as the pores are being plugged and is reasonable to assume that it is a decreasing function, because a more fully saturated wall provides less chance for a daughter cell to find space on it. For example, Freter [18,19] employs

G(X ) = 1 − X

1 − X + γ ,

with γ typically small. One common assumption in Monod’s growth rate for bacteria is that there is only one nutrient that limits the

growth [20]. That is, there is an excess of the other nutrients. In this case the rate of substrate utilization has been usually described by [16]:

µ j (Si ) = µ

j,i max

Si

K j Si

+ Si ,

j = B, K , where Si is the concentration of the limiting nutrient. For simplification, we assume that the effects of the adsorption process, which is aided by the extracellular polymer substances (EPS), are implicitly incorporated into the growth rate function through the maximum specific bacterial growth rate µ

j,i max. To account for the experimental

B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656 649

observation that biofilms do not degrade a substrate that is present below a certain threshold concentration Smi n [21, 22], we modify the Monod expression for the growth rate to force µ(S) to zero when S − Smi n < 0. In this case, the following modified Monod equation [23] is used:

µ j (Si ) =

µ j,i max

2 Si

K j Si

+ Si

 1 + Si − S ji,min∣∣∣Si − S ji,min∣∣∣

  , (8)

where S j i,min represents the minimum/threshold concentration of substrate i below which the energy that species j

gains from the metabolism is less than the energy required to metabolize the substrate i .

2.2. Dual-species biobarrier model II

We assume both nutrients species are present in the system, for example methane, phenol, or toluene, and also another nutrient that is used to grow the biofilm-barrier. Direct interactions in the system are assumed to occur only between the strong biofilm-forming microbial and the biofilm-growth nutrients species, N , and also between the biodegradation microbial and the other nutrients, M , and the contaminant species.

Incorporation of the above simplifying assumptions yields the following governing system of differential equations:

∂ X B ∂t

= µ B (SM )X B G(X B + X K ) − k B X B

∂ X K ∂t

= µ K (SN )X K G(X B + X K ) − k K X K

∂ SN ∂t

+ v ∂ SN ∂ x

− ∂

∂ x

( D

∂ SN ∂ x

) = −

1 YK

µ K (SN )X K G(X B + X K )

∂ SM ∂t

+ v ∂ SM ∂ x

− ∂

∂ x

( D

∂ SM ∂ x

) = −

1 YB

µ B (SM )X B G(X B + X K ),

∂ ST ∂t

+ v ∂ ST ∂ x

− ∂

∂ x

( D

∂ ST ∂ x

) = −

R

YB µ

B (SM )X B G(X B + X K ),

(9)

where SN is the concentration of the biofilm-growth nutrients, and the other scalars/functions are the same as in System (7).

In our mathematical models we use Eq. (8) to describe the kinetics of microbial transformations of the nutrients. In many cases, the limiting nutrient may change with time and position. In this case the Monod kinetics are given by a product of terms of the form (8), as in Section 2.3. Another scenario is when two or more nutrients have a similar role. Then the Monod kinetics are given as the sum of above terms, as in Section 2.4.

2.3. Dual-species biobarrier model III

This case corresponds to the pollutant-reducing bacteria needing both substrates to grow with no unique limiting one, while the biofilm-forming bacteria needing only one of the nutrients, N , for growth. Incorporating the above simplifying assumptions yields the following governing system of differential equations:

∂ X B ∂t

= µ B

(SN , SM ) X B G(X B + X K ) − k B X B

∂ X K ∂t

= µ K (SN )X K G(X B + X K ) − k K X K

∂ SN ∂t

+ v ∂ SN ∂ x

− ∂

∂ x

( D

∂ SN ∂ x

) = −

1 YK

µ K (SN )X K G(X B + X K ) −

F

YB µ

B (SN , SM ) X B G(X B + X K )

∂ SM ∂t

+ v ∂ SM ∂ x

− ∂

∂ x

( D

∂ SM ∂ x

) = −

1 YB

µ B

(SN , SM ) X B G(X B + X K )

∂ ST ∂t

+ v ∂ ST ∂ x

− ∂

∂ x

( D

∂ ST ∂ x

) = −

R

YB µ

B (SN , SM ) X B G(X B + X K ),

(10)

650 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

Table 1 Parameter values used in the numerical experiments

Parameters Values

Initial saturated hydraulic conductivity, K0 0.865 cm/h Initial porosity, φ0 0.35 Hydrodynamic dispersion coefficient, D 0.4 cm2/h Parameter, nk 3 Biofilm microbes decay coefficient, k K 0.0072 /h Degrading microbes decay coefficient, k B 0.005 /h

Biofilm microbes maximum specific growth rate, µK ,Mmax 0.095 /h

Biofilm microbes maximum specific growth rate, µK ,Nmax 0.144 /h

Biofilm microbes maximum specific growth rate, µB,Mmax 0.072 /h

Biofilm microbes maximum specific growth rate, µB,Nmax 0.02 /h Half saturation constant, K KSM

0.3 µg/ml

Half saturation constant, K KSN 1.5 µg/ml

Half saturation constant, K BSM 0.4 µg/ml

Half saturation constant, K BSN 0.5 µg/ml

Threshold substrate concentration, S KM,min 0.1 µg/ml

Threshold substrate concentration, S KN ,min 0.25 µg/ml

Threshold substrate concentration, S BM,min 0.07 µg/ml

Threshold substrate concentration, S BN ,min 0.08 µg/ml Biofilm microbes yield coefficient, YK 0.16 Degrading microbes yield coefficient, YB 0.33 Ratio constant, F 0.5 Degradation constant, R 0.3 Parameter, γ 0.1

where F is the ratio of nutrients N to nutrients M consumed, and µB (SN , SM ) has the following form:

µ B

(SN , SM ) = µBmax

4

∏ i =N ,M

Si K BSi + Si

 1 + Si − S Bi,min∣∣∣Si − S Bi,min∣∣∣

  . (11)

Here it is assumed that both nutrients have a threshold S j i,min that can be different for each species of bacteria.

2.4. Dual-species biobarrier model IV

In this case, we assume that the biofilm-forming bacteria needs only one of the nutrients, N , and that the pollutant- reducing bacteria can use either of the two nutrients for growth. Incorporating the above assumptions yields the following governing system of differential equations:

∂ X B ∂t

=

( µ

B (SN ) + µ

B (SM )

) X B G(X B + X K ) − k B X B

∂ X K ∂t

= µ K (SN )X K G(X B + X K ) − k K X K

∂ SN ∂t

+ v ∂ SN ∂ x

− ∂

∂ x

( D

∂ SN ∂ x

) = −

1 YK

µ K (SN )X K G(X B + X K ) −

F

YB µ

B (SN )X B G(X B + X K )

∂ SM ∂t

+ v ∂ SM ∂ x

− ∂

∂ x

( D

∂ SM ∂ x

) = −

1 YB

µ B (SM )X B G(X B + X K ),

∂ ST ∂t

+ v ∂ ST ∂ x

− ∂

∂ x

( D

∂ ST ∂ x

) = −

1 YB

µ B (SM )X B G(X B + X K ).

(12)

B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656 651

Fig. 1. Log-linear scale plots of the L 1 norms of X B , X K , SM , and ST vs. time at a low-nutrient supply (top) and at a high-nutrient supply (bottom) into the system (7) for the dual-species biobarrier model I.

3. Simulations

Eqs. (4) and (5) represent a coupled system of nonlinear, time-dependent ordinary and partial differential equations that is difficult to solve numerically. A key objective of the numerical simulation is to develop time-stepping procedures that are accurate and computationally stable. We use a sequential solution technique that first solves implicitly for the Darcy velocity v at the current time-level, by solving Eqs. (3) and (4). Then the species transport system (5) is solved implicitly for the concentrations X B , X K , SM , SN , and ST , in a decoupled fashion [24]. New values of the porosity and the hydraulic conductivity are then calculated using Eq. (6) and the cycle is repeated by calculating the new velocities.

For the solution of the ordinary differential Eq. (4) we use a standard finite-difference method to calculate h. Then we numerically differentiate using Darcy’s law (3) to get the velocity field v. The temporal differentiation in the microbial species Eq. (5) uses a forward Euler time integration. The nutrients and contaminants transport equations are solved using a nonstandard finite-difference (NSFD) numerical method [6,25,26]. Since we are mainly interested in the steady states of the system, the low-order accuracy of this method is not an issue.

Values of the parameters used in the numerical experiments are summarized in Table 1.

652 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

Fig. 2. Numerical simulation results for the dual-species biobarrier model II, System (9) — at a steady-state of the system (top) and a log-linear scale plot of the L 1 norms of X B , X K , SM , SN and ST vs. time (bottom).

The initial and boundary conditions used in the numerical simulations are:

h(0, t ) = 0.5 cm, h(1, t ) = 0 cm,

X B (x , 0) = X K (x , 0) = {

0.2, 0.3 ≤ x ≤ 0.4 0, otherwise,

SM (x , 0) = SM (0, t ) = 2.5 µg/ml, SN (x , 0) = SN (0, t ) = 3.5 µg/ml,

ST (x , 0) = ST (0, t ) = 0.2 µg/ml ∂ SM ∂ x

(1, t ) = ∂ SN ∂ x

(1, t ) = ∂ ST ∂ x

(1, t ) = 0 µg/(ml h),

except in the fourth set of simulations, System (12), where

SN (x , 0) = SN (0, t ) = 7 µg/ml.

The boundary and initial conditions considered in the biobarrier models are in agreement with the porous media experiments done by Cunningham et al., 1991 [14] for a 5 cm-long reactor packed with 0.70 mm, in diameter. Most of the reaction parameters are taken from Semprini and McCarty, 1991 and 1992 [27,28], Chang and Alvarez-Cohen,

B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656 653

Fig. 3. Numerical simulation results for the dual-species biobarrier model III, System (10) — at a steady-state of the system (top) and a log-linear scale plot of the L 1 norms X B , X K , SM , SN and ST vs. time (bottom).

1995 [29], Reddy and Ford, 1996 [23], Travis and Rosenberg, 1997 [30], Komlos et al., 2004 and 2006 [31,7], and the parameter γ in the function G is taken from Jones and Smith, 2000 [17]. We assume that the amount of TCE degraded depends directly on the amount of enzymes produced, which depends directly on the amount of methane or phenol consumed, and use and estimate 0.3 for the rate-constant, R, of TCE degradation per substrate consumed. The diffusion that we consider represents the average diffusion in the system, so we are assuming, in our simulations, that the coefficient D is constant. For ease of calculations, the reactor’s length has been scaled to 1. For graphing purposes, the time has been scaled by a factor of 250 and the nutrients concentrations, SM and SN , have been scaled by a factor of 1/5 in all figures, except in Fig. 4 where SN has been scaled by a factor of 1/10.

In the first set of experiments we simulate the case where both species of bacteria use only methane as their carbon source. Both species are competing for methane, so here we have a case where one bacteria is likely to dominate and therefore the degrading biobarrier will disappear as one of the species disappears. In our calculations we found that using SM (x , 0) = SM (0, t ) = 2.5 µg/ml, the biodegrading bacteria X B will eventually disappear (Fig. 1 (top)). By doubling the substrate concentration level (Fig. 1 (bottom)), i.e., using SM (x , 0) = SM (0, t ) = 5 µg/ml, leads to the biofilm-forming bacteria eventually disappearing. In both cases the degrading biobarrier disappears, but for different reasons.

654 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

Fig. 4. Numerical simulation results for the dual-species biobarrier model IV, System (12) — at a steady-state of the system (top) and a log-linear scale plot of the L 1 norms X B , X K , SM , SN and ST vs. time (bottom).

In the second set of numerical experiments we simulate the situation of the two species using different carbon sources. So there is no competition between the two species and both will survive and the biobarrier will persist for a very long time. The results show that all values stabilize as time increases. Fig. 2 shows the results for model II.

In the third and fourth set of numerical experiments we simulate the porous media experiments done by Komlos et al., [7–9] at the Center for Biofilm Engineering, Montana State University. In these experiments the pollutant- reducing bacteria can use either of the two nutrients as a carbon source, while the biofilm bacteria has only one carbon source. Figs. 3 and 4 show the results for model III and IV, respectively.

The dual-species simulation results qualitatively match actual experiment results done by Komlos et al., [7–9]. The presence of a threshold, which limits the amount of bacteria, reduces the amount of TCE eliminated as expected, but does not change the results significantly.

Which of the models is the most realistic, of course, depends on the species of bacteria used. In real life it is better to use bacteria that use different carbon sources since then the two species are more likely to survive and the biodegrading biobarrier will exist and reduce TCE for a long time.

In real life situations it is likely that the TCE will come into the system for a limited amount of time, instead of continuously as in our simulations. Our simulations show more extreme cases but give more information about the interaction of all species.

B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656 655

The TCE degradation potential of the dual-species biobarriers observed in our numerical simulations is as follows:

• 83% of TCE removed via biodegradation using model I; • 71% of TCE removed via biodegradation using model II; • 29% of TCE removed via biodegradation using model III; and • 51% of TCE removed via biodegradation using model IV.

Long-term numerical simulation results of contaminant degradation are in agreement with actual experiments presented in [30–32].

Even though the agreement with experimental results is fairly good, it should be noted that the values of the kinetic parameters reported in the literature have large variations, which could certainly produce different results. The amounts of biodegradation reported also vary very much. By changing the model we were able to obtain much different rates of degradation. Another point is that some parameters that are not given in the literature were estimated from the given ones.

4. Conclusions and future research

The motive for this research was to gain a better understanding of the interactions between two microbial species as part of a single biofilm capable of performing multiple functions (bioremediation, biofilm formation, etc.). We also wanted to study the amount of reduction of TCE and the effectiveness of the degrading biobarrier. We were also interested in observing the persistence of the barrier. We have presented a mathematical model for the flow, the transport of two nutrients and a contaminant, and the growth of two types of microorganisms in porous media, under different hypotheses. The coupled system of equations was solved numerically in a sequential way using mixed- finite elements for the flow and a nonstandard method for the transport equations. Nonstandard methods work well in conditions where there is transport and nonlinear reactions as in this case. And the coupling of the numerical methods produces a method for a very complicated problem that is stable and fairly accurate.

The dual-species biobarrier experiments show that it is possible to construct a degrading biobarrier and therefore enhance its TCE control potential. The experiments also show that these biobarriers can persist for long times. The four models emphasize the importance of choosing the bacteria species and the carbon sources carefully, as the models give very different rates of TCE elimination.

In real life, the biofilm would spread in the porous medium, which our mathematical model does not allow in its current form. In order to include this in the model, spatial transport terms in the equations for the biomass densities would be needed. Our future work will allow the degrading bacteria some movement, will let biofilm bacteria detach and will consider more spatial dimensions.

Acknowledgement

Supported in part by DOE Grant DE-FC26-06NT15568.

References

[1] A. Tabernacka, E. Zborowska, M. Lebkowska, Trichloroethene elimination from air by means of biofoltration, Polish Journal on Environmental Studies 15 (2006) 335–340.

[2] D.M. Bagley, Systematic approach for modeling tetrachloroethene bio-degradation, Journal on Environment Engineering 124 (1998) 1076–1086.

[3] S. Lee, W.M. Moe, K.T. Valsaraj, J.H. Pardue, Effect of sorption and desorption resistance on aerobic trichloroethene biodegradation in soils, Environmental Toxicology and Chemistry 21 (2002) 1609–1617.

[4] D.J. Arp, C.M. Yeager, M.R. Hyman, Molecular and cellular fundamentals of aerobic cometabolism of trichloroethene, Biodegradation 12 (2001) 81–103.

[5] M.G. Parvatyar, R. Govind, D.F. Bishop, Treatment of trichloroethene (TCE) in a membrane biofilter, Biotechnology and Bioengineering 50 (1996) 57–64.

[6] B.M. Chen, H.V. Kojouharov, Non-standard numerical methods applied to subsurface biobarrier formation models in porous media, Bull. Math. Biol. 61 (1999) 779–798.

[7] J. Komlos, A.B. Cunningham, A.K. Camper, R.R. Sharp, Effect of substrate concentration on the dual-species population density of Burkholderia cepatia PR1 − pTOM31c and Klebsiella oxytoca in Porous Media, Biotechnology and Bioengineering 93 (3) (2006) 434–442.

656 B.M. Chen-Charpentier, H.V. Kojouharov / Computers and Mathematics with Applications 56 (2008) 645–656

[8] J. Komlos, A.B. Cunningham, A.K. Camper, R.R. Sharp, Interaction of Klebsiella oxytoca and Burkholderia cepatia in dual-species batch cultures and biofilms as a function of growth rate and substrate concentration, Microbial Ecology (2005) Published online: 28 January.

[9] J. Komlos, A.B. Cunningham, A.K. Camper, R.R. Sharp, Varying substrate concentration to enhance TCE degradation in dual-species bioreactors, In: Sixth International In Situ and on Site Bioremediation Symposium, San Diego, CA, USA, 4–7 June 2001, pp. 117–124.

[10] B.M. Chen-Charpentier, D.T. Dimitrov, H.V. Kojouharov, Numerical simulation of multi-species biofilms in porous media for different kinetics, Mathematics and Computers in Simulation (in press).

[11] B.J. Travis, N.D. Rosenberg, Modeling in situ bioremediation of TCE at Savannah River: Effects of product toxicity and microbial interactions on TCE degradation, Environmental Science & Technology 31 (1997) 3093–3102.

[12] B.J. Travis, Numerical simulation of in situ bioremediation, in: V.N. Burganos, G.P. Karatzas, A.C. Payatakes, C.A. Brebbia, W.G. Gray, G.F. Pinder (Eds.), Computational Methods in Water Resources XII: Computational Methods in Contamination and Remediation of Water Resources, Computational Mechanics Publications, Southampton, UK, 1998, pp. 91–98.

[13] M.B. Allen, Basic mechanics of oil reservoir flows, in: C.A. Brebia, S.A. Orszag (Eds.), Multiphase Flow in Porous Media, M.B. Allen III, G.A. Behie, and J.A. Trangenstein, in: Lecture Notes in Engineering, vol. 34, Springer-Verlag, New York, NY, 1988, pp. 1–81.

[14] A.B. Cunningham, W.G. Characklis, F. Abedeen, D. Crawford, Influence of the biofilm accumulation on porous media hydrodynamics, Environmental Science & Technology 25 (7) (1991) 1305–1311.

[15] T.P. Clement, B.S. Hooker, R.S. Skeen, Microscopic models for the predicting changes in the saturated porous media properties caused by microbial growth, Ground Water 34 (5) (1996) 934–942.

[16] J.E. Bailey, D.F. Ollis, Biochemical Engineering Fundamentals, McGraw-Hill, New York, NY, 1986. [17] D. Jones, H. Smith, Microbial competition for nutrient and wall sites in plug flow, SIAM Journal of Applied Mathematics 60 (5) (2000)

1576–1600. [18] R. Freter, Mechanisms that control the microflora in the large intestine, in: D. Hentges (Ed.), Human Intestinal Microflora in Health and

Disease, Academic Press, New York, 1983, pp. 33–54. [19] R. Freter, H. Brickner, S. Temme, An understanding of colonization resistance of the mammalian large intestine requires mathematical

analysis, Microecology and Therapy 16 (1986) 147–155. [20] J. Monod, The growth of bacterial cultures, Annual Reviews on Microbiology 3 (1949) 371–394. [21] B.E. Rittmann, P.L. McCarty, Model of steady-state-biofilm kinetics, Biotechnology and Bioengineering 22 (1980) 2343–2357. [22] E.J. Bower, P.L. McCarty, Modeling of trace organics biotransformation in the subsurface, Ground Water 22 (1984) 433–440. [23] H.L. Reddy, R.M. Ford, Analysis of biodegradation and bacterial transport: Comparison of models with kinetic and equilibrium bacterial

adsorption, Journal on Contaminant Hydrology 22 (3–4) (1996) 271–287. [24] R.E. Ewing, T.F. Russell, Efficient time-stepping methods for miscible displacement problems in porous media, SIAM Journal on Numerical

Analysis 19 (1982) 1–66. [25] H.V. Kojouharov, B.M. Chen, Non-standard Eulerian–Lagrangian methods for advection–diffusion–reaction equations, in: Applications of

the Nonstandard Finite Difference Schemes, World Sci. Publishing, River Edge, NJ, 2000, pp. 55–108. [26] D.T. Dimitrov, H.V. Kojouharov, Stability-preserving finite-difference methods for general multi-dimensional autonomous dynamical

systems, International Journal of Numerical Analysis and Modelling 4 (2) (2007) 282–292. [27] L. Semprini, P. McCarty, Comparison between model simulations and field results for in-situ biorestoration of chlorinated aliphatics: Part 1.

biostimulation of methanotrophic bacteria, Ground Water 29 (3) (1991) 365–374. [28] L. Semprini, P. McCarty, Comparison between model simulations and field results for in-situ biorestoration of chlorinated aliphatics: Part 2,

Cometabolic Transformations, Ground Water 30 (1) (1992) 37–44. [29] H.-L. Chang, L. Alvarez-Cohen, Model for the cometabolic biodegradation of chlorinated organics, Environmental Science and Technology

29 (9) (1995) 2357–2367. [30] B. Travis, N. Rosenberg, Modeling in situ biodegradation of TCE at savannah River: Effects of product toxicity and microbial interactions on

TCE degradation, Environmental Science and Technology 31 (11) (1997) 3093–3102. [31] J. Komlos, A.B. Cunningham, A.K. Camper, R.R. Sharp, Biofilm barriers to contain and degrade dissolved trichloroethylene, Environmental

Progress 23 (1) (2004) 69–77. [32] J. Komlos, Effect of co-substrate concentration on dual-species Population Distribution, Permeability Reduction and Trichloroethylene

Biodegradation in Porous Media, Ph.D. Dissertation in Civil Engineering, Montana State University; September 2001.