Cellular Autometa Matlab algorithm

cellAutomataEnvironModelChapter.pdf

Home >Mathematics homework help >Statistics homework help > Cellular Autometa Matlab algorithm

and Consequences for Policy Making, © Springer Science + Business Media B.V. 2009

Cellular automata modeling of environmental systems

Salvatore Straface, Giuseppe Mendicino

Dipartimento di Difesa del Suolo, Università della Calabria, Ponte Pietro Bucci, Cubo 41 b 87036, Arcavacata di Rende (CS), Italy.

Abstract

Flow and transport processes in unsaturated soil are analyzed through a simulation environment based on cellular automata (CA). The modeling proposed in this chapter represents an extension of the original computational paradigm of cellular automata, because it uses a macroscopic CA approach where local laws with a clear physical meaning govern interactions among automata. This CA structure, aimed at simulating a large-scale system, is based on functionalities capable of increasing its computational capacity, both in terms of working environment and

taken from benchmarks in literature, showing a good agreement even in the cases where non-linearity is very marked. Furthermore, some analyses have been carried out considering quantization techniques aimed at transforming the CA model into an asynchronous structure. The use of these techniques in a three-dimensional benchmark allowed a considerable reduction in the number of local interactions among adjacent automata without changing the efficiency of the model, especially when simulations are characterized by scarce mass exchanges.

Keywords: flow and transport model, macroscopic cellular automata, parallel simulation, discrete approach, quantization techniques, model validation, environment.

1. Introduction

The capability of understanding and modeling hydrological processes at different spatial scales together with the need for a more detailed knowledge of mechanisms regulating soil interaction between surface and subsurface, has led experts to investigate and develop different forms of modeling.

in terms of the optimal number of processors available for parallel computing. Specifically, the performance of a three-dimensional unsaturated flow model has been verified comparing the results with reference multidimensional solutions

P.C. Baveye et al. (eds.), Uncertainties in Environmental Modelling

S. STRAFACE AND G. MENDICINO

In many cases, the models using fully-coupled system equations to describe the hydrological cycle of a basin completely have shown computational limitations, mainly due to both the reduced dimension of sampling-grids (or mesh), involving a narrow space domain, and smaller time steps which are essential to prevent eventual numerical instability.

These problems become even more evident during the simulation of soil dynamics because of the fast responses of this zone to atmospheric forcing, which limit us to a modeling characterized by very reduced space-time steps (Orlandini, 1999).

It becomes increasingly necessary to find alternative numerical solutions that, always under the aegis of describing physical phenomena, allow us to increase the spatial and temporal domain of simulations with acceptable computational requirements.

Modeling based on cellular automata (CA) represents a valid alternative to analytical-deductive methods based on the analysis of physical equations describ- ing a particular phenomenon and their subsequent resolution carried out through numerical methods (Toffoli and Margolus, 1987; von Neumann, 1966; Wolfram, 1986, 1994). In complex physical phenomena, this modeling allows us to capture the fundamental characteristics of systems whose global behavior is derived from the collective effect of numerous simple components interacting locally. Many CA applications in fluid-dynamics exist, most of these based on microscopic app- roaches: lattice gas automata models were introduced to describe the motion and collision of particles interacting according to appropriate laws on a lattice space (Frisch et al., 1987). These models can simulate the dynamic properties of fluids (Di Pietro et al., 1994; Pot et al., 1996) and their continuum limit leads to the Navier-Stokes equations (Rothman and Zaleski, 1997). Lattice gas models, due to the simplicity of both fluid particles and their interactions allow simulations of a large number of particles, but do not allow to express the velocity in explicit form, because an amount of fluid moves from one cell to another one in a CA step, which is defined as a constant time, involving a constant “velocity” in the CA con- text of discrete space-time. However, velocities can be achieved through the analysis of the global behavior of the system: in the space by considering clusters of cells; and in time by taking into account the average velocity of the advancing flow front in a sequence of CA steps.

Surface and subsurface flow modeling represent complex macroscopic fluid dynamical phenomena. Subsurface solute transport has been traditionally des- cribed by a deterministic advection-dispersion equation based on analogy to Fick’s laws of diffusion. According to this analogy, the spread of a nonreactive contami- nant in a hydrogeologic environment is controlled by a constant directional medium property called dispersivity which, when multiplied by absolute velocity, yields a directional dispersion coefficient. Flow and transport modeling seem dif- ficult to model in these CA frames, because they occur on a large space scale and need, practically, a macroscopic level of description that involves the management of a large amount of data.

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Lattice gas models are only capable of reproducing the birth of macroscopic flow patterns. The drawbacks of the lattice gas models have been overcome by the Lattice Boltzmann Method (LBM), where the state variables can take continuous values (instead of integer variables), as they are supposed to represent the density of fluid particles endowed with certain properties located in each cell (McNamara and Zanetti, 1988; Succi et al., 1991; Chopard and Luthi, 1999). In the LBM space and time are discrete and they represent a very powerful approach which combines numerical efficiency with the advantage of having a model whose microscopic components are intuitive. Also for the Boltzmann models the equivalence with the Navier-Stokes equations for incompressible fluids has been demonstrated (Higuera and Jimenez, 1989; Qian et al., 1992; Chen et al., 1992; He and Luo, 1997).

LBM has been used to model flow in porous media, but the applications so far developed usually take a more microscopic approach and aim at describing phe- nomena which take place at the pore level (Martis and Chen, 1996; Zhang et al., 2000; Knutson et al., 2001; Zhang et al., 2005). Some extensions to porous media and scaling up have been also proposed (Soll et al., 1993; Soll and Birdsell, 1998; Sukop and Or, 2004; Chau et al., 2005), as well as LBM has been used for macro- scopic plots modeling (Deng et al., 2001; Zhang et al., 2002a,b; Ginzburg, 2005), but quantities such as density, velocity and concentration are obtained by taking the moments of the distribution function (Deng et al., 2001).

Empirical CA methods were developed on the macroscopic scale in order to directly deal with the macroscopic variables associated to the representative ele- mentary volume of the physical system under consideration (Di Gregorio and Serra, 1999; Di Gregorio et al., 1999). In these methods an almost unlimited number of states are allowed and, each state is described by a Cartesian product of sub- states, each one describing a specific feature of the space portion related to its own cell. Furthermore, the transition function is split in different parts, each correspond- ing to an elementary process of the macroscopic phenomenon. When an amount of fluid is computed to pass from one cell to another, this particular characteristic of the transition function allows fluid to be added to the cell and, at the same step, to change all the sub-states involved, specifying the evolution of the main quantities of the macroscopic system such as velocity explicitly (Avolio et al., 2003).

But, these empirical CA methods make use of some local laws where automata interactions are based on parameters whose physical meaning is not clear and, as a consequence, heavy calibration phases are necessary to estimate suitable values of the same parameters. The computational effort which results is partially justified by the possibility of making cellular automata directly compatible with parallel programming (Toffoli and Margolus, 1987; Crutchfield et al., 2002).

In this chapter the same extended notion for macroscopic cellular automata has been considered to develop a three-dimensional model which simulates water flux in unsaturated soils, but the local laws governing the automata interactions are based on physically-based rules.

The model, developed in a problem solving CA environment called CAMELOT (Dattilo and Spezzano, 2003) has been used for different multidimensional schemes

S. STRAFACE AND G. MENDICINO

and provides results similar to those of other approaches described in Paniconi et al. (1991), Paniconi and Putti (1994) and Huyakorn et al. (1986).

Moreover, in the second part of the chapter a particular aspect regarding the use of theory of quantized systems (Ziegler, 1998) applied to the proposed CA model has been investigated. This concerns the computational benefits derived from the application of some quantization techniques aimed at decreasing the number of local interactions by neglecting scarce mass exchanges occurring between neigh- boring automata. Through these techniques CA act as an asynchronous system, where according to a common quantum size rule a given automaton evolves or is kept at rest depending on its state and on those of the adjacent cells.

2. Discrete formulation of flow in unsaturated soil

All existing numerical methods for the solution of field equations have a differen- tial formulation as their starting point. A discrete formulation is then obtained by means of one of the many discretization methods, such as Finite Difference Method (FDM), Finite Element Method (FEM) or, in general, a weighted residual or weak solution method. Even Boundary Element Method (BEM) and Finite Volume Method (FVM), which use an integral formulation, have a differential formulation as their substratum (Tonti, 2001). The formation of densities and rates and then the passage to the limit to form the field functions, that is typical of field and continuum theories, deprive physical variables of their geometrical content (for this reason any discretization process restored from the differential formula- tion is not based on geometrical but only on numerical considerations).

So, it follows that to obtain a discrete formulation of the fundamental equation of a physical theory, it is not necessary to go down to the differential form and then go up again to the discrete form. It is enough to apply the elementary physical laws to small regions where the uniformity of the field is attained to a sufficient degree, this being linked to the degree of accuracy of data input and to the obser- vation scale of the physical phenomenon. In this way, we obtain a direct discrete formulation of the field equation.

If we consider a discrete cell system, in which the cells have smaller dimensions where the departure from uniformity is larger, we may use the same constitutive law used in the differential context, as an approximation.

