Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances

ARTICLE IN PRESS

Water Research 38 (2004) 3349–3361

Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances Chrysi S. Laspidou*, Bruce E. Rittmann Department of Civil and Environmental Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3109, USA Received 22 May 2003; received in revised form 13 January 2004; accepted 29 April 2004

Abstract We present the unified multi-component cellular automaton (UMCCA) model, which predicts quantitatively the development of the biofilm’s composite density for three biofilm components: active bacteria, inert or dead biomass, and extracellular polymeric substances. The model also describes the concentrations of three soluble organic components (soluble substrate and two types of soluble microbial products) and oxygen. The UMCCA model is a hybrid discrete-differential mathematical model and introduces the novel feature of biofilm consolidation. Our hypothesis is that the fluid over the biofilm creates pressures and vibrations that cause the biofilm to consolidate, or pack itself to a higher density over time. Each biofilm compartment in the model output consolidates to a different degree that depends on the age of its biomass. The UMCCA model also adds a cellular automaton algorithm that identifies the path of least resistance and directly moves excess biomass along that path, thereby ensuring that the excess biomass is distributed efficiently. A companion paper illustrates the trends that the UMCCA model is able to represent and shows a comparison with experimental results. r 2004 Elsevier Ltd. All rights reserved. Keywords: Biofilm; Cellular automaton; Consolidation; Density; Detachment; EPS; Inert biomass; Modeling; SMP

1. Introduction Because experimental evidence shows that some biofilms have heterogeneous structure and composition (e.g., [3,13,17,28,47]), complex two- and three-dimensional models have been developed to predict different aspects of heterogeneity, e.g., the formation of microcolonies, the development of heterogeneous colonization patterns, and the formation of streamers [15,16,18,21,22,32,34–39,43]. While addressing many complicated features of modeling heterogeneous biofilms in two or three dimensions, these models present a simple view of the components that comprise the biofilm. With the exception of Eberl et al. [15,16], all *Corresponding author. Fax: +1-413-375-3178. E-mail address: [email protected] (C.S. Laspidou).

other advanced models assume a constant biomass density throughout the depth of the biofilm. Furthermore, most models include only one type of biomass and do not account for the formation of extracellular polymeric substances (EPS) and ‘‘inert’’ or dead biomass that microorganisms always produce [25]. Thus, although Kreft and Wimpenny [22] include EPS, none of these models distinguish among the three solid species: active biomass, inert biomass, and EPS. Although Eberl et al. [15,16] predict variable biomass density, only one biomass species is included. In this paper, we present a mechanistic, multicomponent model—the Unified Multiple-Component Cellular Automotan (UMCCA) model—that predicts quantitatively the biofilm’s heterogeneity for many components of a biofilm system: three solid species (active bacteria, inert or dead biomass produced by death and decay, and extracellular polysaccharides

0043-1354/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.watres.2004.04.037

ARTICLE IN PRESS 3350

C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

Nomenclature

qˆ S

b

qˆ UAP

First-order endogenous decay rate coefficient (T
1) B Ratio of secondary consolidation to total consolidation (—) BAP Dimensionless concentration of BAP[ — ] bap Concentration of BAP [MP/L3] bapmax Maximum BAP concentration [MP/L3] bdet Biofilm detachment coefficient [T
1] Bioage ‘‘Age’’ of each biofilm compartment (T) CompDen Composite density of biofilm [MX/L3] d Dimension of each square grid space element [L] D Dimensionless donor substrate [—] dmax Maximum donor substrate concentration [MD/L3] Ds Diffusion coefficient [L2/T] EPS Dimensionless concentration of EPS [—] eps Concentration of EPS [MP/L3] epsmax Maximum EPS packing density [MP/L3] fd Biodegradable fraction of active biomass [—] k1 UAP formation rate coefficient [MP/MS] KBAP Half-maximum-rate concentration for BAP utilization [MP/L3] KD Half-maximum-rate concentration for utilization of donor substrate [MD/L3] kEPS EPS formation coefficient [MP/MS] khyd First-order hydrolysis rate coefficient [T
1] KO Half-maximum-rate concentration for O2 consumption [MO/L3] KS Half-maximum-rate concentration for utilization of original substrate [MS/L3] KUAP Half-maximum-rate concentration for UAP utilization [MP/L3] O2 Oxygen concentration [MO2 /L3] O2,max Maximum oxygen concentration [MO2 /L3] Maximum specific BAP utilization rate [MP/ qˆ BAP MX–T] Maximum specific substrate utilization rate for qˆ D donor substrate [MD/MX–T]

(EPS)) and three soluble components (soluble substrate and two types of soluble microbial products (SMP) [25,31])). This model builds on our unified theory, which reconciles the apparently disparate findings about active and inert biomass, EPS, and SMP [25]. The model presented here is the biofilm adaptation of our multicomponent mathematical model that quantifies the unified theory [26]. Our biofilm model represents a growing biofilm using a cellular automaton (CA) approach in which the biofilm grows in a two-dimensional domain of compartments. One key feature is that

S s smax tc UAP uap uapmax Uc x X Xa xa xa;max Xres xres xres;max YP Ys z Z dtb dts Z l rBAP rD rS rUAP

Maximum specific substrate utilization rate for original substrate [MS/MX–T] Maximum specific UAP utilization rate [MP/ MX–T] Dimensionless concentration of original donor substrate [—] Concentration of original donor substrate [MS/L3] Maximum concentration of original donor substrate [MS/L3] Consolidation time [T] Dimensionless concentration of UAP [—] Concentration of UAP [MP/L3] Maximum concentration of UAP [MP/L3] Consolidation ratio [—] Biofilm width [L] Dimensionless biofilm width [—] Dimensionless density of active biomass [—] Active biomass density [MX/L3] Maximum active biomass packing density [MX/L3] Dimensionless density of true residual inert biomass [—] Density of true residual inert biomass [MX/L3] Maximum true residual inert biomass packing density [MX/L3] True yield for SMP (UAP and BAP) utilization [MX/MP] True yield for substrate utilization [MX/MS] Biofilm depth [L] Dimensionless biofilm depth [—] Time step used for biomass growth [T] Time step used for relaxation of S [T] Creep constant [T
1] 1st order decay coefficient for dissolved oxygen in the biofilm [L
1] Utilization rate of BAP [T
1] Utilization rate of any donor substrate [T
1] Utilization rate of original donor substrate [T
1] Utilization rate of UAP [T
1]

