Cellular Autometa Matlab algorithm

A parallel cellular simulator for bioremediation of contaminated soils S. Di Gregorio*, R. Rongo*, W. Spataro*, G. Spezzano§, D. Talia§ §CRAI, Località S. Stefano, 87030 Rende (CS), Italy *Dept. of Mathematics, Univ. of Calabria, 87036 Rende (CS), Italy

Abstract

The bioremediation of contaminated soils is one of main strategies for site clean-up. The most important principle of bioremediation is that micro- organisms (mainly bacteria) can be used to destroy hazardous contaminants or transform them into a less harmful form. Currently, we are facing this problem in the CABOTO project within the PCI ESPRIT framework. The CABOTO objective concerns the design and implementation of a parallel simulator for the bioremediation of contaminated soils by using models based on the cellular automata (CA) theory. For the parallel implementing of the simulator has been used the CAMEL system, a parallel environment for the simulation and modelling of complex systems based on CA. This paper describes the model used to simulate the contamination and the bioremediation of the soil, the main features of the CAMEL system and the parallel implementation of the simulator by CAMEL. Finally, experimental results are described.

1 Introduction

For the implementation of tools that allow the prediction of the behaviour of complex dynamic phenomena it is necessary to have models that can describe with accuracy the real systems. The implementation of these models often requires massive amounts of computation that can be only offered by parallel processing. Massively parallel computers allow to simulate scientific and engineering problems exploiting high performance. They represent the natural computational support for modeling and simulation of large complex phenomena and systems. The adoption of high performance computing simulation in a wide variety of sectors and activities demonstrates that this

technology can be used to improve complex industrial processes in terms of efficiency, safety and security of service. Soil decontamination is a major environmental problem in all industrial countries: while different decontamination methods are known, the evaluation of their effectiveness requires a number of laboratory or pilot plant tests. However, appropriate forecasting of real field results from these tests could be achieved only by resorting to mathematical models, which are fairly complicated and require high computing power. An important goal in order to face soil pollution problems is the ability to simulate contaminant transport, diffusion and transformation processes for reliably estimating their impact and forecasting the effects of different decontamination strategies. Several factors are involved when the full cycle is considered: contamination by liquid transport, adsorption, precipitation, chemical reactions, biological degradation kinetics and so on. If the classical approach of differential calculus is adopted, particular complexities could arise because of heterogeneous soil composition and because of the intricate physical-chemical-biological interactions. In fact, such a way is very hard to be followed because it is extremely difficult to solve analytically the equations for such phenomena without making simplifications that neglect in many cases substantial characteristics of the phenomenon. Here we propose a CA approach like Karapiperis & Blankleider [1], where the fluid flow with diffusion-transport of contamination-decontamination agents inside the soil and their mutual influences are viewed as a dynamic system based on local interactions with discrete time and space, where the space is represented by cubic cells. Each cell is characterized by specific values (the state) of selected physical parameters, representing physical-chemical specifications, relevant to the evolution of the phenomenon. A transition function is defined in order to describe state variation and chemical and biological interactions. This model has been implemented using CAMEL a parallel software environment based on cellular automata theory developed by Cannataro et al. [2].

2. Bioremediation model

In the specification of a CA for soil contamination-decontamination, we refer to the contamination of a porous soil caused by a phenol rich water solution. This problem involves two phases only (water and air) while of course in other cases it might be necessary to consider additional phases (e.g. an oil phase in the case of hydrocarbon contamination). The CA model that is used to simulate the phenomenon is composed by a large number of cells, where each cell being of mesoscopic size does not describe the phenomena which take place inside a single pore (typical pore dimensions being 10-00 mm) but those inside a portion of soil. Our model is three-dimensional and is necessary for real scale simulations, in order to describe heterogeneous regions.

The cellular automaton specification is given by:

SOIL = (R3, A1, A6a, A6b, X, Q, σ, γ1, γ6a, γ6b)

where

• R3 = {(x, y, z)| x, y, z ∈N, 0 ≤ x ≤ lx, 0 ≤ y ≤ ly, 0 ≤ z ≤ lz} is the set of points with integer co-ordinates in the finite region, where the phenomenon evolves. N is the set of natural numbers.

• A1 ⊂ R3, specifies the water source cells at air-soil contact, where contamination ends; such cells obey to the particular transition function γ1 substitutive, but compatible with σ. In cells of A1 the contaminated water appears.

• A6a, A6b ⊂ R3, specify the two last down layers of cells, which obey respectively the particular transition functions γ6a, γ6b different from σ.

