C h a p t e r X I

overwhelmedblueearthAI and Robotics

Dec 1, 2013 (3 years and 7 months ago)

67 views

Chapter XI
Cellular Automaton Mo dels of
Reaction T ransp ort Pro cesses
Theo KARAPIPERIS

P aul Sc herrer Institut
CH  Villigen PSI
Switzerland
XI  In tro duction
In man y pro cesses o ccurring in geological media mass transp ort and c hemical reactions
are in tricately connected The range of these pro cesses includes natural phenomena suc h
as the c hemical w eathering of ro c ks and the diagenetic c hanges caused b y groundw ater
o w as w ell as pro cedures of tec hnological in terest suc h as oil reco v ery and con tainmen t
of con taminan t migration In case of a failure in a n uclear w aste rep ository  the released
radion uclides will b e transp orted b y o wing w ater while they deca y radioactiv ely  sorb on
solid surfaces and react c hemically with eac h other and the host ro c k Mo delling coupled
transp ort and c hemical reactions in their full complexit y is a formidable task and a great
deal of ingen uit y is required to iden tify the most imp ortan t features and pro cesses in an y
giv en situation
In spatially extended geological systems ho w ev er complexit y at the macroscopic lev el
need not reect particularly in tricate asp ects of the microscopic dynamics On the con
trary  it is w ellkno wn that relativ ely simple microscopic prop erties can lead to v ery ric h
macroscopic b eha viour Th us the motion of uids laminar and turbulen t o w sho c k
w a v es o w in p orous media can b e understo o d in terms of lo cal mass momen tum and
energy conserv ation The random collisions of the molecules of a few c hemical sp ecies
in solution can couple with their c hemical reaction kinetics to pro duce a wide v ariet y of
￿
Presen t address Europ ean P arliamen t Rue Belliard  B Brussels Belgium

Extract from MODELLING IN AQUATIC CHEMISTRY
(OECD Publications, 1997, 724 pp., ISBN 92-64-15569-4)
Scientific Editors: Ingmar Grenthe and Ignasi Puigdomenech Contributors: Bert Allard, Steven A. Banwart,
Jordi Bruno, James H. Ephraim, Rolf Grauer, Ingmar Grenthe, Jörg Hadermann, Wolfgang Hummel,
Andreas Jakob, Theo Karapiperis, Andrey V. Plyasunov, Ignasi Puigdomenech, Joseph A. Rard, Surendra
Saxena, Kastriot Spahiu Secretariat: OECD Nuclear Energy Agency Data Bank: M.C. Amaia Sandino and
Ignasi Puigdomenech Original text processing and layout: OECD Nuclear Energy Agency Data Bank: Cecile
Lotteau. © OECD 1997

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
striking b eha viour oscillations of the concen trations in time stationary spatial patterns
or nonlinear tra v elling w a v es
In theoretical mo dels of suc h phenomena the essen tial microscopic prop erties are usu
ally em b edded in macroscopic p artial di er ential e quations  eg  the con tin uit y equation
and the Na vierStok es equations for uid motion or the reactiondi usion equations for
the concen trations of dissolv ed sp ecies The macroscopic attributes describ ed b y these
equations uid v elo cit y  solute concen tration are related to the underlying microscopic
quan tities b y an a v eraging op eration whic h extends o v er a v olume elemen t small on the
t ypical v olume scale of in terest but large compared with the a v erage v olume p er molecule
In the pro cess of a v eraging suc h information as microscopic uctuations is discarded A
basic premise of the macroscopic approac h is that spatial and temp oral v ariations are
su
cien tly slo w for di eren t parts of the system to b e assumed in a state of lo cal ther
mo dynamic equilibrium at all times The macroscopic equations are usually nonlinear
and ha v e to b e solv ed n umerically When iterating a n umerical algorithm in a computer
roundo errors arising when real n um b ers are truncated to nite computer w ords can
accum ulate exp onen tially  leading to o v ero ws Guaran teeing the stability of n umerical
algorithms is a tedious pro cedure for nonlinear problems In con trast to the macroscopic
approac h mole cular dynamics deals with the motion of individual molecules under the
inuence of realistic in termolec ular p oten tials In this approac h one addresses problems of
a more fundamen tal nature suc h as understanding at a molecular lev el the phenomeno
logical co e
cien ts of macroscopic mo dels BOOYIP Ob viously the length and time
scales relev an t to molecular dynamics are m uc h shorter than those t ypical of a macroscopic
description
There is a large class of problems in v olving natural systems that lie in the in termediate
region b et w een the macroscopic and the microscopic description These naturally include
problems with a c haracteristic length scale that is m uc h larger than that of molecular
ev en ts  

A  
￿ 
cm but also m uc h smaller than macroscopic length scales for
example bulk prop erties of a ro c k suc h as p orosit y  can b e t ypically de ned on a scale
of  cm W e refer to suc h an in termediate scale as mesosc opic Th us pro cesses
taking place at p ore lev el in a ro c k t ypical p ores ha v e linear dimensions of the order
of  m  
￿ 
cm suc h as c hanges in p orosit y due to precipitationdissolution
b elong in trinsically to this in termediate class of problems Di eren tial equations ha v e
b een emplo y ed in mesoscopic mo delling of natural systems but only under con v enien t
assumptions concerning the often v ery irregular b oundaries b y assigning for example the
o w in p orous media to a net w ork of capillary tub es of simple geometry DUL Ev en at
t ypically macroscopic scales some kind of in termediate description ma y b e indisp ensable
when the complexit y of the b oundary conditions or the v ariabilit y of certain quan tities
in space and time call for a degree of detail that cannot b e a orded b y the standard
macroscopic equations
The t yp e of problems addressed at eld scale sev eral metres or more are of a v ery
di eren t nature from those solv ed b y t ypical mesoscopic mo dels The spatial discretisa
tion appropriate for eldscale problems is m uc h coarser than the spatial detail that has
to b e accoun ted for at the mesoscopic scale This is dictated in the rst place b y reasons

In tro duction
of computational feasibilit y  but also and most imp ortan tly for conceptual reasons the
ev olution equations for macroscopic quan tities suc h as solute concen trations can usually
b e expressed in terms of macroscopic bulk parameters whic h are de ned b y a v eraging out
ner details Con tact b et w een mo delling at di eren t scales is established at the in terface
where mesoscopic mo dels pro vide the relations b et w een bulk parameters and their ev olu
tion in time whic h cannot b e deriv ed or are kno wn only empirically at the macroscopic
lev el The same is true ab out di eren t time scales where the times t ypical of geological
phenomena ma y often render irrelev an t the kinetics of sp eci c reactions As w e shall
see later in this c hapter there are notable exceptions where certain phenomena cannot
b e mo delled in terms of purely macroscopic parameters and a mesoscopic description is
indisp ensable irresp ectiv e of the computational resources that this ma y require
In this c hapter w e are going to presen t a new and promising metho d of mo delling
mass transp ort and c hemical reactions at the mesoscopic lev el The basic concept in
this approac h is that of a c el lular automaton i e  a dynamic al system c onsisting of an
inte ger eld dene d on the sites of a r e gular sp atial lattic e and evolving in discr ete time
steps ac c or ding to a lo c al up dating rule the rule determines the new value of the eld at a
c ertain site of the lattic e ac c or ding to the curr ent values of the sites within a neighb ourho o d
surr ounding the site in question W OL In the case of a system of dissolv ed sp ecies
migrating in a p orous or fractured medium and reacting with eac h other as w ell as with
the host ro c k one can think of the lattice as spanning the ro c k v olume of in terest and the
lo cal rule as mo delling the collisions and the c hemical reactions of the solute particles
T ransp ort and c hemical reactions are mo delled in a w a y that preserv es a limited n um b er
of microscopic features most notably  microscopic uctuations and correlations as w ell
as reaction kinetics whic h are exp ected to inuence the macroscopic b eha viour w e wish
to describ e Moreo v er since the b eha viour of the solute particles up on con tact with the
solid b oundaries is mo delled in a ph ysically motiv ated manner the complexit y of the
b oundaries has no b earing on mo delling feasibilit y or ev en computational e
ciency The
sync hronous ev olution of a cellular automaton reects the lo c al and p ar al lel c haracter of
the natural pro cesses it mo dels molecules migrate and react sim ultaneously at di eren t
lo cations inside the ro c k and eac h molecule reacts only with molecules in its imm ediate
neigh b ourho o d The parallel and lo cal prop erties of cellular automata map naturally
on to the arc hitecture of massiv ely parallel computers In the latter a large n um b er of
pro cessors up to sev eral tens of thousands run in parallel manipulating data from their
lo cal memories data are transferred b et w een pro cessors via sophisticated comm uni cations
net w orks whic h are particularly e
cien t for transfers b et w een neigh b ouring pro cessors
Massiv ely parallel arc hitectures are exp ected to dominate the design of new computers in
the coming y ears since they o er the only w a y to impro v e b y orders of magnitude on the
p erformance of curren t scalar and v ector computers it tak es an electromagnetic signal
some nanoseconds to tra v erse a t ypical length of  m and this sets an absolute upp er limit
to the clo c k frequency of serial pro cessors
The remainder of this c hapter is organised as follo ws The dev elopmen t of cellular
automata CA from their initial inception as mo dels of biological systems to their presen t
implem en tati on in the sim ulation of a wide range of ph ysical phenomena is outlined in

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
the next section In the same section w e illustrate the basic concepts with the help of
elemen tary examples F ollo wing that w e presen t CA mo dels for coupled transp ort and
c hemical reactions and describ e their applications to reactiontransp ort systems ranging
from the simplest and b etter understo o d to increasingly more complex and realistic ones
In the nal section w e ev aluate the p oten tial of the CA approac h and indicate p ossible
future dev elopmen ts It should b e emphasised that this c hapter has a rather p edagogical
orien tation The author will allo w a certain p ersonal bias in the selection and presen tation
of the relev an t material and do es not exp ect this w ork to b e seen as a review article
ev en more so since reviews and extensiv e compilations of articles are already a v ailable
F ARTOF MANBOC GUT CHED A W
XI  Cellular automata
XI Historic al development
The theory of automata w as en visaged b y John v on Neumann as a systematic mathemat
ical and logical framew ork that w ould unify the basic principles of natural and arti cial
information pro cessing systems ASPBUR ASP The idea of a uni ed approac h
app eared in the early da ys of the digital computer and w as motiv ated b y the plausible
analogy b et w een computers as w ell as other comm uni cations and con trol devices and
selfregulating biological systems cells the h uman nerv ous system organisms capable of
repro duction and ev olution T ypical questions addressed w ere those of reliabilit y  ie 
the abilit y to detect errors and minim ise their e ect and of the minimal complexit y that
enables an automaton to construct other equally complicated automata A universal c on
structor  ie  an automaton that can construct an y other automaton w ould b e equiv alen t
to the universal T uring machine  the theoretical protot yp e of a generalpurp ose computer
A computer A is called univ ersal if giv en an y other computer B  there is alw a ys an ap
propriate piece of soft w are whic h enables A to em ulate B in the sense that the t w o devices
deliv er the same output for an y giv en input PEN F ollo wing initial di
culties with a
threedimensional kinematic automaton that w ould ph ysically assem ble its duplicate out
of a p o ol of elemen tary parts v on Neumann heeded the advice of Stanisla w Ulam and
considered the abstract setting of c el lular automata  ie  t w odimensional square arra ys
of cells eac h of whic h assumed one out of the same nite set of states