the UMCCA model produces a ‘‘composite density’’ that changes in time and space in the biofilm. We first review past experimental results on biofilm density. This review identifies important phenomena that must be included in the model. Then, we describe all the features of our two-dimensional, multi-component model. A companion manuscript [27] illustrates the trends that the UMCCA model represents concerning how the composite density and the distribution of the biofilm components vary with time and location in the biofilm, as well as the ability of our model to explain

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

experimental results obtained by Zhang and Bishop [53] and Bishop et al. [2].

2. Experimental trends in biofilm density Biofilm density does not remain constant throughout the biofilm, and several factors seem to influence it. Although differences in biofilm density have been explained by assuming the presence of different types of microorganisms [8], pure-culture studies demonstrate that physical forces created by hydrodynamic conditions have a direct influence on biofilm density. Vieira et al. [50] saw that the density of a pure-culture biofilm increased with increasing shear stress on the biofilm, which indicates a physical rather than a species effect. Kwok et al. [23] and Chang et al. [7] conducted experiments with moving-carrier biofilm processes, i.e., a Biofilm Airlift Suspension reactor and a liquid fluidized bed, respectively. Both research groups showed that biofilm density increased with increasing detachment forces related to particle–particle collisions. Kwok et al. [23] hypothesized that a very dense biofilm developed in their airlift reactor because high detachment forces removed low density, newly grown biomass from the outer layers of the carriers, leaving behind only the dense, strongly adhering biofilm. Based on their own observations, Melo and Pinheiro [30] identified inner ‘‘hard’’ and outer ‘‘loose’’ layers in biofilms grown in similar hydrodynamic conditions, with the outer layers being removed more readily by the fluid flow. Melo and Pinheiro [30] determined that the fraction of the loose layer in the biofilm decreased with an increase in the fluid velocity. Droste et al. [14] found a similar pattern for anaerobic biofilms. Observing a similar relationship between biofilm density and shear stress, Vieira et al. [50] expanded the explanation to include a biological response to a large detachment force: The bacteria reinforce the EPS matrix to protect themselves against the more aggressive forces from the surrounding fluid. In an extensive study on biofilm density, Zhang and Bishop [53] showed how biofilm density changes throughout the depth of the biofilm. The data were obtained using a microslicing technique able to measure density at different layers in the biofilm, rather than just an average density of the whole biofilm. They showed that the bottom layers of the biofilm had a higher density of total solids (gTS/cm3). The bottom layers were 5–10 times denser than the top layers. They also measured biofilm porosity, which decreased from 84– 93% in the top layers to 58–67% in the bottom layers. The density and porosity data confirm that the biomass was packed closer together near the substratum. Zhang and Bishop [53] also measured viable biomass. Near the substratum, the ratio of viable biomass to TS

3351

declined, which suggests that the biofilm was enriched in ‘‘dead’’ or ‘‘inert’’ biomass near the substratum. Viable biomass was approximately 30% of total biomass at the bottom of the biofilm, while it was nearly 100% at the outer surface. Furthermore, Ohashi et al. [33] found that the ‘‘live’’ cells, measured by a commercial kit (Live/ Dead Baclight Viability Kit L-7007), were not uniformly distributed throughout the biofilm depth. Again, the deep biofilm portion had a lower live-cell ratio compared to the outer surface. One-dimensional models of multi-species biofilms (e.g., [40,41,51]) also predict that inert biomass accumulates mainly near the attachment surface, while the active bacteria predominate near the outer surface. As stated by Rittmann and Manem [40], long-term operation of biofilm processes results in the accumulation of inert biomass from the endogenous decay of active biomass. Hence, mature biofilms can have a relatively small fraction of the volatile solids being active in substrate metabolism. Precipitated mineral solids are another form of inert solids. Cooke et al. [11,12] used one-dimensional modeling to describe how mineral deposits accumulate near the bottom of the biofilm in which precipitation occurs. According to all these findings, a multi-component biofilm model able to predict heterogeneity must include phenomena that underlie these observations: higher density for biofilms exposed to larger detachment forces, higher density near the attachment substratum, and more inert biomass near the attachment substratum. While the selective loss of ‘‘loose’’ biofilm may explain why low-density biofilm does not survive when detachment forces are strong, it does not explain how some of the biofilm becomes dense in the first place. One idea, expressed by Vieira et al. [50], is that the bacteria somehow ‘‘reinforce’’ the EPS, perhaps by regulating its production. Another hypothesis is that the continual exposure to physical forces from the flowing water causes the biofilm solids to consolidate, or pack more densely as water is squeezed out. We favor the consolidation hypothesis and include it in our model. Consolidation is understood and exploited in the fields of soil mechanics and dewatering of sludges and slurries [10,46,48,49]. Flowing water creates average forces that tend to squeeze water out of the biofilm and consolidate it. In addition, turbulent fluid flow over a biofilm creates fluctuating pressure that can accelerate biofilm consolidation. In simple words, turbulence from the moving fluid vibrates the biofilm, which causes irregularly shaped solids to realign themselves in such a manner that the porosity is minimized [19]. The average pressure load on a biofilm and the pressure fluctuations are increased when the biofilm surface has an irregular shape. Forcing the water to move over clusters increases the pressure loss due to friction, the turbulence, and the normal component of forces on the biofilm. The

ARTICLE IN PRESS 3352

C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

combination of an average normal force and pressure fluctuations results in consolidation. It is the bottom layers of the biofilm, the ones closest to the substratum, that undergo consolidation for the longest time; therefore, the bottom layers should have the highest densities and the lowest porosities. The model we present captures the phenomenon of biofilm consolidation in the context of a growing biofilm.