In this section our aim is to obtain a direct discrete formulation of the unsatu- rated flow, the equation that we use in CA model. This equation, whose solution

Following this approach, starting from the technique suggested by Muzy et al. (2002), dynamic quantization which makes CA structure both asynchronous and non-uniform is developed. In particular, the quantum size rule conditioning the tran- sition function is varied locally according to the sub-state describing the degree of saturation of the soil, and this produces a considerable reduction in the local transitions without adversely affecting the performance of the model.

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

provides the variables configuring the phenomenon, is the final result of the composition of the mass-balance equation and the constitutive equation (Darcy’s law). The discrete mass-balance equation is the following (Straface et al., 2006):

�c� � � �

1 �c

�mc �t

� Sc (1)

K ci � Vc � Vi

Vc Kc�

� Vi

Ki�

(3)

where Vc and Vi [L 3] are the volumes of the cells, while Kc� and Ki� are the ele-

ments on the diagonal of the hydraulic conductivity tensor corresponding to the direction � in the cells c and i, assuming that the xyz Cartesian system is colinear with the principal anisotropy directions.

The remaining terms of equation (2) are the total heads hc and h� [L], the cell dimension l� [L], the hydraulic gradient g� [–] and the surface area A� [L

2] where the flux �c� passes through (Figure 1b). Equation (2) is a particular form of a more general direct discrete formulation of the unsaturated flow equation. For more details see Appendix A.

where c is the cell where the water mass balance is performed, � is the generic direction obtained linking the center of the mass of cell c with that adjacent (see Figure 1a, for a two-dimensional case), �c� is the mass flux [L

3T–1], �c is the

fluid density [ML–3], cm is the mass inside the single cell c [M], �t is the time step

[T], �mc �t

[MT–1] is the mass-time variation in the cell, and Sc is the mass source

term [L3T–1]. Equation (1) is valid both for internal nodes of the discrete domain and for boundary nodes, avoiding the typical differential approach where differen- tial equations and boundary conditions are separated. The mass flux is given by Darcy’s equation, which in discrete terms for a cubic regular mesh becomes:

�c� � K c� h� hc

l�

�

��

�� A� � K c� g� A� (2)

where �cK is the hydraulic conductivity averaged between cell c and that adja- cent along the direction � (e.g., cell i in Figure 1), obtained assuming the energy dissipated between two adjacent cells equivalent to the energy dissipated in a ficti- tious cell containing them (Indelman and Dagan, 1993). Then, if we consider the hydraulic conductivity averaged between cell c and i (Figure 1b), we assume:

S. STRAFACE AND G. MENDICINO

k c i

p c

i �

A�

a) b)

�

k c i

p c

i �

A�

a) b)

�

Fig. 1. (a) Scheme of a squared two-dimensional cells complex; (b) Representation of the flux through the face of two adjacent cells.

The time variation of mass content changes according to the flow type consid- ered. For unsaturated flow, assuming the fluid is incompressible and density- independent, the mass time-variation becomes:

�mc �t

� Vc �c �� c �t

(4)

If the ratio between capillary pressure pc [ML –1T–2] and specific weight of

liquid �w [ML –2T–2] represents the pressure head �c [L] and the specific retention

capacity Cc(�) [L –1] is considered (Bear, 1972), then applying the chain rule to

equation (4) the following equation is achieved:

�� c �t

� �� c �hc

�hc �t

� Cc �hc �t

(5)

from which the complete form of the unsaturated soil flux equation is obtained:

� � K

___ c�

h� hc l�

�

��

�� A� � VcCc

�hc �t

� Sc (6)

For the solution of unsaturated flow it is necessary to specify the non-linear dependencies among the assumed independent variable, total head hc, and terms

The discrete unsaturated flow equation (6) is applied to each cell of the domain, with hc the only unknown term to be estimated. The gravitational effects are taken into account considering that is h = ��+ z, and consequently � chh� � �cc zz � �� .

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

characterizing the hydraulic properties of soil represented by water content �c, specific retention capacity Cc and hydraulic conductivity Kc = K��c). Such rela- tionships can be expressed in tabular form or more usually through empirical equations fitting experimental data using theoretical models (Mualem, 1974, van Genuchten and Nielsen, 1985).

In the discrete formulation that we present the hydraulic head is assumed to be continuous and need not to be differentiable. Each cell may have different consti- tutive properties: this allows composite porous or fractured media to be dealt with. The boundary conditions are of two kinds: on some parts of the boundary the hydraulic head hc can be assigned, while on the remaining parts, the flow rate �

c� can be assigned. The sources can be continuous (drainage trench or superficial re- charge) and also may be concentrated (wells or sinks), by means of the mass source term Sc of equation (1). Our purpose is to find the hydraulic head of all cell nodes (barycenters) where it is not assigned: these can be internal as well as boundary nodes. Equation (6) is valid for interior or boundary cells: in this way, we can avoid the unnatural separation of the differential equations and the bound- ary conditions, which is typical of a differential formulation.

The discrete formulation of physical phenomena gives several computational advantages, such as discussed in Appendix A. This approach proves to be particu- larly suitable to the use of a macroscopic cellular automata environment and to be developed in a parallel computing system.

3. Discrete formulation of solute transport

The discrete form of the mass solute equation may be written as follows (Straface, 1998):

� s c �

c �I (h ) �

1 �s

c �ms

�t � 0 (7)

where � s c is the solute mass flux and �ms

c is the cell mass change. This equation

indicates that, the amount produced in the volume during a time interval is equal to the sum of the outgoing solute flow of the same quantity across the boundary of the volume during the time interval, and of the quantity stored inside the volume in the same interval. The equation is valid either for the internal nodes or for the lateral nodes. The mass time-variation becomes:

�ms

�t � n c � V

c � �C c

�t (8)

where Cc is the solute mass concentration in cell c.

S. STRAFACE AND G. MENDICINO

The mass transport phenomenon in an aquifer is subdivided into three separated mechanisms: (1) convection, (2) molecular diffusion and (3) cinematic dispersion.

Convection is the phenomenon where dissolved substances are carried along by the movement of fluid displacement. In this case we assume:

� s conv � � s

c � C c � 1 �s

c (9)

Molecular diffusion is a physical phenomenon linked to the molecular agita- tion. In a fluid at rest the Brownian motion projects particles in all directions of space. If the concentration of the solution is not uniform in space the point with the highest concentration sends out, on the average, more particles in all directions than the point with a lower concentration. The result of this molecular agitation is then that particles are transferred from zones of high concentration to those of low concentration (de Marsily, 1986).

Fick has found that the mass flux of particles in a fluid at rest is proportionate to the concentration gradient:

�s diff � A c �

1 �s

c � nc � d0 � g c (10)

The coefficient d 0 , knows as the molecular diffusion coefficient, is isotropic and can be expressed by the following equation:

d 0 � RT N

1 6��r

where R is the constant of perfect gases, N the Avogrado’s number, T absolute temperature, � fluid viscosity and r the mean radius of diffusing molecular aggre- gates.

In fine, kinematic dispersion is a mixing phenomenon linked mainly to the het- erogeneity of the microscopic velocities inside the porous medium on whatever scale they are observed:

� Inside a pore the velocities in mobile fraction are not uniformly distributed. � The differences of aperture and travel distance from one pore to another create

a difference in mean velocities. � A stratification or any features of large scale heterogeneity.

The suggested mathematical formula adopts a transport law through dispersion similar to the Fick’s law which accounts for the phenomena of mixing applied to the whole section of the medium, like the Darcy velocity, but with a dispersion coefficient D:

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

�sh disp � A c �

1 �s

c D � g c (11)

Here, in explicit form, the balance equation becomes:

1 �s

c �s c � C c n c � d 0 � D c� �� A c g c� ��

c �I (h ) � 1

�s c

�ms c

�t (12)

The function gc is the concentration gradient, whose formulation depends on type of the used interpolation. Employing a linear interpolation it can be expressed by:

C (x, y, z) � a � g x x � g y y � g z z

According to the analysis showed in Appendix A for the head gradient, we can write

gc � 1

3V c A cC c

f1f3

f2 x

e2 � �

If we use a cubic cell, the reference system XYZ is co-linear with cubic faces

and so the dispersion tensor D c assumes the form:

D c � DL 0 0 0 DT 0 0 0 DT

��

�

!� !� !�

"�

#�

$� $� $�

Because of more general flow directions, we have to express dispersion tensor on a new reference system X�Y�Z�, where X� is the flow direction, Y� and Z� his orthogonal. The change of the basis must be effected:

S. STRAFACE AND G. MENDICINO

D f c � T 1 D c

where f indicates the new reference basis, where tensor will assume the form:

D f c �

Dx' x' Dx' y' Dx' z' Dy' x' Dy' y' Dy' z' Dz' x' Dz' y' Dz' z'

��

�

!� !� !�

"�

#�

$� $� $�

4. Macroscopic cellular automata model

Cellular automata (CA) are a dynamic system where space, time and states are hypothesized as being discrete. They are based on a division of space in regular cells, each one having an identical computational device embedded in it: the finite automaton (fa). The physical quantities (or fa state) take only a finite set of values. The fa input is given by the states of the neighboring cells, including the cell which contains the fa.