Although cel
lular automata could b e directly translated in to parallel computer structures in theory 
the size and p o w er consumption of hardw are comp onen ts a v ailable at the time rendered
a ph ysical impleme n tation of the cellular automaton paradigm impractical In terestingly 
v on Neumanns name is usually asso ciated only with the serial c omputers dev elop ed in his
lifetime essen tially  he also laid the theoretical foundations for the massively p ar al lel c om
puters built later when in tegrated circuit tec hnology drastically enhanced the pro cessing
and memory storage capacit y of single c hips
Although v on Neumanns program of a comprehensiv e theory of automata w as nev er
￿
It is p ossible that other p eople arriv ed at similar ideas indep enden tly
see e g  ZUS  but v on
Neumanns w ork is b y far the most extensiv e alb eit itself incomplete on published record

Cellular automata
completed cellular automata CA dev elop ed their o wn dynamic in the s and  s
Information concerning CA implem en tations in sp ecialpurp ose hardw are and general
purp ose computers can b e found in Refs HIL PREDUF TOFMAR BOG
Largely as a result of the ongoing dramatic rise in a v ailable computing p o w er in general
and the adv en t of massiv ely parallel computers in particular CA algorithms are increas
ingly emplo y ed to sim ulate natural phenomena A wide eld of applications op ened
when CA w ere applied to the sim ulation of uid dynamics DOOFRI DOO fol
lo wing the realisation that a certain class of CA appro ximate the Na vierStok es equations
FRIHAS

The CA of this class ha v e come to b e kno wn as lattic e gas automata
LGA A wide range of o w phenomena ha v e b een mo delled in the mean time  including
the o w of m ultiphase mixtures  eg  oil and w ater R OTKEL SOMREM and
phase transitions  eg  liquidgas APPZAL APPZAL turbulence SUCSAN
magnetoh ydro dynamics  ie  the motion of electrically conducting uids suc h as plas
mas in the presence of a magnetic eld CHEMA T MONDOO colloidal sus
p ensions LADCOL and o w in p orous media BALHA Y R OT Due to the
simplicit y of their dynamics LGA ha v e b een used as mo dels to address fundamen tal
problems in kinetic theory  suc h as the longtime b eha viour of time correlation func
tions HARP AZ FRE FREERN and the relation of transp ort co e
cien ts  eg 
viscosit y  thermal conductivit y to the time in tegral of correlation functions RIV
ERNDUF ERND AS CA ha v e also b een used to mo del v arious ph ysico c hemical
pro cesses of engineering in terest carb onation of concrete BRIBON catalytic COg
o xidation on platin um surfaces WUKAP p olymer c hains in solution VIAK OE
c hargedparticle transp ort in semiconductors K OMZAN biological functions  eg 
the imm une resp onse KA UURB or the bio c hemistry of the cell HASKAP and
astroph ysical phenomena PERLEJ
What in terests us here are CA sim ulations of systems of molecules that mo v e b y means
of di usion andor adv ection and react c hemicall y with eac h other and with the surround
ing medium  eg  reactiv e c hemicals di using in a gel or solutes migrating in a p orous
medium CA mo dels of mass transp ort with c hemical reactions w ere dev elop ed either
along lines similar to the LGA mo dels of uid dynamics D ABBOO D ABLA W
LA WD AB or indep enden tly BLA KAR KARBLA Before w e describ e these
mo dels in detail w e shall illustrate the basic concepts b y means of a few elemen t ary exam
ples whic h will also demonstrate the capacit y of simple CA rules to pro duce in teresting
complex b eha viour
XI Elementary examples
Life is a mathematical game in whic h a t w odimensional square arra y of cells displa ys
a lot of the liv ely activit y t ypical of an assem bly of living organisms GAR Beginning
with some initial con guration of liv e cells state  one follo ws the ev olution of the
￿
Although there is no general prescription for nding the microscopic rules that will sim ulate a giv en
set of macroscopic equations the CA approac h has b een applied to the n umerical solution of v arious
partial dieren tial equations of ph ysical in terest BOGLEV CHECHE CHEMA T

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
p opulation as previously dead sites state  come to life if there are exactly  liv e cells
among the neigh b ours on the square surrounding them and liv e cells p ersist if  or 
of their neigh b ours are aliv e or die otherwise This ev olution rule w as c hosen b y John
Con w a y  the in v en tor of the game so as to a v oid the more common scenarios of rapid
uncon trolled gro wth or extinction of the p opulation T ypical snapshots from the ev olution
of an originally random p opulation of  liv e cells are sho wn in Fig XI  F ollo wing an
initial burst of activit y and o ccasional later ares most con gurations settle do wn to a
sparse distribution of stationary and oscillating patterns last frame in Fig XI  Prop
agating patterns gliders and con gurations that gro w without limit glider guns
are less common but con tribute to the unpredictabilit y of Life When a glider meets
another group of liv e cells a new pattern is in general created Th us glider crashes can b e
designed to pro duce glider guns whic h in turn emit streams of fresh gliders An eater is
a pattern that annihilates other patterns up on colliding with them One can build logical
gates NOT AND OR out of glider guns and eaters It turns out that one can construct
an arbitrary computer out of suc h gates so that Life is capable of univ ersal compu
tation in the sense of the univ ersal computer de ned in Section XI   BERCON
It is a remark able fact that a CA with no more than t w o p ossible states p er site and a
v ery simple ev olution rule can p erform all the complicated tasks that a generalpurp ose
computer can
As a second example w e men tion the CA implem en tati on of the billiardball mo del of
computation In the latter the elastic collisions of a gas of hard spheres with eac h other
and their reections b y a set of mirrors are used to construct a univ ersal computer
FRETOF W e consider a square lattice of cells as b efore eac h with t w o p ossible
states and w e divide them in    blo c ks Giv en suc h a partition of the lattice w e de ne
the up dating rule of a    blo c k as sp eci ed in Fig XI a complem en te d with all p ossible
rotations of the transitions indicated MAR Bet w een successiv e applications of the
rule w e shift the partitioning grid along the diagonal b y
p
 times the lattice spacing from
the solid to the dashed lines in Fig XI b where the legend solid or dotted indicates
the partitioning for the next application of the rule If w e let s and s denote particles
and empt y sites resp ectiv ely  iteration of the ab o v e rule mak es particles propagate and
collide elastically with eac h other Fig XI b Tw o adjacen t blo c ks lled with particles
form a stable con guration that acts as a mirror It has b een sho wn MAR that the
CA w e ha v e just describ ed is equiv alen t to the billiardball mo del and is consequen tly
a univ ersal computer This is an example of a rev ersible CA ie  one whose ev olution
can b e traced bac kw ards in time in direct analogy with the microscopic rev ersibilit y of
ph ysical la ws The realisation that rev ersible univ ersal computers are feasible has far
reac hing implications for the ph ysical limitations of computers since it app ears to remo v e
the lo w er b oundary on dissipation asso ciated with standard irrev ersible logic elemen ts
BEN MAR
As a last example w e sho w in Fig XI  di eren t stages in the ev olution of a CA whic h
dev elops spatial structures similar to the concen tration patterns observ ed in certain c hem
ical exp erimen ts The sim ulation is p erformed with a triangular lattice of cells that exist
in three di eren t states activ ated white receptiv e blac k and quiescen t grey