3. The UMCCA model The unified multi-component cellular automaton (UMCCA) model combines a discrete representation of the solid phase by CA with classical continuous methods for soluble components. The suitability of this combination for representing the structure of a heterogeneous biofilm has been demonstrated by recent works [18,32,34,35,39,52]. An advantage of this approach is that it can provide a fast and accurate model solution with readily available computing resources, such as a high-capacity personal computer. The UMCCA model begins with a differentialdiscrete CA approach similar to that used by Picioreanu et al. [34,35]. It computes properties that are associated with biofilm heterogeneity, such as biofilm density, porosity, and surface shape. For example, ‘‘mushroom’’ shaped clusters and finger-like projections can be generated, something not possible with one-dimensional models. One main goal of this model is to compute a ‘‘composite density,’’ or the density that includes active biomass, inert biomass, and EPS. This composite density corresponds to what is measured experimentally by total solids, volatile solids, or dry weight. To do this, the UMCCA model incorporates not only active biomass—as has been done so far by Picioreanu et al. [34,35], Noguera et al. [32], Pizarro et al. [39] and Eberl et al. [15,16]—but also EPS and residual inert biomass formed as a result of cell decay and lysis. The composite density can vary in time and in space, according to how its three components vary. Because each of the components is represented separately, the UMCCA model also computes the distribution of the different components in time and space. This makes it possible to compare the results of UMCCA model runs to experimental data on each component, such as those existing for viable biomass. The second goal is to incorporate the consolidation phenomenon into the model in a realistic, yet simple manner that is based on the well-established theory from the consolidation field. Consolidation allows the UMCCA model to describe the increases in biomass density that occur over time and deeper in the biofilm. A third goal of the UMCCA model is to find an improved system for distributing newly formed biomass.

Specifically, excess biomass in a compartment must be distributed in a way that is optimum in terms of the distance the biomass has to move until it finds an empty compartment in which it can be placed. Because biomass does not advect or diffuse, the CA algorithm has rules that allow the excess biomass to move out of a full compartment and into a nearby compartment that has room to accept it. Although previous CA algorithms have such rules, the UMCCA model incorporates a more efficient and realistic set of rules. 3.1. System definition The physical space is represented by a rectangular uniform grid with square compartments used to fill the 2D space. The model has Nx ¼ 150 grid points across the x Cartesian direction (parallel to the substratum) and Nz ¼ 70 square elements across the z direction (perpendicular to the substratum). The dimension of each square compartment is d ¼ 4 mm. 3.1.1. Variables and mass-balance equations The variables chosen to represent the status of each compartment are the dimensionless concentrations of the soluble organic species, densities of each of the solid species, and concentration of the electron acceptor, i.e., dissolved oxygen (O2). Soluble organic species are the bacteria’s limiting donor substrate (S) and two types of SMP, i.e., substrate utilization-associated products (UAP) and biomass-associated products (BAP). Solid species are the active biomass (Xa ), EPS, and the inert biomass (Xres ) produced by the decay/lysis of active biomass. For all species, dimensionless quantities are used to enhance the uniformity and stability of the algorithm, while time is not dimensionless, as was illustrated by Picioreanu et al. [34]. The mass balances used in the UMCCA model come directly from Laspidou and Rittmann [26], which quantified the unified theory presented in Laspidou and Rittmann [25]. The only difference is the dual Monod kinetics used for the donor substrate utilization rate, which includes donor substrate and oxygen (Eq. (1) in Table 1). According to the unified theory, substrate utilization generates four electron flows. The first flow is the synthesis of new active biomass, the second flow is the release of UAP, the third flow is the production of EPS, and the fourth flow results in consumption of the electron acceptor for energy generation. Active biomass decays in two ways. Part of the active biomass is oxidized by endogenous respiration to yield energy for maintenance. The rate is proportional to the biodegradable fraction of the active biomass (fd ), the endogenous decay coefficient (b; T
1), and the active biomass concentration (Xa ). Electron acceptor (O2) is consumed in proportion. Biomass decay also produces residual

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

3353

Rittmann [26] and the UMCCA model, utilization of UAP or BAP does not generate more UAP or EPS. Eqs. (2)–(7) in Table 1 are the mass balances for all components except oxygen. All variables are dimensionless and take real values between 0 and 1. All parameters and their units are defined in the Nomenclature list. Each mass balance has the non-steady-state accumulation term (T
1) on the left. Phenomena represented by

inert biomass (Xres ) in proportion to the rate of endogenous decay and the fraction of the active biomass that is not biodegradable (1
fd ). EPS dissolution/ hydrolysis produces BAP with a rate proportional to EPS and with khyd being the first-order rate coefficient (T
1). Hydrolysis of EPS is the only source of BAP. UAP and BAP also can be used as electron-donor substrates. According to the model of Laspidou and Table 1 Equations for all components in the UMCCA model Component

Equation

Utilization rate of any donor substrate rD ¼ Mass balance for original donor substrate

No.



 xa;max D O2 qˆ D Xa KD =dmax þ D dmax KO2 =O2;max þ O2

(1)


 qS Ds q2 S q2 S ¼ 2 þ
rS |{z} qt d qX 2 qZ 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} utilization

(2)


 qUAP Ds q2 UAP q2 UAP smax ¼ 2 þ þ k1 rS
rUAP |ffl{zffl} qt qX 2 qZ2 d uapmax |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} UAP utilization

(3)


 qBAP Ds q2 BAP q2 BAP ¼ 2 þ þ khyd EPS
2 2 |fflfflfflfflffl{zfflfflfflfflffl} qt qX qZ d

(4)

2-D diffusion

Mass balance for UAP

UAP formation

Mass balance for BAP

EPS formation

Mass balance for active biomass (Xa )

rBAP |ffl{zffl}

BAP utilization

qXa smax ¼ YS ð1
k1
kEPS Þ r þ qt xa;max S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

(5)

synthesis from S

uapmax bapmax YP r þ YP r
bX
bdet Xa |{z}a |fflffl{zfflffl} xa;max UAP xa;max BAP |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflffl} endogenous decay detachment synthesis from UAP

Mass balance for residual inert biomass (Xres )

qXres ¼ qt

synthesis from BAP

(6) bð1
fd ÞXa |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflffl ffl}

formation of inert biomass

Mass balance for EPS


bdet Xres |fflfflfflffl{zfflfflfflffl}

detachment