The fa states vary according to a local rule (transition function); i.e., in a given time step a fa state depends on its state and on those of neighboring cells at the previous time step. Finite automata have identical transition functions, which are simultaneously applied to each cell (synchronous CA). At the beginning all the fa are in arbitrary states representing the initial conditions of the system. Then the CA evolve by changing the state of all fa simultaneously at discrete time steps according to the fa transition function. Finally, the global evolution of the CA sys- tem is derived from the evolution of all the cells.

Since every cell has the same neighborhood structure, even the cell at the boundary of a physical domain has neighboring cells that are outside the domain. Conventionally, border cells are assumed to be connected to the cells on the oppo- site boundary like neighbors forming a closed domain. For example, for a two- dimensional rectangular domain, a site on the left border has the site in the same row on the right border as its left neighbor. With the same update rule applied to all the cells, this yields what is called a periodic boundary condition which is rep- resentative of an infinite system. Certainly, the type of the boundary condition to be used in a simulation depends on the physical application under consideration. Other types of boundary conditions may be modeled using preset values of the cell for the boundary nodes or writing suitable update rules for the cells at the boundary.

The previous definition is not sufficient for modeling spatially extended natural macroscopic phenomena such as soil infiltration. More detailed conditions need to permit a correspondence between the system with its evolution in the physical space-time and the model with the simulations in the cellular space-time. Thus, for many cases the complexity of macroscopic natural phenomena requires an exten- sion of the original CA computational paradigm.

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Firstly, the dimension of the cell and the time correspondence to a CA step must be fixed. These are defined as global parameters, as their values are equal for all the cellular space. They represent the set P together with other global para- meters, which are frequently necessary for the simulation of complex phenomena involving time and/or space heterogeneity.

The state of the cell must account for all the characteristics (referring to the amount of space represented by the same cell), which are assumed to be relevant to the evolution of the system. Each characteristic corresponds to a sub-state, where the permitted values for a sub-state have to form a finite set. The set Q of possible states of a cell is given by the Cartesian product of the sets of the values of sub-states: Q = Q1 × Q2 × … × Qn. When a characteristic representing a physi- cal quantity is expressed in terms of a continuous variable referring to a space point, then the cell size chosen must be small enough so that the approximation of considering a single value for the whole cell extension may be adequate to the fea- tures of the phenomenon. Actually, the continuity of the variable is not a problem because a finite, but adequate, number of significant digits are utilized, so that the set of values allowed is extremely large but finite. Furthermore, the cell size must be large enough to describe a macroscopic approach, even though it must be capable of catching the smallest variations in the sub-states, so that they can be hypothe- sized as constant within each cell.

As well as the state of the cell that can be partitioned in sub-states, the transi- tion function % may also be split into p elementary processes, defined by the func- tions %1, %2, … %p. Elementary processes are applied sequentially according to a defined order compatible with the phenomenon under consideration. Different elementary processes may involve diverse neighborhoods, each one given by the union of all the neighborhoods associated to each process. If the neighborhood of an elementary process is limited to a single cell, such a process is defined as an internal transformation.

Specifically, according to the empirical approach suggested by Di Gregorio and Serra (1999), an elementary process is given by %�: Qn

m � Qc, where Qn and Qc are Cartesian products of the elements of sub-sets of Q, m is the number of cells of the neighborhood, involved in the elementary process, Qn describes the sub-states in the neighborhood that effect the change in the sub-state value and Qc identifies the cell sub-states that change their value.

5. Subsurface flow modeling through CA

Using the same extended notion for macroscopic CA, the three-dimensional model simulating water flux in unsaturated soils consists of a three-dimensional domain, regularly subdivided into cubic cells described by the following functional struc- ture (Mendicino et al., 2006):

S. STRAFACE AND G. MENDICINO

A � E d , X, Q, P, %� � (13)

where E d � x, y, z� � | x, y, z � N , 0 & x & 1x , 0 & y & 1y , 0 & z & 1z' ( is the set of cells identified by points with integer co-ordinates in the finite region, where the phenomenon evolves; N is the set of natural numbers; lx, ly and lz represent the limits of the region;

X � 0, 0, 0� �, 1, 0, 0� �, 0,1, 0� �, 0, 0, 1� �, 0, 0,1� �, 0, 1, 0� �, 1, 0, 0� � ' (

identifies the von Neumann neighborhood, which influences the change in the state of the central cell; Q is the finite set of the fa states, given by the Cartesian product of the following sub-states:

Q � Qh ) Q� ) Q� ) Qk ) Q�

where Qh is the sub-state describing the total head of the cell, Q� is the sub-state describing the pressure head, Q� is the sub-state describing the water content (therefore, it may indicate the value of moisture content in volume � or, according to the characteristic equation adopted, the water saturation Sw or the effective satu- ration Se), QK describes the hydraulic conductivity sub-state (if the medium is ani- sotropic, this sub-state is further subdivided according to the principal anisotropic directions) and, for transient condition analyses Q� indicates the sub-state cor- responding to a parameter value necessary to guarantee the convergence of the system;

P is the finite set of CA global parameters which affects the transition function and is made up of some parameters associated with the characteristic equations, the saturated hydraulic conductivity, the automaton dimension and, the time step. Specifically, we have the residual water content �r, the saturation water content �s, the capillary air entry pressure �a, the pore-dimension distribution index n, a con- tinuity parameter for pressure head �0, the specific storage Ss, the saturated hydraulic conductivity Ks, the cell dimension l� and the time step �t.

%: Qn 7 � Qc is the deterministic transition function. Once the initial conditions

(total head, pressure head, conductivity and water content values) and the bound- ary conditions are fixed, it is based on two elementary steps, %1 and %2:

1. The elementary process %1 consists in the update of the hydraulic characteris- tics of the soil (i.e., the hydraulic conductivity Kc, the water content �c and the specific retention capacity Cc), depending on the pressure head through the characteristic equations. Many theoretical models for the constitutive equations � �� and Kr = Kr(�) are available in literature. Those suggested by van Genuchten and Nielsen (1985) are commonly expressed as follows:

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

� �� r � � s � r� � 1� ��

m � * 0 � �� s � + 0

(14)

K r �� 1 � �� 5m / 2 1 � �� m � m� �

2 � * 0

K r �� 1 � + 0 (15)

where ��,�a�� n, m is a parameter given by m = 1 – 1/n, and Kr(�)[–] is the

relative hydraulic conductivity (so that K(�)=Ks Kr(�), with Ks [LT –1] saturated

conductivity); 2. The elementary process %2 is the application of the unsaturated soil flux equa-

tion (6), to update the values of the total head hc and the pressure head �c of the cell.

Starting the simulation, the sub-states condition depend on the initial values assigned to the total head hc, while the boundary conditions can be assigned either in terms of mass flow coming in (infiltration) or out (exfiltration) from the system (Neumann conditions), or fixing the total or pressure head values on some cells of the system (Dirichlet conditions).

6. Cellular automata quantization

Different conditions affect the nature of the CA. Among these the one based on quantized systems theory (Ziegler, 1998) was applied to the proposed model in order to reduce the state update transmission process. Such a process, called quan- tization, starting from the discretization of the states of a continuous process, fixes their evolution only through multiples of a specific value, called quantum size. Specifically, in a CA system the local interactions within the neighborhood, car- ried out through the transition function, involve each fa in temporal changes of the state: if during a time step the application of the transition function does not allow the cell to evolve from state D at least to state D - 1 (where the difference bet- ween two states is given by the quantum size), it maintains its current state, and does not exchange information with the neighborhood. The quantization makes cellular automata asynchronous, because for each iteration every fa decides, with respect to the value of its state and those of adjacent cells, whether to be updated or to remain frozen in its state during the previous time step.

For the proposed model the quantization procedure has been applied in a dif- ferent way to the original approach. In Ziegler (1998) the cell evolves through fixed values from state D to states D - 1, D - 2 , D - 3, etc. (like a step function). Instead we assume that the automaton state evolution, depending on a set of local interactions of the transition function, is only admitted if a fixed threshold is

S. STRAFACE AND G. MENDICINO

exceeded. The local interactions of the transition function are conditioned by the difference between the values of a chosen sub-state of the cell (the driving sub- state) before and after a single time step; if this difference is greater than a fixed threshold value then interactions are allowed and the automaton state can assume any numerical value in the field of real numbers. Obviously, when all local inter- actions are not significant (all the differences are smaller than the fixed threshold) then the automaton is assumed to be at rest.

In the unsaturated flow problems, the mass exchanges among the cells are due to the differences (in the space) among the total heads. Mass exchanges can be considered as not significant when they lead to a total head variation (in time) in the analyzed cell lower than a fixed threshold. Exceeding the threshold has to be checked for every neighbor with the aim of allowing the interaction along the generic ��direction. Therefore rewriting equation (6) it should be verified that (Mendicino et al., 2006):

�h�c � �t

VcCc Sc � K c�

h� hc l�

�

��

�� A�

��

� !�

"�

#� $� . quantumh (16)

where the term quantumh is the threshold, representing the minimum admissible total head variation in the time step �t. This quantumh determines a constant threshold within the CA domain.

On the CA domain a different threshold can also be hypothesized based on variations in time of water content, defined as quantum�. This quantum�, accord- ing to equation (16) has different effects on each cell of the automata, depending on the water content value of the same cell at the beginning of the time step and on the characteristic equation �� adopted.

Specifically, given the water content �c for the cell analyzed at the beginning of the time step, it is:

quantumh � � � c� � � � c quantum�� exfiltration quantumh � � � c � quantum�� c� � infiltration

(17)

In this way, the CA system becomes both asynchronous and non-uniform: the transition function changes in each fa, because in each automaton the threshold depends on a local factor represented by the degree of saturation at the beginning of the time step.

For both thresholds depending on total head and water content, increasing values produce two opposite effects: on the one hand the number of interactions decreases and on the other the model is less accurate. Then, a suitable compromise has to be reached, such as that shown in the next section.

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

7. CA high-performance environment

The three-dimensional unsaturated flow model has been developed using a high- performance environment called CAMELOT, specifically developed for macro- scopic cellular automata simulations (Dattilo and Spezzano, 2003; Folino et al., 2006). In CAMELOT each transition function uses the same local rule, even though it is possible to define some cells characterized by different transient func- tions (heterogeneous CA). In contrast to the classical cellular approaches which make reference to the basic model, where the state of each fa is defined as a single bit or a set of bits, in this case the fa state is defined like a set of typed sub-states. This increases the range of applications which can be simulated through cellular algorithms. Furthermore, in CAMELOT a logical neighborhood representing a wide range of different neighborhoods within the same radius, which can also be time dependent, has been introduced.

The CAMELOT components also include load balancing procedures based on a scatter-type decomposition technique which allows the computation to be equally distributed among processors of the parallel computing system. This is carried out using the standard Message Passing Interface (MPI) library, which allows CAMELOT to run on different hardware platforms for all communications. More specifically, the working environment is capable of determining the best number of nodes to be used in parallel grid computing for a given CA system, making the overheads due to the remote communications among the single processors very small. Such a process, which identifies the computational scalability of the CA system, becomes more efficient as the spatial dimension of the problem to be ana- lyzed increases. Therefore, the up-scaling from micro to meso-scale occurs with- out changing the local rule of the CA and, at the same time, the computational efficiency is improved by the working environment.

8. Results

The CA model was validated considering a wide range of benchmarks proposed in literature. Regarding the flow model, results were compared with numerical simu- lations of one-dimensional cases (benchmarks 1D1 and 1D2), a two-dimensional case (benchmarks 2D) and a three-dimensional case (benchmark 3D), while for the transport model a contaminant transport problem has been solved, starting from field data based on a real site where experimental measurements were available.

8.1. One-dimensional benchmarks

Two one-dimensional benchmarks were considered along two different soil col- umns (Paniconi et al., 1991): the former refers to an infiltration problem and the

S. STRAFACE AND G. MENDICINO

latter to an evaporation one. Both the cases are based on the unmodified Richards’ equation which does not allow general closed form exact solutions. Then, for the estimate of errors produced by the CA model some reference numerical solutions were considered, assuming a very dense grid and small time step (Mendicino et al., 2006).

The characteristic equations taken into account for these simulations are the (14) and (15), with a modified retention curve �� in a form like this (Paniconi et al., 1991):

� �� r � � s � r� � 1� ��

m � * �0 � �� r � � s � r� � 1� �0� �

m � Ss � �0� � � + �0 (18)

where �0 ��(�0) = (��0/�a) n, �0 being a continuity parameter and Ss the specific

storage. The general storage term of the Richards’ equation in this case is given by /(�)=d�/d�. The variation in the original equations is necessary to avoid numeri- cal problems, such as those indicated by Paniconi et al. (1991). The accuracy of the CA model was estimated by using both a first and a second order norm error:

01(t) � ˆ ��(z i , t) �ref (zi , t)

ni�1

n � (19)

02(t) � ˆ ��(zi , t) �ref (z i , t)� �

ni�1

n � (20)

where �̂ and �ref are the simulated and reference pressure head at time t respec- tively, for the i=1..n soil profile points at level zi.

Benchmark 1D1 consists of an infiltration and redistribution simulation into a soil column initially at hydrostatic equilibrium. The boundary condition at the sur- face is a time-varying specified Darcy flux q which increases linearly with time, while the boundary at the base is maintained at a fixed pressure head value of��= 0, allowing drainage of moisture through the water table. The space grid and time step have been chosen to guarantee the convergence of the system, following a cri- terion described in Appendix B (the same criterion has been applied for all simula- tions considered).

Figure 2 shows the comparison between CA simulations and reference numeri- cal solutions: the differences are very small for all the times analyzed.

Benchmark 1D2 simulates evaporation from an initially wet soil with a fixed water table boundary condition at the base of the column. The boundary condition at the surface is a specified and constant Darcy flux, until the pressure head at the surface reaches its air dry value �min, after which the surface becomes a fixed head boundary. The comparison between simulations and numerical reference solutions is shown in Figure 3.

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Pressure head(m)

E le va tio n (m )

-10 -8 -6 -4 -2 0

0 2

4 6

8 10

32Time0 1 2 4 10

Fig. 2. Comparison between CA simulation (solid lines) and reference numerical solutions (points) for the bench-mark 1D1.

Pressurehead(m)

E le va tio n (m )

-2 -1.5 -1 -0.5 0

0 1

2 3

4 5

Time45 024 12 3 1

Fig. 3. Comparison between CA simulations (solid lines) and reference numerical solutions (points) for bench-mark 1D2.

S. STRAFACE AND G. MENDICINO

The CA model shows a good agreement with the reference solutions, although small differences in pressure head at the soil surface were observed. This is mainly due to the constitutive equations in the conversion of surface boundary condition, which is expressed in terms of water content � in the reference solution and in terms of total head h in the CA model.

8.2. Two-dimensional benchmarks

The 2D benchmark is partly based on textural, structural and hydraulic properties of Jornada Test Site soil at Las Cruces (New Mexico) which are described by

in a very-arid heterogeneous soil. The test area consists of three soil layers with different characteristic proper-

ties. Within the bottom-layer, there is a very permeable small zone. The domain extends about 8.00 m horizontally and 6.50 m vertically. The experiment was car- ried out again numerically by Magnuson et al. (1990), utilizing FLASH (FEM) and PORFLOW (FDM) calculus codes, therefore the results (supplied by such codes) were assumed as standard numerical solutions for the comparison with the proposed AC model. Regarding the initial conditions, the soil piezometric head is uniform and equal to –7.34 m, corresponding to a saturation degree variable from 0.31 into the top-layer to 0.37 into the bottom-one. The boundary conditions are given by an infiltration zone of 0.02 m d–1 which extends on the area of 2.25 m starting from the left side of the domain. The remaining boundary is characterized by no flow conditions. Finally, soil properties are based on van Genuchten (1978) relationship according to the Mualem (1974) model.

The AC model which simulates soil moisture evolution consists of square cells with a 0.05 m side. The simulation time-step goes from a maximum value of 100 s to a minimum ones of about 28 s, according to the limits established by the fol- lowing criterion of convergence (Mendicino et al., 2006):

�t & l 2Cc 4 Kc

(21)

As shown in Figure 4 the saturation degree distribution along the analyzed plot is obtained through AC model after 30 days from the starting of the simulation, and compared with the results obtained by means of FLASH and PORFLOW codes (Figure 5). Thirty days step is enough to reach the conditions comparable with the steady ones. The effect of the high-saturation zone is clearly visible, because the curve indicating the saturation degree is strongly influenced by it. In particular, the AC model results shown a good agreement with PORFLOW ones, while the output of the FLASH code shows some differences (a greater horizontal extension of the curve at 0.4 saturation degree, and a less face along the vertical).

Smyth et al. (1989). The test considers the transient phase of a 2D infiltration flow

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Fig. 4. Saturation degree change respect to the initial condition after 30 days. The saturation degree shown is greater than 0.0001.

Fig. 5. Comparison between CA simulations and reference numerical solutions (PORFLOW and FLASH) for benchmark 2D after 30 days from the simulation start.

S. STRAFACE AND G. MENDICINO

8.3. Three-dimensional benchmark

The benchmark 3D is taken from Huyakorn et al. (1986) and analyses a transient flow in a rectangular soil column of dimensions 0.50 × 0.50 × 2.00 m subjected to infiltration and subsequent evaporation. The bottom and the top faces of this column correspond to the water table and the soil surface respectively. The initial pressure

and then subjected to evaporation for another 10 days; for both phenomena the maximum flux rate was assumed as equal to 0.05 m d–1. The minimum allowable pressure head at the top of soil column is –0.90 m, saturated hydraulic conductivity Ks is equal to 0.1 m d

–1, porosity 1 is 0.45, residual water saturation Swr is 0.333 and air entry value �a is equal to 0.0 m. The constitutive equations are given by the equation:

� �� 100 �� 1 Sw� � 1 Swr� � (22)