• The set X identifies the geometrical pattern of cells which influence the cell state change. They are the cell itself and the "up", "north", "east", "west", "south" and "down" neighbouring cells, which are individuated by the indexes 0, 1, 2, 3, 4, 5, 6

X = {(0, 0, 0), (0, 0, 1), (0, 1, 0), (1, 0, 0), (-1, 0, 0), (0, -1, 0), (0, 0, -1)}

The finite set Q of states of the elementary automaton:

Q = Qw x Qp x Qoutf6 x QK_sat x Qcap_thr x Qc_ph_s x Qc_ph_w x Qb

where the substates are:

∗ Qw is the water content in the cell.

∗ Qp is the effective porosity parameter. It depends on physical- chemical characteristics of the soil in the cell and individuates the quantity of water drainable from the cell in condition of full saturation, i.e. the porous volume, which can be filled by water. For example, a cell with completely impermeable soil inside has the effective porosity parameter equal to zero, while the effective porosity of a cell is greater if the soil is more permeable.

∗ Qoutf represents the water flow toward the six neighbourhood directions from the central cell. Note that the inflows are not explicitly considered, they are obtained trivially by the outflows.

∗ QK_sat is the hydraulic conductivity at the saturation, considering that, inside the cell, the soil is homogeneous.

∗ Qcap_thr is the capillary water threshold substate, i.e. the maximum cell water which is effected by the capillary forces.

∗ Qc_ph_w, accounts for the phenol concentration in the cell water.

∗ Qc_ph_s, accounts for the phenol mass adsorbed in the cell per unit mass of dry solid.

∗ Qb, is the biological mass inside the cell.

• σ:Q7→Q is the deterministic state transition function for the cells in R3, a large outline of its specification will be given in the next sections.

• γ1:Q2 x N→Q is the transition function of source cells, where water and solutes have origin at each SOIL time interval t ∈ N ; according to the assigned feeding, the outflow to the inferior cell only is calculated.

• γ6a:Q2→Q is the deterministic state transition function for the cells in the first down layer, it specifies the behaviour of ideal cells, which receive water only from the superior cell and transmit it to the inferior one without retention.

• γ6b:Q2→Q is the deterministic state transition function for the cells in the second down layer, it specifies the behaviour of ideal cells, which receive illimitable water without transmit it at all.

Such ideal cells with strange transition functions γ6a, γ6b are considered in order to describe forced conditions in laboratory experimental devices. At the CA step 0 the initial configuration is defined, specifying all the starting values of the cell substates. Usually the water outflows may vanish. At each next step, the function σ is applied to all cells in R3 {A1 ∪ A6a ∪A6b}, while at the same time step the functions γ1, γ6a, γ6b are respectively applied to A1, A6a, A6b, so the configuration is changed and an evolution step of SOIL is obtained. Let us remark that contamination is a non stationary process, so that some usual simplifications cannot be used and really dynamic models are required.

3 CAMEL overview

CAMEL has been implemented on a parallel computer composed of a mesh of 32 Transputers connected to a host node. The current implementation of CAMEL does not limit the number of Transputers which can compose the parallel computer, so no changes should be necessary in the software of CAMEL whether a very large number of Transputer should be used. The system is composed by a set of macrocell processes each running on a single processing element of the parallel machine and by a controller process running on a master processor. Each macrocell process implements several

elementary cells and makes use of a communication system which handles the data exchange among cells. CAMEL offers the computing power of a parallel computer although hiding the architecture issues to a user. It provides an interface which allows a user to dynamically change the parameters of the simulation during its running, and a graphical interface which shows the output of the simulation. The approach taken in CAMEL is to make parallel computers available to application-oriented users hiding the implementation issues coming from their architectural complexity. CAMEL allows the solution of complex problems as described in Di Gregorio et al. [3], which may be represented as discrete across a square or hexagonal grid or lattice. CAMEL implements a cellular automaton as a SPMD (Single Program Multiple Data) program. CA are implemented as a number of processes each one mapped on a distinct PE that executes the same code on different data. According to this approach a user must specify only the transition function of a single cell of the system he wants to simulate by CARPET, a cellular language defined to program CA applications in CAMEL. CARPET (CellulAR Programming EnvironmenT) is a high-level language that extends C with some additional constructs to describe the rules of the transition function of a single cell of a cellular automaton. The main features of CARPET are the possibility to describe the state of a cell as a set of typed substates each one by a user-defined type, and the simple definition of complex neighbourhoods (e.g., hexagonal), that can be also time dependent, in a n-dimensional discrete cartesian space. The CAMEL system allows, by the Graphical Interface (GI), to rapidly and interactively explore and analyse very large amounts of scientific data gathered during the execution of computer simulations. A tool called IVT (Interactive Visualization Tool) has been added to CAMEL to improve data visualization. Utilizing data computed by simulation, IVT provides a variety of functions and services, including 2 and 3-dimensional graphical displays of data, hard copy of graphical displays and text, interactive colour manipulation, animation creation and display, rotation of the images, saving of data in files according to different data formats. Further, the tool allows quick access to the software for the novice as well as the advanced user and provides an intuitive easy-to-learn-and-use interface. One of the most well-known visualization software package is certainly MATLAB. It offers outstanding numerical computation functionality together with sophisticated graphical presentation. By MATLAB we designed IVT provides a window-oriented graphical interface for an interactive manipulation of graphical displays. The integration of IVT with the CAMEL system allows to extract more information from raw data gathered during the simulation using the MATLAB's graphics capabilities (i.e. 3D visualization of 2D images) and the numerical transformation functions. A software interface controls the communications between CAMEL and IVT. The communication occurs connecting dynamically, under Windows, the two applications. Further, this approach will provide the