qEPS smax ¼ kEPS r
khyd EPS
bdet EPS |fflfflfflffl{zfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} qt epsmax S |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} EPS hydrolysis detachment

(7)

EPS formation

Definition of composite density

Definition of consolidation ratio

i;j xres;max CompDeni;j ¼ Xai;j xmax þEPSi;j epsmax þ Xres

Uc ¼ 1
0:9820 |fflfflffl{zfflfflffl} expð
0:0315 |fflfflffl{zfflfflffl} tc Þ B ½-

Profile for oxygen (O2)

O2 ¼ O2;max e
lz

(8)

(9)

Z ½h
1 

(10)

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

3354

the terms on the right are identified by words under each term. Eqs. (8) and (9) in Table 1 define the composite density (CompDenij) and the consolidation ratio (Uc ), both of which are described in detail below. In order to control the computing demands, we compute the oxygen profile based on first-order decay from a maximum value at the outer surface of the biofilm. The relationship is O2 ¼ O2;max e
lz

ð10Þ;

where l is the first-order decay coefficient for dissolved oxygen in the biofilm (L
1). For our standard condition, the maximum oxygen concentration O2,max is 9.2 mg/L (dissolved oxygen saturation at 25 C) and decreases with biofilm depth (z) with a decay coefficient (l) of 0.005 mm
1, which corresponds to the penetration depth of 800 mm experimentally observed by Bishop et al. [2] for a well-aerated heterotrophic biofilm. Imposing an oxygen profile made it possible for us to solve the UMCCA model for a domain grid that had enough compartments to represent realistic heterogeneity. Thus, we captured the key effects of oxygen limitation in the biofilm without exceeding the computing capacity of our computer. 3.2. Using the UMCCA model The solution strategy, described next and illustrated in the flowchart in Figs. 1 and 2, follows Picioreanu et al. [34,35], but it additionally calculates components unique to the UMCCA model: UAP, BAP, true residual inert biomass, and EPS according to our unified theory [25]. The UMCCA model combines the three solid components (active biomass, inert residual biomass, and EPS) to calculate the biofilm’s composite density, which corresponds to what is measured experimentally by total solids, volatile solids, or dry weight. All components and the composite density vary with time and location. The physical space, boundary conditions, and initial conditions are summarized in Fig. 1. 3.2.1. Solution algorithm For every CA cell—referred to as a ‘‘compartment’’ in this work—we keep track of a series of quantities: Sij — substrate concentration, UAPi;j —UAP concentration, BAPi;j —BAP concentration, Oi;j 2 —dissolved oxygen concentration, Xai;j —active biomass density, EPSi;j — i;j EPS concentration, Xres —residual inert biomass density,
3 i;j CompDen [ML ]—composite biofilm density, and Bioagei;j [T]—‘‘age’’ of each biofilm compartment, i.e., elapsed time since the formation of the biomass in the specific compartment. Superscripts i and j indicate the position of the compartment, with i indicating the row up from the substratum boundary and j indicating the column from the left boundary.

The composite density is calculated by multiplying the fraction of each of the solid species with its maximum density. i;j CompDeni;j ¼ Xai;j xmax þEPSi;j epsmax þ Xres xres;max : ð8Þ

Eq. (8) underscores that each compartment contains all three biomass components, each at a unique fraction of its maximum density. Each compartment is homogeneous with regard to the biomass components. Thus, the compartment is the smallest unit of biomass with a homogeneous composition. The way the age of the biomass in each compartment is computed is explained next. Solution begins with initial conditions for all components in all compartments. Any set of initial conditions are possible. However, we set the following standard initial conditions. For all compartments, Si;j ¼ 1:0; indicating that the maximum substrate is initially available throughout the domain. All other quantities, except for active biomass and oxygen, are set to zero. The substratum is then randomly seeded with active biomass. Half of the compartments adjacent to the substratum (i.e., 35 out of 70) are randomly selected and inoculated with an amount of active biomass that is randomly selected between 0.5 and 1.0 for each compartment. This inoculation simulates an initial random and patchy attachment of biomass across the substratum. The zero-flux boundary condition is assumed at the substratum and at the side walls of the modeled space, i.e., qS=qX ¼ 0 and qS=qZ ¼ 0 at these boundaries. This means that the substratum and side walls are completely impermeable to the substrate. The substrate mass balance (Eq. (4)) is solved next for all locations to define the substrate field at time t þ dtb ; while the biomass field is kept constant. This is done using a strategy similar to that described in [34], i.e., keeping the biomass density unchanged, while iterating until the substrate concentration field converges to a steady-state condition, and using a time step dts to solve the substrate field of about 1 s and a time step for changes in biomass, dtb of 1000 s. The substrate field solution by a numerical method also follows the scheme presented by Picioreanu et al. [34] with a finite-difference discretization scheme that is then solved with the alternating-direction implicit method. The evolution of biomass concentration in every compartment, using the new substrate field, is computed next. The time step in this case is dtb ; as opposed to the much smaller time step used in the previous procedure to obtain the steady-state substrate field (dts ). Once the biomass concentration field is calculated, all other quantities can be calculated. At this point, i.e., once the new biomass concentration field is obtained, the matrix Bioagei,j is updated. Initially, this matrix is set to 0 s. After one iteration,

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

3355

Initialize 1. Randomly choose half of 150 compartments adjacent to the substratum and “seed” them with Cin. Cit,=j0( random ) = Cin (random between 0.5 and 1.0) 2. Bioagei , j ( random ) = 0 3.

S it,=j 0 = 1.0

4.

UAPit, =j 0 = 0

5.

BAPit, =j 0 = 0

6.

EPS it,=j 0 = 0

7.

O2 ti ,=j0 = 1.0

Calculate Si,j by solving t=t+δts

T=0

∂S ∂t

Convergence?

T=T+∆ts

NO YES,

quasi-steady-state

Calculate UAPi,j by solving Eqn (3).

• • • • • •

Calculate Xai,j by solving Eqn (5). Calculate Xresi,j by solving Eqn (6). Calculate EPSi,j by solving Eqn (7). Calculate BAPi,j by solving Eqn (4). Calculate CompDeni.j by solving Eqn. (8). Calculate Uci,j by solving Eqn. (9).

