7
Environmental Modelling & Software 15 (2000) 673–679 www.elsevier.com/locate/envsoft Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry q Bruno Sportisse a,* , Guy Bencteux b , Pierre Plion b a ENPC-CEREVE, Centre d’Enseignement et de Recherche en Eau, Ville et Environnement, Ecole Nationale des Ponts et Chausse ´es (ENPC–CEREVE), rue Blaise Pascal, 77455 Champs sur Marne, France b DER, A3UR project, Electricite ´ de France, 78401 Chatou, France Abstract Operator Splitting methods are widely used for the numerical computation of reactive flows. The classical analysis of the error induced by such techniques is based on asymptotic expansions with respect to the splitting timestep Dt (supposed to be a small parameter). Such an assumption is unfortunately rather difficult to meet in real computations since Dt is actually much larger than the smallest physical timescales (defined by fast chemical phenomena). The usual error analysis has then to be replaced by another approach, which is associated with the multi-scale behaviour of the coupled system. We give here some theoretical indications and illustrate them with a monodimensional application in Air Pollution Modelling. 2000 Elsevier Science Ltd. All rights reserved. Keywords: Operator Splitting; Stiffness; Singular perturbation; Reduction of dynamical systems; Reaction–diffusion; Air pollution modelling 1. Introduction Operator Splitting methods are commonly used in air pollution modelling (McRae et al., 1982; Pu Sun, 1996; Kim and Cho, 1997). The basic idea is to integrate the different physical phenomena one by one. There are mainly two advantages. First, specific tailored (on-the- shelf) solvers can be used for each operator to be inte- grated (e.g. advection, turbulent diffusion, chemical kin- etics, emission processes, aso). This is particularly prag- matic since black boxes can be easily incorporated in existing softwares. The second advantage is associated with CPU requirements. A Method of Lines is indeed based on the following approach: first discretize spatial operators and then time integrate (with an ODE solver). As numerical stiffness is mainly induced by chemical kinetics (due to the wide range of chemical timescales), implicit schemes have to be used for the time integration of the resulting ODES. This increases drastically the CPU costs due to the necessary algebraic manipulations since the dimension is given by the product of the num- ber of grid cells by the number of chemical species. Even q Presented at the Conference APMS ’98, 1998, October 26–29, ENPC and INRIA * Corresponding author. E-mail address: [email protected] (B. Sportisse). 1364-8152/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved. PII:S1364-8152(00)00036-0 if some specific structure can be taken into account (in the case of a diagonal diffusivity matrix, there is only a local coupling through chemical kinetics: see for instance Graf and Moussiopoulos, 1991), the compu- tational cost of the Method of Lines (MOL) remains pro- hibitive and using Operator Splitting is recommended. In return the accuracy may be unfortunately lowered by the uncoupling of operators. The usual error analysis is based on asymptotic expansions of the local error with respect to the splitting timestep Dt (Yanenko, 1971). This is easily performed in the linear case and the exten- sion to nonlinear systems can be done with the use of Lie derivatives (see Lanser and Verwer, 1998). Such an approach raises many questions. In which sense can Dt be actually defined as a small parameter ? As it has been pointed out above, the main reason for the use of Operator Splitting is the stiffness of one phenomenon (for instance chemical kinetics). The only practical case is then the choice of a splitting timestep Dt much larger than the smallest physical timescales (let us say e): there is otherwise no computational advantage of using Operator Splitting. This is the coarse case for Operator Splitting. The key point is that such (stiff) systems have a parti- cularly interesting multiscale behaviour (see Larrouturou and Sportisse, 1997; Sportisse and Rouchon, 1998). After a fast transient phase (whose length may be esti-

Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

Embed Size (px)

Citation preview

Page 1: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

Environmental Modelling & Software 15 (2000) 673–679www.elsevier.com/locate/envsoft

Method of Lines versus Operator Splitting for reaction–diffusiuonsystems with fast chemistryq

Bruno Sportissea,*, Guy Bencteuxb, Pierre Plionb

a ENPC-CEREVE, Centre d’Enseignement et de Recherche en Eau, Ville et Environnement, Ecole Nationale des Ponts et Chausse´es(ENPC–CEREVE), rue Blaise Pascal, 77455 Champs sur Marne, France

b DER, A3UR project, Electricite´ de France, 78401 Chatou, France

Abstract

Operator Splitting methods are widely used for the numerical computation of reactive flows. The classical analysis of the errorinduced by such techniques is based on asymptotic expansions with respect to the splitting timestepDt (supposed to be a smallparameter). Such an assumption is unfortunately rather difficult to meet in real computations sinceDt is actually much larger thanthe smallest physical timescales (defined by fast chemical phenomena). The usual error analysis has then to be replaced by anotherapproach, which is associated with the multi-scale behaviour of the coupled system. We give here some theoretical indications andillustrate them with a monodimensional application in Air Pollution Modelling. 2000 Elsevier Science Ltd. All rights reserved.

Keywords:Operator Splitting; Stiffness; Singular perturbation; Reduction of dynamical systems; Reaction–diffusion; Air pollution modelling

1. Introduction

Operator Splitting methods are commonly used in airpollution modelling (McRae et al., 1982; Pu Sun, 1996;Kim and Cho, 1997). The basic idea is to integrate thedifferent physical phenomena one by one. There aremainly two advantages. First, specific tailored (on-the-shelf) solvers can be used for each operator to be inte-grated (e.g. advection, turbulent diffusion, chemical kin-etics, emission processes, aso). This is particularly prag-matic since black boxes can be easily incorporated inexisting softwares. The second advantage is associatedwith CPU requirements. A Method of Lines is indeedbased on the following approach: first discretize spatialoperators and then time integrate (with an ODE solver).As numerical stiffness is mainly induced by chemicalkinetics (due to the wide range of chemical timescales),implicit schemes have to be used for the time integrationof the resulting ODES. This increases drastically theCPU costs due to the necessary algebraic manipulationssince the dimension is given by the product of the num-ber of grid cells by the number of chemical species. Even

q Presented at the Conference APMS ’98, 1998, October 26–29,ENPC and INRIA

* Corresponding author.E-mail address:[email protected] (B. Sportisse).

1364-8152/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved.PII: S1364-8152 (00)00036-0

if some specific structure can be taken into account (inthe case of a diagonal diffusivity matrix, there is only alocal coupling through chemical kinetics: see forinstance Graf and Moussiopoulos, 1991), the compu-tational cost of the Method of Lines (MOL) remains pro-hibitive and using Operator Splitting is recommended.

In return the accuracy may be unfortunately loweredby the uncoupling of operators. The usual error analysisis based on asymptotic expansions of the local error withrespect to the splitting timestepDt (Yanenko, 1971).This is easily performed in the linear case and the exten-sion to nonlinear systems can be done with the use ofLie derivatives (see Lanser and Verwer, 1998).

Such an approach raises many questions. In whichsense canDt be actually defined as a small parameter ?As it has been pointed out above, the main reason forthe use of Operator Splitting is the stiffness of onephenomenon (for instance chemical kinetics). The onlypractical case is then the choice of a splitting timestepDt much larger than the smallest physical timescales (letus saye): there is otherwise no computational advantageof using Operator Splitting. This isthe coarse caseforOperator Splitting.

The key point is that such (stiff) systems have a parti-cularly interesting multiscale behaviour (see Larrouturouand Sportisse, 1997; Sportisse and Rouchon, 1998).After a fast transient phase (whose length may be esti-

Page 2: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

674 B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

mated bye), the evolution proceeds with a slow rate inthe vicinity of a simplified reduced model (up to first-order in e). We will then check that Operator Splittingmethods respect such a dynamical behaviour by comput-ing the reduced solutions associated with Operator Split-ting methods and comparing them with the reduced sol-ution associated with the MOL. On the technical sidewe use a double-limit approach (see Sportisse, 1998):the analysis of the local error is performed withe tendingfirst to zero and thenDt. Let us mention that suchapproaches can be related to other works for the studyof hyperbolic systems comprising a stiff relaxation term(Shi Jin, 1995; LeVeque and Yee, 1990).

This leads to some surprising conclusions for the prac-tical use of Operator Splitting methods. The main pointis that the stiff operator (chemical kinetics) has always tofinish the splitting process (due to its stabilizing effect).

We briefly review Operator Splitting methods in thefirst section. We give some indications illustrating thetheoretical background of the coarse case in the secondsection. A benchmarking study is presented in the lastsection. The numerical tests have been performed in amonodimensional case with diffusion, chemical kinetics,dry deposition and emission.

2. Operator Splitting methods

2.1. Some classical methods

We consider the following linear evolution system:

dzdt

5Az1Bz, z(0)5z0, , zPRn

whereA andB are linear operators. LetDt be the split-ting timestep. Fist-order methods (with respect toDt) areusually defined in the following way:

5dz∗

dt=Az∗, z∗(0)=z0 on [0,Dt]

dz∗∗

dt=Bz∗∗, z∗∗(0)=z∗(Dt) on [0,Dt]

The final value is then taken asz** (Dt). This scheme willbe denoted by (A–B) splitting. (B–A) splitting is ofcourse defined by reversingA andB.

The local error for the (A–B) splitting is easily com-puted through

le5(exp(BDt)exp(ADt)2exp((A1B)Dt))z0

An asymptotic expansion leads straightforward to

le=BA−AB

2Dt2+O(Dt3). We have a first-order global error

with respect toDt unlessA andB commute.Strang (1968) has proposed to symmetrize the split-

ting scheme for improving the accuracy. We will name(A–B–A) the following scheme:

5dz∗

dt=Az∗, z∗(0)=z0 on [0,

Dt2

]

dz∗∗

dt=Bz∗∗, z∗∗(0)=z∗(

Dt2

) on [0,Dt]

dz∗∗∗

dt=Az∗∗∗, z∗∗∗(0)=z∗∗(Dt) on [0,

Dt2

]

where the final value isz∗∗∗(Dt2

). The local error is then

easily computed and an asymptotic expansion leads to athird-oder local error.

2.2. No-time splitting

Another less usual splitting scheme has been recentlyproposed, mainly for Air Pollution Modelling (Pu Sun,1996; Knoth and Wolke, 1994; Sportisse et al., 1997).It is a slight modification of a first-order scheme as oneoperator (let sayA) is supposed to be non-stiff. In orderto avoid transient phases due to stiffness, the initial con-ditions for the second integration step are not modifiedbut an artificial source term is added:

5dz∗

dt=Az∗,z∗(0)=z0 on [0,Dt]

dz∗∗

dt=Bz∗∗+

z∗(Dt)−z0

Dt, z∗∗(0)=z0 on [0,Dt]

We will call this scheme (NTS) The final value is thenz** (Dt). It is equivalent with an explicit integration ofthe (non-stiff) operatorA, which justifies that we do notreverseA andB.

Asymptotic expansions of the local error confirm that(NTS) is a first-order scheme.

3. Coarse time integration

3.1. Why the classical analysis may fail in the stiffcase

In many applications (such as Air Pollution) splittingis advocated since one operator is stiff. From now on

we will make the assumption thatB=ce

is a stiff operator

(for instance chemical kinetics) whileA is a non stiffone (for instance diffusion).e is (as expected) a smallpositive parameter. The above computed local error offirst-order (A–B) or (B–A) schemes is then:

|le|5Dt2

eu(cA2Ac)z0|

Page 3: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

675B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

There are obviously many reasons to suspect that thisanalysis is highly doubful. First, higher order terms maybe of course much larger and cannot be neglected. Sucha local error is an increasing function in the stiffnessratio e (with a fixedDt): this contradicts the widespreadargument that splitting schemes have to be used for wellseparated timescales. The most hazardeous point is asso-ciated with the symmetry between both operators: theschemes (A–B) and (B–A) should have the same behav-iour. The time evolution of aL2 error (Fig. 1: see nextsection for more details) for a 1D Reaction–DiffusionAir Pollution Model illustrates however the importanceof ending with chemistry.

The main difficulty with the previous analysis is the

presence of terms of magnitudeDte

while thecoarse case

DtÀe is commonly met in practice. An alternative isthen to compare the solutions computed with the Oper-ator Splitting methods in the asymptotic casee tendingto zero (which is satisfied up to first-order ine). This isequivalent with comparing the reduced solutions asso-ciated with each methods and checking that OperatorSplitting methods preserve the dynamical behaviour ofthe system.

We refer to Larrouturou and Sportisse (1997), Sport-isse and Rouchon (1998) and Sportisse (1997) for thetheoretical background. An extensive study can be foundin Sportisse (1998) but is out of the scope of this article.We will then only mention the key points through anenlightening example.

3.2. An example

Let us consider the following system:

edxdt

52x1y2ex, edydt

5x2y, x(0)5x0, y(0)5y0

Fig. 1. Evolution of aL2 error (R2) for Dt=900 s.

The identification of A and B (with the previousnotations) is an easy task left to the reader. We computethe reduced solutions (defined up to first order ine) asso-ciated with the MOL and the first-order schemes. Forconvenience we will sytematically solve the stiff step(B) in the lumped basis (Sportisse, 1997) (u=x+y,y). LetDtÀe be the splitting timescale.

O We have for the MOL (in the lumped basis):

dudt

52u1y, edydt

5u22y, u(0)5x01y0, y(0)5y0

The solution can then be approximated (up to first-order in e) by:

dudt

5u2, y5

u2

which leads easily tox(Dt)=y(Dt)=x0+y0

2exp(2

Dt2

).

The algebraic constraintx=y defines the reduced

model whose slow exponential rate, is12

O The first step for the (A–B) scheme is:

dx∗

dt52x∗,

dy∗

dt50, x∗(0)5x0, y∗(0)5y0

which provides the initial conditions for the secondstep:

du∗∗

dt50, e

dy∗∗

dt5u∗∗22y∗∗, u∗∗(0)5x0·e−Dt

1y0, y∗∗(0)5y0

The reduced solution can then be computed as:

x(Dt)5y(Dt)5x0·e−Dt+y0

2

O In the same way we have for the scheme (B–A)

du∗

dt50, e

dy∗

dt5u∗22y∗, u∗(0)5x01y0, y∗(0)5y0

This gives directly the reduced solution

x∗(Dt)=y∗(Dt)=x0+y0

2to be taken as initial conditions

for the second step:

dx∗∗

dt52x∗∗,

dy∗∗

dt5x∗∗(0)5y∗∗(0), 5

x0+y0

2

leading to:

x(Dt)5x0+y0

2e−Dt, y(Dt)5

x0+y0

2

Page 4: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

676 B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

O No-Time Splitting is performed through:

edx∗∗

dt52x∗∗1y∗∗1e

x0(e−Dt−1)Dt

, edy∗∗

dt5x∗∗2y∗∗

by keeping x∗∗(0)=x0 and y∗∗(0)=y0. In a lumpedform:

du∗∗

dt5

x0(e−Dt−1)Dt

, edy∗∗

dt5u∗∗22y∗∗, u∗∗(0)5x0

1y0, y∗∗(0)5y0

which leads easily to

x(Dt)5y(Dt)5x0+y0+x0(e−Dt−1)

2

We perform now an asymptotic expansion of the local(reduced) errors with respect toDt. Easy calculationsgive:

le(A−B)5Dt4

(x02y0)1O(Dt2), le(B−A)5Dt4

(x01y0)1O(Dt2)

leNTS5Dt4

(x02y0)1O(Dt2)

This seems to be a rather disappointing result sincewe have only first-order local errors (on the contrary tothe usual analysis). The former study has however to beperformed with (x0, y0) being the current value of thesplitting process. The key property of splitting (A–B) and(NTS) is that (x0, y0) (except for the first iteration) satis-fies a constraint: the computed solution is in the vicinityof the reduced modelx=y [on the contrary to the splitting(B–A)]. That is to say that (outer a transient phase) wehave (up to first-order ine) x0=y0 for the splitting (A–B)and (NTS). This justifies that we recover the usual orderfor both methods:

le(A−B)5O(Dt2), leNTS5O(Dt2)

What is the reason why we have lost accuracy for (B–A)? This is rather surprising since the dynamical behav-iour justifies that the initial slow–fast model can beapproximated up to some orders ine by integrating suc-cessively the fast and then the slow operators. Thescheme (B–A) would then seem to be more logical. Letus go back to the computation of the reduced solutionfor the MOL:

O we first integrate the fast operator on [0,Dt]:

edx∗

dt52x∗1y∗, e

dy∗

dt5x∗2y∗

This is exactly the first step of the splitting scheme(B–A) and that explains why we have ‘good’ initialconditions after the first step.

O we integrate then the slow operator on [0,Dt] by forc-ing the algebraic constraint:

du∗∗

dt52

u∗∗

2, 5y∗∗5x∗∗

or in the initial basis:

dx∗∗

dt52

x∗∗

2, 5

dy∗∗

dt5

x∗∗

2

This has to be compared with the second step of the(B–A) scheme:

dx∗∗

dt52x∗∗, 5

dy∗∗

dt50

which for the solution does not actually satisfy thealgebraic constraint.

In few words it is better to loose accuracy at the begin-ning and then project the evolution onto the simplifiedmodel than to be close to the exact solution after thefirst step but without being close to the reduced modelthereafter.

Such results can be easily extended to second-ordermethods (Sportisse, 1998). The method ending with thestiff operator is always more accurate than the other one.Moreover they are only first-order methods at best(unless some uncoupling conditions are met).

We will conclude with a rather confusing remark:theaccuracy order does actually depend on the species. Thelumped (slow) speciesu=x+y is indeed:

u(Dt)5x0+y0

2(11e−Dt)

for the poorly accurate scheme (B–A) while the exactsolution is:

u(Dt)5(x01y0)e−Dt2

A direct computation gives a second-order accuracyfor the local error associated withu! This illustrates howdifficult the estimation of the accuracy order is: crudelyspeaking, we can have the usual order for slow species(since there is no effect of stiffness for them) while theorder is lowered for fast species.

Remark 2.1 (Dealing with boundaryconditions).Boundary conditions have also to be handledin practical cases. This raises many numerical difficult-ies. We will only give here some indications and refer toa future work. LetC denote a constant term (representingboundary conditions). We have then to solve by split-ting:

Page 5: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

677B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

dzdt

5Az1Bz1C

The main question is of course to determine to whichoperatorC has to be associated. A rigorous answer canbe done (see Hunsdorfer, 1996 for instance) but is unfor-tunately difficult to use in practice. We focus once moreon the case whereB is a stiff operator. A similar studyto the previous one indicates that it is better to associateC with the non stiff operator rather than with the stiffone.

3.3. Summary of the main results

We summarize now the main conclusions induced bythis example:

O the usual accuracy order of Operator Splittingmethods is not systematically preserved,

O the stiff operator (chemical kinetics) has always tofinish the splitting process since it has a stabilizingeffect (it projects the evolution onto the reducedmodel),

O it is recommended to include boundary conditions(e.g. dry deposition and emissions) with diffusion(which is physically meaningful),

O No Time Splitting methods are robust since the initialconditions always satisfy the reduced model,

O accuracy is better for slow species than for fast ones.

4. Some numerical tests in air pollution modelling

We will now focus on coupling between chemistryand diffusion. One field of application is air pollutionmodelling. The main errors induced by Operator Split-ting are indeed issued from the uncoupling betweenchemistry and diffusion (see Lanser and Verwer, 1998for a systematic study) and from the treatment of bound-ary conditions.

We have simulated a monodimensional reaction–dif-fusion system with data given in Kim and Cho (1997).The boundary conditions are dry deposition and emis-sions at the ground and no flux at the top. We have usedLSODE in order to perform a 3-h time integration.

A reference solution (ZMOL) has been computed withthe Method of Lines. IfzS is the solution given by asplitting method, we define the usual error norms as fol-lows:

ERRi2(n,m)5

ziMOL(n·m)−zi

S(n·m)zi

MOL(n·m)+atol

Fig. 2. Final relative errorR2 (in %).

Fig. 3. Final relative errorR2 (in %). Dt=900 s.

Ri2(n)!1

HS

m5M

m51ERRi

2(n,m)2.dx(m), R2(n)51I

Si5I

i51Ri

2(n)

wherei, n andm stand respectively for the species, iter-ation and grid cell index (whose whole numbers areI,N and M), atol is a tolerance parameter, dx(j) is thelength of thej grid cell andH is the whole vertical height(|1200 m).

There are 10 vertical cells (which is a realisticnumber). The turbulent diffusion coefficient is constant(k=5 m2 s-1).

We have studied the influence of the sequential orderof operators and of the treatment of boundary conditions(BC) (included either in the diffusion term or in thechemical term) by performing an exhaustive comparativestudy between the following seven Operator Splittingmethods: First-order Splitting (3), Strang Splitting (3)and No Time Splitting (1) (the first-order and second-order schemes which for boundary conditions are solvedwith chemistry and diffusion ends the process have notbeen considered) (Figs. 2 and 3).

The splitting timesteps have been successively takenas 100, 300 and 900 (the usual case). Figs. 4 and 5 Fig.6 give the time evolution of the global errorR2 for eachtimestep. The final values ofR2 can be found in Fig. 2(boundary conditions have been included in the operator

Fig. 4. Dt=100 s.

Page 6: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

678 B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

Fig. 5. Dt=300 s.

Fig. 6. Dt=900 s.

with a star). These results confirm the previous con-clusions: reversing the order of operators can inducehigh errors while boundary conditions have to beincluded with the diffusion term. Strang splitting is onlya first-order scheme (even if it is the most accuratescheme) and No Time Splitting has a good accuracy.

The smaller the splitting timestep is, the more theusual analysis may be recovered. Fig. 3 confirms the lossof accuracy for the first-order splitting ending withchemistry (C+D(BC)). Errors are low for (slow) lumpedspecies, such as NOx=NO+NO2 and Ox=NO2+O3, whilethey are large for fast species (OH, HO2 and HO2NO2).

We have performed another sensitivity analysis withrespect to the grid refinement. The numerical stiffnessassociated with diffusion increases indeed with gridrefinement. We expect that the more the grid is refinedthe more the usual analysis is valid (since the timescalesof diffusion and chemical kinetics become similar). This

Fig. 7. Grid refinement for splitting chemistry+diffusion withDt=100 s.

is confirmed by the results shown in Fig. 7: three gridshave been used (the initial one, a second grid with 20cells and a third grid with 40 cells) and the first-orderschemes have been tested.

5. Conclusion

The usual analysis for splitting errors fails is the stiffcase (met in practice). We have proposed an alternativeanalysis based on the dynamical behaviour of each split-ting scheme. This leads to advocate some special usesof classical methods. Such an analysis is confirmed bynumerical tests applied to air pollution modelling(reaction–diffusion). We now investigate high-ordersplitting methods in this framework.

Acknowledgements

This work has partially been supported by Electricite´de France at CERMICS. It has been done at EdF in theframework of the A3UR project.

References

Graf, J., Moussiopoulos, N., 1991. Intercomparison of two models forthe dispersion of chemically reacting pollutants. Beitrage PhysikAtmosphare 64 (1), 13–25.

Hundsdorfer, W.H., 1996. Numerical solution of advection–diffusion–reaction equations. Technical Report NM-N9603, CWI.

Kim, J., Cho, S.Y., 1997. Computation accuracy and efficiency of thetime-splitting method in soloving atmospheric transport-chemistryequations. Atmosphere and Environment 31 (15).

Knoth, O., Wolke, R., 1994. A comparison of fast chemical kineticsolvers in a simple vertical diffusion model. In: Air Pollution Mod-elling and its Applications. Plenum Press, New York.

Page 7: Method of Lines versus Operator Splitting for reaction–diffusiuon systems with fast chemistry

679B. Sportisse et al. / Environmental Modelling & Software 15 (2000) 673–679

Lanser D., Verwer, J.G., 1998. Analysis of operator splitting for advec-tion–diffusion–reaction problems from air pollution modelling. In:Proceedings 2nd Meeting on Numerical methods for differentialequations. Coimbra, Portugal, February.

Larrouturou, B., Sportisse, B., 1997. Some mathematical and numeri-cal aspects of reduction in chemical kinetics. In: ComputationalScience for the 21st Century, Conference in Honor of Professor R.Glowinski. J. Wiley, pp. 2215–2224.

LeVeque, R.J., Yee, H.C., 1990. A study of numerical methods forhyperbolic conservation laws with stiff source terms. Journal ofComparative Physics 86, 187–210.

McRae, G.J., Goodin, W.R., Seinfeld, J.H., 1982. Numerical solutionof the atmospheric diffusion equation for chemically reacting flows.J.C.P. 45, 1–42.

Pu Sun, 1996. A pseudo non-time splitting method in air quality mode-ling. Journal of Comparative Physics 127, 152–157.

Shi Jin, 1995. Runge–Kutta methods for hyperbolic conservation lawswith stiff relaxation terms. Journal of Comparative Physics 122,51–67.

Sportisse, B., 1997. Reducing chemical kinetics: a mathematical pointof view. In: Proceedings of the Workshop Numerical Aspects ofReduction in Chemical Kinetics, ENPC–CERMICS, September.

Sportisse, B., 1998. An analysis of operator splitting techniques in thestiff case. Technical Report, CERMICS, July 1998, 98–127.

Sportisse, B., Rouchon, P., 1998. Reduction of slow–fast chemistrywith slow processes. Technical Report, CERMICS, July 1998,98–129.

Sportisse, B., Jaubertie, A., Plion, P., 1997. Reducing mechanisms inchemical kinetics for the simulation of reactive transport; an appli-cation to air pollution modelling. In: Proceedings of the St VenantSymposium, August.

Strang, G., 1968. On the construction and comparison of differenceschemes. SIAM Journal of Numerical Analysis 5, 506–517.

Yanenko, N.N., 1971. The Method of Fractional Steps. Springer,New York.