remote visualization of data and the possibility to use IVT with other parallel applications.

4 Model implementation

The bioremediation model has been implemented using the CARPET language. The transition function σ is composed of three sequential steps one for each of phenomena that constitute the model. An evolution step, performed by each cell, of the transition function concerns:

1. computing the water outflows from a cell toward the adjacent cells;

2. computing the solute adsorption/desorption inside a cell and the solute dispersion from a cell toward the adjacent cells;

3. updating the cell substates considering the chemical-biological effects of biomass and its growth/reduction.

During the first step the new water substates values are computed according to the following calculations, which must be executed sequentially in each cell:

• first of all, the new water content in the cell is calculated, considering the water balancing between current content cell_w, inflows and outflows; inflows and outflows are determined by the contributions of capillary (infc and outfc) and gravitational (infg and outfg) inflows and outflows;

• the new outflows determined by the capillary forces (outfc) are calculated;

• the outflows determined by the gravitational forces (outfg) and the new global outflows (outf) are calculated.

The water outflow from a cell (or alternatively part of the inflows into the adjacent cell) depends on the capillary and gravity forces; we will account on the main features of the transition function for the water flow.

Capillary forces involve a suction effect, which covers the porous surface with a water layer. The maximum possible thickness of the water layer depends on the attractive inter-molecular forces between water, air and soil.

We can characterise the equilibrium by a suction threshold (cell_cap_thr), which is the minimum water quantity inside the cell above, when the suction effects are practically negligible. Part of outflows (outfc) from the central cell to its neighbours are originated because of balancing the capillary forces, when the value of cell_cap_thr is not reached in the neighbours. The capillary forces can be considered approximately proportional to the suction index suct_ind, here defined as ((cell_cap_thr - cell_w ) /cell_cap_thr), when cell_w < cell_cap_thr (0 otherwise). So a balance can be given by outflows, which minimise the total difference of the suction index between a cell (the central cell) and its six neighbours after each time interval.

Computation involves two sequential actions :

• Selection of the neighbouring cells (the so-called "eliminated cells"), toward which capillary water outflow is not possible, because their suction index is larger than the suction index of the central cell or because other neighbouring cells have by comparison so low suction index that the capillary water flows necessarily only toward such cells.

• Computation of the values of the outflows toward the not eliminated cells, in order that the same value of the suction index is assumed by the central cell and the not eliminated cells.

The gravitational water in a cell is the part of water which is practically not sensible to capillary forces. Its flow is ruled by the hydraulic conductivity, which depends upon the hydraulic conductivity at the saturation K_sat, the porosity p and the water content w according to a function depending on the type of soil. In order to calculate the flow between the cell and its neighbour down, the conductivity value, which is considered, is the minimum of the conductivities of the two cells.

0

2,5

5

7,5

10

0 100 200 300 400 500

minutes

litres

simul. exper.

Figure 1. The simulation data of water flow.

The simulation of the first hours of an experiment is reported in figure 1. The model captures a fluid-dynamical behaviour, difficult to be simulated. The experimental results reported here for comparison have been performed in a pilot plant, under carefully controlled conditions by Banti et al. [4]. The model involves a number of parameters, which can only partly be directly determined. The remaining parameter values have been found either by tuning them at hand, or by applying the novel optimization technique of genetic algorithms.

4.1 The solute flow through the soil