T=T+δtb

∀i, j NO

Xai,j + Xresi,j + EPSi,j >Uc?

Terminate for all i,j

YES

Cellular Automaton Algorithm (see Fig. 2)

Calculate O2i,j by solving Eqn (10). Fig. 1. Flowchart of the solution strategy for the UMCCA model. The flowchart of the CA Algorithm is shown in Fig. 2.

and since the elapsed time is dtb ¼ 1000 s, all compartments that have some biomass in them increase their bioage by 1000 s. The empty compartments are still at a bioage of 0 s. When the second iteration is completed, all occupied compartments increase their bioage to 2000 s,

and then 3000 s with the third iteration, and so on. When an unoccupied compartment becomes occupied at some point, its bioage changes from 0 to 1000 s and increases accordingly thereafter. Thus, at any time, the biofilm is made of compartments that have different

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

3356

X aXS = 0.5 . X ai , j EPS XS = 0.5 EPS i , j i =1

Find all neighbors of i,j that are I compartments away

Is any of the neighbors unoccupied?

NO I=I+1

Find all neighbors of i,j that are I compartments away

Xam,n = X aXS YES, only compartment m,n

EPSm,n = EPSaXS

YES, 2 or more compartments Randomly select ONE of the unoccupied compartments, call it m,n

Store path from i,j to the current neighboring compartment

NO

Is any of the neighbors unoccupied?

YES, 2 or more compartments

Randomly select ONE of the unoccupied compartments, call it m,n

YES, only compartment m,n

Backtrace path (shortest path) from m,n to i,j. Copy XS XS X iXS , j and EPS to nearest neighbor, its X i , j and XS EPSXS to next-nearest neighbor, its X iXS , j and EPS

Terminate CA module

to next-next nearest neighbor, etc., until compartment m,n is reached. Or, make X iXS , j and EPSXS push biomass in all occupied compartments along the path of least resistance (shortest path) to make room for itself. Fig. 2. Flowchart of the CA algorithm of the UMCCA model.

bioages: most often, the top layers have ‘‘younger’’ biofilm, while the bottom layers have ‘‘older’’ biofilm. The specific compartments that were randomly chosen to be seeded with biomass first (at time=0) have the

‘‘oldest’’ biofilm, since their bioage is the same as the model’s elapsed time. To solve the biomass mass balance (Eq. (5)) and, following that, all the rest of the mass balances, we

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

employ the strategy described in [34], i.e., we use the discretized form of the mass-balance equation dA ¼ B ð11Þ dtb in which A is a solid-phase component, B is the time derivative of A; and dtb is the time step. The solution marches forward in time with a time step dtb Atþdtb ¼ At þ Bt dtb :

ð12Þ

The superscript t indicates the value of the variable A at time t; and t þ dtb indicates the value of that variable at the next time step, i.e., t þ dtb : The active biomass concentration is calculated first, and then all other mass balances for solids—Eqs. (6) and (7)—are calculated in a similar manner, using the substrate field and the computed active biomass field. 3.2.2. The cellular automaton algorithm, which includes consolidation, moving excess biomass, and detachment Active biomass and EPS are distributed according to the CA algorithm. Such biomass distribution simulates the movement of the newly formed ‘‘daughter’’ cells that result from the division of ‘‘mother’’ cells. We expect such cells to be ‘‘wrapped’’ in EPS, often referred to as ‘‘glycocalyx;’’ therefore, we expect EPS to move and get distributed together with the new cells. However, we do not expect the same behavior for the residual inert biomass, which accumulates at the bottom of the biofilm; therefore, we do not distribute the residual inert biomass, but only the active biomass and EPS. Once the decision is made that active biomass and EPS have to be distributed, these quantities are divided in two parts: 50% of each quantity stays at the same site, and 50% is the excess mass that is placed in another compartment. Biomass distribution is performed for a compartment when its composite density exceeds the maximum composite density for that compartment. The maximum composite density—the sum of the individual fractions of the three biomass species: active, residual and EPS— depends on the age of the biofilm in the compartment. As biofilm in a compartment is exposed to pressure and pressure fluctuations for a longer time, its maximum density for each component increases. Thus, increasing the maximum density over time is the simple way by which we incorporate biofilm consolidation in the UMCCA model. 3.2.2.1. Implementing consolidation. Mathematically, consolidation is implemented with the consolidation ratio Uc, which defines the degree to which the maximum composite density of a compartment has been consolidated towards an absolute maximum for a biofilm that has reached its limit of consolidation. Hence, the maximum value of the consolidation factor is

3357

1.0, signifying that the biofilm has reached an absolute maximum packing density that corresponds to being as completely dewatered as is possible. Consolidation in the UMCCA model is represented by Eqs. (9) or (13), which are adapted from the wellestablished consolidation theory [1,4–6,9,10,19,29,44– 46,48,49]. Laspidou [24] provides a review of consolidation theory and how it applies to the biofilm setting. Key to the application here is secondary consolidation, which is an irreversible, non-linear ‘‘densification’’ or ‘‘creep.’’ Secondary consolidation is represented by the creep rate constant Z (T
1) and the fraction of total consolidation that is secondary, B. Eqs. (9) and (13) show the exponential and logarithmic forms with our standard values of B and Z Uc ¼ 1
0:9820 |fflfflffl{zfflfflffl} expð
0:0315 |fflfflffl{zfflfflffl} tc Þ;

ð9Þ

ln ð1
Uc Þ ¼ ln 0:9820 |fflfflffl{zfflfflffl}
0:0315 |fflfflffl{zfflfflffl} tc :

ð13Þ

B ½-

Z ½h
1 

B ½-

Z ½h
1 