Cellular automata
Figure XI  Initial random distribution of  liv e cells and distribution after  
and  iterations of Life
During one up date all receptiv e cells among the nearest neigh b ours of an activ ated cell
b ecome activ e while the activ ated cell itself turns quiescen t A quiescen t cell cannot b e
excited un til it b ecomes receptiv e after t w o time steps MADFRE Di eren t patterns
suc h as circular w a v es and spirals dev elop dep ending on the initial con guration The
one and t w oarmed spirals in Fig XI  arise from isolated lines of activ e cells bu ered b y
lines of quiescen t cells It can b e seen that when the spirals collide they halt eac h other
along their line of con tact but con tin ue to spread unhindered elsewhere Similar pat
terns ha v e b een observ ed in exp erimen ts with auto catalytic c hemical reactions WIN
ZHAZAI W e are going to return to the sub ject of auto catalytic reactions in the next
section
Ha ving seen ho w simple CA can sim ulate the motion and collisions of particles as w ell
as e ects reminiscen t of complex c hemical systems w e can no w pro ceed to describ e in

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
Figure XI  a Up dating rule for CA impleme n tation of billiardball mo del of com
putation b Propagation and collision of particles according to this rule reprin ted with
p ermission from Ref MAR Figs  and 

Cellular automata for transp ort with c hemical reactions
Figure XI  One and t w oarmed spirals ev olving out of isolated lines of activ ated cells
white bu ered b y lines of quiescen t cells grey
detail more elab orate CA mo dels that couple these pro cesses
XI  Cellular automata for transp ort with c hemical reactions
XI
 Mo dels
XI
 T r ansp ort
F rom a macroscopic p oin t of view solutes migrate in a p orous medium b y means of adve c
tion  mole cular di usion and kinematic disp ersion MAR Adv ection tak es place when
a solute is carried along b y the o w of the solv en t it is c haracterised b y a v elo cit y  whic h
for su
cien t dilution coincides with the Darcy v elo cit y of the pure solv en t Molecular
di usion is the macroscopic manifestation of the Bro wnian motion of the solute parti
cles caused b y random impacts with uid molecules in an ideal solution it is describ ed
b y Fic ks la w whic h in tro duces the c o e cient of mole cular di usion Finally  kinematic
disp ersion is the spreading of the solutes around the mean adv ectiv e displacemen t due
to microscopic v ariations in o w v elo cit y inside individual p ores but also b et w een p ores
of di eren t ap erture and length kinematic disp ersion is commonly accoun ted for in the
transp ort equation b y a di usion term whic h can b e com bined with the analogous term
for molecular di usion to de ne an e ectiv e o v erall di usion co e
cien t In what follo ws
w e shall use di usion to describ e the com bined e ect
The breakdo wn of transp ort in to adv ection and di usion is rather arbitrary and b ecomes

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
less meaningful as w e fo cus on the ner structure of the p orous medium A t the mesoscopic
lev el particle transp ort can b e though t of as a single pro cess whic h consists of a collection
of particles p erforming indep enden t random w alks

W e assume that particles mo v e
b et w een the sites of a regular spatial lattice In discrete time steps particles jump to
one of their nearestneigh b our sites with a certain probabilit y If the particles jump
in all directions with the same probabilit y  then the random w alk manifests itself at
the macroscopic lev el as pure di usion If on the other hand particles tend to jump
preferen tially in a giv en direction this leads macroscopically to com bined adv ection in the
fa v oured direction and di usion FEL W e consider for simplicit y a onedimensional
lattice since extensions to higher dimensions are straigh tforw ard Let the lattice spacing
and the time step b e  and  resp ectiv ely W e refer to sites on the lattice b y their lo cation
in space eg  site x During one time step eac h particle mo v es to the next site on the
righ t with a certain probabilit y p  to the next site on the left with probabilit y q of course
p  q   or remains stationary  with probabilit y   p  q It is useful to see the
system as a xed arra y of sites the state of site x at time t is giv en b y the in teger eld
N  x t  the n um b er of particles at that site at the giv en time the ev olution of N  x t  is
determined b y the rule de ning the random w alk The arra y of sites together with the
in teger eld and the ev olution rule constitute a CA that mo dels the transp ort of solute
particles
W e run the sim ulation on an ensem ble of lattices sub ject to the same mesoscopic la ws
 ie  same p and q  with macroscopically iden tical initial and b oundary conditions Let
us rst assume that p and q are constan t in space and time If w e de ne the particle
densit y   x t  as the ensem ble a v erage of N  x t  it can b e easily seen that  ev olv es in
time according to the di erence equation
  x t     p  x   t     p  q    x t   q   x   t  XI 
whic h can b e readily rearranged to
  x t       x t 

  V
  x   t     x   t 
 
 D
  x   t      x t     x   t 


 XI 
with the parameters V and D de ned as
V   p  q 


 D   p  q 


 
 XI 
￿
W e in ten tionally talk ab out solute particles and their densit y rather than ab out solute
molecules and concen tration as it is clear that in the foreseeable future the n um b er of par
ticles in all CA sim ulations of exp erimen tally or naturally realisable pro cesses m ust for practical
reasons fall short of the n um b er of molecules actually in v olv ed b y sev eral orders of magnitude Con
tact with the real pro cess is established b y rescaling the densit y at the lev el of the output KAR
KARBLA

Cellular automata for transp ort with c hemical reactions
If w e tak e the con tin uum limit       and p  q   while k eeping 

 and
 p  q   nite Eq XI  go es o v er to the transp ort equation
   x t 
 t
  V
   x t 
 x
 D


  x t 
 x

 XI 
The righ t hand side of Eq XI  con tains an adv ectiv e and a di usiv e term W e iden tify
V as the adve ction velo city and D as the di usion c o e cient The con tin uum limit is
tak en in suc h a w a y that V and D  as de ned b y Eq XI  remain nite Eq XI  is
the forw ardtime cen tredspace nite di erence appro ximation to the transp ort equation
XI  It is clear that b y making p and q sp ecies dep enden t w e can ha v e sp ecies with
di eren t transp ort prop erties mo v e on the same lattice If V and D are not homogeneous
the system can b e mo delled b y making p and q p osition dep enden t then if V and D
are de ned according to Eq XI  the adv ection v elo cit y has to b e de ned b y

V  x  
V  x   dD  x  dx 
Boundary conditions are mo delled b y mimic ki ng the b eha viour of actual molecules
when they reac h b oundaries with the sp eci ed prop erties Th us when particles reac h an
imp ermeable b oundary  they b ounce bac k or remain stationary The freeo w b ound
ary condition according to whic h only adv ection carries particles across the b oundary
    x t   t   for x at the b oundary is sim ulated b y cancelling the outgoing di usiv e
ux with an equal and opp osite di usiv e ux coming from a ctitious reserv oir on the
other side of the b oundary A xed densit y at the b oundary is mo delled b y c ho osing the
n um b er of particles at b oundary sites from a distribution consisten t with the sp eci ed
a v erage v alue this sim ulates the situation where the b oundary sites are in con tact with
a particle reserv oir of the appropriate densit y A p erio dic b oundary condition is often
computationally con v enien t and is justi ed when the precise b eha viour of the b oundary
sites is not imp ortan t eg  with translationally in v arian t homogeneous systems or when
the b oundary is to o far a w a y to inuence the region of in terest
In a t ypical sim ulation of the kind describ ed ab o v e it cannot b e excluded that at a
certain p oin t in the sim ulation a large n um b er of particles are found at the same lattice
site Giv en the total n um b er of particles initially presen t and the rate at whic h particles
en ter and lea v e the system an upp er limit can b e calculated in adv ance of the sim ulation
This limit is ho w ev er of little v alue in practice since it is highly unlik ely to ev er b e reac hed
in an actual sim ulation It is often desirable in computer sim ulations to allo w a priori 
only a limited n um b er of particles p er site This mak es it p ossible to program at bit
lev el and construct lo okup tables of all conceiv able inputoutput com binations that ma y
o ccur at a site One is then able to pro duce v ery e
cien t CA sim ulations on dedicated
hardw are TOFMAR but also to substan tially impro v e p erformance on serial and
v ector computers SHIDOO A restriction to the the n um b er of particles p er site is
imp osed in LGA mo dels of uid dynamics in the form of an exclusion principle  a small
n um b er of v elo cities are assigned to eac h lattice site  eg  in a singlesp eed mo del the
v ectors of magnitude  p oin ting to the nearest neigh b ours and at most one particle
p er site is allo w ed to ha v e a giv en v elo cit y As a result the state of an LGA is de ned b y

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
a Bo olean eld that giv es the o ccupation n um b ers  or  of the v elo cit y states at eac h
site
One can also mo del di usion with v ery e
cien t CA algorithms whic h ob ey an exclusion
principle CHODR O W e consider again a onedimensional lattice and assign t w o
p ossible v elo cities to eac h site c

     e

p oin ting to the righ t and c

     e

to
the left   e

and  e

are unit v ectors One ev olution step of the automaton consists of t w o
elemen tary op erations First a rotation op erator acts indep enden tly on di eren t no des
and causes the v elo cit y con guration at eac h of them to remain in tact or rotate b y 
with probabilit y p and   p  resp ectiv ely follo wing this particles propagate to the next
lattice site along the direction of their v elo cit y Let the Bo olean v ariable N
i
 x t  denote
the n um b er of particles with v elo cit y c
i
at site x and time t Then the ensem ble a v erage