A dynamic description of the adsorption/desorption of solute onto the pore wall was found necessary, and two different techniques to achieve such dynamic description have been tested:

1. The solute distribution between bulk water and pore walls described by the Freundlich isotherm was considered as the goal, and a linear dynamics towards that goal was defined: at each time step, only a fraction of the difference between the actual and the goal values was allowed to move from the bulk to the wall or vice versa.

2. A mass balance equation for the solute was used, with a term, describing the flow from bulk to walls, proportional to the solute concentration in the bulk, and a term, describing the opposite flow, proportional to the concentration of the adsorbed solute.

Moreover, a further aspect that has been taken into account is the water moves from one cell to another, carrying with it its solute content, a mixing takes place with the immobile water (i.e. residual + irreducible) which is present in the target cell. Therefore the solute concentration (at least initially) decreases, as it is diluted in a larger amount of water. Since the immobile water may take fairly high values, this dilution phenomenon may be quite important. Moreover, it would be unrealistic to suppose that an instantaneous dilution in all the immobile water may take place, so it could be limited to a fraction of the overall immobile water. Here again a set of models has been considered: in all simulations the mixing with the immobile water has been introduced in the model, but in to two different ways, explicit or implicit. In the former case a parameter has been used, which describes the fraction of immobile water involved in mixing, while in the latter case the phenomenon has been implicitly taken into account by assuming a lower value of immobile water than the experimental one..

0

1 0 0 0

2 0 0 0

3 0 0 0

4 0 0 0

0 2 0 0 0 0 4 0 0 0 0 6 0 0 0 0 8 0 0 0 0

m in u tes

m g /l

sim u l.

exper.

Figure 2. The simulation data of drained water contamination.

Figure 2 compares the experimental data of the phenol concentration in solution with the simulation results

4.2 The biomass transition function

Strictly biological effects are here considered, i.e. metabolic processes, which alter the contaminant or the nutrient as well as bacterial reproduction processes. Bacterial colonies may be either mobile or substantially fixed. The motion of bacteria could be described by the same techniques useful for solute transport, adding possibly chemiotropism. So far not enough data are available to decide about the necessity to take into account bacteria transport in the case of interest here. We will therefore limit ourselves to describe the algorithm for bacterial metabolism and growth within each cell: this will suffice if bacterial transport can be neglected, while, if it is not the case, the appropriate transport equations should be also included. Looking to the Monod formulae, we propose to adopt any solution method suitable for a finite cubic closed volume, considering a time interval Δt, equal to the computation step of the CA.

5 Experimental results

We measured the water flow and contamination model execution times (seconds) on our Transputer-based machine for different number of cells per processor, iterations and processor configurations. Figure 3 shows the speed-up results for a 32x17x13 (7072 cells) grid with 10 substates for the water flow model and 12 substates for the contamination event model. Comparing the results of the simulation with the experimental data we deduce that 100 steps of the simulation correspond to about 12 minutes of the real phenomenon.

0

4

8

12

16

20

24

28

32

0 4 8 12 16 20 24 28 32

Processors

S p

e e d

u p Ideal

Water flow simul.

Contamination simul.

Figure 11. Speed-up measures.

Furthermore, we point out that the execution time of a single step is about 0,5 seconds. This result is very important because many iterations are necessary to complete our simulation on a long period basis (64 days). Consequently, it is

worth noting that 24 hours of the real phenomenon can be simulated in about 1 hour and 40 minutes. Figure 4 shows a snapshot of the water flow simulation by IVT tool described in section 3.

Figure 4. A snapshot of the water flow simulation by the IVT tool.

References

1. Karapiperis T & Blankleider B. Cellular Automaton Model of Mass Transport with Chemical Reactions, Bericht Nr.93-06 Paul Scherrer Institut, Villingen, Schweitz,1993.

2. Cannataro M., Di Gregorio S., Rongo R., Spataro W., Spezzano G. & Talia D. A Parallel Cellular Automata Environment on Multicomputers for Computational Science, Parallel Computing, 1995, 21, 803-824.

3. Di Gregorio S., Rongo R., Spataro W., Spezzano G. & Talia D. A Parallel Cellular Tool for Interactive Modeling and Simulation, IEEE Computational Science & Engineering,1996, 4.

4. Banti A., Mazzullo S., Giorgini M. & Serra R. A Detailed Analysis of in-situ Bioremediation: Pilot Plant Experiments with a Phenol-Contaminated Soils, in Proc. of the Conf. Water Resources and the Environment, Como, 1996.

Keywords: Cellular Automata, Parallel Processing, Simulation, Bioremediation.