K r � Sw Swr� � 1 Swr� � (23)

Simulations did not show particular convergence or stability problems because of the linearity of equations (22) and (23). The time step used for simulations was fixed as equal to 60 s while the automaton dimension was assumed as equal to 0.05 m. The results obtained with the CA model are compared with reference solutions in Figures 6a (infiltration) and 6b (evaporation).

Fig. 6. Comparison between CA simulations (solid lines) and reference numerical solutions (points) for benchmark 3D, for infiltration (a) and evaporation (b).

– 0.97 m elsewhere. The soil column was first subjected to infiltration for 10 days head was assumed to be zero at the water table, –0.90 m at the soil surface, and

Pressure head (m) Pressure head (m)

E le

v a

ti o

n (m

)

-1 -0.8 -0.6 -0.4 -0.2 0 -1 -0.8 -0.6 -0.4 -0.2 0

0 0 .5

1 1 .5

E le

v a

ti o

n (m

)

0 0 .5

1 1 .5

Time 0d 0.1 1 2.6 10

Time 20d

Time 10d

(a) (b)

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

8.4. Transport benchmark

In this section we report an example of a validation process using a two dimen- sional contaminant transport problem (Troisi et al., 2000). The two PS sets have been defined starting from field data based on a real site where experimental measurements were available. The study area (Figure 7) is located near the town of Montalto Uffugo in Calabria, Southern Italy. It is a valley at the confluence of the Settimo River on the south, the Mavigliano River on the north and the Crati River on the east. The geology is formed by three layers: sand (0–7 m), clay (7–11 m), and silt (11 to about 40 m) (Troisi, 1995). A basal clay underlies the silty sand. A local perched water table is in the alluvium above the clay layer, and the silty sand layer constitutes a confined aquifer above the basal clay. Measurements of the hydraulic conductivity and its distribution have been made over the past few

Fig. 7. Location of the study area used for the transport benchmark. The boundary of the numerical model is shown as a dotted line.

S e rra

S e

tim o G

ra n

d e

S . A

n to

n e llo

F o

n t.n

a S

u g

h e ro

S . B

ia g

7 7

5 1

0 0

2 3

0 8

0 0

1 8

632 1270 420

M A

V IG

L IA

N O

R IV

E R

HIGHWAY A3

S E TTIM

O R

IV E R

CRATI RIVER

S. STRAFACE AND G. MENDICINO

years using a kriging approach with external drift methodology applied to electrical-

application is defined as the area enclosed by the surface water system as shown in Figure 7. In the flow model a heterogeneous isotropic porous medium is assumed. Dirichlet boundary conditions are according to the water levels in the surface water system. Two permeability zones are assumed in the flow fields, and are indicated in Figure 7. Zone 1 is characterized by K11 = K22 = 10

–5, Ss = 10 –2m while in zone

2, K11 = K22 = 10 –6, Ss = 10

–3m. One extracting well, located at coordinates X = 1640 and Y = 2525 and pumping at a rate of qout = 0.003 m, is coupled with an injecting well, located at X = 114 and Y = 225, where the head is kept constant at 138 m. In the transport test case we simulate a contaminant injection from the Mavigliano River, just east of Highway A3, where a contaminant source of 100 m of length is kept at constant concentration. A second contaminant point source is located at X = 146 and Y = 241. The parameters used in the simulations are: n = 0.2, �L � 5, and �T � 1.

The benchmark code for this test case has been run on a number of successively refined computational grids so that a reliable solution is obtained. The results for the test case at time t = 2.3 days are shown in Figure 8. The solutions obtained with the TRAN2D (Gambolati et al., 1993) and CA models are not significantly different so that the CA model can be considered verified for this benchmark. How- ever, more meaningful and objective comparison procedures need to be devised, as the contour plots give only a qualitative picture. Appropriate simulation quanti- ties useful for these purposes could be, for example, mass balance measures (both local and global), stability measures (e.g., negative concentrations in transport pro- blems, existence of negative transmissivities in stiffness matrices), accuracy measures (e.g., outgoing fluxes from zero Neumann boundaries), and so on.

resistivity data (Troisi et al., 2000). The domain for the flow and transport model

Fig. 8. Results of the transport problem at t = 2.3 days: benchmark solution (left), CA model solution (right).

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

9. Cellular automata quantization effects

The quantization effects considering both total head and water content thresholds were investigated using the benchmark 1D1 in a three-dimensional configuration. Figures 9a, 9b, 10a, 10b, 11a and 11b show the reduction in the number of auto- mata interactions along a column of cells together with the corresponding norm error 01 values (assuming as reference values �ref those obtained from the simulations without quantization), when the values of thresholds depending on total head (quantumh) and water content (quantum�) are changed, respectively at times of 4, 10 and 32 h. Obviously, messages exchanged for quantumh = quantum� = 0 are the same. For the particular structure of the 3D problem, it is clear that considering the whole automata instead of only a column leads to a number of exchanges equal to the values shown in the figures multiplied by the number of automata columns, not varying the 01 error values.

Fig. 9. Reduction in exchanged messages in a column of cells and variation of error norm 01 values at a time of 4 hours, varying (a) quantumh and (b) quantum� threshold values. Points drawn with two concentric circles represent the number of exchanged messages with quantumh = quantum� = 0.

Fig. 10. Reduction in exchanged messages in a column of cells and variation of error norm 01 values at a time of 10 hours, varying (a) quantumh and (b) quantum� threshold values. Points drawn with two concentric circles represent the number of exchanged messages with quantumh = quantum� = 0.

N u

m b

e r

o f

m e

s s a

g e s

(m il li o

n s )

10-7 10-6 10-5 10-4 10-3 10-2 10-10 10-9 10-8 10-7 10-6 10-5

0 3

6 9

1 2

1 5

N u

m b

e r

o f

m e

s s a

g e s

(m il li o

n s )

0 3

6 9

1 2

1 5

0 0

.0 0 5

0 .0

0 0

.0 0 5

0 .0

Number of messages

(m )

quantum h =0

quantum h (m) quantum q (-)

Number of messages

quantum = 0

(a) (b)

(a)

quantum h (m)

N u

m b

e r

o f

m e

s s

a g

e s

( m

il li

o n

s )

(m )

10-7 10-6 10-5 10-4 10-3 10-2

(a)

quantumq (-)

10-7 10-6 10-5 10-4 10-3 10-2

0 1

0 2

0 3

0 4

N u

m b

e r

o f

m e

s s

a g

e s

( m

il li

o n

s )

0 1

0 2

0 3

0 4

0 0

.0 0

5 0

.0 1

(m )

0 0

.0 0

5 0

.0 1

Number of messages

quantum h

= 0 quantum = 0

Number of messages

S. STRAFACE AND G. MENDICINO

Fig. 11. Reduction in exchanged messages in a column of cells and variation of error norm 01 values at a time of 32 hours, varying (a) quantumh and (b) quantum� threshold values. Points drawn with two concentric circles represent the number of exchanged messages with quantumh = quantum� = 0.

For all the times analyzed, starting from a very small threshold value (quantumh equal to 10–7 m and quantum� equal to 10

–10), an appreciable reduction in the number of automata interactions is obtained. This however appears less significant for increasing simulation times.

Once a certain threshold value is exceeded, for both kinds of thresholds the increase in the norm error is not tolerable, even if the reduction in exchanges still continues. Considering quantumh the critical threshold value seems to be between 10–5 and 10–4 m.

For quantumh equal to 10 –5 m reductions in exchanged messages of 82.3%

(with 01 = 9.1 × 10–4 m), 42.7% (with 01 = 2.3 × 10–4 m) and 13.3% (with 01 = 1.0 × 10–3 m) were obtained at times of 4, 10 and 32 hours respectively. Instead, for quantumh equal to 10

–4 m reductions in exchanged messages of 87.3% (with 01 = 9.1 × 10–3 m), 51.8% (with 01 = 5.2 × 10–4 m) and 16.1% (with 01 = 1.0 × 10–3 m) were obtained at the same times respectively. For increasing simulation times it was observed that higher threshold values still provided acceptable results.

No essential differences were observed using quantum� threshold, even though for high values some numerical divergence problems occurred, especially for increasing simulation times.

For the same infiltration case 1D1, characterized by retention curves (18), the relationship between quantumh and quantum� is given by the following equation (Mendicino et al., 2006):

quantumh � � � � quantum��

�a � � quantum� � r

� s � r

�

��

�� 1 m

1 ��

� !� !�

"�

#� $� $�

1 n

�a � � r � s � r

�

��

�� 1 m

1 ��

� !� !�

"�

#� $� $�

1 n

(24)

(a) (b)

N u

m b

e r

o f

m e s s a

g e

s (m

il li o

n s )

0 2 0

4 0

6 0

8 0

1 0 0

1 2 0

N u

m b

e r

o f

m e s s a

g e

s (m

il li o

n s )

0 2 0

4 0

6 0

8 0

1 0 0

1 2 0

0 0

.0 0 5

0 .0

Number of messages

(m )

0 0

.0 0 5

0 .0

(m )

quantum h = 0

quantum h (m)

Number of messages