i
 x t  of N
i
 x t  satis es the discrete ev olution equation



 x t   


 x t   

 p



 x   t 


 x   t 

   p 



 x   t 


 x   t 

 XI 
The rst second term on the righ t hand side of Eq XI  accoun ts for particles that
propagate to the next site without after b eing rotated W e de ne the particle densit y
  x t  
P
i

i
 x t 
Since the n um b er of particles is conserv ed a con tin uit y equation holds in the form
  x t     x t    

X
i 
J
i
 x t     XI 
where J
i
 x t   
i
 x    e
i
 t   
 ￿ i
 x t  is the net ux of particles tra v ersing the lattice
link from x to x    e
i
in the time b et w een t   and t T aking the con tin uum limit of
Eqs XI  and XI         k eeping 

 nite w e arriv e after some simple
algebraic manipulations at the di usion equation Eq XI  without the adv ection term
with
D 
p
  p



 
 XI 
Sp ecies with di eren t di usion co e
cien ts can b e sim ulated on the same lattice b y ad
justing p Di usion in t w o and three dimensions is treated in Ref CHODR O The
adv an tage of this algorithm is that it can b e readily implem en ted on sp ecialpurp ose
computers and its p erformance on v ector sup ercomputers exceeds that of the previous al
gorithm b y more than an order of magnitude On the other hand the exclusion principle
induces correlations b et w een particles and these in turn giv e rise to nonlinear terms in
the transp ort equations deriv ed from mo dels with biased random w alks BOGLEV
As a result mo dels with an exclusion principle are not able to mo del adv ection without
in tro ducing at the same time un w an ted terms in the transp ort equation
Mo delling di usion and adv ection with CA w ould b e an exercise of little practical
in terest w ere it not to b e coupled with c hemical reactions a task that w e immedi ately
pro ceed to ful l

Cellular automata for transp ort with c hemical reactions
XI
 Chemic al r e actions
In the full coupled problem an ev olution step consists of a transp ort op eration follo w ed
b y a c hemical reaction op eration W e rst consider a system with no restriction on the
n um b er of particles p er site T ransp ort tak es place as describ ed ab o v e for suc h a system
P articles of v arious solute sp ecies mo v e on the lattice and when su
cien t n um b ers of
them meet at a site they react with a certain probabilit y A reaction is more lik ely to
o ccur the more the particles of eac h reactan t that are presen t If a reaction tak es place
then the appropriate n um b er of pro ducts reactan ts are added to subtracted from the
lo cal in v en tory of particles
W e consider a system of sp ecies s

and reactions r W e write all reactions in the one
w a y form b y separating in the case of rev ersible reactions the forw ard from the rev erse
pro cess
X


r i
s


X


r f
s

XI 
where
r i
and
r f
are the stoic hiometric co e
cien ts of the initial and nal state of the
reaction resp ectiv ely The summations run o v er all sp ecies but the
r i

r f
 v anish
for sp ecies not presen t among the reactan ts pro ducts If reaction r in v olv es only sp ecies
in solution w e let it pro ceed with probabilit y
P
r
Y


r i
Y
m 
 N

 x t   m    XI 
where P
r
is a constan t and N

 x t  is the n um b er of particles of sp ecies s

at site x
and time t F or example if r is the reaction a   b  c  it tak es place with probabilit y
P
r
N
a
N
b
 N
b
  Of course the reaction probabilit y has to b e   since N
a
and N
b
are nite in a sim ulation w e can alw a ys ensure that P
r
N
a
N
b
 N
b
    b y c ho osing
P
r
su
cien tly small W e shall see shortly that it is the ratio P
r
 that is ph ysically
meaningful so  has to b e c hosen corresp ondingly small with a concomitan t increase in
the n um b er of iterations In the case of reactions b et w een solutes and mineral surfaces
the reaction rule has to b e mo di ed so that precipitation tak es place with a probabilit y
de ned as b efore but only if solid is already presen t or the solution is saturated
W e wish to deriv e an ev olution equation for the particle densit y 

 x t  de ned as an
ensem ble a v erage of N

 x t  T o this end w e ha v e to mak e t w o imp ortan t assumptions
In ev aluating the a v erage of the pro duct in Eq XI  w e assume that there are no
correlations b et w een the N

for di eren t sp ecies  mole cular chaos  The resulting simple
a v erages  eg 
N
b
 N
b
  in the example ab o v e can b e expressed in terms of the densities
if the N

do not v ary appreciably o v er sev eral lattice spacings  smo othness assumption 
then w e can ev aluate 

 x t  b y a v eraging lo cally around x o v er a uniform distribution
of particles the crucial fact is that for a uniform particle distribution N

 x t  ob eys a
P oisson distribution parametrised with 

 x t  FEL W e th us arriv e at the discrete
equation KARBLA


 x t     

 x t  
X
r

r f

r i
 P
r
Y



 r i

 x t   XI 

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
where the sum is o v er all reactions r  the pro ducts o v er all sp ecies
and


 x t   p



 x   t     p

 q

 

 x t   q



 x   t   XI 
is the result of the transp ort op eration on the particle distribution at time t  here p

and
q

are the sp eciesdep enden t displacemen t probabilities
Up on taking the con tin uum limit as b efore with the additional limiting pro cedure
P
r
  so that P
r
 remains nite w e obtain the standard transp ort equations for
reacting solutes
 

 x t 
 t
  V

 

 x t 
 x
 D





 x t 
 x


X
r

r f

r i
 k
r
Y



 r i

 x t  XI 
where
V

  p

 q




 D

  p

 q




 
 k
r
 P
r
  XI 
k
r
is the r ate c onstant of reaction r
W e emphasize that the CA sim ulations are not sub ject to the assumptions whic h w ere
in tro duced ab o v e with the purp ose of deriving the macroscopic equations The sim ulations
accoun t for microscopic e ects of correlations and uctuations that are t ypically a v eraged
out b y the macroscopic approac h The results of the discrete mo del w ere compared
carefully to the solution of the reactiontransp ort PDEs in Ref KARBLA Tw o
essen tial sources of discrepancy w ere encoun tered i correlations b et w een the particles
that arise as a result of c hemical reactions and ii statistical uctuations Systematic
study sho ws that the e ect of c hemicall y induced correlations disapp ears as the con tin uum
limit is approac hed for details see KARBLA If one sees CA algorithms only as
n umerical to ols for the solution of macroscopic equations correlations of this t yp e are no
more than n umerical artefacts asso ciated with the particular metho d of solution This is
to b e con trasted with the p oten tially observ able e ect that statistical uctuations can ha v e
on the b eha viour of reactiondi usion systems as it will b e illustrated in Section XI   Of
course ph ysical systems are after all discrete and so correlations will manifest themselv es
at some lev el Here the discrete mo del could pro vide a more accurate mo del than a orded
b y the PDEs
In reactiondi usion mo dels with an exclusion principle ones needs to ensure careful
b o okk eeping of the o ccupied v elo cit y states b efore and after a reaction In general there
are sev eral sp ecies on the lattice and the rotation and propagation op erations describ ed
previously are applied to eac h of them indep enden tly The dynamics of di eren t sp ecies
is coupled b y c hemical reactions of the t yp e giv en b y Eq XI  W e concen trate for
illustration purp oses on an e ectiv e singlesp ecies system assuming the concen trations of
other sp ecies to b e k ept constan t b y con tact with an appropriate reserv oir D ABBOO

Cellular automata for transp ort with c hemical reactions
D ABLA W LA WD AB In the notation of Eq XI  there is one activ e sp ecies
s

 and the reactions r are iden ti ed b y the stoic hiometric co e
cien ts
r i
and
r f

W e denote for simplicit y  s

b y X and
r i

r f
b y 
 resp ectiv ely The reaction
probabilit y P 
 is indep enden t of the initial and nal v elo cit y states if for a giv en
initial con guration there are n

p ossible nal v elo cit y con gurations of the reaction
X 
X  eac h o ccurs with the same probabilit y P 
 n


In tro ducing indices i
k
to di eren tiate b et w een the m di eren t v elo cit y states  i
k

      m  m   for a onedimensional lattice  for a t w odimensional square lattice etc 
w e note that a reaction con guration is fully de ned in terms of 
and the v elo cit y
states i

     i
q
 q  j 
j  of the particles that are created or destro y ed dep ending on
whether
or
 resp ectiv ely With eac h con guration th us de ned w e asso ciate
a random Bo olean v ariable 

 i

   i
q
 whic h tak es the v alue  when the reaction tak es
place with probabilit y P 
 n

 and is  otherwise If w e then de ne N
i
as b efore w e
can write the e ect of the c hemical reaction op eration as follo ws
N
i
 N
i

X
   
X
￿
i
￿
 i

￿ i
Q  f N
i
g  i

   i


X
￿
j
￿
 j
q
￿ i


 ij

   j
q


X
   
X
￿
i
￿
 i

￿ i
Q  f N
i
g  ii

   i


i

X
￿
j
￿
 j
q
 i
￿


 ij

   j
q
 
XI 
where
Q  f N
i
g  i

   i

 

Y
k 
N
i
k
Y
￿
l ￿
  
  N
i
l
 XI 
p opulates the initial state with particles of v elo cities i

     i

The primes constrain the
indices of the summation or m ultiplic ation to b e distinct In Eq XI  N
i
increases
decreases b y  if a particle with v elo cit y i is among the j 
j particles pro duced
destro y ed
W e tak e the ensem ble a v erage of the Eq XI  de ning 
i
as the a v erage of N
i
and
 