When using the principles and mathematics of consolidation in the context of this work, we are adapting concepts from the field of mechanical dewatering of sludges, slurries, and soils. Therefore, the pressures applied on the biofilm by the moving fluid usually are much smaller than those applied in a mechanical system. The secondary-consolidation effects on the biofilm happen because the small pressures are coupled with vibration due to pressure fluctuation. Thus, the time scale for biofilm consolidation is longer than typical for mechanical consolidation of sludges and slurries. This longer time scale is represented with a relatively small value of Z (0.0315/h), which indicates a high-viscosity dashpot (which is used to represent irreversible creep), and thus a slow consolidation. However, almost all of the consolidation process is secondary consolidation, which is represented by the large value of B (0.982). B ¼ 0:982 means that only about 2% of the total consolidation occurs as primary consolidation, which has the linear ‘‘bounce back’’ behavior. The remaining 98% occurs as secondary consolidation, which is non-linear creep. The small Z value, which indicates slow consolidation, is consistent with the large B value, which indicates an inelastic sludge. In the UMCCA model, the consolidation time (tc ) in Eqs. (9) and (13) is the bioage of each compartment. Therefore, the biofilm is composed of compartments that correspond to different degrees of consolidation, i.e., have experienced consolidation pressure for different times and, therefore, have different Uc values. The bottom compartments, because they normally contain the most mature biofilm with the longest bioages, are consolidated more than the top compartments, which have young biofilms that have had minimum

ARTICLE IN PRESS 3358

C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

consolidation times. Thus, the UMCCA model naturally leads to a situation in which the outer layers of biofilm are relatively loose and fluffy, while the bottom layers are dense and tightly packed. The way we implement consolidation in UMCCA is by using Uc to define when a compartment has excess biomass. Specifically, excess biomass is redistributed when i;j Xai;j þ EPSi;j þ Xres XUc :

ð14Þ

3.2.2.2. Distributing excess biomass. Once a compartment is defined as ‘‘overflowing’’ by Eq. (15), then the 50% of the active biomass and EPS is moved to one of its eight neighboring compartments that has room to receive it, i.e., is totally unoccupied by any solid species. If more than one compartment is unoccupied, then a random choice is made. If all eight neighboring cells are occupied with biomass, then the algorithm identifies the shortest path from the overfilled compartment to the nearest empty compartment. Such a compartment could exist either inside the biofilm, or it may be at the biomass interface with the overlying fluid. This approach represents the ‘‘path of least resistance.’’ The new cell pushes a line of other cells in the direction of the shortest distance and displaces a row or a column or a combination of rows and columns to make space for itself. If two or more paths are of equal distance, a random choice is made. A similar idea has been implemented by Hermanowicz [18]. Laspidou [24] provides a more detailed description of how the UMCCA model distributes biomass, including graphical representation of the grid and examples of distributing excess biomass. The term ‘‘distance’’ in this context does not correspond exactly to physical distance, but to the number of neighboring cells that lie between the overflowing compartment and the candidate empty compartment. For example, the path that is two compartments away along a diagonal is considered equivalent to the path that is two compartments away along a horizontal or vertical line. Naturally, the actual physical distance is shorter in the latter case, but the CA algorithm does not make such a differentiation. Although there is little insight regarding the path of bacteria in a growing biofilm, it is a fact that cells growing near the bottom of a biofilm must displace other components as the biofilm thickness increases. Thus, it seems reasonable to assume that cells will displace the minimum number of cells, therefore exerting the smallest energy possible. This is the concept underlying our ‘‘path of least resistance’’ approach. The bioage of each compartment is not affected by the movement of the biomass. Any newly occupied compartment has its bioage reset to zero and starts increasing from there. Any compartment that remains

occupied increases its bioage with every time step. In other words, bioage is time- and place- (or compartment-) specific. If a compartment has been occupied by biomass for 1 h, then it has accepted the consolidation pressure and vibration for 1 h, which is its bioage, independent of the fact that biomass might be pushed in and out of that compartment. This way, bottom compartments that have the oldest biofilm have the highest bioages and vice versa. 3.2.2.3. Detachment. In the UMCCA model, biomass decays to form residual inert biomass. In addition, detachment with first-order coefficient bdet (T
1) describes the constant erosion caused by the fluid moving by the biofilm [42]. In other models [18,39], decay is handled similar to detachment, i.e., by making microbial cells disappear from the compartments with some preset probability. In the UMCCA model, the outermost layer of the biofilm is identified, and first-order detachment is applied only on those compartments; thus, the UMCCA model makes a clear distinction between detachment and decay. In addition to having detachment by erosion only at the outer layer, the UMCCA model also detaches all of the solid components, not just active biomass. 3.2.3. Computer environment and parameter values The UMCCA is programmed in C++ and is executed on a personal computer. The computer we used was a standard, commercially available Pentium with 4 CPU, 2.4 GHz, with 512 MB of RAM. The standard-case values of the variables in the UMCCA model are listed in Table 2. The unit system for the parameters uses g or mg chemical oxygen demand (COD) for all organic species, liters for volume, meters for distance, and days for time. We chose the parameter values based on the following criteria: *

*

* *

We took values for Ds ; xa;max ; epsmax, Ys ; Yp ; qˆ s ; fd ; k1 ; kEPS ; khyd ; KUAP ; qˆ UAP ; KBAP ; and qˆ BAP directly from the literature. Ds and xa;max were taken from [34]. epsmax was taken from Kreft and Wimpenny [22], who quote a maximum EPS density as 75 g C/L. Assuming that a mole of EPS is similar to that of a cell (C5H7O2N, as in [42]), we can obtain epsmax, using the ratios of 60 g C per mole of EPS and 1.42 g CODx per formula weight. The other variables are identical to the ones used in [26], in which we evaluated our unified theory [25] by comparing it with experimental data found in [20]. The values for Ks ; b; and bdet are typical for aerobic heterotrophs in environmental biotechnology [42]. The smax value is typical for wastewater influent [42]. Values for uapmax and bapmax are needed to have dimensionless concentrations of UAP and BAP. We chose uapmax to be the same as smax ; since it cannot

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

3359

Table 2 Standard parameter values for the UMCCA model Ds smax uapmax bapmax xa;max epsmax xres;max O2,max l Ys Yp b Z

*

1.38 10
4 m2/day 500 mgs/L 500 mg CODp/L 50 mg CODp/L 70 g CODx/L 200 g CODx/L 220 g CODx/L 9.2 mgO2 /L 0.005/mm 0.34 mgx/mgs 0.45 mgx/mgp 0.4/day 0.0315/h

