- Email: [email protected]
Diffusion of reacting pollutants
Diffusion of reacting pollutants K. H. Miiller Commission of the European Communities, Joint Research Centre, Ispra Establishment, I-21 020 Ispra (Va), Italy (Received November 19 79; revised March 1980)
The Boltzmann equation including reactive scattering terms is the correct description of the combined transport-reaction process. A suitable parameterization transforms it into a set of differential equations governing the macroscopic mass balance. An operator splitting integration method composed of the stochastic formulation of chemical kinetics and the Green’s representation of the transport proves to be very efficient.
Problem
description
Thousands of organic and inorganic substances are continuously emitted into the atmosphere by natural and man-made sources. They stay there for some time. Unstable molecules decay, and some of their fragments join other components. The photochemical transformations are of particular importance. Most of the gaseous species released become ‘digested’ and therefore play only a minor role in atmospheric pollution. Main sinks are wash-out, biological interactions and decomposition after absorption of radiation. Unfortunately, there are also chemical reactions with pollutants, surface materials and ‘natural’ gases leading to noxious reaction products. Pollutant transport and chemistry in the atmosphere is a domain where one assumes that the basic principles are known. But theoretical treatment is possible only in the simplest cases. It is obvious that a decision maker and one who is engaged in solving practical problems cannot wait until a theoretical solution of the problem is published. Reliable predictions in air hygienics can only be given by a simulation of the physical and chemical processes governing pollutant dispersion and transformation. Their full complexity, however, cannot be simulated experimentally on a laboratory scale. Here, a numerical simulation, i.e. modelling, using an electronic computer, is inevitable. The aim of atmospheric modelling consists both of obtaining insight into the interactions between radiation, chemistry and dynamics and of predicting possible environmental modifications caused by intentional or unintentional human interventions in nature. Since a detailed simultaneous modelling of pollutant chemistry and dispersion proves rather troublesome, one tends to purchase a detailed description of one phenomenon with a rough approximation of another. The extent to which this is justified has to be checked carefully from case to case. Only a simulation containing very detailed mechanisms can decide whether one phenomenon is more important than another. The bottleneck here is the scarcity of the
268
Appl.
Math.
Modelling,
1980,
Vol 4, August
experimental data necessary both to uniquely characterize a given situation and to verify the adequacy of the model. First, one checks the usefulness and completeness of the available experimental data by applying simply structured models (sub-models) and, after that, if necessary, asks for additional measurements. A more complete data set will then represent a more reliable basis suitable for updating and generalizing the previous model version, etc. Such a procedure is somewhat like that of successive approximations. Finally, one will be able both to distinguish between relevant and unimportant events and to extrapolate to situations which are difficult or impossible to reach experimentally. The vitality of a model lies in the availability of an extensive parameterization which allows a wide span of response characteristics of the physical system in question to be studied. Basically, it is not difficult to construct a model and to fix its parameters by fitting to measured data, if only a few species are to be treated; it becomes difficult, however, if several species and all possible atmospheric conditions occurring in a region, have to be taken into account. Through a sensitivity analysis which varies the characteristic parameters and relates these variations to their consequences, a reaction of minor significance can be traced and subsequently eliminated from the model. Input parameters are: dimensions of reaction-space and boundary conditions, transport coefficients, reaction rates, photodissociation coefficients, etc. It is clear that a specific parameterization of the model input determines the overall structure of the model itself. One could, for instance, say that the realism of a model is seen in the accuracy of the characteristic response time of the individual processes; and one could take this as a guide-line for a ‘realistic’ description of the interplay of the different processes under consideration. Certainly, with rather simple models one can already get a good survey of the whole situation and the trend of an evolution. The crucial thing is to find a suitable parameterization of the characteristic processes and of their causes and consequences. Also, one has to ensure that the fitting approach, which mainly hinges on variations of 0307-904X/80/040268-1 0 1980 IPC Business
l/$02.00 Press
Diffusion
parameters, gives reasonable orders of magnitude for each process under study. Air quality models will usually be judged by a comparison of their output with measured concentrations of atmospheric trace species. The corresponding validation experiments mostly rely on a small number of reactants. It must be remembered that a model for only a few reactants may not be unique, e.g. with respect to smog formation, since several different mechanisms can be constructed and adapted to given reaction-chamber experiments just by a suitable choice of parameters.
Mathematical
formulation
The most general principles in physics, such as the conservation of mass, momentum and energy, form its basis. The mathematical expresssion for a conservation principle is a continuity equation. Such equations must be established for each chemical species with a space-time varying particle density. The resulting equations, coupled to each other, due to all the interactions occurring in the atmosphere, have to be solved simultaneously to obtain the real space-time behaviour of the concentrations. It is clear that this set of equations, a model of the chemical kinetics transport, can only be solved numerically, and that only such a numerical simulation of the true interplay of the different atmospheric phenomena enables us to study it and to infer trends and consequences from the actual situation. The kinetic theory of gases’ represents a uniform basis for models formulating the full interaction between radiative, photochemical and aerodynamic processes, and is, therefore, a suitable starting point for the mathematics. Continuous functions of space and time have to be assigned to the bulk matter properties, such as mass density, energy density, velocity field, etc. Only when it is possible to define a volume element of microscopic dimensions which still contains a sufficient number of molecules, can one speak of local density, local velocity, etc., in the sense of average values. Simultaneously, a time interval has to be defined, macroscopically small but large enough so that a sufficient number of molecular events, such as collisions, fluctuations, etc., can be processed to reasonable average values. The main task of the molecular theory is to show that such mean values can reasonably be defined and that they satisfy the continuum equations. One is interested in the particle distribution function fi(?, t, 6) counting the number of molecules of species Si within the phase-space volume element. The functions fi are considered to be continuous in all their variables. In this sense, the normalization holds:
the acceleration field acting on the molecules. In case of simultaneously occurring collisions, however:
#O
holds. Obviously, one has to know the formal structure of the right-hand ‘collision’ term which is a measure for the rate of change offi due to collisions. The basic assumption of the kinetic theory is: the gas is dilute; i.e. the mean distance between the gas molecules is much larger than the range of their intermolecular forces. Occasionally, one molecule hits another molecule; and, much more infrequently, three molecules collide simultaneously. Within a mathematical analysis, this fact can be expressed by the statement that there are only two-body or binary collisions. Having constructed a corresponding binary theory, one will try to generalize it to three-body collisions, etc. A gas is said to be in non-equilibrium, when its distribution function fi deviates from the Maxwell-Boltzmann distribution. Causes are heterogeneous temperature-, density- or velocity-fields. The system of molecules tries to equalize these irregularities by transport of energy, momentum and mass. The corresponding transport mechanism is molecular collision. One usually assumes, as a working hypothesis, that the molecules fly freely between two subsequent collisions, i.e. a’= 0, and that velocity 5 and position i of a molecule are uncorrelated (‘molecular chaos’). The number of pairs of molecules in d3r of two different substances, Si and Sk, whose velocities are 6i and Gx-,is given by: [fi(;, t, 4) d3r d3vi]. [fk(f, t, &) d3r d3vk] The molecules of substance Si will be ‘attacked’ by the current: fk(?, t, Gk) d3V, IGj - 6, I of Sk molecules. The number of collisions occurring during the time interval dt is: fi(Z’, t, 6j) dtojk($,
iTjk)fk(?, t, Gk) d3Vkl$ - 6kl
where oik denotes the ‘cross section’ characterizing the collision mechanism. During such collisions, chemical transformations may occur. The corresponding changes in particle densities are assumed to be proportional to the number of collisions. The proportionality factors, p and 4, are called (positive or negative) production rates. The velocity offi-changes, due to collisions, is therefore given by: = -hG,
t, G) d3vk s
cPik(j, k
jk)
X fk(?,t,Zk).IG-ZkI + d3Vj
The aim of the kinetic theory is to find the distribution functionsfi corresponding to given molecular interactions. The density fi changes in time, since molecules enter and leave the volume element dV= d3r d3v all the time. This transfer can be split into a ‘streaming’ and a ‘scattering’, interpreting the particle motion as the consequences of collisions between molecules. For collision-free motion, the density function h satisfies the ‘streaming’ equation: fi=O
(2)
collision
i at 1 collision
= number of particles in (d3r dt)
K. H. Miller
a’ denotes
afi
fi(;, t, 6)d371 =Ni(g, t)
of reacting pollutants:
s
X c
qijk(G> 4, ;k)
o.ik(j,
6,)
d3Vk .i
ajk(Gj, Gk)
j,k Xfi(‘:
(1)
t,
zj)fk(i, t,;,)I'$ - Gkl
(3)
The balance equation resulting from equations (2) and (3) and including reactive scattering is a correct dynamical description of a combined transport-transformation process in a dilute gas; it is equivalent to the uncorrelated binary collision theory of chemical kinetics.
Appl.
Math.
Modelling,
1980,
Vol 4, August
269
Diffusion of reacting pollutants: K. H. Miiller
The continuity equations of the mass-, momentum- and energy-transport can be derived from this Boltzmann equation by applying a weighted averaging of the molecular velocity 5. Our interest lies in deriving a compact model for the pollutant-mass transport including the simultaneously occurring chemical transformations in an atmosphere characterized by a given velocity- and temperature-field. Using the already defined density Nj and the particle current: ii = r Gfi d3V J
and the following mean values: d3Vk =@ik(q
pikujk.I;-6k(fk
fk d3Vk =&jk(qNk I
s qijkojk
’
1I$ - ?k 1fi fk d3Vi d3Vk
= ,!&@) s the Boltzmann
fi d3vi fk d3Vk = &jk($ NiNk I
equation
can be rewritten
1
CYijNk h=
k
(4)
as follows:
C PijkNjNk
(5)
i>k
An integration of this equation over the whole velocity space and the introduction of:
.i
oik fi d3v = aik
wi
fid3V = UikNi
leads to the parameterized i3Nj z t div2i = -LiNi
Pijk
.r
d3V = bijk
form (with a’ = 0):
dNi dt=
k
Pi = C biikN;Nk
(10)
k = ATB.exp(-C/T)
‘
denote the ‘loss’ and ‘production’ of Si respectively. Before trying to parameterize the current ii one introduces the mean velocity ii of all molecules actually present in the volume element (d3r dt); i.e.: ii = $fi dav = GNi - (; - ;) fi d3v z ;Ni - fi s I The space-time dependent function ?j expresses a sort of mean-dispersion behaviour of substance Si. According to classical macroscopic diffusion theory one could perform a parameterization by a linear ansatz of the type: ?i = lK(ZNi + grad Ni)
(8)
The tensor 1K serves to quantify the dispersion behaviour. One tries to include the influence of a given temperature field Ton the chemical kinetics transport by the ansatz*: c’ = V (1nT) + constant Inserting the above I?i into equation (6) and putting V . li = 0, a nonlinear diffusion equation results: aNi ~+~.VNi-V.(IKVNj)+(Li-V.lK~)N,=Pi (9)
270
+ Pi
valid for homogeneous mixtures.3 Experiments of this kind show a strong temperature dependence of the reaction rates, &k and bijk. One, therefore, tries to find at least their approximate T-structure. The temperature behaviour of most of the reactions relevant to air pollution follows Arrhenius’ formula:
processes.
Li = CfZikNk
’
~LiNi
t Pi
of the macroscopic mass balance of combined The abbreviations:
j,k
Such diffusion equations may be completely wrong, since they are based on an ad hoc assumption, the linear gradient hypothesis (8). There is, indeed, a series of cases where this hypothesis cannot be justified. Here, there is no other alternative but to integrate equation (6) or even (5). Sometimes, the velocity field ;(;, t) and temperature field T(;, t) are derived from so-called circulation models. In most practical cases, however, the velocity- and temperature-fields are experimentally determined, i.e. by measurements. The advection of pollutants might be described afterwards by correlating meteorological data. Since the topography of an urban region is far too complex to be described analytically, one renounces exact solutions of the basic aerodynamical equations and contents oneself with stochastic numerical solutions. Molecular and turbulent diffusion can be quantified by the above mentioned (K-tensor which is usually calculated from measured temperature-, velocity- and trace-gas distributions, via a more or less complex formalism. It seems best, however, to take 1K as a free parameter fixed indirectly by fitting Nj to measured pollutant concentrations. To evaluate Li and Pi both the solar radiation and the molecular scattering- and absorption cross sections must be considered; i.e. one needs a detailed knowledge of the reaction mechanisms governing the chemical processes in the atmosphere. Because of the multiplicity of those phenomena, the theory faces unsurmountable difficulties. With regard to the determination of the characteristic parameters, aik and bijk, there is no alternative to calculating them from reaction-chamber concentration measurements by an inversion of the kinetic equations:
Appl. Math. Modelling,
1980,
Vol 4, August
where k = reaction rate. The parameter A is related to the collision frequency and to the probability that a reaction occurs during a collision, B depends on the temperature of this collision frequency and C is proportional to the activation energy. In practice, the numerical values of A, B, and C are determined by fitting the Arrhenius equation to experimental data. Photolytical reactions, e.g. Si + hv + Sj, are usually parameterized by the so-called photodissociation coefficients which depend in a rather complex way on wavelength X, attenuation factor s, absorption cross section ai, temperature T and ‘quantum yield’ R of the corresponding process. The attenuation factor (optical depth) is an integral over all absorptions and scatterings of the quanta of a given wave length caused by all the chemical species present in the irradiated atmosphere. It is clear that here also, simplifying assumptions have to be made to obtain, at least approximately, useful relations. The radiation also causes secondary effects, such as diffusion, heating, etc., which in their turn create further consequences. Part of the incoming radiation becomes absorbed and excluded from further radiochemical reactions. When including the radiation effects proportional to the intensity
Diffusion
J of the radiation, the loss- and production operators defined by equation (7) have to be modified as follows: Li = C aikNi + ViJ k
(11) Pi = C bijk NiNk f 1 /.likNk J k i,k The photo dissociation respectively the association parameters, Vi and &k, also have to be determined by fitting to the corresponding experimental data. The hypothesis that the impermeable part of space is proportional to the absorbed radiation, gives a complementary equation (a, Ki = absorption rates):
of reacting pollutants:
K. H. Mtiller
Later one should use a time step At comparable to rmaX = maxr,. Unfortunately, such a large time step causes badly diverging numerical errors for most of the classical integration techniques, i.e. these methods are no longer numerically stable in the presence of stiffness. To integrate the set of ordinary differential equations (14), one replaces the first-order derivatives by a finite difference form, e.g.:
(15) The indexj counts the time steps. The resulting set of algebraic equations can be solved in three different ways. The simplest follows the so-called explicit technique:
(16)
(12) which can be integrated formally. Using the resulting expression, the intensity J can be eliminated from the diffusion equation governing the mass balance, an integro-differential equation results. Most authors, however, assume J to be constant.
Numerical treatment A simulation of combined transport-transformation processes occurring in the atmosphere deals with systems of partial differential equations, symbolically written as: (13) where hi is the combined transport-transformation operator. The corresponding physical problem is an initial-boundary value problem with widely differing time constants. Experimentally the chemical reaction rates span a range of time constants from about 10e6 s to 10+5 s; the atmospheric transport processes, sometimes, have characteristic times up to lo’-109s. Because of these facts, the application of procedures based on the usual techniques of operatorsplitting and fractional steps seem to be adequate. HOWever, it soon becomes clear that each time the transport operator is used, the chemical (equation) system is rudely perturbed out of equilibrium, and that subsequent smoothing requires an excessive effort. In connection with this, the question arises of how the simulation reacts to such perturbations. models An approach frequently used in solving system (13) is to replace the spatial derivatives by the corresponding finite differences. The pollution distribution fi(?, t) will be replaced by the vector A(t) whose components ni(t) satisfy nonlinear differential equations:
Eulerian
which works sufficiently well for small time steps.4 The much stabler implicit techniques5-8 evaluate Fi at time step (j + 1) or at some weighted average ofj and (j + 1); their basic form is: (17) The semi-implicit technique9 takes only part of Fi at time point (j + l), e.g. Li. Such a procedure has the advantage of separating the equations in view of the speciesindex i. A completely different structured integration formalism discretizes equation (13) only with respect to time. Using the abbreviations: D=-ii.0
tV.(lKV),
(14)
which possess the unpleasant property of ‘stiffness’, i.e. the representative time constants r, differ in orders of magnitude. The initial value problem of equation (14) is usually solved by a stepwise method. In the beginning, when the state variables ni are still changing very rapidly due to the fast reactions, one has to use a time step At of the order Of 7,rn = min 7,.
Mi=eA2”Ni
and opting for the explicit technique, equation approximately be replaced by the cascade: DM/+l
(9) can
=(L~-V,IKc’-X2)M~-eip~
The introduction satisfying:
of the corresponding
(18) Green’s function
DC = 4716(s) for the same geometry, allows a formal integration equation (18) resulting in the recursion formula: M!+’ = G(j-;)[(Lj I Ir
of
- V JKZ ~ A’) Mi - e’e’] d3p (19)
which is easy to evaluate starting from the initial concentration M/c;). One can show” that the approximate solution N#) constructed in such a way converges uniformly in i and j with At + 0 to the desired density distribution Ni(i, t). The corresponding implicit form can formally be integrated to the Fredholm integral equation: 1
M;+l
J
= G(r’-;)[(LF+‘-V.IKc’)M;+’ I _ ei+lpj+l I
dn. -=Fi(t,nk/k=1,2,...) dt
h2= l/At
-
X’M;] d3p
Inserting one of the extrapolation
(20)
formulae of the type:
(21) in Lij+’ and8+1, the integral equation (20) can be linearized. There are a series of methods for solving such linear equations. The integral structure of the recurrence expressed by the formulae (19) and (20) plus (21) guarantees an excellent numerical stability.
Appl.
Math.
Modelling,
1980,
Vol 4, August
271
Diffusion of reacting pollutants: K. H. Miller
The key-point of the procedure is the knowledge of the Green’s function, G, representing the pollutant dispersion as realistically as possible. A class of rather realistic dispersions are covered by the following operator. For:
and the infinite ;-space we get the following Green’s function. (For an evaluation of equation (19) with respect to equations (20) and (21) see the Appendix.). G(?) =
uoy2
K;~‘~K;’exp
- 4~ i
uozc
- __
2
(22)
C2K3
where u,, measures the wind velocity and:
. I x
Kj = &($) .
dt
b=--
8+1
c=fl-y+2
c
0
The corresponding Green’s function for the half space z > 0 results from a superposition according to the classical mirror method. Finally, one should mention that numerical Green’s functions, i.e. K-models, can also be inserted into equations (18), (19) and (20). Mixed Euler-Lagrangean models The so-called moving-cell methods12y l3 consider homogeneous air parcels or vertical air columns travelling along their trajectories with the velocity 6 of the wind field and taking up pollution during their crossing of source regions. In nature one observes quasi-stable air parcels, eddies, representing temporarily closed systems. They have a certain life-time, at the end of which they decay or coagulate with other eddies. This fact suggests a generalizationll of the above-mentioned cell models which takes into account the phenomenon of ‘life-time’. The ‘free’ flights of the eddies between subsequent collisions are only another interpretation of this life-time. The pollutant transport can be interpreted as a superposition of a continuous ‘streaming’ (convection) and a statistical motion of air parcels (eddy diffusion).” One can, in addition, imagine the streaming to be the average of stochastic particle motions. By such reasoning one returns to the basic idea of the kinetic gas theory. The most obvious difference between a molecular and an eddy-gas consists in the dimensions of the space time volume elements they refer to; i.e. an eddy-gas requires an eddy-scale volume element. Relating eddy-mass to eddy-velocity one can get, finally, a full conformity in modelling of both molecular and eddy-gas. Neither the mechanisms of the molecular nor the eddy collisions are uniquely determinable; they must be suitably parameterized and indirectly quantified by a fitting to pollutant density measurements. Such an eddy-gas model has, because of the temporarily closed-system structure of the eddies, the advantage that the chemical transformations occurring during the free flight can be treated independently of the transport mechanism itself, as each eddy possesses its ‘own’ chemistry.” At the collision points, however, the different chemistries interact. The modeller has a choice between alternative interacting mechanisms. One possible approach assumes that the chemical transformations occur inside the closed eddies during their free flight and that the collisions only represent instantaneous mixing. An alternative assumes
272
Appl.
Math.
Modelling,
1980,
Vol 4, August
that there are no transformations during the flight, only during the mixing process, this time of finite duration. In nature, of course, these two mechanisms take place simultaneously. The number of collisions per Eulerian box and time element is a measure of the eddy number density, and so for the pollutant densities. The local substance composition, Nj(?, t), in a given box changes in the subsequent time interval At, as a result of both the chemical transformations and the continuously occurring eddy transport. Since most of the chemical reactions possess a time constant AT, far smaller than that of the transport processes, one will quantify the incoming and outgoing flow according to the AT-steps; i.e. assume that, during a AT-step, only chemical transformations occur, and that the modifications of the chemical composition caused by the transport need only be taken into account at the end of the different time steps, (t +k.A7) wherek=1,2,... _ A possible approximation accumulates the incoming and outgoing pollutant masses over the whole transport time step At, before combining them with the chemical transformations. The changing density rates are measured, then, via the difference of the corresponding collision densities Ci, at subsequent time points, i.e. ACj = Ci (;, t + A t) - Cj (+, t)
(23)
AC’, will be related to the rate ANi. Starting from the actual density, Ni(?, t), one integrates, in time steps AT, the kinetics equation (10) over the time interval (t, t t At). An ‘intermediate’ density NF(?, t t At) results. The ‘final’ density including transport is given by: Ni(r’, t + At) = NT t ANi
(24)
Intermittent interventions use fragments of ANj calculated by an interpolation based on the previous time behaviour of Nj. Stochastic models To construct as realistic a model as possible, one also has to pay attention to the results of the theory of chemical reactions in nonhomogeneous mixtures, usually subsumed under the key words: stochastic formulation of reaction kinetics.14-l6 It is an improvement of the classical deterministic formulation, in so far as the former is able to simulate the fluctuations of reactant densities and reaction products within small subvolumes. The stochastic formulation assumes that the molecules follow Markovian paths and are activated during collision with other molecules or radiation quanta. After a series of such energy exchanges, a molecule accumulates enough energy to overcome the potential barrier of a collision partner which previously prevented chemical interaction. For a gas the rate of barrier-crossing equals the rate of the reaction. The solution of the stochastic formulation of chemical kinetics becomes asymptotically identical with the solution of the deterministic formulation - in the thermodynamic limit.” One could try to simulate the random walks of the single molecules, in analogy to the real events. The corresponding calculation time would be directly proportional to the number of individual molecular reactions taken into account. The excessive calculation effort forces us to a simplification which has, at any rate, to cover ail the relevant phenomena. Since the collisions without chemical
Diffusion of reacting pollutants:
consequences are much more frequent than the so-called active reactions, one decides to ignore the inactive ones. Between two subsequent active collisions, as in this picture, a free flight represents the vector sum of the intermediate paths, interrupted by the inactive collisions. This simplification is in harmony with the deterministic formulation of the kinetic equations which obtain their reaction coeffcients from chemical reaction rates alone. The stochasticity is not lost during this Monte Carlo procedure. The reaction volume I’, in addition to a series of inert substances, includes homogeneously distributed molecules of I chemically active species, Si, which react with each other in M different ways. z, distinct molecular reactant combinations for reaction m are possible. Each reaction is characterized by a numerical parameter, c, . Reversible reactions are considered to be different. The c, are defined as follows: c, dt = mean probability that an m-reaction occurs in the subsequent time interval dt. They are related to the corresponding reaction constants, k,, of the deterministic formulation by: Reaction
k,
S1 + hu s1+ s1
cm c, . v/2
Sr+S, s,+s,+s,
Cm .v c, . V’, etc.
The stochastic approach starts from the reaction probability density function, ~(7, m), defined by: p(r, m).dr = probability that the next reaction occurring in the time interval (t + 7, t + 7 + dr) is an m-reaction. With this, one can construct a sort of Monte Carlo procedure simulating the stochastic reaction process. The corresponding formula: p(~, m) = a,.exp(-a.7) M
with a,,, = z, .c,,, and a = 1
a,
can easily be factorized:
m=1
~(7, m) = : C
!
.a.exp(-a.7)
=pl(m).Pz(T)
(25)
Using a random number generator, one calculates, from p1 and pZ, a pair (7, m) of random numbers. With their aid one can proceed from t to (t + 7) and find the index m of the presently most probable reaction. Now one identifies all species, Si, involved in reaction m, and calculates the corresponding reaction rates according to the appropriate stoichiometric ratios, and then updates the density distributions. One repeats the above step sequence until a given time limit is reached or no more reactions occur due to the lack of substance. In this way, one constructs a single possible realization of the stochastic transformation process taking place in V. To obtain a statistically complete picture of the whole evolution, one has to construct, for the same initial conditions, numerous realizations independent from each other, and to calculate averages from these different ‘runs’. Such mean values converge for an increasing number of runs to the real time behaviour of the reactions. In air chemistry, the probability pr = a,/a corresponding to the so-called ‘slow’ reactions, is usually extremely small. One, therefore, separates all reactions into different classes, e.g. fast reactions and slow reactions. Sometimes such a binary classification is not sufficient, and it is advisable to
K. H. Miiller
introduce three or more classes. For an illustration of the procedure to be applied, one contents oneself with a binary classification; a generalization to more classes does not create any difficulty. One divides into: <
m
fast s reactions > slow M
a=
s: a, m=1
+a’
a’ = c a, S+l
and
(26)
and uses the following scheme: aria I
(a, + a’ )/a
a& T-1
0
a,/@, + a’ ) II
ON
III
0 -
a’/(a, + a’)
-1 1 aM/a’
4+r la’
A random number generator produces numbers & equally distributed on the interval (0; 1). & thrown on scale I picks out a reaction m < s. If, however, m = s one calculates a new random number, g2, and throws it onto scale II, which now decides whether a slow reaction is meant or not. If yes, a third random number, g3, thrown on scale III, determines the index of a slow reaction. In the case of reaction velocities differing in several orders of magnitude, as is usual in practice, one can dispense with scale II. Form =s, however, one has to take into account both the reaction m = s and a slow one to be identified via scale III. Additionally, an ‘importance weight’, w = al/a, has to be attached to the slow reaction, while w = 1 holds for the fast one. In this way one is able, with only a limited quantity of random numbers, to pay sufficient regard to slow reactions as well. The corresponding life-times are to be calculated from the adequately modified pZ(r)-distributions. An extension of the presented stochastic approach to spatial inhomogeneities can easily be performed; one has to subdivide the reaction volume V into a system of boxes, each homogeneously structured, and to couple the lifetimes with the distances covered by the particles and calculated from the given particle density distributions. Most of the calculation time used by the stochastic approach requires the determination of the reaction index m. A ‘bundling’ of runs is, therefore, recommended, i.e. one assumes that after each time step several elementary reactions occur simultaneously. The use of (artificial) particles disposing of a complete am/a-spectral composition reduces calculation time considerably. Finally, it should be noted that the above algorithm offers a convenient way of by-passing the stiffness difficulties which generally dominate the integration of the coupled kinetics equations. Conclusions The transport of reacting substances in the gas-cover of our globe can be formulated by reactive Boltzmann equations. A suitable parameterization of the single phenomenon replaces the above representation by a set of parabolic differential equations governing the space-time dependent densitiesNi of the interacting substances Si. Since practical problems, in general, do not allow a purely analytical solution of that nonlinear parabolic
Appl.
Math.
Modelling,
1980,
Vol 4, August
273
Diffusion of reacting pollutants: K. H. Miller
system, one is forced to apply a partially or even a purely numerical integration method. For solving such systems, the literature offers a series of numerical procedures differing significantly with regard to structure and calculation effort. Those approaches discretize the spatial differential operator and result in a system of stiff kinetic equations, the solution of which has to overcome the classical difficulties created by stiffness. One can, however, limit or even avoid those difficulties by discretizing only the time-operator of the parabolic system and by transforming the resulting, now purely spatial, formulation using the so-called Green’s technique, which usually does not impose numerical troubles. like instability, etc. For a rather realistic class of Green’s functions, i.e. one of practical significance, a mainly analytical approach has been elaborated (see Appendix) which reduces the remaining numerical effort drastically. Our mathematical model simulating the combined process profits from the fact that the time constants characterizing transformation differ significantly from those characterizing transport. Thus, one can assume some transformation steps occur during the single transport time step At. This assumption could be justified by the so-called eddy-gas mechanism of the transport process which interprets the eddies as closed and, during statistically varying but in any case finite time intervals, free flying units. Each of these eddies carries a certain volume of mixed gases interacting in a proper way. The eddy structure of the gas transport and the molecular collision mechanism behind the chemical transformations suggested the development of a combined particle model interpreting both transport and transformation as a sequence of particle collisions and random paths in between; e.g. eddies as (macro)particles bearing inside a population of interacting (micro)particles. A publication dealing with a particle model simulating the pollutant transport is in preparation. A concluding remark concerns ozone depletion in the stratosphere.r8,r9 It has long been known that the ozone layer covering our planet protects us from the sun’s harmful ultraviolet rays; but only a few years ago it was shown that fluor-carbons may pose a serious threat to the integrity of the ozone layers. The time between the use of CFMs at ground level and the maximum depletion they subsequently produce in the stratospheric ozone layer is of the order of a decade, and it requires decades for the stratosphere to recover fully. This long time lag, together with the fact that there are appreciable long-term natural fluctuations in the ozone layer, means that if no action were taken until the CFMs produced a directly measured reduction in the ozone layer, it would then be too late to avoid further, more drastic long-term effects. As a consequence, mathematical models must be used to analyse the sequence of events and predict the outcome decades later. The accuracy of such predictions depends not only on many directly measurable molecular properties, such as photochemical decomposition and chemical reaction rates, but also on the manner in which the models approximate the motions of the atmosphere and its interactions with the solar radiation. References 1 2
274
Huang, K. ‘Statistical mechanics’, Colgrove, F. D. et al. J. Geophys.
Appl.
Math.
Modelling,
1980,
N.Y., 1963 Rex, 1965, 70,493l
Vol 4, August
3
10 11 12 13 14 15 16 17 18 i9
Leiehton. P. A. ‘Photochemistrv of air oollution’. Academic Press, N.Y., 1961 Lomax, H. and Bailey, H. E. NASA-Techn. Nore D-4109, 1967 Gear.C.N. Commun. ACM, 1971, 14.176 Char& J. S. et al. UCRL-74823, 1973 Crutzen, P. J. Can. J. Chem., 1974,52, 1569 Graedel. T. E. ef al. Atm. Em.. 1976. 16. 1095 Ames, $. F. ‘Numerical methods for partial differential equations’, Academic Press, N.Y., 1969 Grobner-Hofreither, N. ‘Integraltafel’, Springer, Vienna, Austria, 1958 Miiller, K. H. 2. Angew. Math. Phys., 1961, 12, 242 EUR-5564e (1976) and EUR-5633e (1977) Eschenroeder, A. Q. and Martinez, J. R. Adv. Chem. Ser., 1973, 113,101 Sklarew, R. C. etal. J.A.P.C.A., 1972, 22, 865 McQuarrie, D. A. J. Appl. Prob., 1967, 4,413 Bunker, D. L. et al. Comhust. Flame, 1974, 23, 37 Gillespie, D. T. J. Comp. Phys., 1976, 22,403 Oppenheim, I. et al. J. Chem. Phys., 1969,50,460 Hard, T. M. and Bradwick, A. J. DOT-TSC-OST-75-38, 1975 McCracken, M. C. UCRL-51336, 1975
Appendix The recursive calculation of the pollutant density Nj becomes rather simple when applying class (22) as Green’s function, G. Developing then the content of the rectangular brackets, [. . .], of formula (19), with respect to (20) plus (21) in series of Hermitean polynomials Hk:
[...I = C
Qifl($)Hrn(7))HnK)
(AlI
m,n where
6 = (L v,_Q,
one obtains the recurrence
relation:
+-
dSew-udz
- ~?/(c~dl H,(C)
(A2)
-m
Applying the following H,(T/&)
formulaelo:
= 2-n’2
H,_,(t)
H,(T-
t)
+m
I
e-(T-f)ZHm(t&?)
dt = &r(T&)”
-m
0
teCtz PH,(t)
s
-m
dt =
n!fi __ k!4k
n
for
n-m=2k
the evaluation of both the n- and the f-integral can be done analytically. For complex-structured functions K~(x) and KS(x), however, the subsequent &integration has to be performed numerically. Finally, we obtain a representation of the type: M/+l = c C;,(X)y*Z” (A3) m,n characterizing the pollutant density Ni(?, t) of substance Si at time t=(i+l)At.
The Boltzmann equation including reactive scattering terms is the correct description of the combined transport-reaction process. A suitable parameterization transforms it into a set of differential equations governing the macroscopic mass balance. An operator splitting integration method composed of the stochastic formulation of chemical kinetics and the Green’s representation of the transport proves to be very efficient.
Problem
description
Thousands of organic and inorganic substances are continuously emitted into the atmosphere by natural and man-made sources. They stay there for some time. Unstable molecules decay, and some of their fragments join other components. The photochemical transformations are of particular importance. Most of the gaseous species released become ‘digested’ and therefore play only a minor role in atmospheric pollution. Main sinks are wash-out, biological interactions and decomposition after absorption of radiation. Unfortunately, there are also chemical reactions with pollutants, surface materials and ‘natural’ gases leading to noxious reaction products. Pollutant transport and chemistry in the atmosphere is a domain where one assumes that the basic principles are known. But theoretical treatment is possible only in the simplest cases. It is obvious that a decision maker and one who is engaged in solving practical problems cannot wait until a theoretical solution of the problem is published. Reliable predictions in air hygienics can only be given by a simulation of the physical and chemical processes governing pollutant dispersion and transformation. Their full complexity, however, cannot be simulated experimentally on a laboratory scale. Here, a numerical simulation, i.e. modelling, using an electronic computer, is inevitable. The aim of atmospheric modelling consists both of obtaining insight into the interactions between radiation, chemistry and dynamics and of predicting possible environmental modifications caused by intentional or unintentional human interventions in nature. Since a detailed simultaneous modelling of pollutant chemistry and dispersion proves rather troublesome, one tends to purchase a detailed description of one phenomenon with a rough approximation of another. The extent to which this is justified has to be checked carefully from case to case. Only a simulation containing very detailed mechanisms can decide whether one phenomenon is more important than another. The bottleneck here is the scarcity of the
268
Appl.
Math.
Modelling,
1980,
Vol 4, August
experimental data necessary both to uniquely characterize a given situation and to verify the adequacy of the model. First, one checks the usefulness and completeness of the available experimental data by applying simply structured models (sub-models) and, after that, if necessary, asks for additional measurements. A more complete data set will then represent a more reliable basis suitable for updating and generalizing the previous model version, etc. Such a procedure is somewhat like that of successive approximations. Finally, one will be able both to distinguish between relevant and unimportant events and to extrapolate to situations which are difficult or impossible to reach experimentally. The vitality of a model lies in the availability of an extensive parameterization which allows a wide span of response characteristics of the physical system in question to be studied. Basically, it is not difficult to construct a model and to fix its parameters by fitting to measured data, if only a few species are to be treated; it becomes difficult, however, if several species and all possible atmospheric conditions occurring in a region, have to be taken into account. Through a sensitivity analysis which varies the characteristic parameters and relates these variations to their consequences, a reaction of minor significance can be traced and subsequently eliminated from the model. Input parameters are: dimensions of reaction-space and boundary conditions, transport coefficients, reaction rates, photodissociation coefficients, etc. It is clear that a specific parameterization of the model input determines the overall structure of the model itself. One could, for instance, say that the realism of a model is seen in the accuracy of the characteristic response time of the individual processes; and one could take this as a guide-line for a ‘realistic’ description of the interplay of the different processes under consideration. Certainly, with rather simple models one can already get a good survey of the whole situation and the trend of an evolution. The crucial thing is to find a suitable parameterization of the characteristic processes and of their causes and consequences. Also, one has to ensure that the fitting approach, which mainly hinges on variations of 0307-904X/80/040268-1 0 1980 IPC Business
l/$02.00 Press
Diffusion
parameters, gives reasonable orders of magnitude for each process under study. Air quality models will usually be judged by a comparison of their output with measured concentrations of atmospheric trace species. The corresponding validation experiments mostly rely on a small number of reactants. It must be remembered that a model for only a few reactants may not be unique, e.g. with respect to smog formation, since several different mechanisms can be constructed and adapted to given reaction-chamber experiments just by a suitable choice of parameters.
Mathematical
formulation
The most general principles in physics, such as the conservation of mass, momentum and energy, form its basis. The mathematical expresssion for a conservation principle is a continuity equation. Such equations must be established for each chemical species with a space-time varying particle density. The resulting equations, coupled to each other, due to all the interactions occurring in the atmosphere, have to be solved simultaneously to obtain the real space-time behaviour of the concentrations. It is clear that this set of equations, a model of the chemical kinetics transport, can only be solved numerically, and that only such a numerical simulation of the true interplay of the different atmospheric phenomena enables us to study it and to infer trends and consequences from the actual situation. The kinetic theory of gases’ represents a uniform basis for models formulating the full interaction between radiative, photochemical and aerodynamic processes, and is, therefore, a suitable starting point for the mathematics. Continuous functions of space and time have to be assigned to the bulk matter properties, such as mass density, energy density, velocity field, etc. Only when it is possible to define a volume element of microscopic dimensions which still contains a sufficient number of molecules, can one speak of local density, local velocity, etc., in the sense of average values. Simultaneously, a time interval has to be defined, macroscopically small but large enough so that a sufficient number of molecular events, such as collisions, fluctuations, etc., can be processed to reasonable average values. The main task of the molecular theory is to show that such mean values can reasonably be defined and that they satisfy the continuum equations. One is interested in the particle distribution function fi(?, t, 6) counting the number of molecules of species Si within the phase-space volume element. The functions fi are considered to be continuous in all their variables. In this sense, the normalization holds:
the acceleration field acting on the molecules. In case of simultaneously occurring collisions, however:
#O
holds. Obviously, one has to know the formal structure of the right-hand ‘collision’ term which is a measure for the rate of change offi due to collisions. The basic assumption of the kinetic theory is: the gas is dilute; i.e. the mean distance between the gas molecules is much larger than the range of their intermolecular forces. Occasionally, one molecule hits another molecule; and, much more infrequently, three molecules collide simultaneously. Within a mathematical analysis, this fact can be expressed by the statement that there are only two-body or binary collisions. Having constructed a corresponding binary theory, one will try to generalize it to three-body collisions, etc. A gas is said to be in non-equilibrium, when its distribution function fi deviates from the Maxwell-Boltzmann distribution. Causes are heterogeneous temperature-, density- or velocity-fields. The system of molecules tries to equalize these irregularities by transport of energy, momentum and mass. The corresponding transport mechanism is molecular collision. One usually assumes, as a working hypothesis, that the molecules fly freely between two subsequent collisions, i.e. a’= 0, and that velocity 5 and position i of a molecule are uncorrelated (‘molecular chaos’). The number of pairs of molecules in d3r of two different substances, Si and Sk, whose velocities are 6i and Gx-,is given by: [fi(;, t, 4) d3r d3vi]. [fk(f, t, &) d3r d3vk] The molecules of substance Si will be ‘attacked’ by the current: fk(?, t, Gk) d3V, IGj - 6, I of Sk molecules. The number of collisions occurring during the time interval dt is: fi(Z’, t, 6j) dtojk($,
iTjk)fk(?, t, Gk) d3Vkl$ - 6kl
where oik denotes the ‘cross section’ characterizing the collision mechanism. During such collisions, chemical transformations may occur. The corresponding changes in particle densities are assumed to be proportional to the number of collisions. The proportionality factors, p and 4, are called (positive or negative) production rates. The velocity offi-changes, due to collisions, is therefore given by: = -hG,
t, G) d3vk s
cPik(j, k
jk)
X fk(?,t,Zk).IG-ZkI + d3Vj
The aim of the kinetic theory is to find the distribution functionsfi corresponding to given molecular interactions. The density fi changes in time, since molecules enter and leave the volume element dV= d3r d3v all the time. This transfer can be split into a ‘streaming’ and a ‘scattering’, interpreting the particle motion as the consequences of collisions between molecules. For collision-free motion, the density function h satisfies the ‘streaming’ equation: fi=O
(2)
collision
i at 1 collision
= number of particles in (d3r dt)
K. H. Miller
a’ denotes
afi
fi(;, t, 6)d371 =Ni(g, t)
of reacting pollutants:
s
X c
qijk(G> 4, ;k)
o.ik(j,
6,)
d3Vk .i
ajk(Gj, Gk)
j,k Xfi(‘:
(1)
t,
zj)fk(i, t,;,)I'$ - Gkl
(3)
The balance equation resulting from equations (2) and (3) and including reactive scattering is a correct dynamical description of a combined transport-transformation process in a dilute gas; it is equivalent to the uncorrelated binary collision theory of chemical kinetics.
Appl.
Math.
Modelling,
1980,
Vol 4, August
269
Diffusion of reacting pollutants: K. H. Miiller
The continuity equations of the mass-, momentum- and energy-transport can be derived from this Boltzmann equation by applying a weighted averaging of the molecular velocity 5. Our interest lies in deriving a compact model for the pollutant-mass transport including the simultaneously occurring chemical transformations in an atmosphere characterized by a given velocity- and temperature-field. Using the already defined density Nj and the particle current: ii = r Gfi d3V J
and the following mean values: d3Vk =@ik(q
pikujk.I;-6k(fk
fk d3Vk =&jk(qNk I
s qijkojk
’
1I$ - ?k 1fi fk d3Vi d3Vk
= ,!&@) s the Boltzmann
fi d3vi fk d3Vk = &jk($ NiNk I
equation
can be rewritten
1
CYijNk h=
k
(4)
as follows:
C PijkNjNk
(5)
i>k
An integration of this equation over the whole velocity space and the introduction of:
.i
oik fi d3v = aik
wi
fid3V = UikNi
leads to the parameterized i3Nj z t div2i = -LiNi
Pijk
.r
d3V = bijk
form (with a’ = 0):
dNi dt=
k
Pi = C biikN;Nk
(10)
k = ATB.exp(-C/T)
‘
denote the ‘loss’ and ‘production’ of Si respectively. Before trying to parameterize the current ii one introduces the mean velocity ii of all molecules actually present in the volume element (d3r dt); i.e.: ii = $fi dav = GNi - (; - ;) fi d3v z ;Ni - fi s I The space-time dependent function ?j expresses a sort of mean-dispersion behaviour of substance Si. According to classical macroscopic diffusion theory one could perform a parameterization by a linear ansatz of the type: ?i = lK(ZNi + grad Ni)
(8)
The tensor 1K serves to quantify the dispersion behaviour. One tries to include the influence of a given temperature field Ton the chemical kinetics transport by the ansatz*: c’ = V (1nT) + constant Inserting the above I?i into equation (6) and putting V . li = 0, a nonlinear diffusion equation results: aNi ~+~.VNi-V.(IKVNj)+(Li-V.lK~)N,=Pi (9)
270
+ Pi
valid for homogeneous mixtures.3 Experiments of this kind show a strong temperature dependence of the reaction rates, &k and bijk. One, therefore, tries to find at least their approximate T-structure. The temperature behaviour of most of the reactions relevant to air pollution follows Arrhenius’ formula:
processes.
Li = CfZikNk
’
~LiNi
t Pi
of the macroscopic mass balance of combined The abbreviations:
j,k
Such diffusion equations may be completely wrong, since they are based on an ad hoc assumption, the linear gradient hypothesis (8). There is, indeed, a series of cases where this hypothesis cannot be justified. Here, there is no other alternative but to integrate equation (6) or even (5). Sometimes, the velocity field ;(;, t) and temperature field T(;, t) are derived from so-called circulation models. In most practical cases, however, the velocity- and temperature-fields are experimentally determined, i.e. by measurements. The advection of pollutants might be described afterwards by correlating meteorological data. Since the topography of an urban region is far too complex to be described analytically, one renounces exact solutions of the basic aerodynamical equations and contents oneself with stochastic numerical solutions. Molecular and turbulent diffusion can be quantified by the above mentioned (K-tensor which is usually calculated from measured temperature-, velocity- and trace-gas distributions, via a more or less complex formalism. It seems best, however, to take 1K as a free parameter fixed indirectly by fitting Nj to measured pollutant concentrations. To evaluate Li and Pi both the solar radiation and the molecular scattering- and absorption cross sections must be considered; i.e. one needs a detailed knowledge of the reaction mechanisms governing the chemical processes in the atmosphere. Because of the multiplicity of those phenomena, the theory faces unsurmountable difficulties. With regard to the determination of the characteristic parameters, aik and bijk, there is no alternative to calculating them from reaction-chamber concentration measurements by an inversion of the kinetic equations:
Appl. Math. Modelling,
1980,
Vol 4, August
where k = reaction rate. The parameter A is related to the collision frequency and to the probability that a reaction occurs during a collision, B depends on the temperature of this collision frequency and C is proportional to the activation energy. In practice, the numerical values of A, B, and C are determined by fitting the Arrhenius equation to experimental data. Photolytical reactions, e.g. Si + hv + Sj, are usually parameterized by the so-called photodissociation coefficients which depend in a rather complex way on wavelength X, attenuation factor s, absorption cross section ai, temperature T and ‘quantum yield’ R of the corresponding process. The attenuation factor (optical depth) is an integral over all absorptions and scatterings of the quanta of a given wave length caused by all the chemical species present in the irradiated atmosphere. It is clear that here also, simplifying assumptions have to be made to obtain, at least approximately, useful relations. The radiation also causes secondary effects, such as diffusion, heating, etc., which in their turn create further consequences. Part of the incoming radiation becomes absorbed and excluded from further radiochemical reactions. When including the radiation effects proportional to the intensity
Diffusion
J of the radiation, the loss- and production operators defined by equation (7) have to be modified as follows: Li = C aikNi + ViJ k
(11) Pi = C bijk NiNk f 1 /.likNk J k i,k The photo dissociation respectively the association parameters, Vi and &k, also have to be determined by fitting to the corresponding experimental data. The hypothesis that the impermeable part of space is proportional to the absorbed radiation, gives a complementary equation (a, Ki = absorption rates):
of reacting pollutants:
K. H. Mtiller
Later one should use a time step At comparable to rmaX = maxr,. Unfortunately, such a large time step causes badly diverging numerical errors for most of the classical integration techniques, i.e. these methods are no longer numerically stable in the presence of stiffness. To integrate the set of ordinary differential equations (14), one replaces the first-order derivatives by a finite difference form, e.g.:
(15) The indexj counts the time steps. The resulting set of algebraic equations can be solved in three different ways. The simplest follows the so-called explicit technique:
(16)
(12) which can be integrated formally. Using the resulting expression, the intensity J can be eliminated from the diffusion equation governing the mass balance, an integro-differential equation results. Most authors, however, assume J to be constant.
Numerical treatment A simulation of combined transport-transformation processes occurring in the atmosphere deals with systems of partial differential equations, symbolically written as: (13) where hi is the combined transport-transformation operator. The corresponding physical problem is an initial-boundary value problem with widely differing time constants. Experimentally the chemical reaction rates span a range of time constants from about 10e6 s to 10+5 s; the atmospheric transport processes, sometimes, have characteristic times up to lo’-109s. Because of these facts, the application of procedures based on the usual techniques of operatorsplitting and fractional steps seem to be adequate. HOWever, it soon becomes clear that each time the transport operator is used, the chemical (equation) system is rudely perturbed out of equilibrium, and that subsequent smoothing requires an excessive effort. In connection with this, the question arises of how the simulation reacts to such perturbations. models An approach frequently used in solving system (13) is to replace the spatial derivatives by the corresponding finite differences. The pollution distribution fi(?, t) will be replaced by the vector A(t) whose components ni(t) satisfy nonlinear differential equations:
Eulerian
which works sufficiently well for small time steps.4 The much stabler implicit techniques5-8 evaluate Fi at time step (j + 1) or at some weighted average ofj and (j + 1); their basic form is: (17) The semi-implicit technique9 takes only part of Fi at time point (j + l), e.g. Li. Such a procedure has the advantage of separating the equations in view of the speciesindex i. A completely different structured integration formalism discretizes equation (13) only with respect to time. Using the abbreviations: D=-ii.0
tV.(lKV),
(14)
which possess the unpleasant property of ‘stiffness’, i.e. the representative time constants r, differ in orders of magnitude. The initial value problem of equation (14) is usually solved by a stepwise method. In the beginning, when the state variables ni are still changing very rapidly due to the fast reactions, one has to use a time step At of the order Of 7,rn = min 7,.
Mi=eA2”Ni
and opting for the explicit technique, equation approximately be replaced by the cascade: DM/+l
(9) can
=(L~-V,IKc’-X2)M~-eip~
The introduction satisfying:
of the corresponding
(18) Green’s function
DC = 4716(s) for the same geometry, allows a formal integration equation (18) resulting in the recursion formula: M!+’ = G(j-;)[(Lj I Ir
of
- V JKZ ~ A’) Mi - e’e’] d3p (19)
which is easy to evaluate starting from the initial concentration M/c;). One can show” that the approximate solution N#) constructed in such a way converges uniformly in i and j with At + 0 to the desired density distribution Ni(i, t). The corresponding implicit form can formally be integrated to the Fredholm integral equation: 1
M;+l
J
= G(r’-;)[(LF+‘-V.IKc’)M;+’ I _ ei+lpj+l I
dn. -=Fi(t,nk/k=1,2,...) dt
h2= l/At
-
X’M;] d3p
Inserting one of the extrapolation
(20)
formulae of the type:
(21) in Lij+’ and8+1, the integral equation (20) can be linearized. There are a series of methods for solving such linear equations. The integral structure of the recurrence expressed by the formulae (19) and (20) plus (21) guarantees an excellent numerical stability.
Appl.
Math.
Modelling,
1980,
Vol 4, August
271
Diffusion of reacting pollutants: K. H. Miller
The key-point of the procedure is the knowledge of the Green’s function, G, representing the pollutant dispersion as realistically as possible. A class of rather realistic dispersions are covered by the following operator. For:
and the infinite ;-space we get the following Green’s function. (For an evaluation of equation (19) with respect to equations (20) and (21) see the Appendix.). G(?) =
uoy2
K;~‘~K;’exp
- 4~ i
uozc
- __
2
(22)
C2K3
where u,, measures the wind velocity and:
. I x
Kj = &($) .
dt
b=--
8+1
c=fl-y+2
c
0
The corresponding Green’s function for the half space z > 0 results from a superposition according to the classical mirror method. Finally, one should mention that numerical Green’s functions, i.e. K-models, can also be inserted into equations (18), (19) and (20). Mixed Euler-Lagrangean models The so-called moving-cell methods12y l3 consider homogeneous air parcels or vertical air columns travelling along their trajectories with the velocity 6 of the wind field and taking up pollution during their crossing of source regions. In nature one observes quasi-stable air parcels, eddies, representing temporarily closed systems. They have a certain life-time, at the end of which they decay or coagulate with other eddies. This fact suggests a generalizationll of the above-mentioned cell models which takes into account the phenomenon of ‘life-time’. The ‘free’ flights of the eddies between subsequent collisions are only another interpretation of this life-time. The pollutant transport can be interpreted as a superposition of a continuous ‘streaming’ (convection) and a statistical motion of air parcels (eddy diffusion).” One can, in addition, imagine the streaming to be the average of stochastic particle motions. By such reasoning one returns to the basic idea of the kinetic gas theory. The most obvious difference between a molecular and an eddy-gas consists in the dimensions of the space time volume elements they refer to; i.e. an eddy-gas requires an eddy-scale volume element. Relating eddy-mass to eddy-velocity one can get, finally, a full conformity in modelling of both molecular and eddy-gas. Neither the mechanisms of the molecular nor the eddy collisions are uniquely determinable; they must be suitably parameterized and indirectly quantified by a fitting to pollutant density measurements. Such an eddy-gas model has, because of the temporarily closed-system structure of the eddies, the advantage that the chemical transformations occurring during the free flight can be treated independently of the transport mechanism itself, as each eddy possesses its ‘own’ chemistry.” At the collision points, however, the different chemistries interact. The modeller has a choice between alternative interacting mechanisms. One possible approach assumes that the chemical transformations occur inside the closed eddies during their free flight and that the collisions only represent instantaneous mixing. An alternative assumes
272
Appl.
Math.
Modelling,
1980,
Vol 4, August
that there are no transformations during the flight, only during the mixing process, this time of finite duration. In nature, of course, these two mechanisms take place simultaneously. The number of collisions per Eulerian box and time element is a measure of the eddy number density, and so for the pollutant densities. The local substance composition, Nj(?, t), in a given box changes in the subsequent time interval At, as a result of both the chemical transformations and the continuously occurring eddy transport. Since most of the chemical reactions possess a time constant AT, far smaller than that of the transport processes, one will quantify the incoming and outgoing flow according to the AT-steps; i.e. assume that, during a AT-step, only chemical transformations occur, and that the modifications of the chemical composition caused by the transport need only be taken into account at the end of the different time steps, (t +k.A7) wherek=1,2,... _ A possible approximation accumulates the incoming and outgoing pollutant masses over the whole transport time step At, before combining them with the chemical transformations. The changing density rates are measured, then, via the difference of the corresponding collision densities Ci, at subsequent time points, i.e. ACj = Ci (;, t + A t) - Cj (+, t)
(23)
AC’, will be related to the rate ANi. Starting from the actual density, Ni(?, t), one integrates, in time steps AT, the kinetics equation (10) over the time interval (t, t t At). An ‘intermediate’ density NF(?, t t At) results. The ‘final’ density including transport is given by: Ni(r’, t + At) = NT t ANi
(24)
Intermittent interventions use fragments of ANj calculated by an interpolation based on the previous time behaviour of Nj. Stochastic models To construct as realistic a model as possible, one also has to pay attention to the results of the theory of chemical reactions in nonhomogeneous mixtures, usually subsumed under the key words: stochastic formulation of reaction kinetics.14-l6 It is an improvement of the classical deterministic formulation, in so far as the former is able to simulate the fluctuations of reactant densities and reaction products within small subvolumes. The stochastic formulation assumes that the molecules follow Markovian paths and are activated during collision with other molecules or radiation quanta. After a series of such energy exchanges, a molecule accumulates enough energy to overcome the potential barrier of a collision partner which previously prevented chemical interaction. For a gas the rate of barrier-crossing equals the rate of the reaction. The solution of the stochastic formulation of chemical kinetics becomes asymptotically identical with the solution of the deterministic formulation - in the thermodynamic limit.” One could try to simulate the random walks of the single molecules, in analogy to the real events. The corresponding calculation time would be directly proportional to the number of individual molecular reactions taken into account. The excessive calculation effort forces us to a simplification which has, at any rate, to cover ail the relevant phenomena. Since the collisions without chemical
Diffusion of reacting pollutants:
consequences are much more frequent than the so-called active reactions, one decides to ignore the inactive ones. Between two subsequent active collisions, as in this picture, a free flight represents the vector sum of the intermediate paths, interrupted by the inactive collisions. This simplification is in harmony with the deterministic formulation of the kinetic equations which obtain their reaction coeffcients from chemical reaction rates alone. The stochasticity is not lost during this Monte Carlo procedure. The reaction volume I’, in addition to a series of inert substances, includes homogeneously distributed molecules of I chemically active species, Si, which react with each other in M different ways. z, distinct molecular reactant combinations for reaction m are possible. Each reaction is characterized by a numerical parameter, c, . Reversible reactions are considered to be different. The c, are defined as follows: c, dt = mean probability that an m-reaction occurs in the subsequent time interval dt. They are related to the corresponding reaction constants, k,, of the deterministic formulation by: Reaction
k,
S1 + hu s1+ s1
cm c, . v/2
Sr+S, s,+s,+s,
Cm .v c, . V’, etc.
The stochastic approach starts from the reaction probability density function, ~(7, m), defined by: p(r, m).dr = probability that the next reaction occurring in the time interval (t + 7, t + 7 + dr) is an m-reaction. With this, one can construct a sort of Monte Carlo procedure simulating the stochastic reaction process. The corresponding formula: p(~, m) = a,.exp(-a.7) M
with a,,, = z, .c,,, and a = 1
a,
can easily be factorized:
m=1
~(7, m) = : C
!
.a.exp(-a.7)
=pl(m).Pz(T)
(25)
Using a random number generator, one calculates, from p1 and pZ, a pair (7, m) of random numbers. With their aid one can proceed from t to (t + 7) and find the index m of the presently most probable reaction. Now one identifies all species, Si, involved in reaction m, and calculates the corresponding reaction rates according to the appropriate stoichiometric ratios, and then updates the density distributions. One repeats the above step sequence until a given time limit is reached or no more reactions occur due to the lack of substance. In this way, one constructs a single possible realization of the stochastic transformation process taking place in V. To obtain a statistically complete picture of the whole evolution, one has to construct, for the same initial conditions, numerous realizations independent from each other, and to calculate averages from these different ‘runs’. Such mean values converge for an increasing number of runs to the real time behaviour of the reactions. In air chemistry, the probability pr = a,/a corresponding to the so-called ‘slow’ reactions, is usually extremely small. One, therefore, separates all reactions into different classes, e.g. fast reactions and slow reactions. Sometimes such a binary classification is not sufficient, and it is advisable to
K. H. Miiller
introduce three or more classes. For an illustration of the procedure to be applied, one contents oneself with a binary classification; a generalization to more classes does not create any difficulty. One divides into: <
m
fast s reactions > slow M
a=
s: a, m=1
+a’
a’ = c a, S+l
and
(26)
and uses the following scheme: aria I
(a, + a’ )/a
a& T-1
0
a,/@, + a’ ) II
ON
III
0 -
a’/(a, + a’)
-1 1 aM/a’
4+r la’
A random number generator produces numbers & equally distributed on the interval (0; 1). & thrown on scale I picks out a reaction m < s. If, however, m = s one calculates a new random number, g2, and throws it onto scale II, which now decides whether a slow reaction is meant or not. If yes, a third random number, g3, thrown on scale III, determines the index of a slow reaction. In the case of reaction velocities differing in several orders of magnitude, as is usual in practice, one can dispense with scale II. Form =s, however, one has to take into account both the reaction m = s and a slow one to be identified via scale III. Additionally, an ‘importance weight’, w = al/a, has to be attached to the slow reaction, while w = 1 holds for the fast one. In this way one is able, with only a limited quantity of random numbers, to pay sufficient regard to slow reactions as well. The corresponding life-times are to be calculated from the adequately modified pZ(r)-distributions. An extension of the presented stochastic approach to spatial inhomogeneities can easily be performed; one has to subdivide the reaction volume V into a system of boxes, each homogeneously structured, and to couple the lifetimes with the distances covered by the particles and calculated from the given particle density distributions. Most of the calculation time used by the stochastic approach requires the determination of the reaction index m. A ‘bundling’ of runs is, therefore, recommended, i.e. one assumes that after each time step several elementary reactions occur simultaneously. The use of (artificial) particles disposing of a complete am/a-spectral composition reduces calculation time considerably. Finally, it should be noted that the above algorithm offers a convenient way of by-passing the stiffness difficulties which generally dominate the integration of the coupled kinetics equations. Conclusions The transport of reacting substances in the gas-cover of our globe can be formulated by reactive Boltzmann equations. A suitable parameterization of the single phenomenon replaces the above representation by a set of parabolic differential equations governing the space-time dependent densitiesNi of the interacting substances Si. Since practical problems, in general, do not allow a purely analytical solution of that nonlinear parabolic
Appl.
Math.
Modelling,
1980,
Vol 4, August
273
Diffusion of reacting pollutants: K. H. Miller
system, one is forced to apply a partially or even a purely numerical integration method. For solving such systems, the literature offers a series of numerical procedures differing significantly with regard to structure and calculation effort. Those approaches discretize the spatial differential operator and result in a system of stiff kinetic equations, the solution of which has to overcome the classical difficulties created by stiffness. One can, however, limit or even avoid those difficulties by discretizing only the time-operator of the parabolic system and by transforming the resulting, now purely spatial, formulation using the so-called Green’s technique, which usually does not impose numerical troubles. like instability, etc. For a rather realistic class of Green’s functions, i.e. one of practical significance, a mainly analytical approach has been elaborated (see Appendix) which reduces the remaining numerical effort drastically. Our mathematical model simulating the combined process profits from the fact that the time constants characterizing transformation differ significantly from those characterizing transport. Thus, one can assume some transformation steps occur during the single transport time step At. This assumption could be justified by the so-called eddy-gas mechanism of the transport process which interprets the eddies as closed and, during statistically varying but in any case finite time intervals, free flying units. Each of these eddies carries a certain volume of mixed gases interacting in a proper way. The eddy structure of the gas transport and the molecular collision mechanism behind the chemical transformations suggested the development of a combined particle model interpreting both transport and transformation as a sequence of particle collisions and random paths in between; e.g. eddies as (macro)particles bearing inside a population of interacting (micro)particles. A publication dealing with a particle model simulating the pollutant transport is in preparation. A concluding remark concerns ozone depletion in the stratosphere.r8,r9 It has long been known that the ozone layer covering our planet protects us from the sun’s harmful ultraviolet rays; but only a few years ago it was shown that fluor-carbons may pose a serious threat to the integrity of the ozone layers. The time between the use of CFMs at ground level and the maximum depletion they subsequently produce in the stratospheric ozone layer is of the order of a decade, and it requires decades for the stratosphere to recover fully. This long time lag, together with the fact that there are appreciable long-term natural fluctuations in the ozone layer, means that if no action were taken until the CFMs produced a directly measured reduction in the ozone layer, it would then be too late to avoid further, more drastic long-term effects. As a consequence, mathematical models must be used to analyse the sequence of events and predict the outcome decades later. The accuracy of such predictions depends not only on many directly measurable molecular properties, such as photochemical decomposition and chemical reaction rates, but also on the manner in which the models approximate the motions of the atmosphere and its interactions with the solar radiation. References 1 2
274
Huang, K. ‘Statistical mechanics’, Colgrove, F. D. et al. J. Geophys.
Appl.
Math.
Modelling,
1980,
N.Y., 1963 Rex, 1965, 70,493l
Vol 4, August
3
10 11 12 13 14 15 16 17 18 i9
Leiehton. P. A. ‘Photochemistrv of air oollution’. Academic Press, N.Y., 1961 Lomax, H. and Bailey, H. E. NASA-Techn. Nore D-4109, 1967 Gear.C.N. Commun. ACM, 1971, 14.176 Char& J. S. et al. UCRL-74823, 1973 Crutzen, P. J. Can. J. Chem., 1974,52, 1569 Graedel. T. E. ef al. Atm. Em.. 1976. 16. 1095 Ames, $. F. ‘Numerical methods for partial differential equations’, Academic Press, N.Y., 1969 Grobner-Hofreither, N. ‘Integraltafel’, Springer, Vienna, Austria, 1958 Miiller, K. H. 2. Angew. Math. Phys., 1961, 12, 242 EUR-5564e (1976) and EUR-5633e (1977) Eschenroeder, A. Q. and Martinez, J. R. Adv. Chem. Ser., 1973, 113,101 Sklarew, R. C. etal. J.A.P.C.A., 1972, 22, 865 McQuarrie, D. A. J. Appl. Prob., 1967, 4,413 Bunker, D. L. et al. Comhust. Flame, 1974, 23, 37 Gillespie, D. T. J. Comp. Phys., 1976, 22,403 Oppenheim, I. et al. J. Chem. Phys., 1969,50,460 Hard, T. M. and Bradwick, A. J. DOT-TSC-OST-75-38, 1975 McCracken, M. C. UCRL-51336, 1975
Appendix The recursive calculation of the pollutant density Nj becomes rather simple when applying class (22) as Green’s function, G. Developing then the content of the rectangular brackets, [. . .], of formula (19), with respect to (20) plus (21) in series of Hermitean polynomials Hk:
[...I = C
Qifl($)Hrn(7))HnK)
(AlI
m,n where
6 = (L v,_Q,
one obtains the recurrence
relation:
+-
dSew-udz
- ~?/(c~dl H,(C)
(A2)
-m
Applying the following H,(T/&)
formulaelo:
= 2-n’2
H,_,(t)
H,(T-
t)
+m
I
e-(T-f)ZHm(t&?)
dt = &r(T&)”
-m
0
teCtz PH,(t)
s
-m
dt =
n!fi __ k!4k
n
for
n-m=2k
the evaluation of both the n- and the f-integral can be done analytically. For complex-structured functions K~(x) and KS(x), however, the subsequent &integration has to be performed numerically. Finally, we obtain a representation of the type: M/+l = c C;,(X)y*Z” (A3) m,n characterizing the pollutant density Ni(?, t) of substance Si at time t=(i+l)At.