P
i

i
Assuming that w e can factorise the exp ectation v alue of pro ducts of random
v ariables and taking the con tin uum limit for a homogeneous system     P 
  
k eeping P
￿

  P 
  nite w e deduce the rate la w
d  t 
dt

m
X
n 

n

n
 t   XI 
where

n
   
n

m
n

m
￿ n
n
X

  


n


X
 ￿


  P
￿

  XI 
Iden tifying the 
n
with the macroscopic rate constan ts one can solv e the resulting system
of linear equations for the reaction probabilities P 
 F or an inhomogeneous system

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
w e obtain the standard reactiondi usion equations The generalisation to m ultisp ecies
systems is straigh tforw ard D ABBOO KAPLA W KAPLA W The exclusion
principle p ermits the utilisation of v ery e
cien t programming tec hniques as describ ed
earlier In the mo del describ ed in Refs D ABBOO D ABLA W LA WD AB
there are di eren t sets of mesoscopic reaction rules whic h are consisten t with the same
macroscopic kinetics It is p ossible ho w ev er to construct the reaction probabilities so
as to allo w only those particle n um b er c hanges that are sp eci ed b y the macroscopic
mec hanism WUKAP  GR UKAP
XI
 Applic ations
In this subsection CA mo delling will b e illustrated b y means of applications to v arious
reactiontransp ort pro cesses Our goal is to demonstrate the v ersatilit y of the metho d
as w ell as its abilit y to describ e asp ects of the microscopic dynamics esp ecially when the
latter ha v e macroscopic consequences that lie b ey ond the scop e of the standard a v eraged
equations Th us when microscopic uctuations determine under certain conditions the
macroscopic b eha viour of ph ysical systems the mesoscopic lev el of description a orded
b y CA mo dels b ecomes indisp ensable
XI
 a  b  c
W e rst consider a homogeneous system consisting of equal amoun ts of t w o mobile sp ecies
whic h com bine to form an inert sp ecies There is no external supply of reactan ts so the
t w o activ e sp ecies will deplete eac h other un til there are only inert particles left W e
are in terested in the time dep endence of the reactan t particle densit y as the nal inert
state is approac hed A system of this kind w as analysed in connection with the m utual
annihilation of magnetic monop oles and an timonop oles pro duced in the early univ erse
TOUWIL but the analysis w as relev an t to a m uc h wider range of phenomena Other
relev an t ph ysical pro cesses include the recom bination of electronhole pairs in semicon
ductors V AROCO and conceiv ably  c hemical reactions of the t yp e a  b


c  if the
rev erse reaction can b e reduced to a negligible rate b y v arying a ph ysical parameter suc h
as the temp erature  O V CZEL
F or a homogeneous system w e write the macroscopic densities of a and b as 
a
 t  and

b
 t  If the t w o densities are initially equal the particular kinetics implies that they will
remain so at later times Then the macroscopic rate equation assumes the simple form
d
a
 t 
dt
  k



a
 t   XI  
k

b eing the rate constan t According to Eq XI   
a
has the longtime b eha viour t
￿

W e consider a lattice p opulated b y a  and b particles sub ject to the irrev ersible reaction
a  b  c The sp ecies c pla ys no activ e r!ole in the ev olution of the system Initially the
lattice is o ccupied b y a uniform random distribution of equal n um b ers of a and b W e
de ne the densities 
a
 t  and 
b
 t  b y dividing the total n um b ers of a  and b particles b y

Cellular automata for transp ort with c hemical reactions
the n um b er of sites in the lattice When N
a
a particles and N
b
b particles meet at a site
they react with probabilit y P

N
a
N
b
 where P

  k

Eac h time a reaction tak es place
one particle of eac h reactan t is subtracted so that the o v erall n um b ers and densities of
a and b remain equal to eac h other W e consider three di eren t cases i the reactan ts
di use with the same di usion co e
cien t D  ii the reactan ts di use as b efore but one of
them is additionally sub ject to adv ection with a constan t v elo cit y V  and iii transp ort is
replaced b y a random redistribution of the reactan ts on the lattice The last case em ulates
the instan taneous mixing implicitl y assumed b y the macroscopic rate equation
The results of our sim ulations on a onedimensional lattice are sho wn in Fig XI  The
asymptotic time dep endence of 
a
is t
￿  
 t
￿  
and t
￿
 resp ectiv ely  in the three cases
listed ab o v e In the w ellmixed third case the result coincides naturally with that of the
macroscopic equation When the mixing is less fast as in the other t w o cases a slo w er
depletion mec hanism tak es o v er This happ ens b ecause as the bulk of particles disapp ear
in m utual annihilation the microscopic uctuations presen t already in the initial state
pla y an increasingly imp ortan t r!ole Let us rst consider the case of pure di usion The
time exp onen t     whic h is sp eci c to a onedimensional system for d dimensions one
nds  d  for d  and  for d   can b e understo o d as follo ws TOUWIL A t an y
giv en time t in the ev olution of an in nite system particle transp ort has not had the time
to dissolv e uctuations on a length scale longer than l
D

p
D t F or 
a
  
b
  


a part of the system of length l
D
con tains initially 

l
D
particles of eac h sp ecies and
the di erence b et w een the n um b ers of a  and b particles is of the order of
p


l
D
Th us
the surplus densit y of the sp ecies in excess is of the order of
p


l
D
l
D
 
 

 D t 
￿  

W e see that the deca y of initial uctuations results in a t
￿  
decline of the densit y whic h
is b ound to dominate o v er the t
￿
decline due to bulk depletion for long times F or
v ery long times the system breaks up in segregated regions of a  or b particles formed b y
the remaining excess particles after the bulk has annihilated the reaction then pro ceeds
only along the b oundaries of these regions It turns out that the relativ e drift b et w een
the reactan ts mixes them faster than di usion and the asymptotic decline of the densit y
follo ws a t
￿  
la w This can b e easily seen if the length scale
p
D t sw ept out b y di usion
is replaced b y the adv ectiv e scale V t KANRED It is cle ar that with the exc eption of
the unphysic al c ase of innitely fast tr ansp ort the c orr e ct longtime de cline of the density
c an only b e c aptur e d if micr osc opic uctuations ar e explicitly ac c ounte d for T o the exten t
that a mesoscopic mo del is called for it is ob viously desirable that it should b e able lik e
the mo del used here to accoun t for adv ection whic h w as seen to alter qualitativ ely the
b eha viour of the ph ysical system under consideration
XI
 A uto c atalytic r e actions
As w e sa w in Section XI   CA can b e used to mo del sp ecies of di eren t transp ort
prop erties as w ell as complex c hemical reactions W e are no w going to consider reaction
di usion systems in whic h the di erence b et w een the di usion co e
cien ts of the reactan ts
and a certain degree of complexit y in the reactions pla y an essen tial r!ole in the observ ed
phenomena The reactions of in terest here con tain auto catalytic steps in whic h the

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
Figure XI  a densit y for a onedimensional homogeneous system reacting via a  b  c 
as a function of time a only di usion solid line b di usionrelativ e adv ection
dashed line and c random redistribution dotted line
catalyst is one of the pro ducts  eg   X  Y   X  When auto catalytic systems
are main tained far from c hemical equilibrium  eg  b y suppressing rev erse reactions and
b y remo ving or adding reactan ts they can exp erience phase transitions to states with
remark able spatial or temp oral structures A classic example is the Belouso vZhab otinsky
reaction  ie  the o xidation of citric malonic or bromomalonic acid b y bromate in acidic
aqueous solution with v ariable v alence ions suc h as Ce

 Ce

or F e

 F e

 acting as
catalysts  BEL ZHA ZAIZHA whic h under appropriate conditions displa ys
oscillations in the concen trations of the o xidised and reduced forms of the catalyst as w ell
as p erio dic propagation of concen tration w a v es The emergence of ordered states far from
equilibrium is of extreme imp ortance in c hemistry and biology NICPRI MUR
There is a certain class of steady spatial structures whic h are sustained b y the com
bined e ect of reaction kinetics and di usion The earliest example w as prop osed b y Alan
T uring in  TUR A crucial condition for the formation of suc h patterns is that
t w o sp ecies ha v e di eren t di usion co e
cien ts Then it is p ossible that a uniform steady
state whic h w ould b e stable in the absence of di usion b ecomes unstable and go es o v er to
a stationary spatial concen tration pattern when the ratio of di usion co e
cien ts exceeds
a certain critical v alue The faster sp ecies the inhibitor  hinders b y c hemical action the
pro duction of the slo w er one the activator  and con nes the latter to restricted regions of

Cellular automata for transp ort with c hemical reactions
space giving rise to a pattern of inhomogeneous concen tration The inhibitor itself forms
a complem en tary pattern of high and lo w concen tration One can apply linear stabilit y
analysis to the reactiondi usion equations to deriv e the range of w a v elengths that b ecome
unstable for giv en v alues of the macroscopic parameters The emerging pattern is an ex
ample of sp atial symmetry br e aking since it de nes a new length scale that do es not dep end
on the geometrical prop erties of the system but rather on in trinsic prop erties describ ed
b y the di usion co e
cien ts and the reaction rates including feed and remo v al rates of
reactan ts as w ell as prescrib ed concen trations of external reactan ts Sustained concen
tration patterns w ere obtained recen tly with the c hloriteio didemalonic acid reaction in
op en reactors in whic h the reactan ts di use in to an inert gel from adjacen t reserv oirs of
constan t and uniform concen trations CASDUL KEPCAS QUYSWI
One of the simplest reaction systems that are sub ject to the T uring instabilit y under
appropriate conditions is
X
k
￿