exceed that maximum. We chose a smaller value for bapmax, because BAP is usually smaller than UAP. The values of B; Z; and xres,max are unique to the UMCCA model, because they are new variables associated with biofilm consolidation. The B and Z values represent a high degree of secondary consolidation and relatively slow consolidation kinetics, as described earlier. The values are derived in [24]. The value for xres;max approaches the maximum packing density of biomass solids, because it includes the remains of the active cells, after their lysis when they have released their internal water. In other words, unlike the active biomass that is mostly water, the residual dead biomass includes mostly the partially dehydrated remnants of the cells and all other mineral deposits that can achieve a maximum packing density; therefore, xres;max is higher than the other solid species.

4. Conclusions We present a mechanistic, multi-component model that represents quantitatively the biofilm’s heterogeneity for all components of a biofilm system: three solid species (active bacteria, inert or dead biomass produced, and EPS), and three soluble components (soluble substrate, two types of SMP). As a result, it computes a ‘‘composite density’’ that changes in time and space in the biofilm. The UMCCA model is a hybrid discretedifferential mathematical model that represents soluble components in a continuous field and employs discrete mapping for solid components. In the UMCCA model, the fluid flow over the biofilm creates horizontal and vertical pressures, as well as vibration due to the fluctuating pressure that causes the biofilm to consolidate, or pack itself to a higher density. Consolidation in the UMCCA model is implemented by having the consolidation ratio (Uc ) increase with the age

qˆ s Ks fd bdet k1 kEPS khyd qˆ UAP KUAP qˆ BAP KBAP B

28.5 mg CODs/mg CODx-day 20 mg CODs/L 0.8 0.15/day 0.05 mg CODp/mg CODs 0.18 mg CODp/mg CODs 0.17/day 1.8 mg CODp/mg CODx-day 100 mg CODs/L 0.1 mg CODs/mg CODx-day 85 mg CODs/L 0.9820

of the biomass and, therefore, the time it experiences consolidation. Consolidation is the key feature that allows the UMCCA model to represent the density differences throughout the biofilm depth that have been observed experimentally. The UMCCA model employs a CA algorithm that distributes newly formed biomass in an improved way along the path of least resistance. It also models decay to inert biomass and detachment of all biomass types by erosion from the outer surface as distinct processes. The output of the UMCCA model is a biofilm composed of compartments that have consolidated to different degrees and have different distributions of active biomass, inert biomass, and EPS. Thus, each compartment has its own unique composite density. A companion paper [27] illustrates the trends that the UMCCA model is able to predict. Discrimination among different types of biomass, as well as the prediction of how the composite density varies with location in the biofilm, makes it possible to compare model predictions to experimental observations of these key structural properties. In the comparison paper [27], we make such a comparison with experimental results that quantify the composite density and differentiate among the biomass components.

References [1] Bement RAP, Selby AR. Compaction of granular soils by uniform vibration equivalent to vibrodriving of piles. Geotech Geol Eng 1997;15:121–43. [2] Bishop PL, Zhang TC, Fu Y-C. Effects of biofilm structure, microbial distributions and mass transport on biodegradation processes. Water Sci Technol 1995;31(1):143–52. [3] Caldwell DE, Korber JR, Lawrence DR. Analysis of biofilm formation using 2D vs. 3D digital imaging. J Appl Bacteriol 1993;74:S52–66.

ARTICLE IN PRESS 3360

C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361

[4] Chang IL, Chu CP, Lee DJ. Filtration followed by expression characteristics of polymer flocculated clay sludge. J Colloid Interface Sci 1997;185(2):335–42. [5] Chang IL, Chu CP, Lee DJ, Huang C. Filtration followed by expression characteristics of alum coagulated clay sludge. Environ Sci Technol 1997;31(5):1313–9. [6] Chang IL, Lee DJ. Expression dewatering of mixed sludges. Spectrosc Lett 1996;29(8):1659–72. [7] Chang HT, Rittmann BE, Amar D, Heim R, Ehlinger O, Lesty Y. Biofilm detachment mechanisms in a liquidfluidized bed. Biotechnol Bioeng 1991;38:499–506. [8] Christensen BE, Characklis WG. Physical and chemical properties of biofilms. In: Characklis WG, Marshall KG, editors. Biofilms. New York: Wiley; 1990. p. 93–130. [9] Chu CP, Lee DJ. Three stages of consolidation dewatering of sludges. J Environ Eng 1999;125(10):959–65. [10] Cleveland TG, Tiller FM, Lee JB. Theory of filtration of highly compactible biosolids. Water Sci Technol 1996; 299–306. [11] Cooke AJ, Rowe RK, Rittmann BE, Fleming IR. Modeling biochemically driven mineral precipitation in anaerobic biofilms. Water Sci Technol 1999;39(7):57–64. [12] Cooke AJ, Rowe RK, Rittmann BE, vanGulck J, Millward S. Biofilm growth and mineral precipitation in synthetic leachate columns. J Geotech Geoenviron Eng 2001;127:849–56. [13] Costerton JW, Lewandowski Z, De Beer D, Caldwell D, Korber D, James G. Minireview: biofilms, the customized microniche. J Bacteriol 1994;176(8):2137–42. [14] Droste RL, Andras E, Kennedy KJ. Initial biofilm formation of acetoclastic methanogenic bacteria. Biofouling 1990;2:191–210. [15] Eberl HJ, Parker DF, van Loosdrecht MCM. A new deterministic spatio-temporal continuum model for biofilm development. J Theor Med 2001;3:161–75. [16] Eberl HJ, Picioreanu C, Heijnen JJ, van Loosdrecht MCM. A three-dimensional numerical study on the correlation of spatial structure, hydrodynamic conditions, and mass transfer and conversion in biofilms. Chem Eng Sci 2000;55:6209–22. [17] Gjaltema A, Arts PAM, van Loosdrecht MCM, Kuenen JG, Heijnen JJ. Heterogeneity of biofilms in rotating annular reactor: occurrence, structure and consequences. Biotechnol Bioeng 1994;44:194–204. [18] Hermanowicz SW. Two-Dimensional Simulations of Biofilm Development: effects of external environmental conditions. Water Sci Technol 1999;39(7):107–14. [19] Holtz RD, Kovacs WD. An introduction to geotechnical engineering. New York: Prentice-Hall; 1981. [20] Hsieh KM, Murgel GA, Lion LW, Schuller ML. Interactions of microbial biofilms with toxic trace metals I. Observation and modeling of cell growth attachment, and production of extracellular polymer. Biotechnol Bioeng 1994;44:219–31. [21] Kreft J-U, Booth G, Wimpenny JWT. BacSim, a simulator for individual-based modeling of bacterial colony growth. Microbiology 1998;144:3275–87. [22] Kreft J-U, Wimpenny JWT. Effect of EPS on biofilm structure and function as revealed by an individual-based model of biofilm growth. Water Sci Technol 2001;43(6):135–41.

