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

## Comments 0

Log in to post a comment