k
￿ ￿
A  X  Y
k
￿
  X  B
k
￿
 Y  XI 
whic h w as conceiv ed b y J Sc hnac k en b erg in his searc h for simple systems admitting
p erio dic solutions SCH The concen trations of A and B are held xed for example b y
external supply or remo v al of reactan ts The system XI  can b e analysed analytically
and is kno wn to exhibit a T uring instabilit y for a wide range of parameter v alues MUR
If the rst reaction in XI  is also made irrev ersible the system ma y still oscillate but
it ma y also ev olv e to w ards in nite concen trations  eg  X   Y    with X

Y
nite the system X  A   X  Y   X  B  Y had b een used earlier b y E E Selk o v
to mo del the phosphofructokinase catalysed pro duction of ADP from an A TP substrate
as a p ossible source of glycolytic oscillations  SEL Eq I I Replacing the last reaction
in Selk o vs sc heme with the pair of reactions B  X  Y  D  E  X  w e obtain
the Brusselator one of the earliest systems emplo y ed to illustrate the p ossibilit y of
symmetry breaking instabilities  PRILEF
Sim ulations of the ab o v e reaction systems ha v e b een p erformed with CA mo dels
KAPLA W KAPLA W DR OFRA KAR KARBLA The pattern in
Fig XI  w as obtained from a random macroscopically homogeneous initial state after
a few thousand iterations of a CA sim ulation of the Sc hnac hen b erg mo del on a t w o
dimensional lattice with p erio dic b oundary conditions KAR KARBLA Giv en the
rest of the macroscopic parameters linear stabilit y analysis predicts in this case that
the T uring instabilit y o ccurs if the ratio of di usion co e
cien ts d  D
Y
D
X
exceeds
the critical v alue d
cr it
    The sim ulation w as p erformed with d   for whic h
the macroscopic analysis predicts a windo w of unstable w a v elengths consisten t with the
observ ed pattern
Whereas reactiondi usion equations describ e the ev olution of the macroscopically a v
eraged particle densities CA sim ulations incorp orate microscopic uctuations at all stages
of the ev olution Fluctuations are exp ected to ha v e an enhanced inuence on the dynam
ics of the system near bifur c ation p oints  suc h as near the onset of the T uring instabilit y in

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
Figure XI  Spatial concen tration pattern obtained from t w odimensional sim ulation of
the Sc hnac k en b erg mo del
parameter space The computational e
ciency of CA mo dels with an exclusion principle
allo ws their application to detailed in v estigations of the r!ole of uctuations in particular
c hemical systems Th us in t w odimensional sim ulations of the Sc hnac k en b erg mo del as
w ell as of the Selk o v mo del transien t patterns app ear also b elo w the critical v alue of the
di usion co e
cien t ratio KAPLA W KAPLA W DR OFRA In that region pat
tern formation is not p ossible according to the macroscopic reactiondi usion equations
Apparen tly  b elo w the critical ratio the homogeneous state is destabilised b y densit y
uctuations whic h e ectiv ely widen the range of unstable w a v en um b ers F or parameter
v alues for whic h the macroscopic equations predict a T uring instabilit y  the patterns ob
tained b y CA sim ulation generally displa y a certain roughness and are less regular than
the ones deriv ed from the reactiondi usion equations Whereas the latter patterns are
strictly stationary after a certain time the former are con tin ually p erturb ed b y uctu
ations whic h o ccasionally destabilise a sustained structure and induce reorganisation of
the pattern Fig XI 
XI

 R e actions with miner al surfac es
W e no w mo v e to c hemical reactions taking place in heterogeneous systems where di eren t
phases are in con tact with eac h other Th us solutes often react with the surfaces of the
w aterconducting c hannels along whic h they are transp orted The reaction ma y b e in the
form of mineral pr e cipitation and dissolution  whic h leads to c hanges in ro c k prop erties
 eg  the p orosit y with a direct impact on the transp ort of the reacting solutes Alter
nativ ely  the solid surface ma y act as a c atalyst b y adsorbing one or more reactan ts from
the liquid or gaseous phase and mo difying them so that they react more readily In b oth

Cellular automata for transp ort with c hemical reactions
Figure XI  T uring pattern reorganisation from four to three strip es caused b y densit y
uctuations The densit y pattern sho wn here see grey scale is unpublished w ork b y D
Dab and w as obtained with the Magin u mo del MA G MA G sim ulated on a  
 lattice with parameter v alues k     c     s   D
x
  and D
y
  as
describ ed in Refs D ABBOO D AB

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
cases inhomogeneities in the spatial distribution of reactan ts pla y a crucial r!ole in the
ev olution of the system while mo ving b oundaries ma y further complicate the situation
Due to its abilit y to treat spatially and temp orally v arying parameters with no additional
e ort the CA approac h is particularly appropriate for mo delling systems of the kind de
scrib ed here W e shall illustrate this b y means of t w o simple CA mo dels of heterogeneous
pro cesses
W e rst consider the precipitationdissolution reaction
c  s 
K
￿


K
￿
a  b  XI 
where a and b are aqueous sp ecies c is a mineral and K

 K

are the rate constan ts
F or constan t and uniform p orosit y   the standard reactiontransp ort equations for this
system are

 C
a
 x t 
 t
  V
a
 C
a
 x t 
 x
 D
a


C
a
 x t 
 x

 

K

 

K

C
a
 x t  C
b
 x t 

 C
c
 x t 
 t
  

K

 

K

C
a
 x t  C
b
 x t 

 XI 
where C
a
and C
b
are the concen trations of the solutes in mol p er m

of w ater C
c
is the
concen tration of solid in mol p er m

of ro c k and  is a constan t with
 
  if C
c
 x t   or C
a
 x t  C
b
 x t  K
  otherwise
 XI 
K  K



K

b eing the solubility pr o duct The equation for C
b
is analogous to the one
for C
a

In the sim ulation the a  and b particles are mobile whereas the c particles are sta
tionary The reaction rule dep ends on the lo cal n um b ers of a and b as in the case of
the reaction a  b  c considered ab o v e but it additionally dep ends on the n um b ers of
reactan ts through the conditions for the onset of precipitation Precipitation ma y namely
o ccur only if there is solid already presen t or if the solution is saturated the de nition
of  ab o v e is the mathematical equiv alen t of this statemen t This leads to the follo wing
mesoscopic rule If one or more c particles are presen t at a site one of them go es in to
solution in the form of an a  and a b particle with probabilit y P

If N
a
a particles and
N
b
b particles are presen t then they ma y precipitate and form a c particle only if at least
one c particle is already presen t or if the pro duct of the densities of a and b calculated b y
a v eraging o v er an ensem ble of macroscopically iden tical systems is greater than P

P


If one or b oth of these conditions are ful lled precipitation is allo w ed to tak e place with
probabilit y P

N
a
N
b
P article densities are related to the ph ysical quan tities in Eq XI 
b y C
a
  
a
 C
b
  
b
and C
c
  
c
 where the scaling factor  can b e xed b y the

Cellular automata for transp ort with c hemical reactions
ph ysical data and the initial conditions of the sim ulation The reaction parameters are
related to the rate constan ts b y P

  K

  and P

   K


The CA mo del describ ed ab o v e is no w applied to a system with onedimensional
geometry with translational in v ariance assumed in the other t w o dimensions W e
lo ok at the dissolution of a blo c k of p orous ro c k b y an inert uid o wing through it
KAR KAR W e consider a piece of ro c k of length of  cm with a uniform
p orosit y of    The left half of the ro c k is o ccupied at time t   b y a sol
uble mineral with a densit y of   gcm

and a molecular w eigh t of  The rate
constan ts for dissolution and precipitation are K

     
￿ 
mol  m
￿ 
 s
￿
and
K

    m

 mol
￿
 s
￿
 resp ectiv ely they are compatible with a solubilit y pro d
uct of K  
￿   
mol

 m
￿

 
￿   
mol  l 

The ph ysical and c hemical parameters
are c hosen to b e consisten t with calcite CaCO

 as a t ypical mineral in the geological
systems of in terest The v alue of K

in particular follo ws from the v alue of the dissolution
rate p er unit surface of mineral K
s
     
￿
mol  cm
￿ 
 s
￿
     
￿
mol  m
￿ 
 s
￿
obtained in an in v estigation of the pH dep endence of the dissolution rate of v arious car
b onates at 
￿
C using a uidized b ed reactor CHOGAR T o obtain the dissolution
rate p er unit v olume of mineral w e ha v e m ultiplied K
s
b y a surface area to v olume ratio
of  cm
￿
 whic h corresp onds to grains of size    cm It is clear that K

is in v ersely
prop ortional to the grain size and a go o d kno wledge of the latter is necessary for mo d
elling dissolution kinetics in particular situations The v alue of K

is then calculated from
K

 K



K and is consisten t within  with the v alue of K
s
found indep enden tly
in Ref CHOGAR The dissolv ed sp ecies Ca

and CO
 ￿

mo v e with an adv ection
v elo cit y of V  
￿
ms     mda y and a di usion co e
cien t of D  
￿
m

 s F or