[23] Kwok WK, Picioreanu C, Ong SL, van Loosdrecht MCM, Ng WJ, Heijnen JJ. Influence of biomass production and detachment forces on biofilm structures in a biofilm airlift suspension reactor. Biotechnol Bioeng 1998;58(4): 400–7. [24] Laspidou CS. Modeling heterogeneous biofilms including active biomass, inert biomass and extracellular polymeric substances. PhD dissertation, Department of Civil and Environmental Engineering, Northwestern University, Evanston, IL, USA, 2003. [25] Laspidou CS, Rittmann BE. A unified theory for extracellular polymeric substances, soluble microbial products, and active and inert biomass. Water Res 2002;36:2711–20. [26] Laspidou CS, Rittmann BE. Non-steady state modeling of extracellular polymeric substances, soluble microbial products, and active and inert biomass. Water Res 2002;36:1983–92. [27] Laspidou CS, Rittmann BE. Evaluating trends in biofilm density using the UMCCA model. Water Res 2004, in press, doi:10.1016/j.watres.2004.04.051. [28] Lawrence JR, Korber DR, Hoyle BD, Costerton JW, Caldwell DE. Optical sectioning of microbial biofilms. J Bacteriol 1991;173:6558–67. [29] Lee DJ, Wang CH. Theories of cake filtration and consolidation and implications to sludge dewatering. Water Res 2000;34(1):1–20. [30] Melo LF, Pinheiro JD. Fouling tests: equipment and methods. In: Suitor JW, Pritchard AM, editors. Fouling in heat exchange equipment. New York: American Society for Mechanical Engineering; 1984. p. 43–9. [31] Namkung E, Rittmann BE. Soluble microbial products (SMP) formation kinetics by biofilms. Water Res 1986;20(6):795–806. [32] Noguera DR, Pizarro G, Stahl DA, Rittmann BE. Simulation of multispecies biofilm development in three dimensions. Water Sci Technol 1999;39(7):123–30. [33] Ohashi A, Koyama T, Syutsubo K, Harada H. A novel method for evaluation of biofilm tensile strength resisting to erosion. Water Sci Technol 1999;39(7):261–8. [34] Picioreanu C, van Loosdrecht MCM, Heijnen JJ. A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads. Biotechnol Bioeng 1998;57(6):718–31. [35] Picioreanu C, van Loosdrecht MCM, Heijnen JJ. Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach. Biotechnol Bioeng 1998;58:101–16. [36] Picioreanu C, van Loosdrecht MCM, Heijnen JJ. Discretedifferential modeling of biofilm structure. Water Sci Technol 1999;39(7):115–22. [37] Picioreanu C, van Loosdrecht MCM, Heijnen JJ. A theoretical study on the effect of surface roughness on mass transport and transformation in biofilms. Biotechnol Bioeng 2000;68(4):355–69. [38] Picioreanu C, van Loosdrecht MCM, Heijnen JJ. Twodimensional model of biofilm detachment caused by internal stress from liquid flow. Biotechnol Bioeng 2001;72:205–18. [39] Pizarro G, Griffeath D, Noguera DR. Quantitative cellular automaton model for biofilms. J Environ Eng 2001;127(9):782–9.

ARTICLE IN PRESS C.S. Laspidou, B.E. Rittmann / Water Research 38 (2004) 3349–3361 [40] Rittmann BE, Manem JA. Development and experimental evaluations of a steady-state, multi-species biofilm model. Biotechnol Bioeng 1992;39:914–22. [41] Rittmann BE, Stilwell D, Ohashi A. The transient-state, multiple-species biofilm model for aerobic biofiltration processes. Water Res 2002;36:2342–56. [42] Rittmann BE, McCarty PL. Environmental biotechnology: principles and applications. New York, NY: Mc-Graw Hill; 2001. [43] Rittmann BE, Pettis M, Reeves HW, Stahl DA. How biofilm clusters affect substrate flux and ecological selection. Water Sci Technol 1999;39(7):99–105. [44] Shirato M, Murase T, Kato H, Shibata M. Analysis for constant pressure expression mechanism. J Ferment Eng 1965;255–65 (in Japanese). [45] Shirato M, Murase T, Tokunaga A, Yamada O. Calculations of consolidation period in expression operations. J Chem Eng Japan 1974;7(3):229–31. [46] Sorensen PB, Hansen JAA. Extreme compressibility in biological sludge dewatering. Water Sci Technol 1993; 133–43.

3361

[47] Stoodley P, Lewandowski Z, Boyle JD, Lappin-Scott HM. Structural deformation of bacterial biofilms caused by short-term fluctuations in fluid shear: an in situ investigation of biofilm rheology. Biotechnol Bioeng 1999;65: 83–92. [48] Tiller FM, Kwon JH. The role of porosity in filtration: XII. Filtration with sedimentation. AIChE J 1998; 1153–64. [49] Tiller FM, Yeh CS. The role of porosity in filtration: XI. Filtration followed by expression. AIChE J 1987;1241–56. [50] Vieira MJ, Melo LP, Pinheiro MM. Biofilm formations: hydrodynamic effects on internal diffusion and structure. Biofouling 1993;7:67–80. [51] Wanner O, Gujer W. A multispecies biofilm model. Biotechnol Bioeng 1986;28:314–28. [52] Wimpenny JWT, Colasanti R. A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models. FEMS Microbiol Ecol 1997;22: 1–16. [53] Zhang TC, Bishop PL. Density, porosity and pore structure of biofilms. Water Res 1994;28:2267–77.