Cellular Autometa Matlab algorithm
aseelismail
bioremediationAndCA02.pdf
ELSEVIER Theoretical Computer Science 217 (1999) 131-1.56
Theoretical Computer Science
Applying cellular automata to complex environmental problems: The simulation of the bioremediation
of contaminated soils
S. Di Gregorio a, R. SerrabT*, M. Villani b
a Dipartimento di Matematica, Universitri della Calabria, Arcavacata, I-87036 Rende (CS), Italy b Centro Ricerche Ambientali Montecatini, via C. Menotti 48, I-48023 Marina di Ravenna, Italy
Abstract
Cellular automata can be applied to modelling the dynamics of spatially extended physical
systems, and represent an alternative to the classical PDE approach. In this paper a macroscopic cellular automata model for simulating the bioremediation of contaminated soils is introduced.
The choice of macroscopic automata is motivated by the aim to simulate large-scale systems.
It is suggested that in some cases, where the basic laws of continuum mechanics cannot be
directly applied without adding phenomenological assumptions, and where the equation system
is not amenable to analytical solution, direct discrete modelling may represent a convenient
alternative to the use of continuum models, followed by numerical discretization. This hypothesis
is empirically tested in the bioremediation case.
The model describes the bioremediation of contaminated soils, which relies upon the use of
indigeneous microorganisms to degrade the contaminant: bioremediation models pose particular
challenges as several physical, chemical and biological phenomena interact in a disordered and
partially unknown matrix (the soil). The model is hierarchical, and is composed by a fluid
dynamical layer, a solute description layer and a biological layer. The model has been tested
in a pilot plant, in the case of contamination by phenol. The values of the phenomenological
parameters have been determined by the use of genetic algorithms. The model has proven capable
to carefully describe experimental results in a wide range of experimental conditions. It has also
been run on a MIMD parallel architecture, achieving a high speed-up. It therefore represents
an example of application of cellular automata to a real-world problem which has a very high
social and economic importance, and where progresses in modelling may greatly improve the
effectiveness of the decontamination interventions. @ 1999-Elsevier Science B.V. All rights
reserved
Keywords: Cellular automata; Bioremediation; Genetic algorithms; Multiphase flow; Porous medium; Bacterial dynamics
* Corresponding author. E-mail: [email protected].
0304-3975/99/$-see front matter @ 1999-Elsevier Science B.V. All rights reserved PII: SO304-3975(98)00154-6
132 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
1. Introduction
Cellular automata, introduced in the 1950’s by von Neumann to study self- reproducing systems [lo], have been later used for modelling complex dynamics and physical systems [37,39,40]. In particular, macroscopic phenomena like, e.g. diffusion, which are usually described within a continuum approach, have been actively investi- gated. Most attempts to apply CA to modelling these physical phenomena have dealt with a “microscopic” approach, where the state variables can take only a limited num- ber of different discrete values [22,37]; for example, the state of a cell represents the presence or absence of a particle with a given momentum at a given point. A similar approach has been applied also to modelling percolation phenomena [3].
In the microscopic approach, the system is considered as composed of a number of discrete particles, which move in the lattice space and interact according to approp~ate laws: the laws which rule the system at a macroscopic level are then obtained by averaging over a large number of particles.
This description can be viewed as “more microscopic” than the continuum equations, like the diffusion or Navier-Stokes equations, because the basic tokens of the latter (concen~ations, pressures, etc.) are themselves the result of an ~derlying discrete dynamics, which is directly described in the former approach. This view is close to that of molecular dynamics, but it tries to avoid the computational burden using very simple basic discrete particles and interaction laws, so that the state variables associated to each cell can take only a limited number of discrete values.
We adopt here a different attitude, considering the cells as portions of space which include several pores: it is therefore me~ing~l to assign to each cell a complex state space, which is the Cartesian product of different subspaces, which may refer, e.g. to the water content, the concentration of a given chemical, the density of a given kind of bacteria, etc. While some of these variables are usually considered as continuous, we will admit their use, assuring the formal compliance with the CA definition by limiting the range of allowed values to a discrete set (as it is always the case when digital simulations are performed). However, in order to avoid unnecessary linguistic burdens, we will refer to these variables in the usual (i.e. continuous) way, always assuming that discretization is performed.
This kind of model is reminiscent of the so-called lattice Boltzmann models [35]; however, our model is macroscopic, and makes use of variables such as water satura- tion, concentration of a given chemical, etc., while the applications of lattice Boltzmann models to flow in porous media so far developed take a more microscopic approach and aim at describing phenomena which take place at the pore level [l]. Macroscopic CA models can be regarded as intermediate between classical PDE models and the so-called compa~ment models of flow in soils [9], which provide a simpler description of the processes which take place in soil, by dividing the space into cells of finite size, which are assumed to be in local equilibrium; this latter condition is not required here.
The use of these “macroscopic” CA models seems promising to simulate large sys- tems because they allow one to choose a cell size which is appropriate for the scale
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 133
of the simulation. Moreover, they avoid the need for a limiting operation, followed
by the introduction of a discrete approximation. Encouraging examples come from the
application of this kind of models to landslide and lava flow simulations [6,21].
Ultimately, the validity of this approach has to be judged by its usefulness, which
implies a reasonable ease to develop new models and a good agreement with experi-
mental data. We tried to verify the effectiveness of this approach in the development
of a model of the bioremediation of contaminated soils, which is both a challenging
scientific problem and a very important environmental technology.
Section 2 contains a short description of bioremediation, while in Section 3 the
basic features of the model are outlined. The physical reasons for most model choices
are only sketched, and the reader is referred to [ 15, 17, 191 for a more comprehensive
discussion of these aspects.
The model has a three-layer structure, where the first layer describes the fluid dy-
namical phenomena, the second layer the physico-chemical dynamics of the solutes
and the third layer describes biomass growth and those phenomena where microorgan-
isms play a role (e.g. contaminant transformation). The three layers are described in
Sections 4-6.
The model has a number of parameters which cannot be (or have not been) experi-
mentally measured, and which can be estimated by comparing simulation results with
remediation experiments. This optimization has been performed with genetic algorithms,
and the particular version of GAS used here is described in Section 7.
While this paper concentrates on theoretical aspects, it is important to know how the
model compares with experimental results: the present model has been able to precisely
describe the experimental data obtained on a pilot plant. The model testing which has
been performed is limited but meaningful. The major limitations are due to the size of
the equipment (a pilot plant of about 1 m3) and to the choice of a single contaminant,
namely phenol (which is most frequently encountered in practice). On the other hand,
an unusually high number of data has been recorded, which allow a thorough testing
of the model; moreover, different experimental conditions have been examined, thus
allowing to test the model in cases which had not been used for parameter optimization.
While referring the reader to [ 17, 191 for a more complete analysis, some comparisons
are reported in Section 8.
The differences between CA and PDE modelling are discussed in the final Section,
together with indications for further model improvements.
2. Soil bioremediation
Soil pollution is one of the major environmental problems in industrial countries:
indeed, although present day industrial technologies have achieved high environmental
standards, so that in many countries ongoing pollution has been largely reduced, soil
contamination has a long decay constant, and it is therefore necessary to remediate the
effects of a long period of low environmental care in plant and waste management.
Different technologies are available, including physical and chemical ones [ 12,311 while
134 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
the traditional use of landfills is being discouraged in many countries (landfills displace
the problem in space and time, rather than solving it).
The use of biotechnological methods, and in particular of the so-called in situ biore-
mediation techniques [2,30,32], is particularly attractive both economically and envi-
ronmentally. It is based upon the capability of many microorganisms in the polluted soil
to degrade organic contaminants. In most cases, indigenous bacteria are used, so that
bioremediation can be viewed as a way to accelerate and amplify natural phenomena.
The stimulation of indigenous microorganisms is achieved by providing them appro-
priate nutrients (e.g. oxygen, phosphorous, nitrogen, etc.), pH, moisture content, etc.,
in order to create in the soil environmental conditions which favour the growth and
the metabolic activity of those microorganisms which degrade the target contaminant.
Our model describes the (frequent) case where the nutrients are provided in aqueous
solution.
Summarizing, we will introduce a model for the simulation of in situ bioremediation
processes, using indigeneous bacteria which are stimulated by percolating an aqueous
nutrient solution through the soil.
These models may be very useful to improve the quality of today’s bioremediation
interventions, as they can allow a more reliable forecasting of the field behaviour,
starting from laboratory and pilot plant data. In a typical bioremediation operation, it
is necessary in the laboratory to assess whether there are indigeneous microorganisms
that can degrade the contaminant, to find satisfying conditions for their growth and to
assure that highly dangerous intermediate metabolites are not produced. After that, it
is necessary to resort to pilot plant studies, in order to determine the effects of the
soil matrix, which may profoundly affect the process (let it suffice to consider e.g.
the bioavailability of pollutants and nutrients). The pilot plant study can be performed
either on a constructed soil, i.e. a vessel which is filled with soil coming from the
contaminated site (also called microcosms or mesocosms, depending upon the size) or
on a small portion of the site, or both.
The results of the laboratory and pilot plant studies are then used to design the
full-scale field intervention. This sequence of tests at different scales is typical of many
areas of chemical engineering, and in most of these areas mathematical models are
widely used to achieve a proper scale up. In the bioremediation field it is not yet so,
and in most cases the process choice are based upon experience and heuristics (“black
magic”) rather than models. This is partly due to the complexity of bioremediation
models. Our model can be regarded as a contribution to the development of reliable
models of bioremediation processes, which take into account physical, chemical and
biological phenomena.
3. The model
While there are several models describing multiphase flow in porous media and
the fate of the solutes (for a review see e.g. [26,33]), few models exist which de-
scribe the wealth of interacting phenomena which take place in bioremediation (see
S. Di Gregorio et al. I T~eoret~cul Computer Science 217 (1999) J3I-156 135
[8,11,13,14,26,28] and further references cited in the last two). All these models make use of a continuum approach, while we adopt a discrete model whose general features will now be introduced; further comments on the relationship between the two approaches will be postponed to the final section.
As it was mentioned before, the different phenomena which interact in the system can be grouped into three classes; the model has a layered structure which reflects this division. l The lower layer is the fluid dynamical one, which describes the flow of different
fluid phases in a porous medium (the soil). The motion of one or more liquid or gaseous phases can be completely described at this level.
l A second layer describes the physical and chemical fate of the solutes, including hydrodynamic dispersion, interactions of the solutes with other solutes as well as with chemicals which may be present in the soil, and the processes of adsorption and desorption to and from the pore walls. The contamination of soil by a chemical or a mixture of chemicals can be described at this level.
l The third layer is the biological one, and it describes the interaction of biomass with its environment: the growth of microorganisms which depends upon the availability of nutrients, the disappearance and the appearance of chemicals which depends upon the metabolic activity, the interaction among different bacterial populations. A complete model should take into account also the possible modi~cations of the flow pattern which can be caused by huge bacterial growth leading to a reduction of the pore space available for the fluids. Note that this model can be regarded as “layered” in a specific sense, which differs
from the one used by [4]. In the present case, there is no hierarchy, like e.g. the one which might derive from progressively coarser graining; instead, the fluid dynamical layer is regarded as more fundamental as it is a necessary part of every meaningful model. For example, if we were to consider the flow of pure water in a sterile porous medium composed by water beads, the lower layer would be able to fully describe the phenomenon. If we add chemicals, either as solutes in the fluid phases or desorbed from the soil, which is still sterile, then a model with layers 1 and 2 would be needed. If we relax the sterile soil assumption, then all the three layers are needed. Even in this latter case (which is the one involved in all the experiments shown below), the behaviour of the variable of an “upper” layer may sometimes be itrelevent, and a simpler model suffices to describe the variables of interest (consider e.g. the percolation of a dilute aqueous solution, which is described to a very good approximation by layer 1 alone),
While the above framework is a general one, which should be applicable to any bioremediation process, this work describes a specific case, namely that of a phenol contaminated soil. The major simplifications of the general model structure, which have been applied in this case, are the following: l Two-phase pow: Only the flow of a gaseous phase (air) and a liquid phase (aqueous
solution) has been considered. If the contaminant had been less soluble, a third non- aqueous liquid phase might have been present (like it is often the case e.g. in the case of hy~~~bon con~mination).
136 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
l No chemical reactions: There are no indications of relevant chemical reactions taking
place in the system; anyway, the introduction of these reactions in a CA scheme
would be straightforward, as they take place within each cell.
a No predator-prey interactions: The bacterial population has been divided into three
classes, according to their behaviour with respect to the chosen contaminant. No
attempt has been made to provide a more detailed description of the different species
involved and of their mutual interactions. This simplification seems to be of wide
applicability.
l No pore clogging: The growth of microorganisms does not appreciably change the
fluid dynamical properties of the soil.
Apart from these limitations, it should be stressed that the model is very ambitious, as
it aims at describing a wealth of different interacting phenomena.
As it was mentioned, we will not try to describe the phenomena which take place
at the pore level, but will make use of macroscopic cellular automata, where each cell
comprises several pores.
We are interested in describing phenomena which occur in a finite region, where
some very important phenomena happen on the borders, as is the case for the example
of the upper soil layer in bioremediation by spraying a nutrient solution; the phe-
nomenon is induced by a flow of nutrients through the system boundaries. While the
PDE approach has extensively dealt with the study of boundary conditions, most CA
models involve infinite systems. In order to properly formalize the model for a finite
region, we will introduce a non-homogeneous cellular automata [41], which represents
the simpler and more natural approach. Let us notice that one could also describe the
system within a homogeneous CA formalism, at the expense of a more complicated
state space and transition function. In the following, we will always consider formally
infinite automata, and the inhomogeneities of the boundaries will allow us to consider
phenomena which take place in a finite region.
The specification of a cellular automaton (CA) involves its topology, its state space
and its transition function [27]; formally [20,23] a CA is a quadruple A = (G, I’, Q, f),
where
l G is the cellular space; usually G c Zd, the set of points with integer coordinates
in a d-dimensional Euclidean space (in order to deal with finite cellular spaces, one
sometimes uses toroidal topologies, see e.g. [40])
l V defines the set of neighbours of each cell: V = { 51, 52, . . . , (,}, where all the tk’s
are d-dimensional vectors with integer coordinates; without loss of generality, we will
consistently adopt the convention that (1 = (0, 0, . . . , 0) always; the neighbourhood of
cell i, N(i), is N(i) = {i + 51, i + 52,. . . , i + g,} n G = {i, i + 52,. . . , i + g,} n G; the number of elements in N(i) will be denoted m(i) < m.
l Q, the state space, is a discrete set; in the present case, as it has been stressed, we
will assume that the cardinality of Q can be high, to account for discretized versions
of variables which are considered continuous like, e.g. contaminant concentration.
l f is the transition function; let Y = {y ] y: G + Q}; each function y E Y defines a
possible state assignment in A, and will be called a configuration of A; let y(i) be
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 137
the state of cell i in configuration y, and let (Pi : Qmci) + Q; vi(y) = cpi({y(k) 1 k E
N(i)}) be the transition function for the elementary automaton, i.e. a rule which
determines the state of cell i at time t + 1 from the states of the cells in N(i) at time
t. Then the overall transition function f : Y --) Y is given by f(y) = {vi(y) 1 i E G}.
Whenever vi has the same functional form for every i, and N(i) has the same shape
for every i, the automaton is homogeneous, while we will consider in the following
also inhomogeneous systems, where inhomogeneities come from the existence of
boundaries, so that cells belonging to the boundary, or close to it, have a set of
neighbours which differs from that of the cells in the bulk of the CA.
Let us also define the notion of quiescent state: a state x of an automaton A is quiescent
if go(x, x, . . . ,x) = x, which means that if all the cells of the neighbourhood of cell i are
quiescent at time t, then i will be quiescent at time t + 1.
In our case the CA will be called Asoil =( G, V, Q, f).
Themodelisthree-dimensional: G={~I~EZ~~O<Y,QL,; O<r,<L,; O<r,<L,}.
This simple choice allows one to describe phenomena which take place in a simple
parallelepipedal geometry, like those which have been used for experimental purposes,
while obvious modifications would be needed for different geometries.
G can be partitioned into disjoint regions which reflect the physical conditions of
the system to be modelled:
where
l GUb={YIrEZ3 IO<r,<L,; Odry<L,; r,=L,} is the upper boundary of the sys-
tem (physically: the upper layer of soil cells);
0 Glb={rIrEZ3 10<rX<LX; O<r,<L,; r, = 0) is its lower boundary (physically:
the lower layer of soil cells);
l G,={r/r~Z3~(O=rx) u(rx=Lx)v(O=r,)v(ry=L,); O<r,<L,} denotes the
system lateral boundaries (physically: the soil cells which are close to the lateral
impermeable walls);
0 Gb={rIrEZ3 IO<r,<L,; 0 <ry XL,; 0 -cr, CL,} is the set of the system’s inter-
nal cells.
The neighbourhood V of each cell is composed by itself and by six other cells (up,
north, east, south, west, down); using relative coordinates:
The state space Q is the Cartesian product of different subspaces
Q = Qw x Qcap x Q,ra x Qmph x !&h x Qnrb X &nab X Qrdb
where Qw is the water content in the cell, QcaP is the “capillary water” substrate, i.e. the
water which is held by capillary forces (see Section 4), Qgra is the “gravitational water”
(see Section 4), Qmph accounts for the mass of contaminant adsorbed or precipitated in
138 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
the cell per unit mass of dry solid, Q@, accounts for the contaminant concentration in the cell water, &t, accounts for the concentration of bacteria which are not resistant to the con~min~t, Q&b accounts for the concen~atio~ of those bacteria which are resistant to the contaminant but are not able to degrade it and Q& refers to those bacteria which are able to resist to and to degrade the contaminant. As it has been mentioned, it is assumed that each of these variables can take a (large but) finite number of discrete non negative values.
The main features of the dete~inistic transition function f : &’ -+ Q, i.e. the law which determines the state of a cell at time t + 1 from the state of its neighbouring cells at time t, are described in Sections 4-6.
Updating is synchronous; if Xi(t) E Q is the state of cell i at time t, and N(i) denotes the neighbourhood of cell i, then for each i E G and each time t > 0: Xi(t +
1) = “fxX_ift)h j EN(i). In order to describe the boundaries of the experimental vessel, boundary cells are
given a neighbourhood different from the bulk cells (due to the intersection of G and V) and also a different transition function; for example, the upper layer can take an inflow of water from the outside, while the lower layer describes the collecting apparatus at the bottom of the pilot plant, and the lateral boundary cells do not allow lateral flux through the impermeable walls. More details can be found in [17].
4. The fluid dynamical layer
Let us consider what happens when water sa~tion SW (which is defined as the ratio of volume of water to the volume of voids) increases, in the case where water is the wetting phase. When SW is low, a minimum quantity of water surrounds the grains of soil, there are no connected paths and water is therefore immobile. As water saturation increases, it covers the pore walls, until connected paths are formed, so that water can Aow; in this condition, water is mainly influenced by its micr~nviro~ent (surface forces, etc.). If water saturation increases further, the water motion is easier and gravitational effects become important.
The main phenomenological assumption of the fluid dynamical model, which is based upon several experimental studies of imbibition and drainage [5], is therefore that different kinds of water can be identified in a single cell: l immobile (or irreducible) water, which cannot flow l capillary water (also called diffusion water), which can flow unaffected by gravity l gravitational water, which can flow more rapidly under the action of gravity. Moreover, as in some cases macroscopic fissures can be found, these can also be described in a way similar to gravitational water, but with a faster flow rate.
An important problem concerns the way how the overall water (the only physically measurable quantity) is divided into the three kinds of water: from the above qualitative description, the rule is induced that all the water is regarded as immobile until a certain threshold is met. Above this level, water is capillary until a second threshold is found.
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 139
Cdl: i-l
Cell: i
Cell: i+ 1
n IRR DIFF /I- Z-ILL
Fig. 1. IRR indicates irreducible (or immobile) water, DIFF refers to diffusional (or capillary) water, TRAS
to transport (or gravitational) water and FSR to fissural water.
Above this latter threshold, water is considered gravitational (if appropriate, part of the
water can be regarded as fissural).
The relationship among the different kinds of water is described in Fig. 1.
In the following, unless otherwise stated, fissures will not be taken into account.
According to the model, the motion of capillary water tends to equalize its value
in neighbouring cells. This equalization is achieved with a CA rule which distributes
the excess water of a given cell among those of its neighbours whose water content
is smaller than the local average water content. The procedure for capillary water is
shown below.
Procedure new_capillary_water_flow (...)
/* “0” is the central cell */
/* “1” ) “2”,... are the other cells in the neighbourhood *I
/* “w[i]” is the capillary water content of the i-cell */
/* “max_w[i]” is the maximum quantity of capillary water
content (of the i-cell) */
140 S. Di Gregorio et al. I Theoretical Computer Science 217 (19993 131-156
/* “new_outfc[i]” is the capihary water flow toward the i-cell */
I’* Individuation of the neighbou~ng cells (the so-called “eliminated cells”) */ for i from 1 to 6 do:
eliminated[i] :=I FALSE
do: new-control := FALSE sum-w := w[O] count := 1
for i from 1 to 6: if NOT e~iminated[i]
snm_w := sum-w + w[i] count := count + 1
average-w := sum-w/count
for i from 1 to 6 do: if (w[i] >average_w) AND (NOT eliminated[i])
new-control := TRUE eliminated[i] := TRUE
until new_control = TRUE /* Computing the values of the outflows toward the not eliminated cells */
for i from 1 to 6 do: if NOT eliminated[i]
new_outfc[i] := average-w-w[i]
else new_outfc[i] := 0
f* Correction of the outflow computed values according to the clock performed by the constant capsate */
for i from 1 to 6 do: new_outfc[i] := new_outfc[i]/cap_rate
/* Correction of the outflow computed values according to the ratio of the capillary water inside the cells to the maximum quantity of capillary water */
I* See Fig. 3 for a representation of sigma function *f
for i from 1 to 6 do: new_outfc[i] := new_outfc[i] * sigma(w[i]/max_w[i])
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 141
1 1
30
~
13 17 7
3
Average: 70/5=14 Average: 40/4=10
I
Average: 2113=9
I
l **
=-i-:
13 17 7
3
Fig. 2. The algorithm for capillary water (in a ho-dimensional case, with Ni = 5).
Let us illustrate the main features of the algorithm; first of all, for each cell, a maximum quantity of transmissible water is computed as follows.
Let I$ denote the neighbourhood of cell i, which is composed by A$ cells, and let qi denote the capillary water in cell i; an average (q)i between all the cells of the
neighbourhood is first computed:
The neighbouring cells whose water content is higher than such average are identified and discarded. Let y(t) be the subset of I$ composed by the A$‘) cells which have not been eliminated:
#l) = {k; qR d (4)i).
Then the average between the remaining cells is computed:
(q)!l’=l c q. 1 N,U) J’
1 ,E r”’
The process is iterated until no more cells are found which are higher than the average (so that 5 @+‘) - V”‘). the value of this latter average is attributed to the central cell - i , and to all the cells which have “survived”, as shown in Fig. 2.
The algorithm here described is not limited to the case of capillary water, and can be applied to describe several processes where a conserved quantity tends to a
142 S. Di Gregorio et al. I Theoretical Computer Science 217 11999) 131-156
0
0 W/P 1
Fig. 3. The dependence of the capillary Row rate coefficient k Cap upon the fraction of capillary water present in the cell (capillary water/m~im~ capillary water) (sigma function).
homogeneous eq~lib~~ dis~bution, It has already been applied to describe l~dslide d~~ics, lava flow and heat diffusion [6,21,29]. It has also been shown that in some cases this algorithm provides a good approximation to the solution of the heat diffusion equation [29].
The quantity computed according to the previous algorithm is the maximum (capil- lary) water which can flow from the central cell; however, if we let all this water flow at once, we would implicitly fix the clock of the CA so to match the time constant of the capillary water flow, and we would thus loose the possibility to describe the actual system kinetics. Therefore, at each CA time step, only a fraction k,, of the excess water so computed actually flows: kinetic parameters of this kind allow the desc~ption of phenomena which have di~erent time constants (e.g. capillary water flow is slower than gravi~tional water flow).
It is well known [7] that the flow rate of a given phase depends upon its sa~ation level. Based on phenomenological considerations and experimental data about the so- called relative permeability curves, the dependence of the flow rate coefficient upon saturation has been assumed to be S-shaped, and it has been approximated by a three straight line shape like that depicted in Fig. 3.
The gravitational water simulates the water fall from the superior cell to the inferior cell under the action of gravity (or, more generally, a flow driven by a difference in hy~ulic head), A coefficient of relative permeability for each cell is calculated: the flow of water to the inferior cell is propo~ional to the smaller of the two pe~eabilities and in each case cannot be higher than the receptivity (i.e. the quantity of water that a cell can accept) of the inferior cell.
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
The procedure is listed below.
143
Procedure new-ve~ical_~~tationalflow (. . .)
/* “0” is the central cell */ /* “6” is the inferior cell */ /* “tot_w[i]” is the total water content of the i-cell */ /* “por[i]” is the porosity of the i-cell */ /* “‘perm[i]” is the permeability of the i-cell */ /* “gw[i]” is the ~avitational water content of the i-cell */
f* The correction of the outflow computed values according to the clock is performed by the constant grav_rate */
perm[O] := sigma(tot_w[O]/por[O]) perm[6] := sigma(tot_w[6]/por[6]) if perm[O] < perm[6]
permmin := perm[O] else permmin := perm[6] maxflow := permmin * gw[O] if max_fiow -K receptivity[6]
new_outfg[6] := receptivity[6] / gravrate else
new_outfg[6] := max_flow ,/ gravrate
. . . . . . . . . . . .
For the same reasons just discussed for capillary water, also the gravitational flow rate constant depends upon saturation in a way similar to that of Fig. 3 (a S-shaped function with different parameters - see Fig. 3).
The procedure for fissural water flow will not be given here, as the examples de- scribed later refer to the case without fissures; it can be found in [ 171.
5. The dynamics of chemicals
The behaviour of the different chemicals which are transported by the aqueous phase (i.e. nutrients, contaminants, others) is described at this level. Solutes are transported by the water flowing in the system. However, there are also interactions with the pore walls and with the water which was already present in the cell. The adsorption/desorption process has been modelled by first-order kinetic equations, where the adsorption (re- spectively, deso~tion} rate to (from) the pore waits is propo~ional to the contamin~t bulk concentration (adsorbed concentration).
144 S. Di Greg&o et al. I Theoretical Computer Science 217 (1999) 131-156
Uncontaminated water
Fig. 4. The two kinds of water used to describe the contamination phase.
It should also be recalled that immobile water can play a key role in contamination and remediation. Let us describe only the contamination phase, where an aqueous contaminant solution percolates through the soil (actually, in the test apparatus, the amount of water used for contamination was smaller than the overall immobile water of the soil): according to the spirit of the CA approach, only two kinds of water are considered, ~o~ta~in~t~~ water and not c~~ta~~~~te~ water (Fig. 4). In order to describe the kind of experiments performed, it has been assumed that only contaminated water flows from one cell to another (of course, if it were necessary to describe different phenomena, like the flow of pure water through a contaminated soil, this simplification should be avoided). It is assumed that, as time passes, the contamination front advances within the cell, so that a larger portion of water becomes contaminated; mass balance is assured by determining the contaminant concentration as the ratio between the quantity of contaminant and the quantity of contaminated water in the cell.
The procedure for contamination is listed below.
Procedure new~olute_con~entmtions (. . .)
/* “‘0” is the central cell */ /* “‘tot_w[i]” is the total water content of the i-cell */
I ,
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 145
* “soil[i]” is the soil mass of the i-cell */
* “cont_w[i]” is the contaminated water of the i-cell */
/* “ph_w[i]” is the phenol mass in the water of the i-cell */
/* “ph_s[i]” is the phenol mass in the soil of the i-cell */
:
* “c_ph_w[i]” is the phenol concentration in the water of the i-cell */
* “c_ph_s[i]” is the phenol concentration in the soil of the i-cell */
/* alpha (beta) is the rate of the disappearance of phenol in the water
(in the soil) */
/* w_cont is the rate of the water contamination *J
. . . . .
flow_to_soil := alpha*c_ph_w[O]
flow-from-soil := beta*c-ph_s[O]
ph_w[O] := ph_w[O]-flow-to_soil+flow_from_soil
ph_s[O] := ph_s[O]+flow_to_soil-flow-from-soil
cont_w[O] := cont_w[O]+w_cont*(ph_w[O]/cont_w[O])*( 1-(cont_w[O]/tot_w [0])
c_ph_w[O] := ph_w[O]/cont_w[O]
c_ph_s[O] := ph_s[O]/soil[O]
. . . . . . .
The introduction of chemical reactions is straightforward, as these reactions take
place within each cell and do not involve transport among neighbours. In the present
version of the model we limited to monitor the disappearance of phenol, without
providing a quantitative description of the degradation path and of the intermediate
metabolites.
6. Biological phenomena
In environmental microbiology one is usually faced with microbial consortia, com-
posed by different species, instead of pure cultures which are typical of the laboratory
or engineered bioreactors. Determining the different bacterial strains and their rela-
tive proportion would require long tests; moreover, these proportions undergo major
changes in time, and microorganisms grown in laboratory conditions (like in plate
counting techniques) may lead to results which differ from those which are found
in situ [34]. Therefore, it is necessary to resort to a more lumped description. While
most existing models limit to one single kind of bacteria, whose parameters represent
the collective behaviour of the consortium, we considered three kinds of microorgan-
isms: those which are able to degrade the contaminant, those which are able to survive
146 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
in the presence of the contaminant, but not to degrade it, and those which are killed
by the contaminant.
The time evolution of these different populations can be at least roughly estimated,
from an experimental viewpoint, using plate counting techniques (although these may
underestimate the real number of bacteria sometimes even by two orders of magnitude
[31]) on selected culture media: for example, the degrading population is estimated
by using the medium where the chosen contaminant is the only carbon source, the
resistant population is estimated by using a medium where the contaminant has been
added, etc. (see [5] for further details).
The evolution of degrading bacteria in each cell is ruled by the following law:
Bd(t + l)=Bd(l)(l +KNd) - PdBd(tj2 +&d(t) F(t)
A + F(t) - &dBd(t)F(T), (6.1)
where &j(t) is the concentration of degrading bacteria at time t, F is the contaminant
concentration in the soil, c1 and A are the coefficients of the Monod growth term
(Monod kinetics has been supposed to be appropriate for phenol metabolization in
the concentration ranges which have been experimentally investigated, see also [38]),
/$j is the coefficient of the generic death term for bacteria, &d is the coefficient
of the birth term, &d is the coefficient of the death term caused by contaminant
(the chosen contaminant, phenol, has long been used for disinfection purposes). When
contamination by phenol is simulated, the phenol-related death term is immediately
active, while the Monod growth term does not immediately become active, because
the bacteria need some time to activate the appropriate metabolic pathway (like e.g.
in substrate induced enzyme synthesis). This delay is modelled by a “clock”, so that
the coefficient of the Monod term takes its nonzero value only after a suitable time
interval has elapsed since the appearance of phenol in the cell. Similar remarks apply
to the role of water: when the soil is very dry, only a few “survival forms” resist for
a long time, but when water is added the bacterial population starts to grow. This is
also modelled by a “clock” operating on the birth rate term.
The disappearance of contaminant in a cell, due to bacterial metabolism, is
F(t + 1) = F(t) - d&(t) F(t) A + F(t)'
where u’ is the phenol consumption rate (~‘>a, as not all the carbon coming from
phenol is used to produce new biomass).
The law for the phenol-resistant (but not degrading) bacteria is similar to that of
Eq. (6.1) without the Monod term:
&(t + 1) =&(t)(l + KNr) - /$h(# - K~r&(t)F(t). (6.3)
The form of the equation for the non-resistant bacteria B,(t) is similar to that of the
resistant but not degrading ones B,(t), the only difference being that the phenol-related
death term is larger (Z& <<Km) so to lead to an almost vanishing final concentration
in the case of high contamination.
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 147
The effect of nutrients has been modelled as an influence upon the kinetic parameters
of the growth equation. For example, as the phenol degradation pathway is aerobic, it
is assumed that the coefficient of the Monod growth term CY. be a growing function of
the available oxygen concentration.
As far as bacterial motion is concerned, although it is known that bacteria can
sometimes move along with the fluid flow, either as individuals or colonies (see [26]
for a review), and that more complicated behaviours have also been observed, like
e.g. chemotaxis, which implies a selective movement towards the regions of higher
concentration of a given chemical, we took an “Occam razor” approach and choose to
try to keep the model complexity at a minimum level, compatible with experimental
data: therefore bacteria are assumed to be immobile, bound to the pore walls. Bacterial
transport in the aqueous phase is supposed to be negligible.
7. Optimization methods
As there are in the model some parameters which cannot be directly measured or
derived from first principles, it is necessary to estimate them by comparing simulation
results and experimental data and applying a suitable optimization method.
In this work we resorted to the well known genetic algorithms [24,25], as they are
not bound to get stuck in local minima like, e.g. gradient descent techniques and do
not require that the cost function to be minimized possess any peculiar property (like
continuity, differentiability, etc.). The parameters are coded in binary form, so that
each set of parameter values represents an individual, and a random initial population
is generated. The genetic algorithm then tries to produce new generations of individ-
uals of higher fitness, i.e. able to give a closer agreement between simulations and
experiments, by applying the classical genetic operators (i.e. selection, crossover and
point mutation). In this work the selection probability is proportional to the fitness f of each individual, which measures the “goodness” of the simulation performed with
those parameter values, and is defined as follows. Let S and 0 be two M-dimensional
vectors whose components are the values, at M time steps, obtained in the simulation
(S vector) and in the pilot plant experiment (0 vector). The fitness function is then
defined as the (normalized) square of the euclidean distance between these two vectors:
f = II0 - sl12/llol12. The crossover operator which has been used is the classical single point crossover,
while point mutation is introduced by mutating the value of each single bit with a
fixed probability (mutations are independent of each other).
In order to prevent the risk of premature convergence to a suboptimal solution
which may affect genetic search, competition for reproduction is limited to subsets
of the whole population. For this purpose, the individuals are regarded as distributed
in space, interacting at each time step only with their neighbours and not with distant
individuals. More precisely, the individuals are arranged on a two-dimensional lattice
(with wrap around, thus giving rise to a toroidal topology); at each step an individual
148 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156
is chosen (either in a probabilistic way or in a sequential deterministic way) and its neighbourhood is isolated from the remaining population: selection, crossover and mutation will only act on the selected zone. The children strings that will have been generated from the genetic operators will come in conflict only with the individuals of the subpopulation and (if fairly strong) they will replace a part of it; after this, a new individual will be chosen, and the same procedure will be applied, a.s.o. A more precise description of the algorithm is given below, while further details can be found in 1191.
Procedure genetic (. . .) I * “tag.x”, “tag.y” are the central coordinates of the new subset of the
population *f i* “nei~bours[i]” is the i-individual of the neighborhood *I /* “‘random(x,y)” return a integer value between x and y */ I* "n_identi fy(. . . ),’ individuate the neighbouring individuals *I /* “end_criterion(. . .)” checks the condition of end program */ . . . . . . . . . . ..*...*..*... control:=FALSE do:
tag_x:=random(O,XMax) tag_y:=random(O,YMax) neighbours:=n_identify(population) genetic_action(tag~, tag-y, neighbors) con~ol:=end-c~te~on(population)
until control = TRUE . . . . . . . . . . . . . . . . . . . . . . . End. Procedure genetic_action(. . .) i* “mate 1” “‘mate2” are the two new children strings *I I* ‘“objet&e” is the position of the competitor string of the two new
strings *I /* “neighbours[i]” is the i-individual of the neighbourhood */ /* “selectf.. .).’ returns an individual, selected proportionally to its fitness */ /* “select-bad (. . .>,, returns an individual, selected in inverse ratio to its
fitness *I /* “cross(. . .)” effects the single point crossover */ /* “calculate_fitness(. . .)” return the fitness of an individual */ /* ‘“wimrer(. . .),, select a winner string, proportionally to its fitness */
. . . . . . . . . . . . . . . . . . . . . . .
matel:=select(neighbours)
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 149
mate2:=select(neighbours)
cross(mate1, mate2)
matel:=mutation(matel )
mate2:=mutation(mate2)
fit[ l]:=calculate_fitness(matel )
fit[2]:=calculate_fitness(mate2)
objective:=select_bad(neighbours)
fit[3]:=calculatefitness(neighbours[objective])
neighbours[objective]:=wimrer(fit, matel, mate2, neighbours[objective])
. . . . . . . . . . . . . . . . . . . . . . .
End.
The layered structure of the model also allows one to optimize subsets of parameters
in different phases: in this way a large search space can be broken in more manageable
portions. We first simulated the fluid dynamical layer alone, thus optimizing the values
of the relevant parameters, then those of the contamination experiments, so to calibrate
the parameters of the second layer, and finally those of the bioremediation experiments.
In this way, at each stage a limited number of parameters was involved, and the length
of the corresponding chromosome remained manageable. Further details can be found
in [17,18].
8. Comparison with experimental results
Several experiments have been performed on pilot scale constructed soils, i.e. on two
twin vessels of about 1 m3, filled with soil, which came from an industrial site which
had been contaminated by phenol for several years, and which had already undergone
a successful field scale bioremediation intervention. Two kinds of experiments have
been performed, the first referring to the contamination process and the second to
the (engineered) bioremediation process. In the first kind of experiments, an aqueous
phenol solution was percolated through the soil, while in the subsequent bioremediation
phase an aqueous nutrient solution was percolated in order to stimulate the activity of
indigeneous bacteria. Further details can be found in [5].
The two twin pilot plants were contaminated in two slightly different ways, i.e. in
a case with an instantaneous flooding from the upper side of the container, and in the
other case with a slower flow generated by a spraying system (in this case the same
amount of water was provided in about 45 min).
The bioremediation experiments were based on the use of a nutrient solution which
was sprayed on the soil, while field operations are usually performed in a closed loop,
where percolated water is collected by wells, added with oxygen and nutrients and
re-injected in the field, the tests were performed in an open loop way to ease model
150 S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) i3i-1.56
16 I
0 20000 40000 60000 80000 100000 minutes
Fig. 5. Total percolated water (MRI, flooding).
6
t 4
t ‘- 3
- Simulation
--c-- Real
0 50 100 150 200 250 300 350 400 450 500
minutes
Fig. 6. Total percolated water (MRI, fiocding).
validation. In order to assess the role of nutrients, while one container was sprayed with a nutrient-rich solution, the other was sprayed with water only.
Fig. 5 shows a comparison of model and experimental results concerning the total percolated water on a long time scale (in the contamination experiment). Fig. 6 shows a similar comparison on a short time scale: it can be seen in both cases that the two differ at most by a few per cent. As integrals tend to smooth differences, in order to provide a more demanding test, the flows are shown in Fig. 7, where it can be seen that the model predictions are accurate. Fig. 8 reproduces the short time data for the case of contamination by rain.
s. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-156 151
90
80
70
60
--r Simulation
-Real
0 100 200 300 400 500
minutes
Fig. 7. Flow of percolated water in MRl (flooding)
0 100 iU0 30 400 500
Fig. 8. Total percolated water (ME, rain). The simulation refers to the 3D model with ponds.
While the above data all refer to contamination experiments, it is interesting to note
that the model proved able to describe also the percolation dynamics in the bioreme-
diation phase (Fig. 9), where the amount of water added every day was much smaller
than that of the contamination phase, without fnrther parameter modification.
Also for the phenomena described by the other two layers, the simulation proved able
to match very closely the experimental results. While a detailed comparison is provided
in [17,19], to which the interested reader is referred, the major accomplishments have
been:
l a precise description of the outflow of phenol with drainage water, in different
conditions (flooding vs. rain, different initial values of soil saturation)
l a behaviour of phenol in soil which matches the main observed experimental features,
i.e. (i) phenol density decreases with increasing depth, at each time, and (ii) it has
a faster decay rate in the upper soil layers; the combined effect of the two leads
152 S. W Ckyorio et al. / Theoretical Computer Science 217 llYYY) 131-156
25
0 2oooo4cQoofzOooom 1cKloca
Minutes
Fig. 9. Total percolated water in bioremediation (MRI ),
to an almost simultaneous disappearance of phenol in all the soil layers (a precise
quantitative comparison between model and experimental results was not allowed
by the latter’s high variability in different points)
a behaviour of the concentration of bacterial populations which closely follows the
one observed experimentally, under different conditions including water introduction
in a dry soil, contamination by phenol and remediation by a nutrient solution.
Conclusions
The agreement between experimental data and simulation results is remarkable, tak-
ing into account the complexity of the system under study and the fact that the phe-
nomena are fully non stationary. This agreement provides an a posteriori test of the
validity of the cellular automata approach.
Let us come back to the comparison of CAs and PDEs. In continuum models, a
macroscopic approach is also taken; it is assumed [7] that a representative volume
element exists, which is small enough (with respect to the scale length of the phenom-
ena of interest) to allow a meaningful limiting operation AV + 0, yet large enough to
comprise several pores so to allow the use of average quantities which vary smoothly
in space. The water content of a cell at a given point would then be the average of
the water content on a number of neighbouring pores, and it would therefore be a
smoothly varying quantity.
The PDEs are obtained by performing mass and momentum balance on a volume
element AV for a small time At, and then letting both AV and At tend to zero.
Constitutive equations are also used, like e.g. Darcy law or its multiphase analogous.
If it were possible to integrate the Navier-Stokes equation this would not have been
S. Di Gregorio et al. I Theoretical Computer Science 217 (1999) 131-1.56 153
necessary, but the disordered nature of soil and the lack of information on its internal
structure makes it necessary to resort to empirical laws.
The major difference with respect to the discrete approach described in this paper
lies in the limiting operation. This has of course the advantage of eliminating any
dependence on the space scale (provided that the conditions for a REV are fulfilled)
and on the shape of the cell. If analytical solutions could be found, the advantage of
the continuum approach would be overwhelming. However, for systems as complex
as those which are considered here, this is not the case, so that, in the continuum
approach, a limit for d V, At -+ 0 is followed by a discretization with finite time and
volume increments. It is likely that these combined operations do not provide any
particular advantage with respect to discrete modelling, in the case of soil remediation.
Here, we have only provided empirical arguments in favour of this hypothesis, and
further investigation would be needed to assess its validity and to better understand
which phenomena lend themselves to such a discrete approach.
Moreover, there is a subtle difference between CA-based and PDE-based modelling:
when using a CA model, one is free to adopt phenomenological assumptions which
may be well suited for the case, but uncommon in PDE modelling. Let us consider,
e.g. the case of the bacterial clocks, which are mentioned in Section 6, like the phe-
nol clock: the coefficient of the Monod growth term for degrading bacteria, related
to contaminant concentration, is nonvanishing only after a certain time has elapsed
since the appearance of phenol (even for a bacterial flora which is already adapted
to the contaminant). This behaviour is well known to biologists as “lag phase”. In a
PDE system, this could be simulated by introducing a discontinuous coefficient (e.g. a
Heaviside function), or approximating it with some S-shaped steep function. In both
cases, numerical instabilities might arise, so that modellers tend to avoid these terms
(which are absent in the existing models, reviewed in [ 11,261). While they may not be
necessary under some operating conditions, they are needed to describe some transient
phenomena which have been observed.
The presence of modifiable parameters, although it may be disturbing, it basically
unavoidable in this kind of macroscopic models of complex interacting phenomena. Let
us also recall that the testing has been very demanding, involving different experimental
settings and conditions, so that we think that the good results achieved are not just “a
matter of fitting”. During model development we often could not achieve a satisfactory
performance until we identified the dominating physical phenomena: for example, be-
fore understanding the crucial role of immobile water on the fate of the contaminant, it
had not been possible to provide good results for the phenol concentration in drainage
water, no matter which adsorptionjdesorption coefficients were chosen.
Let us also remark that the model described here has been implemented on a MIMD
parallel computer, achieving a very high efficiency (80-90% with up to 32 processors)
[16]. This result proved that this class of models can be efficiently run on parallel
hardware.
Let us recall that the model described here is well suited for the case which has
been also experimentally investigated, namely contamination by phenol, while further
154 S. Di Gregorb et al. i Tkoretical Computer Science 217 (1999) 131-156
terms should be introduced to deal with a wider set of cases; some of them are listed below. l It is necessary to consider different soil types. l It is necessary to consider different contaminant types. In particular, hydrocarbons
can give rise to a three phase case (air, water, oil) and the model needs to be modified accordingly.
a The bacterial equations are fairly classical, and better models of bacterial dynamics might improve the overall model.
l In some cases bacterial growth can alter the fluid dynamical properties (pore clog- ging), and this phenomenon is not yet present in the model.
l Mobile bacteria (including, if necessary, ehemotaxis) could be introduced. All these model improvements can be inserted in the framework presented here. The major issue which is still open concerns the accuracy of the scaling of the model to the field, which of course requires expensive large scale measurements. The tests which have been performed, albeit on a limited scale, are however very stringent, and have provided a thorough testing set. We believe that these results, which show how even such a complex phenomenon can be precisefy described by a model, can facilitate the diffusion of model based techniques in the sector of in situ bioremediation, which is still largely dominated by a highly empirical approach.
This work has been partially funded by the UE, Esprit project, Capri initiative, subproject Caboto. We gratefully acknowledge the contribution of many friends and colleagues who contributed to the successful completion of the Cabot0 project: Rocco Rongo and William Spataro from the Universita della Calabria, Giandomenico Spezzano and Domenico Talia from CRAI and CNR, Alessandro Ban& Federica Abbonda~i, Antonella Iacondini, Nello Lombardi and Massimo Andretta from the Centro Ricerche Ambientali Montecatini.
References
[I] P.M. Adler, Multiph~e Flow in Porous Media, Kluwer, Dordrecht, 1995.
[Z] K.H. Baker, D.S. Herson, Bioremediation, McGraw-Hill, New York, 1994.
[3] S. Bandini, G. Cattaneo, M. Rigotti, A cellular automation model for the hydrodynamics of percolating
systems, in: S. Di Gregorio, G. Spcuano (Eds.), Proc. ACRJ’94, CRAI, Cosenza, 1994. [4] S. Bandini and G. Mauri, Towards multilayered automata networks, in: S. Bandit& G. Mauri (Eds.),
Proc. ACRE 96, Springer, Berlin, 1997.
[5] A. Ban& R. Sena, F. Abbondanzi, M. Giorgini, A. Iacondini, S. Mazzullo, A detailed analysis of in-situ
bioremediation: pilot plant experiments with a phenol contaminated soil, Proc. 1st Intentat. Conf. - The
Impact of Jndustty on Groundwater Resources, Associazione Georisorse e Ambiente, Torino, 1996.
[6] D. Barca et al., Cellular automata for simulating lava flows: a method and examples from the etnean
eruptions, Transport Theory Statist. Phys. 32 (1994) 195-232.
[7) J. Bear, Hydraulics of Groundwater, McGraw-Hill, New York, 1979.
S. Di Greyorio et al. I Theoretical Computer Science 217 (1999) 13X-156 155
[S] P.B. Bedient, H.S. Rifai, Modeling in situ bioremediation, in: B.E. Rittman (Ed.), In situ Bioremediation, National Academy Press, National Research Council, Washington, 1993.
[9] M. Bonazountas, D. Kallidromitou, Mathematical hydrocarbon fate modelling in soil systems, in:
E.J. Calabrese, P.T. Kostecki (Ed.), Principles and Practices for Petroleum Contaminated Soils, Lewis
Publishers, Boca Raton, FL, 1993.
[IO] A.W. Burks, Essays on Cellular Automata, University of Illinois Press, Urbana, 1970. [I l] C. Dacomc, Modelli di biorisanamento di suoli contaminati, Ph.D. Thesis, Universid di Padova, 1997. [12] E. DC Fraja Frangipanc, G. Andreottola, F. Tatano, Tcrrcni contaminati, C.I.P.A., Milano, 1995.
[13] S. Dhawan, L.E. Erickson, L.T. Fan, Model development and simulation of bioremediation in soil beds
with aggregates, Ground Water 31(2) (1993) 271-284.
[14] S. Dhawan, L.T. Fan, L.E. Erickson, P. Tuitemwong, Modeling, analysis and simulation of
bioremediation of soil aggregates, Environ, Prog. lO(4) (1991) 251-260.
[15] S. Di Gregorio, R. Rongo, R. Serra, W. Spataro, M. Villani, Simulation of water tlow through a porous
medium, in: S. Bandini, G. Mauri (Eds.), Proc. ACR1’96, Springer, Berlin, 1997.
[16] S. Di Grcgorio, R. Rongo, W. Spataro, G. Spezzano, D. Talia: High performance scientific computing
by a parallel cellular environment, Future Generation Computer Systems 12 (1997) 357-369.
[I71 S. Di Gregorio, R. Serra, M. Villani, A cellular automata model of soil bioremediation, Tech. Report
08/96, CRA Montecatini, 1996. [18] S. Di Gregorio, R. Serra, M. Villani, Combining cellular automata and genetic search in complex
environmental modelling, in: E. Pessa et al. (Eds.), Proc. 3rd European Congress on System Science, Edizioni Kappa, Roma, 1996.
[19] S. Di Grcgorio, R. Serra, M. Villani, Applying cellular automata to soil bioremediation modelling, II
Nuova Cimento D, submitted.
[20] S. Di Gregorio, G. Trautteur, On reversibility in cellular automata, J. Computer Systems Sci. 1 I (1975) 382-39 I.
(211 S. Di Gregorio et al., Landslide simulation by cellular automata in a parallel environment, in: Proc. 2nd
Intcmat. Workshop on Massive Parallelism; Hardware, Software and Application, 3-7 October 1994,
Capri, Italy, World Scientific, Singapore, 1994, pp. 392-407.
[22] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice gas automata for the Navier-Stokes equation, Phys. Rev.
Lett. 56 (1986) 1505-1508.
[23] E. Gales, Cellular automata, dynamics and complexity, in: P. Manneville et al. (Eds.), Cellular Automata
and Modelling of Complex Physical Systems, Springer, Heidelberg, 1989. [24] E.D. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wclsey,
Reading, MA, 1989. [25] J.H. Holland, Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor,
1975.
[26] A.A. Jennings, A. Manocha, Modeling soil biorcmcdiation, in: D.L. Wise, D.J. Trantolo (Eds.),
Remediation of Hazardous Waste Contaminated Soils, Marcel Dekker, New York, 1994. [27] C.G. Langton, Computation at the edge of chaos: phase transitions and emergent computation, Physica
42D (1990) 12-37.
[28] F.J. Molz, M.A. Widdowson, L.D. Benefield, Simulation of microbial growth dynamics coupled to
nutrient and oxygen transport in porous media, Water Resources Res. 22(S) ( 1986) 1207-12 16.
[29] C. Murmura, Automi cellulari e metodi di discretizrazione per lo studio della conduzione del calore,
Thesis, Universita della Calabria, Facoltl di Scienze, 1996. [30] R.D. Norris et al., Handbook of Bioremediation, Lewis Publishers, Boca Raton, 1994.
[3 I] J.L. Pepper, C.P. Gerba, M.L. Brusseau, Pollution Science, Academic Press, San Diego, 1996.
[32] B.E. Rittman (Ed.), In situ Bioremediation, National Academy Press, National Research Council, Washington, 1993.
[33] M. Sahimi, Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata
and simulated annealing, Rev. Mod. Phys. 65 (1993) 1393--1534.
[34] J.H. Slater, D. Lovatt, Biodegradation and the significance of microbial communities, in: D.T. Gibson
(Ed.), Microbial Degradation of Organic Compounds, Marcel Dekker, New York, 1984. [35] S. Succi, R. Benzi, F. H&era, The lattice Boltzmann equation: a new tool for computational fluid
dynamics, Physica 47D (1991) 219230.
[36] T. Toffoli, Cellular automata as an alternative to (rather than an approximation of) differential equations
in modeling physics, Physica IOD (1984) 117-127.
156 S. Di Greqorio et (11. I Theoretical Computer Science 217 (1999) 131-156
[37J T. Toffoli, N. Margolus, Cellular Automata Machines, MIT Press, Cambridge, MA, 1987. [38J C. Vipuianand~, S. Wang, S. Krishnan, Biodegradation of phenol, in: D.L. Wise, D.J. Trantolo (Eds.).
Remcdiation of Hazardous Waste Con~minated Soils, Marcel Dckker, New York, 1994.
[39J S. Wolfram, Physica 1OD (1984) 1.
[40] S. Wolfram, Theory and Application of Cellular Automata, World Scientific, Singapore, 1986.
[41] T. Worsch, Programming environments for cellular automata, in: S. Bandini, G. Mauri (Eds.), Proc.
ACR1’96, Springer, Berlin, 1997.