the sim ulation w e c ho ose   
￿
m      s and      mol  m


Fig XI  sho ws the pro les of C
c
and C
a
at time t  

s    h
In the sim ulation w e obtain net precipitation of solid do wnstream from the original solid
edge Due to uctuations in the solute densities the pro duct of these densities o ccasionally
exceeds the saturation threshold and precipitation ensues This e ect is absen t from the
macroscopic description where the edge remains sharp The net precipitation of solid in
the righ t half of the lattice is accompanied b y a reduction in solute concen tration relativ e
to the prediction of the macroscopic equations Of course the amoun t of solid that is
found do wnstream of the edge is tin y on the scale of the total amoun t of solid but the
consequences can b e more serious for the concen tration of solutes In Fig XI  of  
million solid particles initially presen t ab out  ha v e b een remo v ed from the solid blo c k
and ab out  ha v e precipitated in the righ t half of the lattice This su
ces ho w ev er
to pro duce a  suppression in the solute concen tration as a comparison of the dotted
and dashed curv es sho ws Sim ulations sho w that the do wnstream precipitation of solid
and the suppression of solute concen tration diminish as w e reduce the magnitude of the
statistical uctuations b y increasing the size of the ensem ble o v er whic h w e a v erage Since
the particles of our sim ulations do not mo del individual molecules of the actual system
the in terpretation of statistical uctuations and their consequences is not immedi ately
ob vious W e can assume a phenomenological p oin t of view allo wing for microscopic
pro cesses not explicitly accoun ted for to con tribute to the uctuations presen t in our

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
Figure XI  Pro le of solute concen tration horizon tal dotted line and solid line running
on the a v erage parallel to it reaction transp ort equations and CA resp ectiv ely the
shortdashed line is a smo othed v ersion of the solid line axis lab elled on the righ t and
the concen tration of solid p er m

of ro c k v ertical longdashed line and solid line running
partly parallel to it reaction transp ort equations and CA resp ectiv ely axis lab elled on
the left
mo del Then one has to normalise the amplitudes of the uctuations b y comparing the size
of the resulting e ects with real systems The impact of the observ ed e ect will therefore
dep end on the problem at hand F or times long enough for precipitationdissolution to
cause c hanges in the p orosit y  it is lik ely to inuence those c hanges to a degree ho w ev er
that remains to b e in v estigated
Remark ably  the e ect describ ed ab o v e o ccurs only in the presence of adv ection and is
absen t when di usion is the only transp ort mec hanism The dep endence of the e ect on
the ph ysical parameters of the problem w as studied in Ref KAR It w as found there
that the suppression of solute concen tration do wnstream from the solid edge p ersists up
to the longest times t for whic h the sim ulation can b e p erformed without exceeding a w eek
of CPU time If one v aries the relativ e time scales for transp ort esp ecially adv ection and
reactions then the e ect of uctuations is w ashed out when the reaction kinetics is not
fast enough to k eep up with transp ort The e ect app ears clearly only when transp ort is
slo w relativ e to the kinetics These observ ations underline the imp ortan t role of adv ection
in the in terpla y b et w een transp ort and c hemistry The abilit y of the mesoscopic mo del
used here to accoun t for adv ection clearly op ens a wide eld of in teresting phenomena to
b e in v estigated

Cellular automata for transp ort with c hemical reactions
As another example of a heterogeneous pro cess w e tak e up the catalytic o xidation of
COg on platin um Pt surfaces This pro cess w as in v estigated in Ref WUKAP
from the p oin t of view of the in terpla y b et w een sorption kinetics reactions and surface
phase transformations The basic pro cess consists of the follo wing steps gaseous CO and
O

adsorb on sites of the Ptcr surface and O

molecules disso ciate to O atoms whic h
o xidise adsorb ed CO to form CO

The adsorption rates of COg and O

g dep end on
the a v ailabilit y of free surface sites Adsorb ed CO and O ma y of course desorb b efore they
react but once they react CO

desorbs rapidly If 

and 

are the densities of adsorb ed
CO and O resp ectiv ely note that  and  are here sp ecies lab els the summation o v er
v elo cit y states ha ving already b een p erformed and 

the constan t densit y of adsorption
sites free  o ccupied one can write do wn sc hematically the rate equations
d

dt
 k

p
CO
 

 

 

  k
￿


 k





d

dt
 k

p
O
￿
 

 

 



 k
￿ 



 k





 XI 
where k

 k
￿
 and k

 k
￿ 
 are the adsorption desorption rate constan ts of CO and O


p
CO
and p
O
￿
are the resp ectiv e partial pressures and k

is the rate constan t for the surface
reaction CO  O  CO

The rst term on the righ t hand side of b oth rate equations
mak es explicit the dep endence of adsorption rates on the densit y of free sites
If the COO reaction tak es place on the Pt crystallographic plane then the surface
ma y undergo phase transformations dep ending on the fraction of sites co v ered b y CO
Th us if the surface is in the hexagonally ordered phase hex and the CO co v erage
exceeds  a transition o ccurs to the    phase c haracterised b y a square lattice
of surface atoms on whic h CO molecules o ccup y alternate sites When the CO co v erage
falls b elo w  the rev erse surface transformation tak es place Since a reorganisation
of the surface a ects the sorption rates and hence the o xidation kinetics the coupling
b et w een the di eren t pro cesses can induce oscillations in the surface concen tration of the
reactan ts and the fraction of sites in eac h phase These oscillations pro ceed b y means of
the app earance and gro wth of domains of one phase in the midst of the other coupled
with the di usion of adsorb ed particles across the domains Suc h oscillatory b eha viour
has b een observ ed exp erimen tally NORBIN
The mo del of Ref WUKAP applies an exclusion principle as describ ed in Sec
tion XI   making sure ho w ev er that the reactiv e transitions at the mesoscopic lev el
are only those implied b y the rate la w Eq XI  Extending the notation used
earlier for onesp ecies systems w e in tro duce the probabilit y P  ￿￿  for the transition


CO 

O 


CO 


O with ￿  



 and ￿  





 P  ￿￿  is written as
P  ￿￿   p

 ￿  

￿

￿



￿

￿
 p

 ￿  

￿

￿
￿


￿

￿
 p

 ￿  

￿

￿


￿

￿

 p

 ￿  

￿

￿


￿

￿
￿ 
 p

 ￿  

￿

￿
￿


￿

￿
￿
 XI 
In Eq XI  the p
i
 ￿  i         are the probabilities for CO adsorption and

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
desorption O adsorption and desorption and CO o xidation resp ectiv ely
p

 ￿  

"
k

 m 



  if



m 
  if



 m 
p

 ￿  
"
k
￿



p

 ￿  



"
k

 m 



  m   



  if



m   
  if



 m   
XI 
p

 ￿  


"
k
￿ 




  
p

 ￿  
"
k






If one c ho oses the parameters
"
k
i
according to
"
k

 k


"
k

 k

p
CO

"
k
￿
 k
￿
"
k


m
m  
k

p
O
￿

"
k
￿ 

m
m  
k
￿ 
 XI 
one obtains up on a v eraging and taking the con tin uum limit rate equations v ery similar
but not iden tical to Eq XI  There is a clear connection b et w een the di eren t terms
in the rate equations XI  and the corresp onding expressions for the p
i
 ￿  whic h are
reminiscen t of Eq XI 

Eqs XI  XI  and XI  de ne the mesoscopic
reaction dynamics
The phase transformations of the surface can also b e mo delled b y a CA rule Th us
one ascrib es eac h site to one of t w o phases and allo ws it to mak e a transition to the
other phase with a certain probabilit y  if the lo cal CO co v erage crosses a certain critical
v alue The CO co v erage at a site is de ned b y coun ting the n um b er of CO molecules
in a neigh b ourho o d consisting of the site itself and its nearest neigh b ours The phase
transition couples bac k to the adsorption and reaction kinetics since the reaction rates
k
i
ha v e di eren t v alues for eac h phase
Fig XI sho ws the oscillations in the CO surface densit y obtained for a particular
c hoice of the ph ysico c hemical parameters WUKAP Initially most of the sites are
in the hex phase whic h adsorbs O m uc h more w eakly than the    phase As a
result the CO densit y increases Fig XI a the critical co v erage for transition to the
   phase is lo cally exceeded and domains of that phase b egin to app ear and spread
Fig XI b These domains adsorb O at a high rate and cause o xidativ e depletion of
the CO molecules that are originally presen t in their in terior as w ell as of those that
arriv e there b y di usion Th us the new phase wins gradually o v er and the densit y of
adsorb ed CO diminishes Fig XI c When ho w ev er the CO co v erage falls lo cally
b elo w a certain critical v alue domains of the original phase b egin to app ear and spread
Figs XI de and the CO densit y increases again Fig XI f  Th us the relativ ely
￿
The
i
in Eq
XI  cannot exceed ho w ev er the n um b er of p ossible v elo cit y states and ob ey
a binomial distribution at equilibrium whereas the N

in Eq
XI  are unrestricted and ha v e a
P oisson equilibrium distribution