quantum = 0

10-7 10-6 10-5 10-4 10-3 10-2

quantumq (m) 10-10 10-9 10-8 10-7 10-6 10-5

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Critical quantum�� threshold values seem to vary between 10 –7 and 10–6. For

quantum� equal to 10 –7 reductions in exchanged messages of 82.0% (with 01 = 4.8 ×

10–4 m), 39.8% (with 01 = 7.0 × 10–6 m) and 12.4% (with 01 = 1.0 × 10–6 m) were obtained at times of 4, 10 and 32 h respectively. Instead, for quantum� equal to 10–6 m reductions in exchanged messages of 87.0% (with 01 = 5.6 × 10–3 m), 47.9% (with 01 = 1.5 × 10–3 m) and 15.2% (with 01 = 1.2 × 10–3 m) were obtained at the same times respectively.

The reduction in exchanges varying the thresholds as well as specific effects produced by static (quantumh) and dynamic (quantum�) approaches should not be thought of as general results, but are associated with the problem analyzed. In fact, considering benchmark 1D1 we start from a top-perturbed initial condition and, consequently, the gradients between cells are initially very low and so even a small threshold significantly affects the behavior of the system.

10. Conclusions

Three-dimensional unsaturated flow and transport modeling has been developed by means of an extended notion for macroscopic cellular automata where local laws governing the automata interaction are based on physically-based rules. The modeling uses a CA structure based on functionalities capable of increasing its computational capacity, both in terms of working environment and in terms of the optimal number of processors available for parallel computing.

In terms of performance, the model has been verified considering multidimen- sional reference solutions taken from benchmarks in literature. The simulations carried out on one-, two- and three-dimensional benchmarks have shown results similar to the reference solutions, with first and second order norm error values never greater, respectively, than 0.028 m and 0.5 m2. The AC characteristic equations are subjected to parameters which amplify the non-linearity of the phenomenon. Even if not explicitly shown in this chapter, the effects of non-linearity can be reduced using smaller dimension automata, which evolve with lower mass flow values, which better approximate the non-linearity condition.

Furthermore, in the three-dimensional configuration of the benchmark 1D1 the benefits of the use of quantization techniques to reduce the number of interactions occurring locally among adjacent automata have been observed. The CA structure working in an asynchronous and non-uniform manner, allowed a considerable reduction in the local transitions, which becomes very significant in the first steps of simulations characterized by scarce mass exchanges. In this case, in the first four hours of simulation a reduction greater than 80% in the messages exchanged among automata was achieved with negligible error values by considering both thresholds depending on total head and water content.

Therefore, at the end of this analysis, the modeling proposed appears to be: (1) innovative with respect to other macroscopic CA models, because its physically- based structure differs from the others commonly based on empirical methods

S. STRAFACE AND G. MENDICINO

where local laws show some parameters whose physical meaning is not clear and, as a consequence, heavy calibration phases are necessary to estimate suitable values of the same parameters; (2) innovative with respect to other unsaturated flow models, because it solves the problem through a direct discrete approach; (3) reli- able for all the simulations carried out with one-, two- and three-dimensional benchmarks; (4) faster than some other numerical approaches not based on quanti- zation techniques aimed at reducing the local interactions among the elements of the discrete system.

The proposed discrete approach, coupled to the CA environment implementa- tion, allows us to deal with heterogeneous and anisotropic soils and composite porous or fractured media, and with different (continuous or concentrated) sources in a quite simple way, by just varying the parameters of specific regions of the CA and not changing the transition function.

Finally, from a computational point of view the higher efficiency values shown by the model running on parallel machines increase the capability of the CA sys- tem also with respect to the fully-coupled modeling of complex macroscopic phenomena such as those represented in this chapter.

Appendix A. Direct discrete formulation of the Darcy equation

The aim of this appendix is to show that the direct discrete formulation of the unsaturated flow equation on three-dimensional cubic cells is a particular form of a more general direct discrete formulation of the Darcy equation.

Figure A1a shows the tetrahedral element c with a volume Vc whose vertices correspond to the barycenters of four tetrahedral cells belonging to a more com- plex cell system. Referring to this figure, we hypothesize that the variable H (total head) is given by the following equation:

Fig. A1: (a) The elementary tetrahedron of a generic cell system; (b) The elementary cube of a three-dimensional cubic cell system.

Lj Vc

Ahjk

Ahij

Ahik

Lhi

Lhk

Lhj

a) b)

Lj Vc

Ahjk

Ahij

Ahik

Lhi

Lhk

Lhj

a) b)

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

H (x, y, z) � a � gx x � gy y � g z z (A.1)

where the coefficients gx, gy and gz are the hydraulic gradient vector components, and a is a constant. Given h, i, j and k the tetrahedron vertices, with Li, Lj and Lk the three sides starting from the vertex h, we can write:

H i H h � gx xi xh� �� gy yi y h� �� gz zi zh� �

H j H h � gx x j xh� �� gy y j y h� �� gz z j zh� � H k H h � gx xk xh� �� gy yk y h� �� gz zk zh� �

2�

3� 4�

5� 4�

(A.2)

or in matrix form:

Lix Liy Liz L jx L jy L jz Lkx Lky Lkz

��

�

!� !� !�

"�

#�

$� $� $�

gx gy gz

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

� H i H h H j H h H k H h

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

(A.3)

where the matrix of L terms contains the projections of the sides Li, Lj and Lk on the x, y, z axes. The solution of the system, obtained by Cramer’s rule, is the fol- lowing:

gx gy gz

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

� 1

3Vc

Ahx Aix A jx Akx Ahy Aiy A jy Aky Ahz Aiz A jz Akz

��

�

!� !� !�

"�

#�

$� $� $�

H h H i H j H k

2�

3� 4� 4�

5� 4� 4�

6�

7� 4� 4�

8� 4� 4�

(A.4)

where the matrix of A terms contains the projections of the oriented areas Ah, Ai, Aj and Ak of the tetrahedron on the x, y, z axes. System (A.4) can also be written in the following way:

g' (c � 1

3Vc A� �c H' (c (A.5)

The Darcy equation links the mass flux � to the hydraulic gradient g through the following system:

� x (x, y, z) � y (x, y, z) � z (x, y, z)

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

� kxx kxy kxz kyx kyy kyz kzx kzy kzz

��

�

!� !� !�

"�

#�

$� $� $�

gx (x, y, z) gy (x, y, z) gz (x, y, z)

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

(A.6)

S. STRAFACE AND G. MENDICINO

that is:

�' (c � k� �c g' (c (A.7)

If a cubic cell is considered (Figure A1b), the off-diagonal terms of the system (A.2) are null. In matrix form results:

Lhi 0 0 0 Lhj 0 0 0 Lhk

��

�

!� !� !�

"�

#�

$� $� $�

gx gy gz

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

� H i H h H j H h H k H h

2�

3� 4�

5� 4�

6�

7� 4�

8� 4�

(A.8)

and the solution of the system is given by:

gx � H i H h� � Lhi gy � H j H h� � Lhj g z � H k H h� � Lhk

2�

3� 4�

5� 4�

(A.9)

then, assuming that the xyz Cartesian system is colinear with the principal aniso- tropy directions, the Darcy equation becomes:

� x � kx H i H h� � Lhi � y � ky H j H h� � Lhj � z � kz H k H h� � Lhk

2�

3� 4�

5� 4�

(A.10)

The equations of the system (A.10) are equivalent to equation (2), showing that the direct discrete formulation of the unsaturated flow equation on three-dimensional cubic cells is a particular form of a more general direct discrete formulation of the Darcy equation.

The three-dimensional cubic cell system is a particular case of a Delaunay tessellation. With this kind of tessellation the discrete governing equation system is similar to the one achieved using Finite Difference or Finite Volume Method schemes (Mattiussi, 1997, Manzini and Ferraris, 2004). If we do not use a Delaunay tessellation (as in the case of an irregular mesh), we obtain a discrete governing equation system by means of an interpolation of the hydraulic head on the cells. It has been shown that, for linear interpolation the discrete governing equation sys- tem coincides with the one of the Finite Element Method (Tonti, 2001). However, for quadratic interpolation the discrete equation system, which is asymmetric, dif- fers from the FEM scheme, achieving a convergence of the fourth order, greater than the one obtained with FEM using the same interpolation (Tonti, 2001).

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

When we use a method with a differential or integral formulation, the choice of the cell or the tessellation type depends on the method selected (FEM, FDM, FVM…). In contrast, if we use a direct discrete formulation, the cell or the tessel- lation type and the time interval can be chosen considering the physical laws gov- erning the problem, the spatial and temporal scale of the phenomenon and, if the cells are considered as finite automata of a CA, the macroscopic CA environment can be chosen, with the opportunity, for example, of using different shapes (e.g., squares and hexagons) for the single elements.

Appendix B. Convergence of discrete unsaturated flow equation

The problem of convergence, in a finite difference method, consists of finding the conditions under which the difference U(i,j,k,t) u(i,j,k,t) between the theoretical solutions of the differential (U) and difference (u) equations at a fixed point (i,j,k,t) tends to zero uniformly, as the net is refined in such a way that �x,�y,�z,�t � 0, and m1,m2,m39n � �, with m1�x ( = i), m2�y ( = j), m3�z ( = k) and n�t ( = t) remains fixed, and m1,m2,m39n being integers, with m1 = m2 = m3��n = 0 the origin.

The fixed point (i,j,k,t) is anywhere within the region under consideration, and it is sometimes necessary in the convergence analysis to assume that �x,�y,�z,�t do not tend to zero independently but according to some relationship (Mitchell and Griffiths, 1980).

If we consider a discrete cell system directly, of course we can assume that a difference exists between the ‘exact’ infinitesimal (differential) solution and the discrete solution. Since the independent variable in the unsaturated soil flux equa- tion (6) is the total head h, we can introduce the error e:

i, j,k te � i, j ,k

tH i, j,k th (B.1)

that is the difference between the exact solution H and the discrete solution h at the grid point i,j,k,t.

Rearranging equation (6) considering �x = �y = �z = l, we obtain:

�hc �t

� 1

Cc l 3 Sc �

1 Cc l

2 Kc� _____

h� hc� � � � (B.2)

Expanding the terms and considering the introduced error term, we have with respect to the grid point i,j,k,t (i.e., the cell c becomes the grid point i,j,k,t):

S. STRAFACE AND G. MENDICINO

H i, j,k t ��t ei, j,k

t ��t H i, j,k t � ei, j,k

�t �

1 Ci, j,k

t l 3 Sc �

1 Ci, j,k

t l 2 [ Ki 1, j,k

t H i 1, j,k t ei 1, j,k

t� �� Ki�1, j,kt H i�1, j,kt ei�1, j,kt� �� Ki, j 1,k

t H i, j 1,k t ei, j 1,k

t� �� Ki, j �1,kt H i, j �1,kt ei, j �1,kt� �� Ki, j,k 1

t H i, j,k 1 t ei 1, j,k

t� �� Ki, j,k �1t H i, j,k �1t ei, j,k �1t� � Ki-1, j -1,k -1t H i, j,kt ei, j,kt� �]

(B.3)

where Ki-1, j -1,k -1 t � Ki 1, j,k

t � Ki�1, j,k t � Ki, j 1,k

t � Ki, j �1,k t � Ki, j,k 1

t � Ki, j,k �1 t is the

sum of the hydraulic conductivities averaged between the grid point i,j,k,t and the adjacent points i ± 1,j ± 1,k ± 1,t, obtained considering in equation (3) the ele- ments on the diagonal of the hydraulic conductivity tensor corresponding to the linking directions.

Isolating the error term in the time step t + �t, we get

ei, j,k t ��t � 1 Ki-1, j -1,k -1

t �t l 2Ci, j,k

��

� !� !�

"�

#� $� $� ei, j,k

t �

� �t

l 2Ci, j,k t (Ki 1, j,k

t ei 1, j,k t � Ki�1, j,k

t ei�1, j,k t � Ki, j 1,k

t ei, j 1,k t � Ki, j �1,k

t ei, j �1,k t �

�Ki, j,k 1 t ei, j,k 1

t � Ki, j,k �1 t ei, j,k �1

t ) �t

Ci, j,k t l 3

Sc � f (H , K)

(B.4)

where f (H,K) is a function of permeability K and of the exact solution H:

f (H , K ) �

H i, j,k t ��t 1 K i-1, j -1,k -1

t �t l 2C i, j,k

��

� !� !�

"�

#� $� $� H i, j,k

t �

�t

l 2C i, j,k t (K i 1, j,k

t H i 1, j,k t � K i�1, j,k

t H i�1, j,k t � K i, j 1,k

t H i, j 1,k t

�K i, j �1,k t H i, j �1,k

t � K i, j,k 1 t H i, j,k 1

t � K i, j,k �1 t H i, j,k �1

t )

(B.5)

From now on, for the sake of simplicity, we will not consider the source term Sc. Let , ,

t i j kE denote the maximum value of ei, j,kt at time t and M the maximum

modulus of f (H,K) for all i,j,k,t. When the term in the square brackets in equation (B.4) is equal or greater than zero, all the coefficients of e in the equation are posi- tive or zero. The term in the square brackets is equal to or greater than zero when:

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

�t & l 2Ci, j,k

Ki-1, j -1,k -1 t

(B.6)

and, working on equation (B.4), we have:

ei, j,k t � �t & Ei, j,k

t � M

As this is true for all values of i,j,k it is true for max i, j,k

ei, j,k t ��t & Ei, j,k

t ��t � M . Hence:

Ei, j,k t ��t & Ei, j,k

t � M & Ei, j,k t �t � M� �� M � Ei, j,kt �t � 2M

etc., from which it follows that:

d� dH

:H :t

:Kx ��

:x :H :x

Kx �� :2H :x 2

:Ky ��

:y :H :y

Ky �� :2H :y 2

:Kz ��

:z :H :z

Kz �� :2H :z 2

�

��

� � � �

��

i, j,k

As H i, j,k t hi, j,k

t & Ei, j,k t , this proves that h converges to the exact solution H

and that H is the solution of the Richards’ equation, as �x,�y,�z � 0, when the condition expressed in equation (B.6) is respected and t is finite.

References

Avolio, M.V., Crisci, G.M., D’Ambrosio, D., Di Gregorio, S., Iovine, G., Rongo, R., Spataro, W., An extended notion of Cellular Automata for surface flows modelling, WSEAS Transactions on Computers, 2, 1080–1085, 2003.

Bear, J., Dynamics of Fluid in Porous Media, American Elsevier, New York, 1972. Chau, J.F., Or, D., Sukop, M.C., Simulation of gaseous diffusion in partially saturated

porous media under variable gravity with lattice Boltzmann methods, Water Resources Research, 41, W08410, 2005.

Ei, j,k & Ei, j,k � n � M � n � M

because the initial values for h and H are the same, i.e., 0, , 0i j kE � . When �x,�y,�z� 0, l � 0 and, for the equation (B.6), also �t � 0, then the retention capacity tends to the general storage term of the Richards’ equation ( C d dH�; ) and M tends to be a solution of the Richards’ equation in a form like this:

S. STRAFACE AND G. MENDICINO

Chen, H., Chen, S., Matthaeus, W.H., Recovery of Navier-Stokes equations using lattice- gas Boltzmann method, Physical Review A, 45, 5339–5342, 1992.

Chopard, B., Luthi, P.O., Lattice Boltzmann computations and application to physics, Theoretical Computer Science, 217, 115–130, 1999.

Crutchfield, J.P., Mitchell, M., Das, R., The evolutionary design of collective computation in cellular automata, in Crutchfeld, J.P. and Schuster, P.K., editors, Evolutionary Dynamics-Exploring the Interplay of Selection, Neutrality, Accident, and Function, Oxford University Press, New York, 2002.

Dattilo, G., Spezzano, G., Simulation of a cellular landslide model with CAMELOT on high performance computers, Parallel Computing, North Holland, 29(10), 1403– 1418, 2003.

de Marsily, G., Quantitative Hydrogeology: Groundwater Hydrology for Engineers, Academic, San Diego, USA, 1986.

Deng, J.Q., Ghidaoui, M.S., Gray, W.G., Xu, K., A Boltzmann-based mesoscopic model for contaminant transport in flow systems, Advances in Water Resources , 24, 531– 550, 2001.

Di Gregorio, S., Serra, R., An empirical method for modelling and simulating some com- plex macroscopic phenomena by cellular automata, Future Generation Computer Sys- tem, 16, 259–271, 1999.

Di Gregorio, S., Serra, R., Villani, M., Applying cellular automata to complex environ-

Di Pietro, L.B., Melayah, A., Zaleski, S., Modeling water infiltration in unsaturated porous media by interacting lattice gas-cellular automata, Water Resources Research, 30(10), 2785–2792, 1994.

Folino, G., Mendicino, G., Senatore, A., Spezzano, G., Straface, S., A model based on cel- lular automata for the parallel simulation of 3D unsaturated flow, Parallel Computing, 32, 357–376, 2006.

Frisch, U., d’Humieres, D., Hasslacher, B., Lallemand, P., Pomeau, Y., Rivet, J.P., Lattice gas hydrodynamics in two and three dimensions, Complex Systems, 1, 649–707, 1987.

Gambolati, G., Paniconi, C., Putti, M., Petruzzelli, D., Helfferich, F.G., Numerical model- ing of contaminant transport in groundwater, Migration and Fate of Pollutants in Soils and Subsoils, NATO ASI Series G: Ecological Sciences, 32, 381–410, 1993.

Ginzburg, I., Equilibrium-type and link-type lattice Boltzmann models for generic advec-

1195, 2005. He, X., Luo, L.S., Lattice Boltzmann model for the incompressible Navier-Stokes equation,

Journal of Statistical Physics, 88, 927–945, 1997. Higuera, F.J., Jimenez, J., Boltzmann approach to lattice gas simulations, Europhysics Letters,

9, 663–668, 1989. Huyakorn, P.S., Springer, E.P., Guvanasen, V., Wadsworth, T.D., A three-dimensional

finite element model for simulating water flow in variably saturated porous media, Water Resources Research, 22(12), 1790–1808, 1986.

Indelman, P., Dagan, G., Upscaling of permeability of anisotropic heterogeneous forma- tions, 1. The general framework, Water Resources Research, 29(4), 917–923, 1993.

Knutson, C.E., Werth, C.J., Valocchi, A.J., Pore-scale modeling of dissolution from varia- bly distributed nonaqueos phase liquid blobs, Water Resources Research, 37, 2951– 2963, 2001.

mental problems: the simulation of the bioremediation of contaminated soils, Theo- retical Computer Science, 217, 131–156, 1999.

tion and anisotropic-dispersion equation, Advances in Water Resources, 28, 1171–

CELLULAR AUTOMATA MODELING OF ENVIRONMENTAL SYSTEMS

Manzini, G., Ferraris, E., Mass-conservative finite volume methods on 2-D unstructured grids for the Richards equation, Advances in Water Resources, 27, 1199–1215, 2004.

Martys, N.S., Chen, H., Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method, Physical Review E, 53, 743–750, 1996.

Mattiussi, C., An analysis of Finite Volume, Finite Element, and Finite Difference Methods using some concepts for algebraic topology, Journal of Computational Physics, 133, 289–309, 1997.

McNamara, G.R., Zanetti, G., Use of the Boltzmann equation to simulate lattice-gas auto- mata, Physical Review Letters, 61, 2332–2335, 1988.

Mendicino, G., Senatore, A., Spezzano, G., Straface, S., Three-dimensional unsaturated flow modeling using cellular automata. Water Resources Research, 2006, 42, W1141, W1– W18.

Mendicino, G., Senatore, A., Spezzano, G., Straface, S., Automi cellulari macroscopici per la simulazione del moto tridimensionale in mezzi porosi non saturi. Proceedings of “XXX Convegno Idraulica e Costruzioni Idrauliche”, Rome, 2006.

Mitchell, A.R., Griffiths, D.F., The Finite Difference Method in Partial Differential Equa- tion, J. Wiley & Sons, Chichester, 1980.

Mualem, Y., A conceptual model of hysteresis, Water Resources Research, 10, 514–520, 1974.

Muzy, A., Wainer, G., Innocenti, E., Aiello A., Santucci J.F., Cell-DEVS Quantization techniques in a fire spreading application, Proceedings of the Winter Simulation Con- ference 2002 – Exploring new frontiers, San Diego, USA, 2002.

Orlandini, S., Two-layer model of near-surface soil drying for time-continuous hydrologic simulations, Journal of Hydrologic Engineering , 4(2), 91–99, 1999.

Paniconi, C., Aldama, A.A., Wood, E.F., Numerical evaluation of iterative and noniterative methods for the solution of the nonlinear Richards equation, Water Resources Research, 27(6), 1147–1163, 1991.

Paniconi, C., Putti, M., A comparison of Picard and Newton iteration in the numerical solu- tion of multidimensional variably saturated flow problems, Water Resources Research, 30(12), 3357–3374, 1994.

Pot, V., Appert, C., Melayah, A., Rothman, D.H., Zaleski, S., Interacting lattice gas automaton study of liquid–gas properties in porous media, Journal de Physique II , 6, 1517–1534, 1996.

Qian, Y., d’Humières, D., Lallemand, P., Lattice BGK models for Navier-Stokes equation, Europhysics Letters, 17, 479–484, 1992.

Rothman, D.H., Zaleski, S., Lattice–Gas Cellular Automata: Simple Models of Complex Hydrodynamics, Cambridge University Press, Cambridge, UK, 1997.

Soll, W.E., Celia, M.A., Wilson, J.L., Micromodel studies of three-fluid porous media sys- tems: pore-scale processes relating to capillary pressure-saturation relationships, Water Resources Research, 29(9), 2963–2974, 1993.

Smyth, J.D., Yabusaki, S.B., Gee, G.W., Infiltration evaluation methodology – letter report 3: Selected tests of infiltration using two dimensional numerical models, Technical Report, Pacific Northwest Laboratory, Richland, WA, 1989.

Magnuson, S.O., Baca, R.G., Sondrup, A.J., Independent verification and benchmark testing of the porflo-3 computer code, version 1.0, Technical Report EGG-BG-9175, Idaho National Engineering Laboratory, Idaho Falls, ID, 1990.

Soll, W.E., Birdsell, K., The influence of coatings and fills on flow in fractured, unsaturated tuff porous media systems, Water Resources Research, 34(2), 193–202, 1998.

S. STRAFACE AND G. MENDICINO

Succi, S., Benzi, R., Higuera, F., The lattice Boltzmann equation: a new tool for computa- tional fluid dynamics, Physica, 47 (D), 219–230, 1991.

Sukop, M.C., Or, D., Lattice Boltzmann method for modeling liquid–vapor interface con- figurations in porous media, Water Resources Research, 40, W01509, 2004.

Toffoli, T., Margolus, N., Cellular Automata Machines, MIT Press, Cambridge, MA, 1987.

Engineering and Science, 1(1), 11, 2001. Troisi, S., Problems of mathematical model validation in groundwater hydrology. Excerpta,

9199–236, 1995. Troisi, S., Fallico, C., Straface, S., Migliari, E., Application of kriging with external drift to

estimate hydraulic conductivity from electrical resistivity data in unconsolidated deposits near Montalto Uffugo, Italy, Hydrogeology Journal, 4, 356–367, 2000.

van Genuchten, M.Th., Nielsen, D.R., On describing and predicting the hydraulic proper- ties of unsaturated soils, Annales Geophysicae , 3(5), 615–628, 1985.

von Neumann, J., Theory of Self-Reproducing Automata, University of Illinois Press, Urbana, Illinois, 1966.

Wolfram, S., Cellular Automata and Complexity, Addison-Wesley, Reading, Mass, 1994. Wolfram, S., Theory and Application of Cellular Automata, World Scientific, Singapore,

1986. Zhang, D., Zhang, R., Chen, S., Soll, W.E., Pore scale study of flow in porous media: scale

dependency, REV, and statistical REV, Geophysical Research Letters, 27, 1195– 1198, 2000.

Zhang, X., Bengough, A.G., Crawford, J.W., Young, I.M., A lattice BGK model for advec- tion and anisotropic dispersion equation, Advances in Water Resources, 25, 1–8, 2002a.

Zhang, X., Bengough, A.G., Deeks, L.K., Crawford, J.W., Young, I.M., A novel three- dimensional lattice Boltzmann model for solute transport in variably saturated porous media, Water Resources Research, 38, 1167–1177, 2002b.

Zhang, X., Deeks, L.K., Bengough, A.G., Crawford, J.W., Young, I.M., Determination of soil hydraulic conductivity with the lattice Boltzmann method and soil thin-section technique, Journal of Hydrology, 306, 59–70, 2005.

Ziegler, B.P., DEVS theory of quantized systems, Advanced simulation technology thrust DARPA contract, 1998.

Tonti, E., A discrete formulation of field laws: the cell method, Computational Methods in

Straface S., A flow and transport numerical model based on the method of cells. PhD (in Italian), University of Calabria, Cosenza, Italy, 1998.

Straface, S., Troisi, S., Gagliardi, V., Application of the cell method to the simulation of unsaturated flow. Computers, Materials, & Continua, 3(3), 155–165, 2006.

van Genuchten, M.Th., Calculating the unsaturated hydraulic conductivity with a new closed- form analytic model, Water Resources Program Report 78-WR-08, Department of Civil Engineering, Princeton University, Princeton, NJ, 1978.

Agent-based modeling of socio-economic processes related to the environment: Example of land-use change

J. Gary Polhill

Macaulay Institute, Craigiebuckler, Aberdeen AB15 8QH, UK.

Abstract

This chapter discusses some of the principles behind multi-agent modeling and shows through the examples from land-use change how they can be applied to deal with socio-economic aspects of environmental issues. An underlying theme is the divide between qualitative and quantitative approaches in the social sciences, though the chapter is also aimed at presenting agent-based modeling to those accustomed to mathematical modeling approaches.

Keywords: agent-based modeling, qualitative-quantitative debate, land-use change.

1. Introduction

Traditionally, there have been two options for describing and analyzing phenomena of interest. If sufficient measured data are available, then various mathematical tech- niques may be employed, essentially fitting a function (or functions) to the data describing the trajectory of the system. The resulting equations, and subsequent analysis conducted with them, particularly in a policy context, amount to a process of developing a narrative of the phenomena concerned. That is, the equations are not an end in themselves but a tool used as part of a supposedly rigorous process. The analytical techniques used to develop and assess the equations have ultimately derived from physics. It is a testament to their utility that such techniques have been applied beyond physics into other sciences, including ecology and the social sciences. In the social sciences (and arguably in ecology), however, the techniques may become somewhat strained and unwieldy in their application. The phenomena are, furthermore, so intricately interrelated, and individuals so varied, that doing justice to the complexity is often not possible with the language of mathematics. The only other language available has been natural language.

61 P.C. Baveye et al. (eds.), Uncertainties in Environmental Modelling