Cellular automata for transp ort with c hemical reactions
Figure XI  Oscillations in CO surface densit y obtained from t w odimensional sim ulation
of CO catalytic o xidation on the Pt plane reprin ted with p ermission from Ref
WUKAP Fig  The frames al are for times t   at equally spaced
time in terv als of  time steps The sim ulation w as carried out on a    hexagonal
lattice with p erio dic b oundary conditions F or the v alues of system parameters used in
the sim ulation the reader is referred to WUKAP

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
simple CA w e are considering here is capable of repro ducing the exp erimen tall y observ ed
oscillatory b eha viour Moreo v er as the Pt surface oscillates in space and time b et w een
the t w o phases whic h are c haracterized b y di eren t reaction rate constan ts the presen t
application illustrates the capacit y of CA mo dels to describ e systems with spatially and
temp orally v arying parameters
It should b e remark ed that b oth mo dels of heterogeneous pro cesses that w e ha v e pre
sen ted are relativ ely simple They already sho w the p oten tial of the mesoscopic metho d
to describ e in teresting ph ysical b eha viour but they should rather b e tak en as the basis
for more elab orate and realistic mo dels
XI  Conclusion
The mo dels of reactiontransp ort phenomena discussed in this c hapter o ccup y a p osition
in termediate b et w een the standard macroscopic approac h of partial di eren tial equations
PDEs and the fully microscopic approac h of molecular dynamics W e ha v e indicated
with sev eral examples that there is a wide range of ph ysical phenomena for whic h the
mesoscopic p oin t of view pro vides the necessary amoun t of detail that PDEs cannot
a ord and optimises the use of resources b y neglecting inessen tial details of the mi
croscopic dynamics Cellular automata CA describ e transp ort and c hemical reactions
at an elemen tary lev el k eeping ho w ev er only those asp ects of the microscopic dynam
ics whic h are imp ortan t for the problem at hand Due to their simplicit y  CA mo dels
are conceptually transparen t and can b e transcrib ed in to simple computer co de This
is a signi can t adv an tage when it comes to mo delling complex geometries and b oundary
conditions or the kinetics of arbitrary c hemical reactions CA algorithms iterate in te
ger v ariables and are therefore free of roundo errors and asso ciated instabilities This
remains true for the sto c hastic CA mo dels considered here since real n um b ers app ear
only as b ounded probabilities By b eing faithful to the lo cal and parallel c haracter of
the ph ysical pro cesses they mo del cellular automata are suitable for implem en tati on on
massiv ely parallel computers It should b e emphasised that it w as largely the dramatic
progress in computer hardw are and arc hitecture in recen t y ears that made CA mo delling
of natural phenomena tec hnically p ossible
The mo dels w e ha v e discussed sim ulate the transp ort of particles as a random w alk on
a regular spatial lattice When the particles meet at a lattice site they react c hemical ly
and are transformed according to a lo cal probabilistic rule Radioactiv e deca y  sorption on
conduit w alls and the dissolution of solids are also mo delled b y simple probabilistic rules
W e ha v e presen ted a mo del whic h places no restriction on the n um b er of particles that
can b e found at a site as w ell as one deriv ed from lattice gas automata for uid dynamics
whic h allo ws only a small set of p ossible v elo cit y states to eac h particle The particles in
b oth mo dels are mathematical abstractions of the actual molecules of the ph ysical system
The n um b er of particles is large enough to mak e statistical concepts meaningful but is
still man y orders of magnitude smaller than the real n um b er of molecules in v olv ed Meso
scopic mo delling relies crucially on the assumption that this is legitimate pro vided that
the macroscopic b eha viour of in terest do es not dep end on the details of the microscopic

Conclusion
dynamics but follo ws from general prop erties of the latter Of the t w o mo dels the former
has the merit of b eing able to accoun t for adv ection b y imparting a directional bias to
the random w alk and of mo delling c hemical reactions in a more transparen t and exible
w a y The lattice gas mo del on the other hand o ers substan tial computational adv an
tages thanks to the limited n um b er of allo w ed particle con gurations at a site it can b e
programmed at bit lev el th us making optimal use of computer memory and pro cessing
p o w er and up dating can b e optimised b y using lo okup tables whic h con tain stored all
elemen tary ev en ts that can tak e place at an y site The abilit y to describ e adv ection and
a relativ e ease in the inclusion of v arious c hemic al reactions are indisp ensable for mo d
elling coupled transp ort and reaction phenomena in geological media but computational
p erformance is an imp ortan t criterion of the usefulness of a mo del as a practical to ol
Despite the recen t adv ances in computer resources computational e
ciency still has to
b e rec k oned with esp ecially for problems of the complexit y w e wish to address It is
therefore desirable to com bine the adv an tages of the simple and more exible mo del with
the exp erience accum ulated with n umerous sp eci c applications of lattice gas mo dels
In the sim ulations discussed here the transp ort of solutes w as mo delled assuming that
their adv ection follo w ed a giv en v elo cit y eld If the solute concen tration is small enough
the v elo cit y eld coincides with that of pure w ater o wing under the giv en mec hanical
conditions The dynamics of w ater o w is go v erned b y the Na vierStok es equations
whic h for a p orous medium go o v er to Darcys la w Although the applications presen ted
here in v olv ed only uniform o ws the transp ort mo del without an exclusion principle w as
form ulated for an arbitrary o w eld In general the problem of solute transp ort can
b e solv ed b y a t w ostep pro cedure where one rst solv es the appropriate h ydro dynamic
equations and then uses the resulting v elo cit y eld as an input parameter to the mo del
describing solute transp ort This approac h is conceptually similar to the Lagrangian
approac h for mo delling tracer transp ort in geological formations  eg  Ref PRINA Y
Chemical reactions can b e added and the full algorithm iterated in time If the p orosit y
and the p ermeabilit y of a p orous medium c hange substan tially as a result of c hemical
reactions the inuence of these c hanges on w ater o w ha v e to b e explicitly considered
T o quan tify computational e
ciency in CA sim ulations one usually in tro duces the
concept of site up dates p er second u p s  The densit y of particles has an inuence on
up dating sp eed but except for sim ulations with v ery di eren t particle densities it is
useful to compare the p erformance of v arious algorithms in terms of u p s W e consider
rst the mo del with no a priori limitation on the n um b er of particles p er site The
sim ulation of the reaction a  b  c that pro duced the solid curv e in Fig XI  had an
a v erage p erformance of     


u p s and required ab out  CPU hours on a massiv ely
parallel Connection Mac hine CM The pattern sho wn in Fig XI  to ok ab out 
CPU hours at a rate of   

u p s on a V AX in scalar mo de the sim ulation of
the reaction a  b  c on the same computer w as roughly a factor of  less e
cien t than
on CM The sim ulations can b e made faster b y utilising a lo wlev el instruction set
a v ailable on the Connection Mac hine but this w ould impair the transparency of the co de
The p oten tial of applicationsp eci c optimisation b ecomes clear if w e compare the ab o v e
rates with those ac hiev ed with sp ecialised lattice gas mo dels Th us the reaction a  b  c

Cellular Automaton Mo dels of Reaction T ransp ort Pro cesses
w as sim ulated with a lattice gas mo del using a lo okup table on a CRA YYMP at a rate
of   

u p s  while a four times higher rate w as attained on a CMa t ypically slo w er
than CM b y a factor of  CHODR O It should b e men tioned ho w ev er that
the co de without limitation on particle n um b er w as applied unmo di ed to the case with
adv ection the dashed curv e in Fig XI  at no cost in p erformance whereas adv ection
lies b ey ond the scop e of the lattice gas mo del
The relativ e merits of di eren t CA mo dels and standard n umerical metho ds for the
solution of the corresp onding PDEs dep ends on the sp eci c applications Th us with
presen tda y computers pattern formation in auto catalytic reaction systems can b e stud
ied in general more e
cien tly with con v en tional metho ds if ho w ev er one is in terested
in the r!ole of uctuations the mesoscopic approac h pro vides a p o w erful alternativ e It is
therefore imp ortan t to seek areas of applications where eac h metho d maximi ses its adv an
tages CA mo dels are v ery w ell suited for describing transp ort with inhomogeneous and
timev arying parameters They also naturally describ e the c hemical kinetics of systems
whic h are not in c hemical equilibrium According to the discussion in the In tro duction
the guaran teed stabilit y of CA algorithms is a v aluable asset when nonlinear pro cesses
come in to pla y Finally  the abilit y of cellular automata to accommo date arbitrary b ound
ary conditions including mo ving b oundaries has not b een su
cien tly exploited so far
Among p oten tial applications suitable for CA mo delling of particular theoretical and
practical in terest are the in teraction of solutes with the solid surfaces of p orous and frac
tured media and the coupling of p orosit y c hanges induced b y precipitationdissolution
with solute transp ort CA mo dels app ear to b e particularly called for when for example
the ph ysical pro cesses of in terest sho w strong v ariations in space SARBRA or when
c hemical reactions are slo w compared to transp ort W OL It is in suc h cases where
con v en tional metho ds are confron ted with serious conceptual and tec hnical problems that
cellular automata can b est utilise their p oten tial for incorp orating microscopic detail in
the macroscopic description of the ph ysical w orld
XI  Ac kno wledgemen ts
A large fraction of the authors original w ork presen ted in this article w as p erformed in
collab oration with B Blankleider no w at Flinders Univ ersit y Australia The author
has also b ene ted a great deal from discussions with U Berner J Hadermann F Neall
PSI Switzerland B Chopard M Droz Univ ersit#e de Gen $ ev e J P Bo on and D Dab
Univ ersit#e Libre de Bruxelles D Dab kindly pro vided one unpublished gure P artial
nancial supp ort b y NA GRA is gratefully ac kno wledged