Bruno

Download as pdf or txt
Download as pdf or txt
You are on page 1of 77

Research Habilitation

University Pierre and Marie Curie

Defended on 4 June 2007 by

Bruno Sportisse
Field : Geophysics

Air Pollution
Modeling and Simulation

Jury composed of

Matthias Beekmann CNRS reviewer


Greg Carmichael University of Iowa, USA reviewer
Katia Laval University Pierre and Marie Curie President
François-Xavier Le Dimet University Joseph Fourier reviewer
Jan Verwer CWI, The Netherlands
2
3

Acknowledgements

I thank Katia Laval for having accepted to chair the jury.

I would like to thank François-Xavier Le Dimet for his constant help: his report about my
work was only a small part of the large amount of sollicitations to him.

Thanks to Matthias Beekmann for having written a report about my work. I hope that an
improved knowledge of my (our) activities will be an opportunity for joint works in the future.

I have been honored by the report of Greg Carmichael and by his presence at my defense in
Paris.

Thanks to Jan Verwer for having attended to my defense, even if my research works have
switched from simulation to modeling.

Thanks to Alain Saliot for his help related to my defense.

I would like to thank my colleagues for many joint works, especially Marc Bocquet, Jaouad
Boutahar, Édouard Debry, Rafik Djouad, Olivier Isnard, Vivien Mallet, Luc Musson-Genon,
Denis Quélo, Yelva Roustan, Karine Sartelet and Marilyne Tombette.
4
Contents

1 Introduction 7

2 Bibliography of my research works 11


2.1 Parameterizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Numerical analysis for reactive dispersion . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Modeling and simulation of aerosol dynamics . . . . . . . . . . . . . . . . . . . . 12
2.4 Data assimilation and inverse modeling . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 Sensitivity analysis, uncertainties, ensemble forecast . . . . . . . . . . . . . . . . 14
2.5.1 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5.2 Incertainties and ensemble forecast . . . . . . . . . . . . . . . . . . . . . . 15
2.6 Model reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.7 CTM and applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.8 Books . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.9 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 Research works 19
3.1 Parameterizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1 Turbulent diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.2 Scavenging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.3 Segregation effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.4 Plume In Grid models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.5 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Numerical analysis for reactive dispersion . . . . . . . . . . . . . . . . . . . . . . 25
3.2.1 Splitting methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.2 Time integration of chemical kinetics . . . . . . . . . . . . . . . . . . . . . 27
3.2.3 Advection and mass consistency . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.4 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3 Aerosol modeling and simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.3.1 The General Dynamics Equation for aerosols . . . . . . . . . . . . . . . . 32
3.3.2 Stochastic methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3.3 Modal methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3.4 Size-resolved methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.5 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4 Inverse modeling and data assimilation . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4.2 Sequential versus variational methods . . . . . . . . . . . . . . . . . . . . 42
3.4.3 Coming issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4.4 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.5 Sensitivity analysis, ucertainty propagation and ensemble forecast . . . . . . . . . 48

5
6 CONTENTS

3.5.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.5.2 Some possible strategies for assessing the impact of uncertainties . . . . . 49
3.5.3 Multi-modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.5.4 Ensemble forecast . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.5.5 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.6 Proxy and reduced models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.6.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.6.2 Slow/fast models and dynamic reduction . . . . . . . . . . . . . . . . . . 52
3.6.3 Building look-up tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6.4 POD/EOF methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6.5 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7 Chemistry-Transport Models and applications . . . . . . . . . . . . . . . . . . . . 55
3.7.1 Development of CTM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7.3 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Chapter 1

Introduction

Chemistry Transport Models


A precise description of the chemical composition of the atmosphere is required for many en-
vironmental issues: at local scale for risk assessment or urban air quality, at regional scale for
transboundary pollution and at global scales for the study of the climate change, the study of
the atmospheric oxidizing or the stratospheric ozone hole.
There is a hierarchy of available models, supposed to applied in many different applications
(Table 1.1):

• in an operational context (for instance for risk assessment of chemical, biological or radio-
logical releases), the model has to be robust and will be used with a small amount of input
data;

• for environmental forecast (urban air quality and photochemistry), the model will be fitted
in order to provide good skills for target species;

• for impact studies (or for process studies), a comprehensive model describing the whole
processes will be used in order to be runned in many conditions (far from the current
condition).

Table 1.1: A model hierarchy.

Gaussian industrial risk


Eulerian box forecast 80s
Lagrangian box impact 80s-90s
3D Eulerian, gas-phase Chemistry-Transport (CTM) early 90s
3D Eulerian, multiphase Chemistry-Transport (CTM) late 90s-2000
Coupled models (meteo/CTM) 2000s ?

Air pollution modeling is based on the so-called Chemistry Transport Models (CTM here-
after). The core of a CTM (see Seinfeld and Pandis [1998]) is the simulation of the time evolu-
tion of chemical concentrations. The forcing conditions (meteorological conditions, emissions)
are fixed. Let ci stand for the concentration of a chemical gas-phase species i (in mol l −1 , ppb,
µg m−3 ...). The evolution is then governed by a reaction-diffusion-advection equation (equation

7
8 Chapitre 1 - Introduction

gaseous chemistry

advection VOC Ozone


aqueous
chemistry
NOx
SVOC
activation

nucleation condensation

H2SO4 evaporation
rainout

HNO3
coagulation

SO2 NOx
turbulent VOC washout
diffusion
VOC VOC
NOx primary dust
particles
NH3

dry deposition
TRAFFIC AGRICULTURE INDUSTRY BIOGENIC

Figure 1.1: Processes described in a CTM.

of reactive dispersion):
 
∂ci ci
+ div (V ci ) = div ρK∇ + χi (c, T, I) − Λi (x, t)ci + Si (x, t) (1.1)
∂t | {z } ρ | {z } | {z } | {z }
advection | {z } gas phase chemistry scavenging volumic emission
turbulent diffusion

~ is the wind velocity, K the eddy diffusion coefficient, T the temperature, I the solar actinic
V
flux and ρ the air density. These variables are supposed to be known. χi stands for the chemical
sources and depends on the solar flux (through photolysis) and on temperature (through thermal
reactions). Λ is a parameterization of scavenging processes (by rain for instance). S i describes
the point sources.
Some boundary conditions are specified, for instance at ground for emission and deposition:

∂ci
−Kz = Ei (~x, t) − videp (~x, t)ci (1.2)
∂z | {z } | {z }
surfacic emission dry deposition

Kz is the vertical eddy coefficient, Ei is the flux of surfacic emission and videp is the dry deposi-
tion velocity.

A key component of these models is the computation of the chemical sources (χ i ). Gas-
phase chemical kinetics is well understood for some chemical species. Current comprehensive
gas-phase mechanisms describe up to hundreds of chemical species. In practice, only surrogate
lumped species (the so-called model species) are taken into account in CTM.
9

Table 1.2: Dispersion of scales.

Horizontal grid cell [1-100] kilometers


Characteristic chemical timescales [10−6 second-year]
Cloud droplet [10-100] micrometers
Rain droplet 0.1 millimeter
Aerosols [1 nanometer-10 micrometer]

The condensed matter (in liquid or solid forms, Table 1.2 for the typical scales) has a key
role for many processes. One usually distinguishes the condensed matter related to water (rain
drop, cloud drop) and the aerosols (particles in suspension composed of a mixture of chemical
species). The gas-phase models have therefore to be extended to multiphase processes by de-
scribing the coupling to cloud droplets (diphasic models: gas/cloud) and the aerosol dynamics.
Source terms have then to be added to Equation (1.1).

A CTM is then based on the following components:

• a preprocessing model (emissions, meteorological fields);

• a description of physical and chemical processes (gas-phase mechanism, multiphase model);

• some appropriate parameterizations for unresolved processes;

• some numerical algorithms.

A CTM can then be used for simulating the chemical composition of the atmosphere (forward
mode). It can be used for forecast by using data assimilation methods (coupling to observational
data). Another use is the identification of some poorly known input data (inverse modeling of
emissions for instance). The assessment of the uncertainties is another key point that can be
performed through multi-model methods or ensemble runs.

Outlines
My research activity has been led in this framework. I have focused on the following topics:

1. Parameterization of subgrid processes:


How to parameterize processes that are not solved at the numerical scale ? A typical
example is given by scavenging.

2. Numerical simulation of the dispersion equation:


Once the physical model is fixed, which are the most appropriate numerical algorithms
required for the simulation of reactive dispersion ? The challenging issues are related to
the system dimension and to the wide range of timescales.

3. Modeling and simulation of the aerosol dynamics:


Which are the parameterizations to be used for aerosol dynamics ? Which are the numeri-
cal algorithms for aerosol dynamics (a system of non-linear integro-differential equations)?
10 Chapitre 1 - Introduction

4. Data assimilation and inverse modeling:


How to apply the usual data assimilation methods to CTM ? Which are the specicifities
of CTM ? What is the robustness (with respect to uncertainties) of inverse modeling of
emissions ?

5. Treatment of uncertainties:
Which are the uncertainties of CTM ? How to extend mono-model runs to multi-model
runs ?

6. Model reduction:
How to deal with the high dimension ? How to build reduced models, especially for
integrated modeling ?

7. Development, validation and application of CTM:


What is a “good” modeling system for air quality ? How to build a system devoted to
different scales (from local to continental scales), different kinetics (photochemistry, heavy
metals, aerosols, radionuclides) and different uses (forecast, data assimilation, inverse mod-
eling, impact studies) ?

This manuscript is organized as follows. Chapter 2 gives a bibliography of my works. Chapter


3 is a summary of my research activities. The main stream is given by the answer to the question
“How a CTM is built and which are the limitations ?”. I have therefore made the choice of paying
attention to the 6 first topics (with a methodological focus), the last point (applications of CTM)
being only briefly illustrated.
Chapter 2

Bibliography of my research works

My research works are detailed in the following articles. I have omitted technical reports (related
to my application works especially with industrialists) except those with methodological points
or interesting case studies.

2.1 Parameterizations
Articles
Sportisse, B. (2001b). Box models versus Eulerian models in air pollution modeling. Atmos.
Env., 35:173–178

Sportisse, B. and Du Bois, L. (2002). Numerical and theoretical investigation of a simplified


model for the parameterization of below-cloud scavenging by falling raindrops. Atmos. Env.,
36:5719–5727

Sportisse, B. and Djouad, R. (2003). Mathematical investigation of mass transfer for atmo-
spheric pollutants into a fixed droplet with aqueous chemistry. J. Geophys. Res., 108(D2):4073

Sportisse, B. (2007d). A review of parameterizations for modeling dry deposition and scavenging
of radionuclides. Atmos. Env., (41):2683–2698

2.2 Numerical analysis for reactive dispersion


Articles
Sportisse, B. (2000). An analysis of operator splitting techniques in the stiff case. J. Comp.
Phys., 161:140–168

Sportisse, B., Bencteux, G., and Plion, P. (2000). Method of Lines versus Operator Splitting
methods for reaction-diffusion systems in air pollution modelling. Env. Mod. Soft., 15(6-7):673–
679

Djouad, R., Sportisse, B., and Audiffren, N. (2002). Numerical simulation of aqueous-phase
atmospheric models : use of a non-autonomous Rosenbrock method. Atmos. Env., 36:873–879

11
12 Chapitre 2 - Bibliography of my research works

Djouad, R. and Sportisse, B. (2003). Solving reduced models in air pollution modelling. Appl.
Numer. Math., 44(1):49:61

Sportisse, B., Quélo, D., and Mallet, V. (2007). Impact of mass consistency errors for atmo-
spheric dispersion. Atmos. Env., 41:6132–6142

Mallet, V., Pourchet, A., Quélo, D., and Sportisse, B. (2007b). Investigation of some numerical
issues in a Chemistry-Transport Model: gas-phase simulations. J. Geophys. Res., 112(D15301).
doi:10.1029/2006JD008373

Others (chapter, teaching notes, proceedings, reports)


Sportisse, B. and Djouad, R. (2002). Some aspects of multi-timescales issues for the numerical
modeling of atmospheric chemistry, volume 130, chapter 4, page 39:61. IMA Volume in Mathe-
matics and its applications. Springer New York

Verwer, J. and Sportisse, B. (1998). A note on operator splitting analysis in the stiff linear case.
Technical Report MAS-R9830, CWI

Sportisse, B. and Mallet, V. (2005). Calcul Scientifique pour l’Environnement. Cours ENSTA.
86 pages

2.3 Modeling and simulation of aerosol dynamics


Articles
Debry, E., Sportisse, B., and Jourdain, B. (2003). A stochastic approach for the numerical sim-
ulation of the General Dynamics Equation for aerosols. J. Comp. Phys., 184:649:689

Sartelet, K. N., Hayami, H., Albriet, B., and Sportisse, B. (2005). Development and preliminary
validation of a Modal Aerosol Model for tropospheric chemistry: MAM. Aerosol Sci. and Tech-
nol., 40(2):118–127

Debry, E. and Sportisse, B. (2007a). Numerical simulation of the General Dynamics Equation
(GDE) for aerosols with two collocation methods. Appl. Numer. Math., 57(8):885–898

Debry, E. and Sportisse, B. (2007b). Solving aerosol coagulation with size-binning methods.
Appl. Numer. Math., 47:1008–1020

Debry, E., Fahey, K., Sartelet, K., Sportisse, B., and Tombette, M. (2007). A new SIze REsolved
Aerosol Model: SIREAM. Atmos. Chem. Phys., 7(6):1537–1547

Tombette, M. and Sportisse, B. (2007). Aerosol modeling at regional scale: Model-to-data com-
parison and sensitivity analysis over Greater Paris. Atmos. Env. In press

Sartelet, K., Debry, E., Fahey, K., Tombette, M., Roustan, Y., and Sportisse, B. (2007a). Simula-
tion of aerosols and gas phase species over Europe with the Polyphemus system. Part I: model-
to-data comparison for year 2001. Atmos. Env., 41:6116–6131. doi:10.1016/j.atmosenv.2007.04.024
Data assimilation and inverse modeling 13

Sartelet, K., Hayami, H., and Sportisse, B. (2007c). MICS-Asia Phase II: sensitivity to the
aerosol module. Atmos. Env. doi:10.1016/j.atmosenv.2007.03.005

Others (chapter, teaching notes, proceedings, reports)


Debry, E., Jourdain, B., and Sportisse, B. (2001). Modeling of aerosol dynamics: a stochastic
algorithm. In Sportisse, B., editor, APMS 2001, Geosciences, pages 308–319. Springer

Fahey, K., Debry, E., Foudhil, H., and Sportisse, B. (2005). Formulation, development and pre-
liminary validation of the SIze REsolved Aerosol Model, SIREAM. In Proceedings GLOREAM
2004, pages 7–17

Sportisse, B., Debry, E., Fahey, K., Roustan, Y., Sartelet, K., and Tombette, M. (2006b). Rap-
port du projet Primequal-Predit PAM (Pollution Atmosphérique Multiphasique). Modélisation.
Technical Report 2006-1, CEREA. 171 pages

Sportisse, B., Sartelet, K., Debry, E., Fahey, K., Roustan, Y., and Tombette, M. (2006c). The
SIze REsolved Aerosol Model (SIREAM) and the Modal Aerosol Model (MAM). Technical
documentation. Technical Report 2006-8, CEREA. 115 pages

2.4 Data assimilation and inverse modeling


Articles
Sportisse, B. and Quélo, D. (2002). Observational network design, adaptive observations and
targeting strategies in atmospheric data assimilation: sketch of a methodological review. In
INRIA-CEA-EDF Summer School on Data Assimilation for Geophysical flows

Sportisse, B. and Quélo, D. (2003). Data assimilation and inverse modeling of atmospheric
chemistry. Proc. of Indian National Science Academy. Part A Physical Sciences, 69

Quélo, D., Sportisse, B., and Isnard, O. (2005b). Data assimilation for short-range dispersion of
radionuclides: a case study for second-order sensitivity. J. Environ. Radioactivity, 84:393–408

Quélo, D., Mallet, V., and Sportisse, B. (2005a). Inverse modeling of NO x emissions at regional
scale over Northern France. Preliminary investigation of the second-order sensitivity. J. Geo-
phys. Res., 110(D24310)

Krysta, M., Bocquet, M., Sportisse, B., and Isnard, O. (2006). Data assimilation for short-range
dispersion of radionuclides: an application to wind tunnel data. Atmos. Env., 40(38):7267–7279

Bocquet, M. and Sportisse, B. (2007). Modélisation inverse pour la qualité de l’air : éléments
de méthodologie et exemples. Pollution Atmosphérique. Accepted for publication

Lacour, S. and Sportisse, B. (2007). Estimation of indoor deposition velocity for ozone with a
simplified reactive box model. Atmos. Env. In revision

Sportisse, B. and Mallet, V. (2007). Data assimilation and inverse modelling of slow/fast atmo-
spheric chemistry. Journal of Computational Geosciences. Submitted
14 Chapitre 2 - Bibliography of my research works

Others (chapter, teaching notes, proceedings, reports)

Quélo, D., Sportisse, B., Berroir, J., and Charpentier, I. (2002). Some remarks concerning in-
verse modeling and data assimilation for slow-fast atmospheric chemical kinetics. In Sportisse,
B., editor, APMS 2001, Geosciences, pages 499–513. Springer

Sportisse, B. (2004). Assimilation de données. I Eléments théoriques. Technical Report 2004-24,


CEREA. Notes de cours ENSTA et DEA M2SAP

Sportisse, B. and Quélo, D. (2004). Assimilation de données. II Implémentation des approches


variationnelles: le cas d’un modèle de Chimie-Transport. Technical Report 2004-25, CEREA.
Notes de cours ENSTA et DEA M2SAP

Krysta, M., Bocquet, M., Isnard, O., Issartel, J.-P., and Sportisse, B. (2004). Data assimilation
of radionuclides at short and regional scales. In Farago, I., Georgiev, K., and Havasi, A., editors,
Proceedings of the NATO Advanced Research Workshop on Advances in Air Pollution Modeling
for Environmental security, pages 275–285. Kluiwer

Isnard, O., Krysta, M., Bocquet, M., Dubiau, P., and Sportisse, B. (2005). Data assimilation of
radionuclides atmospheric dispersion at small scale: a tool to assess the consequences of radio-
logical emergencies. In Proceedings of the IAEA Conference. Rio Conference

Korsakissok, I. and Sportisse, B. (2006). Simulations et évaluation de réseaux pour la détection


de traceurs à l’échelle locale. Technical Report 2006-46, CEREA. Rapport de contrat DGA

2.5 Sensitivity analysis, uncertainties, ensemble forecast

2.5.1 Sensitivity analysis

Articles

Djouad, R., Audiffren, N., and Sportisse, B. (2003a). Sensitivity analysis using automatic dif-
ferentiation applied to a multiphase chemical mechanism. Atmos. Env., 37(22):3029:3038

Mallet, V. and Sportisse, B. (2004). 3D Chemistry-Transport Model Polair3D: numerical issues,


validation and automatic-differentiation strategy. Atmos. Chem. Phys. Discuss., (1):1371:1392

Mallet, V. and Sportisse, B. (2005). A comprehensive study of ozone sensitivity with respect to
emissions over Europe with a chemistry-transport model. J. Geophys. Res., 110(D22)

Others (chapter, teaching notes, proceedings, reports)

Tombette, M., Fahey, K., Sartelet, K., and Sportisse, B. (2005). Aerosol modelling at regional
scale: a sensitivity study with the Polyphemus platform. In Proceedings of GLOREAM 2005.
Also as Report CEREA 2005-50
Model reduction 15

2.5.2 Incertainties and ensemble forecast


Articles
Mallet, V. and Sportisse, B. (2006d). Uncertainty in a chemistry-transport model due to phys-
ical parameterizations and numerical approximations: an ensemble approach applied to ozone
modeling. J. Geophys. Res., 111(D01302)

Mallet, V. and Sportisse, B. (2006b). Ensemble-based air quality forecasts: a multi-model ap-
proach applied to ozone. J. Geophys. Res., 111(D18):18302

Mallet, V. and Sportisse, B. (2006a). Air quality modeling: from deterministic to stochastic
modeling. Computers and Mathematics with Application. In press

Others (chapter, teaching notes, proceedings, reports)


Mallet, V. and Sportisse, B. (2006c). Peut-on modéliser la qualité de l’air de manière déterministe
? In Actes 38ièmes Journées de Statistiques. Société Française de Statistique

2.6 Model reduction


Articles
Sportisse, B. and Djouad, R. (2000). Reduction of chemical kinetics in air pollution modelling.
J.Comp.Phys., 164:354–376

Djouad, R. and Sportisse, B. (2002a). Partitioning techniques for reduction in chemical kinetics.
APLA: an Automatic Partitioning and Lumping Algorithm. Appl. Numer. Math., 43(4):383:398

Djouad, R., Sportisse, B., and Audiffren, N. (2003b). Reduction of multiphase atmospheric
chemistry. Journal of Atmospheric Chemistry, 46:131–157

Debry, E. and Sportisse, B. (2006). Reduction of the condensation/evaporation dynamics for


atmospheric aerosols: theoretical and numerical investigation of hybrid methods. J. Aerosol Sci.,
37(8):950–966

Sportisse, B. and Djouad, R. (2007). Use of Proper Orthogonal Decompositions for the reduction
of atmospheric chemistry. J. Geophys. Res., 112(D06303)

Others (chapter, teaching notes, proceedings, reports)


Rouchon, P. and Sportisse, B. (September 1997). Slow and fast kinetic scheme with slow dif-
fusion. In Proceedings of the Workshop Numerical Aspects of Reduction in Chemical Kinetics,
ENPC-CERMICS.

Sportisse, B. (September 1997). Reducing chemical kinetics: a mathematical point of view.


In Proceedings of the Workshop Numerical Aspects of Reduction in Chemical Kinetics, ENPC-
CERMICS.

Larrouturou, B. and Sportisse, B. (1997). Some mathematical and numerical aspects of reduc-
tion in chemical kinetics. In Bristeau, M. and al, editors, Computational Science for the 21st
16 Chapitre 2 - Bibliography of my research works

century, conference in honor of Professor R.Glowinski. John Wiley and Sons

Sportisse, B., Jaubertie, A., and Plion, P. (1997). Reducing mechanisms in chemical kinetics for
the simulation of reactive transport; an application to air pollution modelling. In Proceedings
of the St Venant Symposium

Jaubertie, A., Plion, P., Sportisse, B., Aumont, B., and Toupance, G. (1998). Reduction of
kinetic schemes for atmospheric chemistry. In Proceedings of the Conference APMS98. INRIA,
France

Djouad, R. and Sportisse, B. (2002b). Some reduction techniques for simplifying air pollution
models. In Sportisse, B., editor, APMS 2001, Geosciences, pages 235–245. Springer

Boutahar, J. and Sportisse, B. (2002). Reduction methods for atmospheric chemistry. In Global
and Regional Atmospheric Modeling Meeting, pages 133–139. Un. of Aveiro, Portugal

Boutahar, J. and Sportisse, B. (2005). Reduction methods and uncertainty propagation: appli-
cation to a Chemistry-Transport-Model. In Proceedings of the TAM-TAM conference

2.7 CTM and applications


Articles
Sportisse, B., Boutahar, J., Debry, E., Quélo, D., and Sartelet, K. (2002). Some tracks in air
pollution modeling. RACSAM Journal of the Spanish Science Academy, Real Academia de Cien-
cias de Espana, 96(3)

Sportisse, B. (2003b). Quelques aspects de la modélisation de la pollution atmosphérique. An-


nales des Ponts et Chaussées, 107-108:44–55

Boutahar, J., Lacour, S., Mallet, V., Quélo, D., Roustan, Y., and Sportisse, B. (2004). De-
velopment and validation of a fully modular platform for numerical modelling of air pollution:
Polair3D. Int. J. Env. and Pollution, 22(1-2)

Roustan, Y., Bocquet, M., Musson-Genon, L., and Sportisse, B. (2006). Modélisation du mer-
cure, du plomb et du cadmium à l’échelle du continent européen. Pollution Atmosphérique,
(191):317–327

Mallet, V., Korssakissok, I., Quélo, D., and Sportisse, B. (2007a). Polyphemus : un système
modulaire de modélisation pour la dispersion atmosphérique et l’évaluation des risques. Pollu-
tion Atmosphérique

Mallet, V., Quélo, D., Sportisse, B., Ahmed de Biasi, M., Debry, É., Korsakissok, I., Wu, L.,
Roustan, Y., Sartelet, K., Tombette, M., and Foudhil, H. (2007c). Technical Note: The air qual-
ity modeling system Polyphemus. Atmos. Chem. Phys. Discuss., 7(3):6,459–6,486

Quélo, D., Krysta, M., Bocquet, M., Isnard, O., Minier, Y., and Sportisse, B. (2007). Validation
of the polyphemus system: the ETEX, Chernobyl and Algeciras cases. Atmos. Env., 41:5300–
Books 17

5315

Sportisse, B. (2007c). A review of current issues in air pollution modeling and simulation. Com-
putational Geosciences, 11(2):159–181

Sartelet, K., Hayami, H., and Sportisse, B. (2007b). Dominant aerosol processes during high-
pollution episodes over Greater Tokyo. J. Geophys. Res. Accepted for publication

Others (chapter, teaching notes, proceedings, reports)


Lacour, S., Megueulle, C., Carlotti, P., Soulhac, L., and Sportisse, B. (2003). Unsteady effects
on pollutant dispersion around a tunnel portal. In Proceedings of Urban Air Quality Conference
UAQ4, Prague

Quélo, D., Lagache, R., and Sportisse, B. (2004). Etude de l’impact qualité de l’air des scénarios
PDUs sur Lille à l’aide du modèle de Chimie-Transport Polair3D. Technical Report 2004-28,
CEREA. Rapport PREDIT

Taghavi, M., Musson-Genon, L., and Sportisse, B. (2004). Impact des rejets de la centrale
thermique de Martigues sur la qualité de l’air. Modélisation avec Polair3D. Technical Report
2004-27, CEREA. Rapport pour la Mission Thermique EDF

Roustan, Y., Bocquet, M., Musson-Genon, L., and Sportisse, B. (2005). Modeling of mercury at
the European scale with the chemistry-transport model Polair3D. In Proceedings GLOREAM
2004, pages 121–130. Also as Technical Report CEREA 2005-04

Tombette, M. and Sportisse, B. (2006). Rapport de contrat préliminaire: évaluation des impacts
du CPT Porcheville. Technical Report 2006-55, CEREA. Rapport pour la Mission Thermique
EDF

Lagache, R., Declercq, C., Quélo, D., Sportisse, B., Palmier, P., Quetelard, B., and Haziak, F.
(2006). Evaluation du PDU de Lille-Métropole sur le trafic, les concentrations de polluants at-
mosphériques et la mortalité. In Actes de la 15ième conférence sur les transports et la pollution
de l’air

Roustan, Y., Lasry, F., and Sportisse, B. (2007). Modélisation de l’impact des émissions des
centrales EDF et SNET en France sur le transport atmosphérique des PM 10 et PM2.5 . Technical
Report 2007-9, CEREA

2.8 Books
As Editor
Sportisse, B., editor (1998). Air Pollution Modeling and Simulation APMS98, INRIA/ENPC
APMS Conference. INRIA

Sportisse, B., editor (2001a). Air Pollution Modeling and Simulation APMS01. Springer Verlag
18 Chapitre 2 - Bibliography of my research works

Le Dimet, F. and Sportisse, B., editors (2002). Data Assimilation for Geophysical flows, INRIA-
CEA-EDF School on Applied Non-linear Problems. INRIA

As Author
Sportisse, B. (2007b). Pollutions atmosphériques. Enjeux, processus, modélisation. Springer.
350 pages

2.9 Miscellaneous
Dumas, L., Coron, F., Dupuis, D., Pallegoix, J., and Sportisse, B. (1994). Applications of DSMC
methods to aerodynamic problems in an engineering context. In Méthodes de Monte-Carlo et
applications aux gaz raréfiés. Ecole CEA-INRIA-EDF, 1994

Coron, F., Pallegoix, J., and Sportisse, B. (1994). DSMC computations of complex test cases and
industrial configuration. In 19th International Symposium on Rarefied Gas Dynamics, Oxford

Sportisse, B. (2003a). Ozone des villes, ozone des champs: qualité de l’air et pollution pho-
tochimique. PCM Le Pont

Sportisse, B. (2007a). Management de la recherche. Enjeux et perspectives, chapter Partenariat


recherche publique/entreprise : l’exemple du CEREA, Laboratoire Commun ENPC/EDF R&D.
Chapitre 6. De Boeck
Chapter 3

Research works

The following presentation is mainly based on the review article Sportisse [2007c]. My works
are referenced through underlined references. The general framework (atmospheric chemistry
and physics for air pollution) is also detailed in a book written on the basis of a 10 years course
at Ecole Nationale des Ponts et Chaussées (Sportisse [2007b]).
Most of the works have been done with PhD students or post-doctoral fellows I have su-
pervised. For PhD works, I can cite: Rafik Djouad for model reduction and aqueous-phase
chemistry, Jaouad Boutahar for model reduction and 3D modeling, Edouard Debry for aerosol
modeling, Denis Quélo for inverse modeling and dispersion of radionuclides (among many other
applied studies), Yelva Roustan (I have initiated his PhD work with Luc Musson-Genon and
the supervision has then be made jointly with Marc Bocquet) for mercury and heavy metals,
Monika Krysta (I have initiated her PhD work and the supervision has then be made jointly with
Marc Bocquet) for inverse modeling of radionuclides, and Marilyne Tombette (I have initiated
her PhD work and the supervision has then be made jointly with Patrick Chazette) for aerosol
modeling and radiative effects. For postdoctoral works, I can cite the works of Karine Sartelet
(for air quality modeling and modal aerosol modeling) and Kathleen Fahey (for 3D modeling
of aerosols). The main part of the works devoted to radionuclides have been made with my
colleague from the French Institute for Nuclear Safety and Radiological Protection, Olivier Is-
nard. Even if we have not published papers together, many works have been initiated through
discussions with Luc Musson-Genon, especially when applied to topics relevant for Electricité de
France. Many works have also been initiated and discussed with my colleague Stephanie Lacour
for applications relevant for the French Ministry of Transport.
To conclude, many works have been performed with Vivien Mallet, first with his PhD work
devoted to uncertainties and ensemble forecast, and then for many topics related to Chemistry
Transport Models within the Polyphemus project.

3.1 Parameterizations
Subgrid parameterization is a key component of geophysical models. Many unresolved processes
have indeed characteristic scales much finer than the resolved numerical scale (typically the size
of a cell in a Finite Differences framework).
hΨi will stand for the average value of a physical field Ψ (the average has to be precised
0
rigorously). In practice, this corresponds to the numerical value. The local fluctuation Ψ is
0
defined through Ψ = Ψ − hΨi.
For
D non-linear
E processes, given for instance by quadratic terms such as Ψ 1 Ψ2 , a correlation
0 0
term Ψ1 Ψ2 has then to be parameterized as a function of the resolved values hΨ1 i and hΨ2 i.

19
20 Chapitre 3 - Research works

The issue is similar to the one arising for turbulence modeling (closure scheme).
The applications for air pollution modeling are related to the vertical turbulent fluxes (pa-
rameterization of Kz ), the microphysical mass transfer, the coupling between chemistry and
turbulence and Plume-In-Grid modeling.
My works have focused on mass transfer (scavenging by rain drops and cloud drops). Any-
how, I have made the choice to briefly review the whole topics.

3.1.1 Turbulent diffusion


A first example is given by the vertical transport of pollutants. In Eq. (1.1), K is the eddy
coefficient. The molecular diffusion can be neglected far from a thin laminar layer just above
the ground. After having averaged the dispersion equation for the “physical” fields, the
D advec-
E
0 0
tion term generates a mean advection term div (hV i hci) and a correlation term div V c .
Under some assumptions, the K theory (Lamb [1973]) gives the eddy vertical flux (given by the
correlation between the vertical velocity w and the concentration c) as:
D E  
0 0 ∂ hci
wc = − hρi Kz (3.1)
∂z hρi
D 0 0E
hci
I recommand this formulation rather than w c = −K∇ hci, in order to get zero if hρi is
constant. With this formulation, the properties of mass consistency have only to be checked for
the averaged advection term (see Section 3.2.3 and Sportisse et al. [2007]).
The eddy coefficient has then to be parameterized as a function of the meteorological fields
(see for instance Louis [1979] or Troen and Mahrt [1986]). This parameterization is a key
component of the CTM. Notice that the horizontal diffusion is usually neglected due to the
lack of a rigorous framework and due to the numerical diffusion induced by advection schemes
(supposed to be larger than the physical diffusion, Mallet et al. [2007b]).

3.1.2 Scavenging
The mass transfer between gas phase, aqueous phase and condensed matter (cloud droplets, rain
drops and aerosols) requires appropriate parameterizations. The typical microphysical scales are
indeed the following ones: an aerosol diameter range from a few nanometers to a few micrometers,
a cloud droplet diameter from a few tens of micrometers to one hundred micrometers, etc. Such
scales cannot be solved in a numerical model.

Loss terms for radionuclides

The atmospheric residence time of radionuclides is stronly related to wet scavenging and dry
deposition. The processes to be taken into account may vary for the radionuclides (iodine and
cesium). For instance, cesium is mainly bound to aerosols, which emphazises the role of wet
scavenging for its atmospheric residence. A review of the available parameterizations is presented
in Sportisse [2007d]. A key point is the representation of the rain distribution. The scavenging
coefficients are non-linear functions of the aerosol size and of the rain characteristics (see Figure
3.1 for the collision efficiency, a key component for scavenging).
A size resolved aerosol model is usually not available for radionuclides. It is therefore
mandatory to deal with parameterized models (see Sportisse [2007d]). A reasonable model
is Λ = A × p0 B (p0 is the rain intensity in mm/h, A ' 10−5 typically, B ' 0.8).
Parameterizations 21

0.5
Dr=0.1 mm
Dr=0.5 mm
0.0 Dr=1 mm
Dr=2 mm

-0.5

-1.0
log10 E(Dr,dp)

-1.5

-2.0

-2.5

-3.0

-3.5

-4.0
-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
log10 dp (in m)

Figure 3.1: Scavenging of an aerosol. Collision efficiency as a function of the aerosol size for
different representative raindrop diameters (Dr ). dp is the wet aerosol diameter.
22 Chapitre 3 - Research works

Gas-phase diffusion Mass transfer

S(IV)
S(VI)

Aqueous-phase diffusion Aqueous chemistry

Figure 3.2: Mass transfer in a cloud and aqueous-phase production of sulfate.

Mass transfer for cloud drops (in-cloud scavenging)


For gas-phase species, wet scavenging by cloud drops may be related to mass transfer (dissolu-
tion). Inside the cloud droplets, aqueous-phase chemistry may play a key role (for instance for
sulfate production, Figure 3.2).
The underlying models are diphasic models based on reaction-diffusion equations. Let c g be
the gas-phase concentration and ca be the aqueous-phase concentration.
Let Ωa be the cloud droplet (supposed to be a sphere of surface Γ and of radius a) and Ω g
be the interstitial gas. The microphysical model is then:

 ∂cg + Dg ∆cg = χg (cg ) sur Ωg

∂t (3.2)
 ∂ca
 + Da ∆ca = χa (ca ) sur Ωa
∂t
The mass transfer flux at the surface Γ is then given by:
1 ca
Dg ∇cg · ~n = −Da ∇ca · ~n = αv̄( − cg ) (3.3)
4 HRT
χg (respectively χa ) is the source term for gas-phase chemistry (respectively for aqueous-phase
chemistry). Dg (respectively Da ) is the molecular diffusion for the gaseous phase (respectively for
the aqueous phase). ~n is the unitary vector outward oriented from Γ; α is the accommodation
coefficient, v̄ is the thermal velocity for gaseous species, H is the Henry coefficient, R is the
perfect gas constant and T the temperature.
It is not possible to solve these equations coupled to a 3D CTM. A key property for the
parameterization is the wide range of the timescales: the molecular diffusion is much faster than
the chemical reactions at theses scales (this is the opposite of the situation in Section 3.6). In
Parameterizations 23

practice, one can define a lumped model by replacingR the microphysical concentrations c a and
cg by averaged concentrations, for instance c̄a = Ω1a Ωa ca (ω, t)dω. The time evolution of the
averaged concentrations is then given by a system of Ordinary Differential Equations (ODE):

 dc̄g = χg (c̄g ) − L I(c̄a , c̄g )

dt (3.4)
 dc̄a
 = χa (c̄a ) + I(c̄a , c̄g )
dt
L is the liquid water content (ratio of the volume Ωa to the volume Ωg ). The mass transfer
c̄a
flux I is parameterized as I(c̄a , c̄g ) = kmt (c̄g − ), with kmt a mass transfer coefficient. A
HRT
classical formula is:  2 −1
a 4a
kmt = + (3.5)
3Dg 3v̄α
Such derivations are rigorously presented in Sportisse and Djouad [2003] with the use of
asymptotic expansions with respect to “small” parameters (defined as Peclet-Dahmköhler num-
bers). The mathematical framework can be found in Fife [1979] for reaction-diffusion systems
and Vasileva et al. [1994] for boundary layer methods.
It is for instance proven that the double film parameterization has to be preferred under
some conditions:  2 −1
a a2 4a
kmt = + + (3.6)
3Dg 3HRT Da 3vα
There are still many other challenging issues: how to extend this approach to polydisperse
models ? How to take into account surface reactions (not described in Eq. 3.3) ?

Below-cloud scavenging
A second example is given by below-cloud scavenging. In this case, the scavenging is related to
falling raindrops (of falling velocity U and of radius a). With notations similar to below, the
microphysical model is:
 
 ∂cg L ca 

 =− cg −
∂t τtransfer HRT
∂ca ∂ca 1  ca  (3.7)

 + U = c −
 g
∂t ∂z τtransfer HRT

τtransfer is a timescale associated to mass transfer (not detailed here), given as the inverse value
of a mass transfer coefficient.
This model can be replaced by the following parameterization after having applied an asymp-
totic expansion (the “small” parameter is in this case related either to a large solubility of to a
small liquid water content):
dcg
= − Λ g cg (3.8)
dt
with Λg (p0 , a, HRT, τtransf er ) a scavenging coefficient. One refers to Sportisse and Du Bois [2002]
for more details.

3.1.3 Segregation effects


Another example is provided by the coupling between chemistry and
D turbulence.E After having
0
averaged the dispersion equation, the chemical term appears as χi (hci + c ) . One usually
24 Chapitre 3 - Research works

0.1

-0

-0.1

-0.2

-0.3
Is

-0.4

-0.5

-0.6

-0.7
0 10 20 30 40 50
Time (UTC)

Figure 3.3: Evolution of surface segregation intensity (15-16 September 1989), from observations
in the Netherlands (from Vila-Guerau de Arellano et al. [1993]).

assumes that this term is equal to χi (hci), which may be a crude approximation for the fastest
reactions. The underlying assumption is indeed that the averaging is quicker than any chemical
reactions. We refer for instance to the review Vila-Guerau de Arellano et al. [2004].
A typical example is given by the so-called segregation effect between an updraft of emitted
NO and a downdraft of O3 (entrainment) in the atmospheric boundary layer. The chemical
reaction of O3 with NO (titration) is:

O3 + NO → NO2 + O2

The associated reaction rate is ω = k cO3 cNODwith k the


E kinetic rate. After averaging, the true
0 0
averaged reaction rate is hωi = k(hcO3 cNO i + cO3 cNO ), where we suppose that the kinetic rate
has a constant value (it actually depends on temperature). This can also be written as:

hωi = k hcO3 cNO i (1 + Is ) (3.9)


D 0 0
E
cO cNO
3
with Is = the segregation intensity. Mathematically speaking, Is ≥ −1; physically
hcO3 cNO i
speaking, Is ≤ 0, because O3 is associated to downdrafts while NO is associated to updrafts.
Measured values of Is can be found in Figure 3.3.
In the “laminar” assumption, Is is taken equal to 0. A parameterization is then required
for Is (Molemaker and Vila-Guerau de Arellano [1998]; Vinuesa and Vila-Guerau de Arellano
[2003, 2005]). Once more, the key parameters are the so-called Dahmköhler numbers defined as
the ratio of a transport timescale to a chemical timescale. We also refer to Sportisse [2001b] for
a simple estimation of the segregation errors in a 1D case.
To date, there are no available parameterizations of Is to be used in 3D models with com-
prehensive chemical mechanisms.
Numerical analysis for reactive dispersion 25

3.1.4 Plume In Grid models


The last example is provided by the so-called Plume-In-Grid models. Emissions of point sources
in a large grid cell induce a numerical artefact by artificially creating diffusion. One way to
reduce this diffusion is to parameterize short-term dispersion of freshly emitted pollutants as
an alternative to the use of refined meshes (see Section 3.2). The parameterizations are usually
based on Gaussian models (stationary solutions of the dispersion equation). We refer for instance
to Karamchandani et al. [2000]. An issue is the rigorous coupling between long-range Eulerian
models and short-range Gaussian models for reactive flows.

3.1.5 Future works


These examples illustrate the need for rigorous approaches in order to derive the appropriate
parameterizations. An alternative approach is based on numerical sampling. In Larson et al.
[2005], a numerical sampling technique (a Latin Hypercube method) is for instance used in order
to define the parameterization. A Probability Density Function is then computed on the basis
of repeated calls to a microphysical function.
This method is an illustration of “micro/macro” methods that are used in other fields.
The principle is to couple a deterministic “macro” model to a stochastic “micro” model. The
application of such methods to air quality modeling should be investigated.

3.2 Numerical analysis for reactive dispersion


We now assume that the model is fixed. The numerical simulation may provide a large num-
ber of challenging numerical points and we pay attention to four numerical issues: operator
splitting (the processes are usually not solved simultaneously), time integration of stiff ODEs
(resulting from chemistry or aerosol dynamics), advection (with the specific problem of mass
consistency) and grid resolution. The most challenging issues are related to aerosol dynamics
and are described in the next section.
A detailed review can be found in Verwer et al. [1998] for some of these topics (excluding
aerosol dynamics). We can also refer to the book Hundsdorfer and Verwer [2003].
The numerical issues are investigated in Mallet et al. [2007b] for some real cases (air quality
modeling over Europe for summer 2001 and simulation of the Chernobyl accident).
These topics are also presented in the course textbook Sportisse and Mallet [2005].

3.2.1 Splitting methods


Chemistry-Transport models describe many processes related to dispersion, microphysics and
chemistry. Although all these processes are coupled, operator splitting methods (for instance
Marchuk [1986]; Yanenko [1971]) are widely used for solving the dispersion equation. There are
at least three motivations for this numerical strategy.
First, models can be used as black boxes and are easy to insert in an existing comprehensive
three-dimensional model. Second, appropriate specific numerical methods can be used for each
process. Third, coupling dispersion and chemistry may result in a huge computational burden
because implicit methods are used for chemistry (due to the dispersion of timescales, see below).
In fully coupled methods, the size of the state vector is then typically equal to the product of
the number of species by the number of grid cells (for Finite Differences solvers). The resulting
Jacobian matrices are then difficult to invert.
26 Chapitre 3 - Research works

For the sake of clarity we consider two processes that do not depend on time. The evolution
equation is therefore:
dc
= A1 (c) + A2 (c) , c(0) = c0 (3.10)
dt
A1 (for instance diffusion) and A2 (for instance chemistry) are the source terms associated to
the processes.
Let ∆t be the splitting timestep (for a CTM, it typically ranges from 100 to 1000 seconds
at continental scales). We formally write LA1 +A2 (c0 , t) the solution of (3.10) at time t.
The simplest method is the first-order splitting: one first solve the process A 1 over [0, ∆t]
and then (after having modified the initial condition) the process A2 over [0, ∆t]. The solution
is then LA2 (LA1 (c0 , ∆t), ∆t). Notice that we can also reverse the sequence between A1 and A2 .
The so-called Strang splitting (Strang [1968]) eliminates the lack of symmetry of the previous
method. There are three steps in the sequence: A1 is integrated over [0, ∆t/2], then A2 is inte-
grated over [0, ∆t] and finally A1 over [0, ∆t/2]. The solution is therefore LA1 (LA2 (LA1 (c0 , ∆t/2), ∆t), ∆t/2).
Operator splitting methods induce a loss of accuracy. This is usually estimated with the
computation of the local splitting error in the linear case for Ai (c) = Ai c (Ai is a matrix).
In this case, LA (c0 , ∆t) = eA∆t c0 . The splitting error is related to the commutation default:
eA1 +A2 = eA2 eA1 if and only if A1 A2 = A2 A1 .
Asymptotic analysis with respect to ∆t → 0 prove that the Strang splitting method is of
second-order (and that the first-order method is well named!).

Some results are available for the reaction/diffusion/advection PDE with the extension to
the nonlinear case by using the Lie formalism (Lanser and Verwer [1998] or Sportisse [1999] for
direct calculations). The key results are the following ones (boundary conditions are not taken
into account):

• advection and chemistry commute if chemistry does not depend on space and if the wind
field is divergence free;

• advection and diffusion commute if the wind field and the diffusion matrix do not depend
on space;

• diffusion and chemical kinetics commute if chemical kinetics is linear and does not depend
on space.

It is clear that the assumption concerning the independence with respect to time and space is
not rigorously met for the eddy diffusivity matrix and the temperature field. This is however a
rather good local assumption. On the contrary, the assumption of linearity for chemistry is clearly
wrong (chemical reactions are supposed to describe... reactions between reactants!). However,
these results give an indication for the main source of errors related to operator splitting methods
in air pollution models: chemistry and vertical diffusion have to be handled with a special care
(Graf and Moussiopoulos [1991]).
A difficult point (not addressed here) is also related to the way the boundary conditions have
to be taken into account (see for instance Sommeijer et al. [1981] or Sportisse et al. [2000] for
an illustration).

For atmospheric chemistry, the previous results (especially the order analysis) cannot be
applied as such. A major point is indeed that one process (for instance chemistry) contain fast
dynamics, that quickly reach an equilibrium (see below for more details and the relation to
stiffness). In practice, this means that there exist processes whose timescales τ are much lower
Numerical analysis for reactive dispersion 27

than the splitting timestep ∆t. The asymptotic analysis ∆t → 0 is therefore no more justified.
We refer for instance to Sportisse [2000]; Verwer and Sportisse [1998]; Sportisse et al. [2000] for
a deeper understanding.
One key result is the impact of the splitting sequence: it is advocated to end the sequence
with the stiff processes (chemistry for instance in the case of gas-phase chemistry, mass transfer
for aerosol dynamics) in order to ensure that the system is still along the equilibrium manifold
at each step (see Figure 3.4). A second result is related to the so-called order reduction of the
numerical algorithms. Similar results have been obtained in other fields, for instance in combus-
tion (Yang and Pope [1998]). See also the discussion for oceanography in Higdon and Bennett
[1996] or for climate models in Williamson [2001].

Methods that avoid splitting are also used. Among them, one can cite implicit-explicit
methods that are based on a partitioning between the stiff processes (with an implicit treatment)
and the nonstiff processes (with an explicit treatment). There are many variations of this method
(for instance Verwer et al. [1996a]; Sun [1996]; Wolke and Knoth [1998]; Nguyen and Dabdub
[2003]).
For instance, let us suppose that A1 and A2 are the nonstiff process and the stiff process,
respectively. A possible implementation (usually referred as Source Splitting) is to first solve
LA1 (c0 , ∆t) − c0
A1 and then to solve A2 after having added a correction term Ã1 = without
∆t
modifying the initial condition. The result is then: LA2 +Ã1 (c0 , ∆t). The key advantage is that
the integration is always performed in the vicinity of the equilibrium manifold of the stiff process
(Sportisse [2000]).
Another approaches are the so-called Approximate Matrix Factorizations (also referred as
“internal splitting”, see Berkvens et al. [2002]). The idea is to perform splitting at the algebraic
level. Suppose that an implicit algorithm (let say Backward Euler) is used for solving equation
(3.10) in the linear case.
We have therefore (I − (A1 + A2 )∆t) cn+1 = cn . The method is based on approximations
such as (I − (A1 + A2 )∆t) ' (I − A2 ∆t)(I − A1 ∆t) + O(∆t2 ).

3.2.2 Time integration of chemical kinetics


In this subsection, ∆t is the subcycling timestep.
The CTM are usually based on the so-called Method of Lines: the space discretization is
first performed (especially with Finite Differences or Finite Volumes methods) and then the time
integration of the resulting ODE. Let us consider the case of reactive dispersion (the processes
are advection, diffusion and chemistry). The time integration of the advection step is not an
issue (or specified by the advection scheme); the time integration of (linear) diffusion is usually
performed with an implicit solver; the most challenging point is then the time integration of
chemical kinetics.
Let the system of ODE describing chemical kinetics be:
dci
= χi (c, t) , 1 ≤ i ≤ n (3.11)
dt
The first difficulty is the large dimension (n  1): the dimension of a comprehensive gas-phase
chemical mechanism is typically of magnitude 100 (for instance, RACM Stockwell et al. [1997]).
The second difficulty is related to the wide range of timescales, which induces the so-called
stiffness (Byrne and Hindmarsh [1987], see Figure 3.5). Benchmarks of explicit and implicit
methods have been performed in many works. One refers for instance to Sandu et al. [1997b,
1996, 1997a]; Verwer et al. [1996b].
28 Chapitre 3 - Research works

y
y
1 2 y = h(x) 1-2

M OL NT S

x x

y
y 2 2
1
1

(T − χ) (χ − T )
x
x

y
y 3 3
2
2

1
1

(χ − T − χ) (T − χ − T )

x
x

Figure 3.4: Impact of the splitting sequence. M OL stands for the Method of Lines, N T S for
No Time Splitting. T stands for transport and χ for chemical kinetics. x and y are the solw
and fast components of concentrations, respectively. The equilibrium manifold defined by fast
processes is defined by the constraint y = h(x).
Numerical analysis for reactive dispersion 29

10

−2

−4

−6

−8

−10

−12
0 10 20 30 40 50 60

Figure 3.5: Typical distribution of the eigenvalues for the Jacobian matrix associated to the
chemical mechanism RACM (logarithmic scale).

Among many candidates, the second-order Rosenbrock solver (Verwer et al. [1999] for gas-
phase chemistry, Djouad et al. [2002] for aqueous-phase chemistry) is a promising algorithm as it
realizes a good trade-off between a reasonable accuracy and a rather low computational burden.
The Rosenbrock method is defined in the following way from iteration n to iteration n + 1:

k1 + k 2
 cn+1 = cn +

2 (3.12)
 (I − γ∆t J)k1 = ∆tχ(cn )
 n
(I − γ∆t J)k2 = ∆tχ(c + ∆t k1 ) − 2γ∆t J k1

J = ∂χ∂c is the Jacobian matrix for chemical kinetics. γ = 1 + 1/ 2 ensures the L-stability
of the method (see Verwer et al. [1998] and also Mallet et al. [2007b] for the non-autonomous
case). Preprocessing the computation of J with appropriate techniques that take into account
the sparse nature (a chemical species is actually involved in a small set of reactions) is a powerful
technique (Sandu et al. [2003]; Daescu et al. [2003] and Mallet and Sportisse [2004]).

Another candidates are provided by explicit and quasi-explicit methods. The so-called QSSA
(Quasi Steady State Assumption) methods (Young and Boris [1977] for the original paper and
then for instance Jay et al. [1997]; Verwer and Van Loon [1994]; Mott et al. [2000]) are still
widely used even if they do not perform well in fair benchmarks (with error/CPU diagrams).
Such algorithms are strongly related to the so-called production-loss form of (3.11):

χi (c) = Pi (c) − Li (c)ci (3.13)

with Pi and Li two nonnegative terms. This form is easy to get by partitioning the chemical
reactions into two sets: reactions in which the species i is a reactant (loss term), reactions
in which it is a product (production term). The fast species and the slow species are defined
by timescales Li  1 and Li  1, respectively. There are many variations on the basis of
this partitioning of species. I refer for instance to my work Djouad and Sportisse [2002a] for a
30 Chapitre 3 - Research works

precise definition of the partitioning. In the simplest form of the QSSA method, the fast species
are computed from iteration n to iteration n + 1 with the equilibrium cn+1
i =L Pi
i
while the slow
Pi n Pi
species are computed with ci = Li + (ci − Li ) exp(−Li ∆t).
Beyond QSSA methods, quasi-explicit methods are also very popular. The Two-Step method
(Verwer [1994]; Verwer and Simpson [1995]) is an interesting algorithm that can be seen as an
extension of the EBI method (Euler Backward Iterative, Preussner and Brand [1981]; Hertel
et al. [1993]). Using the Backward Euler method leads to the following nonlinear algebraic
system:
cn+1 − cni
i
= χi (cn+1 ) (3.14)
∆t
This can be solved by an iterative algorithm, for instance the Newton algorithm, which re-
quires to compute the Jacobian matrix J. Another approach is to approximate χ i (cn+1 ) '
Pi (cn ) − Li (cn )cn+1
i . The algebraic system can then be solved with:

cni + Pi (cn )∆t


cn+1
i = (3.15)
1 + Li (cn )∆t

Variations of this idea lead to the EBI (first-order) and Two-Step (second-order) algorithms.
Notice that P (.) and L(.) can be evaluated with the components of cn+1
i that have been already
computed (which is the starting point of many tailored techniques).

Let us mention for the sake of completeness the combined explicit/implicit methods (for
instance Gong and Cho [1993]), that are based on the partitioning of chemical species in a set
of slow species (integrated with an explicit solver) and a set of fast species (integrated with an im-
plicit solver). We will go back to this approach in Section 3.6 (see also Djouad and Sportisse [2003]).

To conclude, a specific issue for the time integration of chemical kinetics is related to the
positivity of concentrations. A usual tricky approach is the clipping (set the negative numerical
concentrations to 0 or to a small predefined parameter ). We refer to Sandu [2001b] for a rigorous
method. The impact for real cases is also evaluated in Mallet et al. [2007b].

3.2.3 Advection and mass consistency


Solving advection is the other challenging step for the CTMs. The results strongly depend on
the accuracy of the wind velocity (computed in meso-scale meteorological models). There are
however some remaining issues for CTMs.
The first point is that the advection has to take care of sharp gradients. This is is a key
concern for the dispersion of accidental releases, for instance of radionuclides. There is indeed no
background field and the point accidental emission implies strong gradients. Advection solvers
with flux limiting are therefore often advocated (see a general presentation in Verwer et al.
[1998]).
Another specific point is related to the so-called mass consistency problem (Sportisse et al. [2007])
related to the off-line coupling between atmospheric dispersion models and meso-scale numerical
models that provide meteorological fields such as wind velocity.
∂ρ
The continuity equation for density ρ can be written under the form + div(V ρ) = 0 with
∂t
V the wind field. The advection of a passive tracer of concentration c (defined for instance as
∂c
the ratio of mass per air volume) is given by + div(V c) = 0. The evolution of the mass
∂t
∂q
mixing ratio q = c/ρ (mass of species per air mass) is then given by + V · ∇q = 0.
∂t
Numerical analysis for reactive dispersion 31

The conservation of uniform mixing ratio in the case of passive tracers is then easy to
prove: if q(0) = 1, one gets for any time t, q(t) = 1. This property is usually referred as mass
consistency. The key point is the consistency between the wind field used for ρ and the wind field
used for c. The discrepancy may result from different numerical algorithms and discretization
strategies between the computation of the wind velocity (inside the meteorological model) and
the advection step (inside the CTM).
We briefly present three possible methods in order to circumvent these difficulties.
A first approach is the renormalization technique. Let cn+1 = S(cn , V n ) be the advection
scheme used by the CTM (supposed to be an explicit solver, which is usually the case in practice).
Let us define ρ̃n+1 = S(ρn , V n ) the density that is computed with this scheme (this is not the
“true” density ρn+1 if the meso-scale model does not use the same scheme). The renormalization
approach is then based on the correction:

ρn+1
cn+1 = S(cn , V n ) (3.16)
ρ̃n+1

It is of course easy to check that mass consistency is met. One refers for instance to Byun
and Lee [2002]; Lee et al. [2004]. A major drawback of this approach is that it may artificially
introduce mass.
A second approach is based on the use of mass mixing ratio in the computation of fluxes after
∂c c
having noticed that + div(ρ V ) = 0. Convergence results for this algorithm are however not
∂t ρ
proved. One refer to Petersen et al. [1998] for the special case of the upwind scheme (without
flux limiting).
A third approach is to modify the vertical component of wind velocity, w, in order to satisfy
mass consistency. Suppose that horizontal advection is defined by cH = SH (cn , un , v n ) (even-
tually with directional splitting). One defines ρH = SH (ρn , nn , v n ). Let the vertical advection
scheme be given by cn+1 = SV (cH , wn ). We therefore try to compute w n such that the con-
straint ρn+1 = SV (ρH , wn ) is met. Of course, this implies that this constraint may be solved in
wn , which implies that a nonlinear scheme (for instance with flux limiting) cannot be used for
vertical advection. One example of a scheme easy to “invert” is the upwind scheme (Hu et al.
[2006]; Odman and Russel [2000]).
These approaches are clearly defined and benchmarked in Sportisse et al. [2007] for two case
studies (air quality modeling over Europe and Chernobyl accident). The third approach is an
outlyer, especially due to the strong modification of the vertical wind profile (Figure 3.6).

3.2.4 Future works


There are still many open points in this field, especially for stiff problems and advection.
In the case of point sources, having a fine resolution may be a key issue. We have already
briefly mentioned the Plume-In-Grid approach. An alternative to this approach is of course
to use fine meshes. A first method is based on adaptive gridding (Srivastava et al. [2001]).
Unstructured meshes are more appropriate for this aim (see for instance Tomlin et al. [1997]).
In Van Loon [1996], local refinement is applied with a uniform-grid based cascade technique that
produces nested grids within one timestep. A second method is based on nesting techniques.
The idea is to compute solutions on different fixed meshes of variable resolutions with one-way
or two-way fluxes of informations (for instance Frohn et al. [2002]).
The grid resolution may have a strong impact on the results. The future generation of CTM
will have a fine resolution (let say 0.1◦ for the horizontal dimension, less than 5 kilometers) in
coherence with the high-resolution meteorological models, currently under development.
32 Chapitre 3 - Research works

3000
dst3
upwind

2500

2000
height (m)

1500

1000

500

0
-2.0 -1.5 -1.0 -0.5 0.0 0.5
x1e-3
vertical wind (m/s)

Figure 3.6: Profile of the vertical wind velocity: reference case (DST3) and vertical adjustment
method (from Sportisse et al. [2007]).

One cas also wonder if the next generation of CTM will not be based on unstructured
methods, similar to those used in CFD (Computational Fluid Dynamics).

3.3 Aerosol modeling and simulation


The 3D applications (including model-to-data comparisons) can be found in Section 3.7. This
Section details the parameterizations and the numerical algorithms.

3.3.1 The General Dynamics Equation for aerosols


The numerical simulation of gaseous species is the first component of a CTM. For many ap-
plications, the modeling of condensed matter is also required. This includes the description of
chemical composition and of size distribution of aerosols (particulate matter). Both have an
influence for the radiative behavior of aerosols, for microphysical processes (through the activa-
tion of the finest aerosols to generate cloud droplets) and for the assessment of health impact
(the most adverse particles are the finest ones).
Simulating a polydisperse distribution of aerosols is a challenging task for many reasons (at
the numerical level):

• there is a wide range of scales (from a few nanometers to at least 10 µm),

• one has to deal with high-dimensional models (the number of chemical species inside the
aerosol mixture has to be multiplied by the number of degrees of freedom chosen to describe
the size),
Aerosol modeling and simulation 33

Hysteresis for ammonium nitrate

500

deliquescence branch
metastable branch
400

liquid water content ( microgram )


300

200

100

0
0 20 40 60 80 100
relative humidity

Figure 3.7: Hysteresis for a salt of ammonium nitrate (computed with Isorropia: T =300 K,
[HNO3 ]=[NH3 ]=2 µg).

• the processes that affect the size and chemical composition of aerosols are nonlinear and
discontinuous (thresholds for phase transition: see for instance Figure 3.7 for the hysteresis
of deliquescence/cristallisation).

The evolution of size distribution and of chemical composition for aerosols is governed by the
so-called General Dynamics Equation for aerosols (GDE in the sequel). It describes the impact
of some processes (Figure 3.8) such as nucleation, coagulation, and mass transfer (condensation
and evaporation).

Nucleation The formation of the smallest aerosols is induced by the agregation of gaseous
molecules through thermodynamically stable clusters. One usually takes into account binary
nucleation (H2 O-H2 SO4 ) or ternary nucleation (H2 O-H2 SO4 -NH3 ).

Coagulation Coagulation is driven for instance by the agitation of the fluid. In practice,
one only describes Brownian coagulation due to thermal agitation. Other effects (the so-called
phoretic effects due to gradients in fields (such as temperature, electric fields, etc) are not taken
into account. Coagulation may be neglected for the evolution of particles above a few µm.

Condensation/evaporation Many gaseous species with a low saturation vapour pressure may
condense onto or evaporate from existing particles. This mass transfer is governed by the
gradient between the concentration at the aerosol surface (supposed to be at local equilibrium
with the internal mixture) and the concentration far from the aerosol.
Condensation is a fast process for submicronic aerosols.

The composition of the continental inorganic aerosol is mainly governed by the interaction
between sulfate, nitrate and ammonium. The corresponding gaseous species are:
34 Chapitre 3 - Research works

primary particle emission


primary gas emission

heterogeneous condensation
reactions coalescence

evaporation

homogeneous nucleation coagulation coagulation hygroscopic activation


reactions

activation and scavenging

wet
scavenging
dry deposition

D < 0.01 0.01 < D < 0.1 0.1 < D < 1 1 < D < 10
gas molecules aerosols aerosols aerosols aerosols cloud droplets

Figure 3.8: Aerosol dynamics (the aerosol diameter, D, is in µm).

• sulfur dioxide (SO2 ) to be oxided to sulfuric acid (H2 SO4 ) whose saturation vapour pressure
can be neglected (H2 SO4 is only under particulate form);

• nitric acid (HNO3 ) that results from the oxidation of NOx ;

• ammoniac (NH3 ), mainly related to emissions from agriculture.

Nitrate and sulfate are in competition inside aerosols for the formation of ammonium nitrate
(NH4 NO3 =NH+ – + 2–
4 +NO3 ) and ammonium sulfate ((NH4 )2 SO4 =2NH4 +SO4 ). A key point is
that the second salt is preferred.

Mass transfer to cloud droplets Part of the aerosol distribution acts as Cloud Condensation
Nuclei. The chemical composition of the activated aerosols is then governed by aqueous-phase
chemistry.

Heterogeneous reactions Surface reactions at the surface of condensed matter (aerosols and
cloud droplets) may have a significant impact for some species like ozone. These processes are
usually described by first-order kinetics.

Loss processes Scavenging and dry deposition have also to be taken into accound. For the
coarsest aersols, gravitational settling (sedimentation) has also to be described.

The aerosols are complicated mixtures of organic species, inorganic species, soot (black car-
bon) and liquid water. A typical composition of the continental aerosol over Europe can be
found in Figure 3.9. Notice that part of the composition is not known and that there are many
Aerosol modeling and simulation 35

Nitrate
Organics

15.0%
20.0%
Ammonium

7.0%
Black Carbon
5.0%

13.0% 9.0%
Sulfate
Mineral Dust
4.0%

27.0%
Sea salt

Other (27 %)

Figure 3.9: Composition of the “European” aerosol (PM10 ) in urban conditions, from Putaud
et al. [2003]). ’Other’ refers to the unknown part.

issues for having accurate measurements.

Let n(m, t) be the number distribution (for instance with respect to aerosol mass) and q i (m, t)
the mass distribution for species i (we use the same notation as for the mass mixing ratio in
the previous section). m0 is the smallest aerosol dry mass (in practice given by the nucleation
threshold). The GDE is composed of the following equations:
• For the number distribution:
Z
∂n 1 m−m0
(m, t) = K(u, m − u)n(u, t)n(m − u, t) du
∂t 2 m0
| {z }
coagulation gain
Z ∞ (3.17)
∂(I0 n)
− n(m, t) K(m, u)n(u, t) du − + δ(m0 , m)J0 (t)
| m0
{z } | ∂m{z } | {z
nucleation
}
mass transfer
coagulation loss

• For the mass distributions:


Z m−m0
∂qi
= K(u, m − u)qi (u, t)n(m − u, t) du
∂t m0Z

−qi (m, t) K(m, u)n(u, t) du (3.18)
m0
∂(I0 qi )
− + Ii (m1 , . . . , mnc , t)n(m, t) + χaer (m1 , . . . , mnc )n(m, t)
∂m |i {z }
Internal chemistry

P c
An aerosol of mass m is composed by nc components of mass mi (m), such that i=ni=1 mi (m) =
m. Notice that qi (m) = mi (m)n(m). An aerosol with a given mass is assumed to have a unique
composition (internal mixing).
36 Chapitre 3 - Research works

K(., .) is the coagulation kernel (due to Brownian motion, not detailed here). I i is the
condensation/evaporation rate for volatile species i:

Ii = ai (m) (ci − ηi ceq


i (m1 (m), . . . , mnc (m))) (3.19)

where ai (m) stands for a parameterization of a mass transfer coefficient (a function of accommo-
dation, of molecular diffusion and of particle size), ci is the gas-phase concentration of species i
and ceq
i is the surface concentration of species i (supposed to be at equilibriumP with the aerosol
mixture). ηi is the Kelvin coefficient that describes the curvature effect. I0 = i Ii is the total
rate.
J0 is the nucleation rate, given by parameterizations as a function of thermodynamic con-
ditions (temperature and relative humidity) and of some gaseous concentrations (water, ammo-
nium, sulfate).
The internal chemistry is usually neglected.

This set of Partial Integro-Differential equations is strongly coupled, especially through the
condensation/evaporation term. One key component is the availability of a thermodynamic
model for computing ceq i (it is the most expensive part of the model, as far as CPU is concerned,
Zhang et al. [1998b]).
Notice that the GDE is usually not solved simultaneously with the other processes. We re-
fer to Müller [2001] for the investigation of splitting errors related to microphysical-multiphase
chemical systems.

One refers to Sportisse et al. [2006a,b] for a description of the parameterizations used in the
aerosol models that have been developed under my leadership: the Model Aerosol Model MAM
(Sartelet et al. [2005]) and the SIze REsolved Aerosol Model SIREAM (Debry et al. [2007]). I
will focus on numerical points thereafter.

3.3.2 Stochastic methods


Analytical solutions of the GDE only exist for some academic cases (Gelbard and Seinfeld [1979];
Ramabhadran et al. [1976]). This means that it is difficult to benchmark numerical algorithms
with respect to a reference (exact) solution. The lack of reference solution has already been
underlined in Tzivion et al. [1999]. Searching for a numerical reference solution is therefore
a key issue. Such a solution is not necessary computationally efficient but has to be highly
accurate.
The stochastic methods are good candidates for being such reference solutions. They have
a rigorous theoretical basis (through the stochastic formulation of the GDE) and are easy to
implement. In stochastic methods, one deals with numerical particles that represent, in a way
to specify (see below), the aerosol distribution.
We refer to Domilovskii et al. [1978] for the application to the coagulation equation (the so-
called Smoluchowski equation). A numerical particle is associated to a fixed number of physical
aerosols of a given size. The drawback is that the number of numerical particles decreases as
coagulation makes the number of aerosol decrease. This is a strong limitation as the method
converges in the limit of the number of particles going to infinity. Babovsky [1999]; Eibeck and
Wagner [2000] have associated to one numerical particle a fixed mass of aerosols. As coagulation
preserves mass, the number of numerical particles can remain constant. This kind of algorithms
is usually referred as mass flow algorithms.
We refer to Debry et al. [2003] for the extension to the GDE. The key idea is to use mass
flow algorithms with a varying mass due to condensation/evaporation.
Aerosol modeling and simulation 37

6000
Mode 1
Mode 2
Mode 3
Total
Distribution en nombre (n) 5000

4000

3000

2000

1000

0
2 1 0 1
10 10 10 10



DiameÁtre ( m) 

Figure 3.10: Modal aerosol distribution (measured at Leipzig, Germany, from Putaud et al.
[2003]).

3.3.3 Modal methods


For 3D modeling, a popular approach is based on the modal approach. The aerosol distribution
is represented by a given number of modes (typically log-normal functions of aerosol mass or
volume, Whitby and McMurry [1997], Figure 3.10).
For instance (the Modal Aerosol Model, MAM, Sartelet et al. [2005]), one usually defines
three modes in the fine range (nucleation i, Aitken j and accumulation k), and one mode
in the coarse range (c). Each of these modes correspond to a distinct range of particle dry
diameters D: D < 0.01µm for the nucleation mode, 0.01µm < D < 0.1µm for the Aitken mode,
0.1µm < D < 2.5µm for the accumulation mode and D > 2.5µm for the coarse mode. Each of
the four modes has its own chemical composition.
The aerosol distribution is then given by summing the mode contributions, for example for
the number distribution: n(D, t) = ni (D, t) + nj (D, t) + nk (D, t) + nc (D, t), where ni is the
log-normal distribution for mode i, etc. For l = i, j, k, c, one has:
Nl 1 ln2 (D/dl )
nl (D) = √ exp[− ( )] (3.20)
2π ln(σl ) dl 2 ln2 (σl )
where Nl is the total number of particles within mode l, dl is the ”dry” median diameter and
σl is the standard deviation. The moment of order h of the distribution is defined as
Z ∞
l h2
Mh = Dh nl (D) d(D) = Nl dhl exp( ln2 σl ). (3.21)
−∞ 2
The modal distribution is known once the three parameters N , dg and σg are. For each
mode, an evolution equation is derived from the GDE for three moments (M0 , M3 and M6 ),
38 Chapitre 3 - Research works

from which the three parameters may be computed:


s
M34 1 1 M0 M6
N = M 0 , dg = ( ) 6 , σg = exp( ln( )) (3.22)
M6 M03 9 M32

The third moments are computed for each chemical species in order to follow the chemical
composition in each mode.
The resulting system of (stiff) ODEs can be solved with appropriate time integration solvers
(see previous section).
A specific issue for the 3D implementation of modal methods is related to the management
of mode merging and splitting. Merging between two modes is required to force modes to be
of distinct size ranges throughout the simulations (Binkowski and Roselle [2003]; Whitby et al.
[2002]). When the combined effect of nucleation and condensation and the effect of coagulation
are of same order of magnitude but act in opposite direction, a mode may split into two modes.
This is particularly likely to happen for the nucleation mode i during high nucleation episodes.
A splitting scheme has then to be applied (Binkowski and Roselle [2003]).

3.3.4 Size-resolved methods


Modal methods are based on an a priori analytical form of the aerosol distribution. Another
approach is of course to discretize the GDE, viewed as a system of integro-differential equations.
We briefly present variational methods (based on a weak formulation of the GDE) and the more
classical sectional method (based on Finite Differences).

Variational methods
Let us write the GDE under the following form:

dn dqi
= f (n, x, t) , = fi (n, qi , x, t) (3.23)
dt dt
x = ln m is the logarithm of the dry mass (the variable that is classically used for numerical
algorithms). If V is a functional space (of functions of x) with an appropriate scalar product
(., .), one defines the weak formulation of the GDE (for instance for the number distribution) as:

dn
∀v ∈ V, ( , v) = (f (n, x, t), v) (3.24)
dt
A numerical method is then based on a projection onto a subspace Vk of V , of finite dimension
k:
j=k
X
Πn(x, t) = nj (t)Lj (x) (3.25)
j=1

where the functions Lj (x) are the basis functions of Vk .


In Pilinis [1990], the basis functions are the so-called “hat functions” and one gets a Finite
Element Method. Another family of functions with a bounded domain of definition is given by
cubic splines (Gelbard and Seinfeld [1977]; Sandu [2001a]).
When the basis functions do not have a bounded domain of definition, for instance in [x 0 , ∞[,
one usually refers the method as a spectral method. In Sandu and Borden [2003], the basis
functions are polynomial functions of Lagrange. In Debry [2004], they are polynomial functions
of Laguerre and of Legendre.
Aerosol modeling and simulation 39

Apart from the choice of the basis functions Lj , the test functions v have also to be chosen.
In the case of a collocation method, the test functions are Dirac functions at collocation points
(xi )i=1,..., k : vi (x) = δ(x − xi ). For a basis of Lagrange polynomials, one easily gets for the
number distribution (Sandu and Borden [2003] and Debry and Sportisse [2007a]):

d~n
= A(~n)~n − diag(B~n)~n − |{z}
C~n (3.26)
dt | {z }
coagulation c/e

where (~n)j = nj . The matrices A(n) (the line i of A(n) is under the form ~nT Ai~n with Ai a
matrix) and B are related to coagulation. The matrix C is related to condensation/evaporation.
The resulting system of ODEs has then to be integrated with appropriate solvers. We refer
to Debry and Sportisse [2007a] for the benchmark of two collocation methods.
The advantage of such methods is that one can solve the GDE without using splitting
methods. However, such methods may suffer of numerical instabilities and are computationally
demanding.

Sectional methods
Another approach is based on Finite Differences. The aerosol distribution is discretized in a
fixed number of bins (let say nb ) or sections. This justifies the name of sectional method. The
bins correspond to the cells in the classical terminology of numerical analysis.
The objective is then to compute Qji and N j , the integrated values over the bin j for the
mass distribution of species i and the number distribution, respectively.
Splitting of coagulation and condensation/evaporation is often used. For the sake of clarity,
we use splitting in our presentation.

Coagulation Here, we only present the equations for the number distribution. After integra-
tion of the coagulation equation over bins, one typically gets:
j j nb
dN j 1X X j j1 j2 j
X
(t) = X j1 j2 N N − N Xjj1 N j1 (3.27)
dt 2
j1 =1 j2 =1 j1 =1

The coefficients Xjj1 j2 are the key components of the method and describe the fraction of the
aerosol generated by the coagulation of bins j1 and j2 located in bin j. Xj1 j2 is defined by
Pj=nb j
summing over all possible resulting bins: Xj1 j2 = j=1 X j1 j2 .
A popular method is based on Jacobson et al. [1994]; Stratton et al. [1994]. We refer to
Debry and Sportisse [2007b] for a review and a rigorous derivation of such equations. The key
point is to choose the closure scheme inside each bin (the a priori local form of the aerosol
distribution).
The resulting equation is then integrated with an ODE solver. A usual approach is based
on a semi-implicit approach, for instance:
∆t Pj−1 Pj j
j n+1
(N j )n + 2 j1 =1 j2 =1 Xj1 j2 (N )
j1 n+1 (N j2 )n
(N ) = P nb j
(3.28)
1+ ∆t j1 =1 (Xjj1 − Xjj 1
)(N j1 )n

Condensation/evaporation If the bins are fixed, the method is referred as an Eulerian


method. A crucial point is to lower the numerical diffusion (condensation/evaporation may be
40 Chapitre 3 - Research works

viewed as an advection process). Classical methods are used in this aim (for instance, semi-
Lagrangian solvers Dhaniyala and Wexler [1995]; Dabdub and Nguyen [2002]).
An alternative approach is to use Lagrangian methods by following the evolution of bins
according to condensation/evaporation. This method is also referred as a moving sectional
algorithm (Kim and Seinfeld [1990]; Jacobson and Turco [1995]). A motivation is that one
can only use a small number of bins in 3D applications. A crucial point is then to map the
resulting moving distribution onto a fixed grid after having solved the GDE because the aerosol
distribution is given on a fixed grid for 3D processes, such as advection or diffusion. One refers
to Debry et al. [2007] for a method that ensures mass and number conservation.
From (3.19), one gets a system of ODEs that gives the time evolution of Lagrangian variables
(Debry and Sportisse [2006]). For instance for mass distributions:
 
nb
dQji X 0
= N j ai (mj ) Ki − Qji − ηij ceq j j 
i (Q1 , . . . , Qnc ) (3.29)
dt 0
j =1

The mass conservation law for the total mass (gas and aerosol) is written under the form c i +
Pn b 0
j
0
j =1
Q i = Ki .
This system is highly stiff due to the wide range of timescales for condensation/evaporation.
We refer to the previous section for the issue of time integration. In the SIREAM model (SIze
REsolved Aerosol Model, Debry et al. [2007]), a Rosenbrock method is used. Some specific flux
limitations can also be applied in order to reduce the stiffness (especially this due to the fluxes
of H+ , Pilinis et al. [2000]).
We refer to Zhang et al. [1999] for a comparative review of algorithms.

3.3.5 Future works


There are of course many challenging issues in aerosol modeling. I can cite on the basis of
research projects that I have already initiated:

• the modeling of Secondary Organic Aerosol, a strong weakness of current aerosol models;

• the modeling of ultrafine particles at local scales (with the competition between nucleation
and condensation, Figure 3.12);

• the modeling of non-spherical soot aerosols (see for instance Naumann [2003]);

• the extension from internal mixing to external mixing, a key point for the assessment of
radiative effects (Tombette et al. [2007], Figure 3.11).

The numerical simulation of the GDE is also still a difficult problem, especially with the increase
in complexity to arise in near future. See for instance Mallet [2007] for some strange behaviours
of the GDE.

3.4 Inverse modeling and data assimilation


3.4.1 Background
We now assume that the model has been defined (that is to say that parameterizations have
been proposed, Section 3.1) and that appropriate numerical solvers are used (Sections 3.2 and
3.3).
Inverse modeling and data assimilation 41

0.8
55

La_Hague 0.7
Cabauw Belsk
Oostende IFT-Leipzig

Lille
50 Mainz

Karlsruhe
Palaiseau Paris 0.6

Laegeren

Ispra
Venise
45 0.5
Modena
Carpentras
Villefranche
Toulon

Palencia
Barcelona 0.4

40

Cabo_da_Roca
Evora
Messina
0.3

El_Arenosillo Granada
Blida

Lampedusa
35
-10 -5 0 5 10 15 20

Figure 3.11: Average value of the AOT (Aerosol Optical Thickness) at 440 nm, simulated with
the Polyphemus system for year 2005, Tombette et al. [2007].

A field of growing interest is the use of data assimilation techniques. The objective is to couple
model outputs and observational data provided by monitoring networks in order to improve the
forecasts (data assimilation) and/or to reduce the uncertainties of some input parameters (inverse
modeling, for instance of emissions). We refer for instance to Tarantola [1987]; Bennett [2002];
Kalnay [2003] for a general presentation in geophysics.
Let use write the evolution system (1.1) under the form:
dc
= f (c, Ψ) , c(0) = c0 (3.30)
dt
Here, Ψ stands for some uncertain forcing parameters (for instance emissions or boundary con-
ditions). The uncertainties are related to f (through the parameterizations), to the initial
conditions c0 and to the forcing parameters Ψ. Notice that in meteorology, the key point is
the sensitivity with respect to initial conditions, due to the chaotic behavior of atmospheric
dynamics (for instance Kalnay [2003]). In atmospheric chemistry, the sensitivity with respect to
the initial conditions is much lower due to the dissipative nature of chemical kinetics. A typical
spin-up is of magnitude one week at continental scale. Notice that in atmospheric chemistry
the eigenvalues of J = ∂f ∂c have strictly negative real parts. This is not a wave problem as in
meteorology.
Some observational data may be provided by ground stations (for instance the local mon-
itoring networks) or satellital data related to earth observation. Let us assume that some
observations obs1 , . . . , obsno are available at discrete times t1 , . . . , tno . At time tn , the model
output is c(tn ) = F n (Ψ) (we use this notation in order to underline the input/output relation).
The observation is not necessarily the model output and we write H n the observation operator
at time tn : if there were no errors, one should get obsn = H n F n (Ψ).
Data assimilation techniques are based on the reduction of discrepancies between the model
outputs and the observational data, measured by the cost function of the model parameters Ψ
(when required, it could also depend on the initial conditions):
n=n
Xo
Jo (Ψ) = (obsn − H n F n (Ψ))T Rn−1 (obsn − H n F n (Ψ)) (3.31)
| {z }
n=1
Jn
42 Chapitre 3 - Research works

Number Distribution at 10 m Number Distribution at 10 m


16 Winter ; case II ; strong wind 16 Summer ; case II ; strong wind
10 10
14 14
10 10
12 12
Nombre par m d air

Nombre par m d air


10 10
10 10
10 10
3

3
8 8
10 10
6 6
10 10
4 4
10 10
Reference (250 ppm) Reference (250 ppm)
2 2
10 Sulfur / 5 (50 ppm) 10 Sulfur / 5 (50 ppm)
Sulfur / 10 (25 ppm) Sulfur / 10 (25 ppm)
0 0
10 

3 

2 

1 0 1 10 

3 

2 

1 0 1
10 10 10 10 10 10 10 10 10 10
D i am e t r e ( m i c r om e t r e s ) Diametre (micrometres)

Figure 3.12: Simulation of the impact of the reduction in sulfur composition (gasoline) in the
near exhaust of a car. Computed with the CFD model Mercure Saturne and the MAM aerosol
model. Left: winter; right: summer. In Albriet [2007].

Rn is the covariance matrix for observation error (to be specified). A second part of the cost
function can be defined by taking into account the so-called background term (an a priori
knowledge of Ψ):
Jb (Ψ) = (Ψ − Ψb )T B −1 (Ψ − Ψb ) (3.32)

with Ψb the background value and B the background error covariance matrix. Adding J b to Jo
is a penalization technique (Tikonov regularization) that may be necessary for ill-posed (high-
dimensional) problems.
The problem is then to minimize J = Jo + Jb with respect to Ψ ∈ Vadm . This provides
an optimized value Ψ? = argminJ, that can be used to compute, when required, an optimal
forecast c? = F (Ψ? ). Vadm is the space of admissible parameters (for instance, some parameters
have nonnegative values).
There are many applications of this approach in inverse modeling of emissions: for CH 4
(for instance Bergamaschi et al. [2005]), for carbon (Chang et al. [1997] at regional scale and
Bergamaschi et al. [2000] at global scale), for NOx (Quélo et al. [2005a], Figure 3.13) and VOC
(Chang et al. [1996]). Another related problem is the so-called localization problem for point
sources of accidental releases (for instance, radionuclides: see the works of my colleague Bocquet
[2005a,b]).
A brief review of inverse modeling in air quality can also be found in Bocquet and Sportisse [2007].

3.4.2 Sequential versus variational methods

Two families of methods are widely used: the variational algorithms, based on a minimization
procedure (3D-Var when time evolution is not taken into account, 4D-Var when it is), and the
sequential methods, based on the Kalman filter theory.
My works have not been devoted to the development of new methods but to the application
of the existing approaches to CTM, especially the 4D-Var approach (Sportisse and Quélo [2003]
and Sportisse and Quélo [2004]).
Inverse modeling and data assimilation 43

Figure 3.13: Optimized hourly emission factors for NOx (simulation over Lille in Northern
France, Quélo et al. [2005a]).

Variational methods
The variational methods are based on gradient-like minimization methods. They require to
compute ∇Ψ J. The Finite Differences techniques are not useable due to the computational
burden and one often advocates the use of the so-called adjoint models. The CPU cost required
for the adjoint model ranges typically from 3 to 7 times the CPU cost required for the model.
Let us detail the algorithm for a CTM. Here, we assume that the CTM (actually the com-
putation of the cost function) is given by the sequence:
1. initialization of the cost function: J = 0 (or a background term);

2. time loop labelled by 1 ≤ n ≤ no :

• get the forcing fields φn−1 ;


• compute the new state cn = F (cn−1 , φn−1 );
• update the cost function: J = J + Jn (cn ).

The forcing fields φn include the meteorological data, the results of parameterizations, etc.
A direct calculation of ∇J could be performed but the storage requirement would be too
high because (φ0 , . .P
. , φno −1 ) and the trajectory (c0 , . . . , cno ) should be stored. The idea is to
notice that ∇Ψ J = n ∇Ψ Jn , and:
X ∂cm T ∂cn ∂Jn X ∂cn ∂Jn
∇Ψ J = ( )|tm−1 . . . ( n−1 )T|tn−1 ( n )T|tn + ( )T|tn−1 ( n )T|tn (3.33)
∂Ψ ∂c ∂c ∂Ψ ∂c
1≤m<n≤no 1≤n≤no
n n
The adjoint model gives, for any vector v, ( ∂c T ∂c T
∂Ψ )|tn−1 v and ( ∂cn−1 )|tn−1 v. If the differentiation is
performed with respect to the initial conditions (Ψ = c0 ), the formula is:
P 1 n
∇c0 J = 1≤n≤no ( ∂c )T . . . ( ∂c∂cn−1 )T|tn−1 ( ∂J
∂c0 |t0
n T
∂cn )|tn
h
∂J −2 T no −1
T [( ∂Jno −1 )T + ( ∂cno )T ( ∂Jno )T ]
i (3.34)
= . . . × ( ∂cnnoo−2 ) + ( ∂c∂c on −2 ) n
∂c o −1 n
∂c o −1 n
∂c o
44 Chapitre 3 - Research works

8e-06
Computation

7e-06
(d) Measurements

6e-06

Activity [Bq] 5e-06

4e-06

3e-06

2e-06

1e-06

0
-500 -400 -300 -200 -100 0 100 200 300 400
Y [m]

Figure 3.14: Model-to-data comparison for a wind tunnel experiment (Bugey nuclear power
plant). The model outputs are simulated with data assimilation; the controlled inputs are either
emissions and/or “physical” parameters. Observations given by ’x’ and model outputs by ’+’.
In Krysta et al. [2006].

In this case, the algorithm for computing ∇c0 J, is then:

1. initialization: (cno )? = ( ∂J no T
∂cno ) ,

2. reverse time loop labelled by no ≥ n ≥ 1:

• get the forcing fields φn−1 ;


• get the stored trajectory cn−1 ;
• update (cn−1 )? = ( ∂c∂cn−1 )T|tn−1 (cn )? + ( ∂J
n n−1 T
∂cn−1
) .

3. conclude with: ∇c0 J = (c0 )? .

A specific issue is related to the availability of the adjoint model. There are two solutions:
the adjoint model may be computed on the basis of the discretization of the continuous adjoint
model or it may be given by the adjoint code of the discretized model. We refer for instance to
Sirkes and Tziperman [1997].
In many cases, automatic differentiation with appropriate tools is performed (Giering and
Kaminski [1998]; Griewank and Walther [1999] and Mallet and Sportisse [2004]). Notice that
discontinuities exist in the CTM (especially related to aerosol thermodynamics with thresholds)
and may require a special care.
We refer for instance to Elbern et al. [1997]; Elbern and Schmidt [1999]; Elbern et al. [2000];
Elbern and Schmidt [2001] and Quélo et al. [2005a] for the application of variational methods to
air quality modeling at regional scale. Another example is given by the inverse modeling of ra-
dionuclides with a Puff model (in the context of a wind tunnel experiment, Krysta et al. [2006],
Figure 3.14).
Inverse modeling and data assimilation 45

In the linear case, when the chemical production term is given by a linear term (χ(c) = χ c,
with χ a matrix), a simplification can be used. For the sake of clarity, we consider the continuous
RT
case. The cost function is then given over an assimilation window [0, T ] as J(c) = 0 j(c(t)).
For instance, we suppose that the control parameter is the source term S. Let us define the
adjoint variable c? governed by (we omit density for the sake of clarity):

∂c?i
− − V div(c?i ) = div(K∇c?i ) + (χT c? )i − Λi c? i + ∇c j , c?i (T ) = 0 (3.35)
∂t
When the wind field is divergence free, one gets:

∂c?i
− − div(V c?i ) = div(K∇c?i ) + (χT c? )i − Λi c? i + ∇c j , c?i (T ) = 0 (3.36)
∂t
The key result is that ∇S(t) J = c? (t). The initial model can be used in order to compute the
sensitivity by modifying the input parameters and performing a backward integration (compare
(1.1) and (3.36)). The so-called retroplume methods (for instance Hourdin and Issartel [2000]
and Krysta et al. [2004]) are based on this approach.

Sequential methods
The sequential methods are based on the Kalman filter theory (for instance Tarantola [1987]).
We distinguish the forecast and the analysis obtained by taking into account the observational
data. Both are supposed to be stochastic variables that follow normal laws, N(c f , Pf ) and
N(ca , Pa ), respectively. P is a covariance error matrix. The model is supposed to be given by
cn+1 = F (cn ) from tn to tn+1 , with an unbiased model error following N(0, Q).
The most basic algorithm (BLUE for Best Linear Unbiased Estimator and the extension to
the nonlinear case with the Extended Kalman Filter) is given by the following sequence:

1. Forecast from tn to tn+1 :


cn+1
f = F (cna ) (3.37)

2. Propagation of the error covariance matrix:

Pfn+1 = Mn,n+1 Pan Mn,n+1


T
+ Qn (3.38)

3. Analysis at tn+1 :  
cn+1
a = cn+1
f + K n+1
obs n+1
− H n+1 n+1
c f (3.39)

where the gain matrix is:


 −1
K n+1 = Pfn+1 (H n+1 )T H n+1 Pfn+1 (H n+1 )T + Rn+1 (3.40)

4. Computation of the analysis covariance matrix at tn+1 :



Pan+1 = I − K n+1 H n+1 Pfn+1 (3.41)

Mn,n+1 is the linear tangent model from time tn to tn+1 . normal PDF.
This method is easy to implement but a key issue is the large dimension for reactive problems
and the computational burden associated to the step (3.38). A powerful technique is given by
the use of reduced filters.
46 Chapitre 3 - Research works

Notice that the covariance error matrices can be written as P f = Sf SfT , Pa = Sa SaT , Q =
Sq SqT and R = Sr SrT . It is straightforward to check that Sf = (M Sa , Sq ), K = Sf A(AT A +
R)−1 and Sa = Sf (I − A(R + AT A)−1 AT )1/2 with AT = HSf . The advantage of using the
square matrices S rather than the covariance matrices P is that these matrices are symmetric
and strictly positive. Moreover, with S = (s1 . . . sns ) (ns depends on the matrix, si is a column),
one gets:
Xns
SS T = si si T (3.42)
i=1
Once the sTi si values have been ranked in a decreasing order, it is easy to define truncated
matrices related to S.
For the RRSQRT filter (Verlaan [1998]), the modes are defined from the leading eigenvalues
of SaT Sa (that has the same eigenvalues as Pa ). For the SEEK (Singular Evolutive Extended
Kalman, Pham et al. [1998]) and the SEIK (Singular Evolutive Interpolated Kalman, Pham
[1996]) filters, the truncation is based on the Empirical Orthogonal Functions (EOF) of a com-
puted trajectory.
An alternative approach is to use an ensemble filter by computing the covariance matrix on
the basis of a Monte Carlo method:
j=N
1 X (j)
Pf ' (c − c̄)(c(j) − c̄)T (3.43)
N −1
j=1

with an ensemble of N possible states c(j) of average c̄.
We refer for instance to Segers [2002] for a benchmark of Kalman filters applied to air
pollution modeling and also to Constantinescu et al. [2006] for the application of Ensemble
Kalman filters.
The Polyphemus system (Mallet et al. [2007c]) includes 4D-Var, optimal interpolation, the
RRSQRT and EnKF methods, applied to the generic class of models available in the platform.
These methods are currently benchmarked for the simulation of air quality over Europe (Figure
3.15).

3.4.3 Coming issues


The sensitivity of the results produced by a data assimilation procedure is a key question. The
sensitivity can be computed with respect to three items:
• the physical behaviour of the model (for instance its slow/fast behaviour): in this case,
this is related to the conditioning of the optimization system (is it ill-posed ?);
• uncertain parameters that are not optimized: one usually refers to the second-order sen-
sitivity (the sensitivity of the optimized system, not the one of the forward model);
• the observational network: this is related to the topics of network design (how to optimize
the observational network ?).
These topics require the study of the Hessian matrix related to the cost function.

Slow/fast behaviour
The impact of the slow/fast behaviour of chemical kinetics (Section 3.6) is investigated in
Sportisse and Mallet [2007]. A theoretical framework is given in order to study the dynami-
cal issues for data assimilation. It is proven that the data assimilation problem is ill-posed for
the fast components (fast species). One also refers to Section 3.6.
Inverse modeling and data assimilation 47

Montandon
140
reference
observation
Optimal Interpolation
EnKF
120
4DVar
RRSQRT

100

Concentration
80

60

40

20
0 10 20 30 40 50
time

Figure 3.15: Comparaison of different data assimilation methods for the ozone concentration:
without assimilation or with Optimal Interpolation, RRSQRT, EnKF or 4D-Var methods. Sta-
tion of Montandon (continental simulation over Europe with Polyphemus in summer 2001).
Credit: Lin Wu, CEREA/CLIME Project, INRIA.

Second-order sensitivity
The first issue is related to the so-called second-order sensitivity. Once the optimized parameter
Ψ? (for instance emissions) has been computed, one can investigate its robustness with respect
to other uncertainties. For instance, if Φ stands for uncertain parameters that are not optimized
(due the large dimension for instance), Ψ? is defined as a function of Φ:

Ψ? (Φ) = argminJ(Ψ, Φ) (3.44)

The second-order sensitivity (Wang et al. [1992]) is defined as:

∂Ψ ∂∇Ψ J
= − (Hess J)−1 (3.45)
∂Φ ∂Φ
with the use of the implicit function theorem. We refer for instance to Quélo et al. [2005b] with
the application to a simple Gaussian model for the atmospheric dispersion of radionuclides.
Computing the Hessian matrix of J may be a difficult task for high-dimensional systems. A
brute-force approach is then to investigate the sensitivity of Ψ? with respect to modified model
configurations or model parameters. An example is given by the inverse modeling of NO x at re-
gional scale and the sensitivity analysis with respect to parameterizations in Quélo et al. [2005a].
Notice that this point is crucial in order to use the optimized parameters for themselves and
not only for forecast.

Network design
An interesting question is related to what is usually called network design: how to optimize a
mobile or a fixed monitoring network (related to the observational operator H) ? This issue is a
classical issue, for instance in meteorology. Mathematically speaking, the quality of a monitoring
48 Chapitre 3 - Research works

network can be estimated by the condition number of the Hessian matrix of J (that depends on
the network through H).
We refer for instance to Sportisse and Quélo [2002]. A preliminary application to air pollu-
tion modeling can also be found in Sandu [2006] (Targeting) and in Krysta et al. [2006] (Gen-
eralized Cross Validation).

3.4.4 Future works


The two last issues are key issues in air quality modeling. The assessment of the quality of
optimized emission fluxes is required for any practical use. Network design may be related to
many applications. Many projects have been already initiated under my leadership: monitoring
of accidental releases (Korsakissok and Sportisse [2006]) or of radionuclides (Vercauteren and
Bocquet [2006]). The economical issues are of course of great interest. The extension to air
quality (photochemistry) is a challenging issue.
The benchmark of data assimilation methods for air quality is still a challenging point. It is
the subject of a Research Action at INRIA that joints 3 projects: ADOQA for Data Assimilation
for Air Quality (with the search for non-linear methods). Another interesting point could be to
optimize with respect to other “appropriate” parameters (for instance eddy coefficient, emissions,
dry deposition velocities, kinetic rates, etc).
The extension to multiphase models (aerosols) is another logical follow-up of these works.
The opportunity of using radiative observations provided by LIDAR (at regional scale) and
satellites (at global scale) are currently investigated. Data assimilation with LIDAR data is the
subject of a research project jointly initiated with Patrick Chazette (LSCE); data assimilation
of satellital data is investigated in a project of the European Space Agency (MetOp platform),
for which I am the Principal Investigator.

3.5 Sensitivity analysis, ucertainty propagation and ensemble


forecast
3.5.1 Background
Air quality modeling and simulation suffer from many uncertainties (for instance Russell and
Dennis [2000]):

• many input data are poorly known (see for instance Table 3.1), especially emissions. More-
over, some forcing fields, such as meteorological fields, are computed with numerical models
that may be uncertain.

• the models are highly parameterized (see Section 3.1).

• the numerical algorithms and the discretization also induce uncertainties (remember that
the numerical resolution may be coarse, due to the computational burden, especially for
aerosols).

It is therefore not relevant to view outputs of CTMs as deterministic values (Mallet and Sportisse [2006a]).
Even if the models are “validated” (it is more rigorous to say that model-to-data comparisons
are performed, when possible), one must keep in mind that there are a large amount of degrees
of freedom (especially in parameterizations) and only a small number of model outputs can
be measured. For instance, most of the existing CTMs have been extensively tuned to meet
acceptable model-to-data error statistics for ozone peaks at ground. It does not ensure that the
Sensitivity analysis, ucertainty propagation and ensemble forecast 49

variable range (95% of data) σ (log-normal)


ozone initial conditions 3 0.549
NOx/VOC initial conditions 5 0.805
major point NOx/VOC emissions 1.5 0.203
rainfall amount 2 0.347
cloud liquid water content 2 0.347
biogenic NOx/VOC emissions 2 0.347

Table 3.1: Uncertainties for some input parameters of CTM (parameters of log-normal proba-
bility density functions, following Hanna et al. [2001].)

results are satisfactory for 3D fields and many other trace species. There is therefore a danger
of using “overtuned” models, especially for impact studies or long-term scenario studies.

3.5.2 Some possible strategies for assessing the impact of uncertainties


Sensitivity analysis
A first strategy is to evaluate the sensitivity of some model outputs, let say c = F (Ψ), with
∂F
respect to uncertain input parameters, let say Ψ, through the partial derivative .
∂Ψ
This approach has been widely used with many variations in the way the sensitivities are
computed. We can cite among many other works Seigneur et al. [1982]; Carmichael and al
[1997]; He and al [2000]; Daescu et al. [2003]; Sandu et al. [2005] for a general presentation,
Pandis and Seinfeld [1989]; Zhang et al. [1998a] and Djouad et al. [2003a] for aqueous-phase
chemistry, Mallet and Sportisse [2005] for the sensitivity of ozone with respect to emissions,
Zhang et al. [2005] for the sensitivity of ozone and Roustan and Bocquet [2006] for the sensitivity
of atmospheric mercury at European scale with respect to the boundary conditions, etc.
Notice that this item is also related to variational data assimilation (see Section 3.4) because
the same techniques can be used to compute gradients and sensitivities.
This approach provides an estimation of the local uncertainty.

Monte Carlo simulations


A second approach is based on a classical Monte Carlo simulation on the basis of Probability
Density Functions (PDF) for the uncertain inputs.
The application to box models (for instance chemical kinetics in Rodriguez and Dabdub
[2003]) is not too computationally demanding. In 3D applications, a major drawback is that
the number of possible runs is usually low (typically 50 or 100). Due to the large dimension
of systems, the brute-force approach is not feasible and some appropriate methods have to be
used. A very attractive method is the so-called DEMM (Deterministic Equivalent Monte-Carlo
Method, Tatang et al. [1997]) which is based on polynomial expansions of the targeted outputs.
For instance, for the simple case of one input (Ψ ∈ R):
i=N
X
c = α0 + αi Pi (Ψ) (3.46)
i=1

where Ψ stands for the uncertain input parameter, Pi are orthogonal polynomial functions
(“chaos” expansion) and N is the order of the approximation.
The coefficients (αi ) are computed at collocation points based on the PDF for Ψ (we refer
to Tatang et al. [1997]). The essence of the algorithm is to simulate the model response surface
50 Chapitre 3 - Research works

Figure 3.16: An ensemble for hourly profiles of ozone (averaged over Europe, summer 2001),
generated by the Polyphemus system. The points indicate the result of the ensemble forecast
for an ensemble algorithm. In Mallet and Sportisse [2006b].

by using orthogonal polynomials of the uncertain input parameters. Using (3.46) is much more
efficient than using the exact model.
The application of this method to simulations over Europe for uncertainties in NO x , NH3
and SO2 emissions can be found in Boutahar [2004]; Boutahar and Sportisse [2002].

3.5.3 Multi-modeling
As pointed out before, a key point is to avoid “overtuning” of models due to the large amount
of output data and the small number of observed data. For instance, a typical version of a
sectional aerosol model may include up to hundreds of species per grid cell (let say 10 bins times
10-20 chemical species per grid cell). In comparison, observational data are routinely available
for aggregated values, such as PM2 .5 or PM1 0 (PMx is the mass of aerosols whose aerodynamic
diameter is less or equal to x micrometers) but not for size distribution or chemical composition.
A promising approach is not to rely on one single model but to use a set (an ensemble) of
models or of model configurations to deliver the forecast (Figures 3.16 et 3.17). We refer for
instance to Mallet and Sportisse [2006d] for the study of model spreads due to uncertainties in
model inputs, numerical algorithms and parameterizations.
The extension to multiphase models can be found in Tombette et al. [2005] (at regional scale,
over Paris) and in Sartelet et al. [2007c] (at continental scale, in the framework of the MICS
Asia exercice).

3.5.4 Ensemble forecast


The next step is to combine the results of the ensemble in order not only to assess the uncer-
tainties but also to improve the forecast. We now detail an interesting approach based on what
is usually called the superensemble technique (Krishnamurti et al. [1999]).
Let M be the set of available models, labelled by m. The objective of ensemble forecast is
to compute a forecast Mt,x on the basis of the model outputs Mm,t,x for model m (at time t and
Sensitivity analysis, ucertainty propagation and ensemble forecast 51

Figure 3.17: Ensemble relative variance (in %) for ozone, summer 2001 (4 months). In
Mallet and Sportisse [2006b].

at position x). The simplest approach is to use the ensemble mean (EM in the sequel):

1 X
EMt,x = Mm,t,x (3.47)
|E|
m∈E

We refer for instance to Galmarini et al. [2004] with the application to the forecast of radionu-
clides. Another similar approach is the use of the median model.
A promising approach is to take into account observations. Let obst,x be the observation
corresponding to the model output Mt,x . An improved forecast may be given as a linear com-
bination of model outputs as: X
ELSt,x = αm Mm,t,x (3.48)
m
P P 2
where α minimize the cost function t,x [obst,x − m αm Mm,t,x ] which measures a model-
to-data discrepancy. The minimization can be performed in the Least Square sense (ELS
stands for Ensemble Least Square). The applications to air quality modeling can be found
in Mallet and Sportisse [2006b] and Pagowski et al. [2005].
A challenging point is to forecast the weights α by only using past observational data.
Some techniques based on machine learning and/or expert selection (some popular tools in
statistics) can be used. One of the simplest methods uses a so-called loss function L t (αt ) =
!2
X
αm,t Mm,t − obst . The weights can be updated with a gradient algorithm (Cesa-Bianchi
m
et al. [1996]) αt = αt−1 − η∇Lt−1 (αt−1 ). We refer for instance to Mallet and Sportisse [2006b].

3.5.5 Future works


These works can be extented to other species and multiphase models. The ensemble calibration
is a key point (through indicators such as the Brier skill or the rank diagram, some usual
indicators used in meteorology).
52 Chapitre 3 - Research works

10000

NH3
NO3
1000 Cl

condensation time scale (sec)


100

10

0
0.001 0.01 0.1 1 10
aerosol dry diameter (µm)

Figure 3.18: Timescales for the condensation of inorganic species (ammonium, nitrate, chloride)
as a function of the aerosol size (Debry and Sportisse [2006]).

Using ensemble forecast for the input meteorological fields could be also relevant.
There are many open issues at the theoretical level: how to define spatialized ensemble
forecast ? How to couple to data assimilation methods ? etc.

3.6 Proxy and reduced models


3.6.1 Motivations
An important application of Chemistry-Transport Models is the long-term estimation of the
impacts of emissions reduction. This is for instance illustrated by the CAFE program in Europe
(Clean Air For Europe). In such applications, the physical model is only a small part of a
modeling chain that include models for economic activities (and then emissions) and monetary
evaluation of impacts due to air pollution. Such a modeling chain is used in an optimization
mode in order to find the optimal partitioning of emissions that minimizes a global economic
cost. This implies that the physical model (the CTM) has to be used many times.
It is of course not possible to use comprehensive detailed models. A classical tool is provided
by the so-called source-receptor matrices, that are based on the linearization of the function
c = F (E) (c are the concentrations, E are the emissions, F is the model viewed as a black box).
We refer to Bartnicki [2000] for the evaluation of the resulting errors.
An alternative technique is to use highly reduced models. We now present some methods
than can be used for the reduction of such high-dimensional systems.

3.6.2 Slow/fast models and dynamic reduction


A first approach is related to the wide range of timescales that characterizes atmospheric chem-
istry (Figure 3.18). Let ε  1 be the ratio of fast timescales to slow timescales. We first assume
that the system can be partitioned in a slow part (cs ) and a fast part (cf ). An ideal system is
Proxy and reduced models 53

ref
no solver
4 bulk fixed
bulk moving
iter fixed

relative error (adim)


iter moving
3 bfgs moving

0
0 0.5 1 1.5
cpu time (sec)

Figure 3.19: CPU/RMSE (relative error) diagram for different reducing methods applied to
condensation/evaporation (Debry and Sportisse [2006]).

then (after scaling):


dcf dcs
ε = f (cf , cs ) , = g(cf , cs ) (3.49)
dt dt
Here, f and g stand for the fast and slow dynamics, respectively. A key result given by the sin-
gular perturbation theory (Tikhonov et al. [1985]; see also Sportisse [1997] for the mathematical
framework) is that this system may be approximated up to first-order in ε after a fast transient
phase by the following differential-algebraic system:

dcs
0 = f (cf , cs ) , = g(cf , cs ) (3.50)
dt

This is the theoretical framework for the QSSA methods (Section 3.2). The key point is of
course to partition the system under the slow/fast form. Simple criteria can be derived from the
production-loss form of chemical kinetics (Djouad and Sportisse [2002a]). The use of reduced
models (that are no more stiff) is investigated in Sportisse and Djouad [2000] and Lowe and
Tomlin [2000]; Heard et al. [1998]; Kalachev and Field [2001]; Neophytou et al. [2004]; Lovas
et al. [2006] for gas-phase chemistry and Djouad et al. [2003b] for multiphase chemistry. We
also refer to Djouad and Sportisse [2003] for the numerical solution of the resulting model.
For aerosol modeling, the fast dynamics is related to condensation/evaporation. It appears
that the finest aerosols quickly reach equilibrium for mass transfer. This leads to the definition
of the so-called hybrid methods (Capaldo et al. [2000]): the smallest aerosols whose diameter
is below a cut-off diameter are supposed to be at equilibrium (they define the fast part) while
the coarsest aerosols are solved with kinetic mass transfer (they define the slow part). We refer
to Debry and Sportisse [2006] for a benchmark of different methods and the application to the
GDE (Figure 3.19).
54 Chapitre 3 - Research works

3.6.3 Building look-up tables


An alternative approach is to replace the most computationally demanding modules of CTMs
by look-up tables. There is of course a problem of dimension. If there are n model inputs and
s sampling points, the table is of dimension O(sn )!
A powerful technique is provided by the multivariate expansions, a tool used in many other
fields. One example is for instance the HDMR method (High Dimensional Model Representation,
Shorter and Precila [1999]; Rabitz and Alis [1999]).
Let c = F (Ψ) be the module to be replaced by a look-up table. The Ansatz is that one can
write the multivariate exact expansion:

i=n
X X
F (Ψ) = F0 + Fi (Ψi ) + Fij (Ψi , Ψj )
i=1 1≤i<j≤n
X (3.51)
+ Fijk (Ψi , Ψj , Ψk ) + . . .
1≤i<j<k≤n
+F1,2,...,n (Ψ1 , Ψ2 , ..., Ψn )

Fi describe the independent action of Ψi , Fij the cooperative action of Ψi and Ψj , etc. A reduced
model is then to use a truncated development (of second order or of third order in practice).
The functions Fi are computed such that the truncation error is minimal for a given norm, over
the space of functions that satisfy the Ansatz. This procedure can be repeated for each order.
A crucial point is the choice of the error metrics. We refer to Wang et al. [1999] and
Boutahar [2004]; Boutahar and Sportisse [2002] for the application to atmospheric chemistry.

3.6.4 POD/EOF methods


Another method is based on the so-called Proper Orthogonal Decomposition (POD, Berkooz
et al. [1993]), also known as Empirical Orthogonal Functions in atmospheric modeling (see for
instance Fiore and Jacob [2003] for an application to the understanding of the surface ozone
patterns over the United States). The idea is to map the initial model, let say given by an
Nr
evolution equation dc dt = f (c) onto a reduced space, span{Φi }i=1 :

i=N
Xr
c= (c, Φi )Φi (3.52)
i=1

where the POD vectors Φi are orthonormal with respect to a given scalar product (., .). The
POD basis is defined such that for each 1 ≤ k ≤ Nr , (Φ1 . . . , Φk ) minimizes the residual Jk for
0
any subspace of dimension k, span{Φi }ki=1 :

Z T i=k 
X  0
0 0 0
Jk (Φ1 , . . . , Φk ) = kc(t) − Φi , c(t) Φi k2 dt (3.53)
0 i=1

along a trajectory c(t). The norm k.k is related to the scalar product (., .). A classical result is
that the POD basis can be computed on the basis of the eigenvalues of the correlation matrix
[(c(ti ), c(tj ))]ij . The choice of the scalar product is once more a challenging point.
In Sportisse and Djouad [2007], the application to comprehensive chemical mechanism for
tropospheric chemistry indicates that there are only 2 or 3 local degrees of freedom for a species
like ozone. This result is promising but the extension to 3D modeling is still an open question.
Chemistry-Transport Models and applications 55

3.6.5 Future works


The reducing methods are a posteriori powerful tools in order to understand the model be-
haviour. The most illustrative example is the slow/fast behaviour of chemical kinetics (with
many impacts for numerics, data assimilation, etc).
It is difficult to use in practice reduced models in 3D models due to the constraints of validity.
The search for other reducing methods is still a challenging point, especially for integrated
modeling (from economic activity to impact assessment). The integrated models have the advan-
tage of having a small set of outputs (typically a cost function) which means that the reduced
model has not to be accurate for 2D or 3D fields. This work is under progress in a joint
project with colleagues from CIRED (Centre International de Recherche sur l’Environnement
et le Développement, Jean-Charles Hourcade).

3.7 Chemistry-Transport Models and applications


Apart from these research activities, my works have also focused on the development of CTM
and on applications to real cases. The focus is then less academic but these works have strongly
structured my activities in the last five years. I briefly summarize them in this Section.

3.7.1 Development of CTM


Polair3D (Boutahar et al. [2004]) is the first CTM that has been developed at CEREA under
my leadership. The motivations were to have a modular model (able to deal with many dif-
ferent applications) and to use automatic differentiation tools for dealing with variational data
assimilation.
This model has first been applied to simple models for acid rains over Europe (on the basis of
the first EMEP models). It has then be extended to photochemistry (Mallet and Sportisse [2006d]
for instance), to aerosol modeling with size-resolved approaches (Debry et al. [2007]) or modal
approches (Sartelet et al. [2005]), to mercury and heavy metals (Roustan et al. [2006]) and to
radionuclides (Quélo et al. [2007]).
The methods related to variational data assimilation have also been developed (Quélo et al. [2005a]
for an application to inverse modeling). Multi-modeling simulations and ensemble forecast have
then completed the system (Mallet and Sportisse [2006b,d]).
These works have been made perennial through the development of the Polyphemus sys-
tem (Mallet et al. [2007a,c]). Polyphemus is a platform hosting different models (Polair3D;
another CTM: Castor based on Chimère, Schmidt et al. [2001]; short-range dispersion models)
for different kinetics (photochemistry, aerosols, heanvy metals, mercury, radionuclides, biological
tracers) and different uses (forecast, ensemble forecast, uncertainties propagation, data assimi-
lation).

3.7.2 Applications
I have been implied in many applications of the developed models.
I refer at regional scale to air pollution modeling over Paris (Sartelet et al. [2002] and
Tombette and Sportisse [2007]), Lille (Quélo et al. [2004] and Quélo et al. [2005a]), Marseille
(Taghavi et al. [2004]) and Tokyo (Sartelet et al. [2007b]). I have already cited many works de-
voted to the continental scale, to which I can add Sartelet et al. [2007c] for Asia and Sartelet et al. [2007a]
for Europe (Figures 3.2 and 3.3). For both cases it was related to a benchmark exercice.
I have also been strongly implied in the model development for the forecast of the dispersion
of radionuclides. This is the subject of a joint project with the Emergency Center of the French
56 Chapitre 3 - Research works

Database Stations Obs. Mod. RMSE corr MFB MFE


mean mean
NO2 EMEP 35 7.5 9.0 5.7 50.0 22 59
AirBase 1049 23.9 15.3 14.2 57.2 -37 62
BDQA 75 22.1 13.4 14.6 55.9 -47 70
NH3 EMEP 3 7.4 6.3 5.5 28.8 12 52
AirBase 7 12.7 6.6 13.7 27.3 -23 101
HNO3 EMEP 7 0.7 1.3 1.4 26.2 40 89
SO2 EMEP 43 2.0 5.3 4.9 46.7 97 106
AirBase 992 6.5 7.3 6.6 44.2 25 70
BDQA 10 7.8 6.8 6.5 35.7 -13 59

Table 3.2: Model-to-data comparison for Polyphemus over Europe in 2001 for gaseous species:
number of stations (from the networks EMEP, Airbase and BDQA), observed mean (µg m −3 ),
simulated mean (µg m−3 ), RMSE (µg m−3 ), correlation (%), MFB (%), MFE (%). In
Sartelet et al. [2007a].

Database Stations Obs. Mod. RMSE corr MFB MFE


mean mean
PM10 EMEP 26 16.9 15.4 12.5 55.1 -9 50
AirBase 537 25.4 15.2 17.0 44.9 -45 59
BDQA 23 19.8 15.5 9.6 57.7 -27 41
PM2.5 EMEP 17 12.6 8.3 8.6 54.4 -40 61
Sulfate EMEP 57 2.5 2.0 1.7 55.6 -6 50
AirBase 11 1.9 2.3 1.7 49.4 39 65
Nitrate EMEP 14 2.6 4.0 3.1 41.3 30 75
AirBase 8 3.5 4.3 2.7 71.7 6 54
Amm. EMEP 9 1.8 2.0 1.3 51.9 19 49
AirBase 8 1.8 2.0 0.9 74.4 14 36
Sodium EMEP 3 1.3 2.4 2.2 62.8 47 68
Chloride AirBase 7 0.9 3.1 3.5 69.8 83 102

Table 3.3: Model-to-data comparison for Polyphemus over Europe in 2001 for aerosols: num-
ber of stations (from the networks EMEP, Airbase and BDQA), observed mean (µg m −3 ),
simulated mean (µg m−3 ), RMSE (µg m−3 ), correlation (%), MFB (%), MFE (%). In
Sartelet et al. [2007a].
Chemistry-Transport Models and applications 57

1986-04-26T13:03:00

Chernobyl

1E-01 1E+00 1E+01 1E+02

Figure 3.20: Simulated cesium 137 for 26 April 1986 at 3:00 (TU) in the framework of the
Chernobyl simulation (Quélo et al. [2007]).
58 Chapitre 3 - Research works

102

Polyphemus [Bq.m^{-3}]
101

0
10

1
10


2
10


3
10


1
10 3
10 2
10 100 101 102


 

Measurements [Bq.m^{-3}]

Figure 3.21: Model-to-data comparison (scatter plot) for cesium 137 for the first two weeks of
the Chernobyl release (Quélo et al. [2007]).

Institute of Radiological Protection and Nuclear Safety (IRSN, France). I refer for instance to
Quélo et al. [2007] (Figures 3.20 and 3.21) for the evaluation of the Polyphemus system that
constitutes the core of the coming forecast model at IRSN.

To date, the developed models are used by EDF R&D for impact studies of power plants
(for instance Tombette and Sportisse [2006] or Roustan et al. [2007], Figure 3.22), by IRSN for
its emergency center and by other partners of CEREA (EDF Polska, Meteo Chile, Technical
Center of the French Ministry for Defence/DGA, some academic teams, etc).

3.7.3 Future works


I refer to the previous topics which are related to CTM. A key point is the extension to transconti-
nental long-range transport (see for instance, among many other references, Jacob et al. [1999];
Hudman and al [2004] for transpacific transport or Li et al. [2002]; Guerova et al. [2006] for
transatlantic transport).
Chemistry-Transport Models and applications 59

PM2.5 mean 1.50

55 1.25

1.00

50

0.75
lat

45
0.50

0.25
40

0.00

35
-10 -5 0 5 10 15 20
-0.25
lon

Figure 3.22: Impact for the year 2005 of the emissions from the French thermal power plants
(EDF and SNET). The target is PM2.5 . The impact is a relative impact computed in % by
comparing a reference simulation and a simulation without the emissions (Roustan et al. [2007]).
60 Chapitre 3 - Research works
Conclusion

This report summarizes my research works related to air quality modeling and simulation for
the period 2000-2006.
Many points are of course still open. Some have been briefly presented and are follow-up of
my previous works. I focus here on five challenging topics for Chemistry-Transport Models.
In near future, a first trend will be the increase of the spatial resolution. At continental scale,
the current model resolution range from 10 to 100 kilometers and will go down to 1-10 kilometers.
This point is of course related to the improvement of the meso-scale meteorological models and
of the emission inventory. This will induce many new open questions at the numerical level, at
the parameterization level and at the data level. One can even wonder if a CTM will not use
meshes similar to those used for CFD (with millions of cells and unstructured adaptative grids).
The second point is related to the improvement of the “chemical” resolution. The complexity
of physical models, especially for aerosols, will increase: a better understanding of Secondary
Organic Aerosol (SOA) implies a better representation of organic species (more variables, Griffin
et al. [2002]), a coupling between organic and inorganic thermodynamics (more computational
burden, Pun et al. [2002]) and new processes (surface heterogeneous chemistry, polymeriza-
tion, etc). Apart from modeling issues, this will require the search for appropriate numerical
algorithms.
A third point is related to the couplings and the feedbacks. This concerns for instance
the coupling between radiative transfer models, CTM and meso-scale models for the radiative
behavior of the atmosphere. On-line coupling will be the standart. Another example is provided
by detailed multi-media models in the estimation of environmental impacts through multi-media
models.
One must keep in mind that the final use of such models always needs a drastic reduction
of the CPU costs. Finding the appropriate trade-off between highly detailed physical models
(the “arms race”) and tailored models that can be used in integrated modeling for impact
assessment will still remain a key issue. The current situation is that one often uses low-level
models and sometimes detailed models for such applications. Models of intermediate complexity
or algorithms for using detailed models are therefore of great interest. The fourth point is
therefore the search for reduced models to be used for integrated modeling and to be derived
with a rigourous basis.
The last point is related to uncertainties. All-in-one models, tuned to a small set of target
species, should be replaced by platforms able to generate a large set of models and model con-
figurations, not only restricted to a small set of models that are based on the same assumptions.
A key issue will be the ability of computing Probability Density Functions for most outputs.

61
62 Chapitre 3 - Research works
Bibliography

Albriet, B. (2007). Modélisation numérique des aérosols à l’échelle locale. PhD thesis, ENPC.

Babovsky, H. (1999). On a Monte Carlo scheme for Smoluchowski’s coagulation equation. Monte
Carlo Methods and Applications, 5(1):1–18.

Bartnicki, J. (2000). Nonlinear effects in the source-receptor matrices computed with the EMEP
Eulerian acid deposition model. Technical Report 4/2000, EMEP/MSC-W.

Bennett, A. (2002). Inverse Modeling of the ocean and atmosphere. Cambridge University Press.

Bergamaschi, P., Hein, R., Heimann, M., and Crutzen, P. (2000). Inverse modeling of the global
CO cycle. 1. Inversion of CO mixing ratios. J. Geophys. Res., 105(D2):1909–1927.

Bergamaschi, P., Krol, M., Dentener, F., Vermeulen, A., Meinhardt, F., Graul, R., Ramonet,
M., Peteres, W., and Dlugokencky, E. (2005). Inverse modelling of national and european
CH4 emissions using the atmospheric zoom model TM5. Atmos. Chem. Phys., 5:2431–2460.

Berkooz, G., Holmes, P., and Lumley, J. (1993). The proper orthogonal decomposition in the
analysis of turbulent flows. Annu. Rev. Fluid Mech., pages 539–575.

Berkvens, P., Botchev, M., Krol, M., Peters, W., and Verwer, J. (2002). Solving vertical trans-
port and chemistry in air pollution models. In Chock, D. and Carmichael, G., editors, Atmo-
spheric Modeling, IMA volumes in Mathematics and its interactions. Springer Verlag.

Binkowski, F. and Roselle, S. (2003). Models-3 Community Multiscale Air Quality


(CMAQ) model aerosol component. 1. Model description. J. Geophys. Res., 108(D6-
4183):doi:10.1029/2001JD001409.

Bocquet, M. (2005a). Reconstruction of an atmospheric tracer source using the principle of


maximum entropy I : Theories. Quart. J. Roy. Meteor. Soc., 131(Part B(610)):2191–2208.

Bocquet, M. (2005b). Reconstruction of an atmospheric tracer source using the principle of


maximum entropy II : Applications. Quart. J. Roy. Meteor. Soc., 131(Part B(610)):2209–
2223.

Bocquet, M. and Sportisse, B. (2007). Modélisation inverse pour la qualité de l’air : éléments
de méthodologie et exemples. Pollution Atmosphérique. Accepted for publication.

Boutahar, J. (2004). Reduction and evaluation of uncertainties in Chemistry-Transport Models.


PhD thesis, Ecole Nationale des Ponts et Chaussées. Sous la direction de Bruno Sportisse.

Boutahar, J., Lacour, S., Mallet, V., Quélo, D., Roustan, Y., and Sportisse, B. (2004). Devel-
opment and validation of a fully modular platform for numerical modelling of air pollution:
Polair3D. Int. J. Env. and Pollution, 22(1-2).

63
64 BIBLIOGRAPHY

Boutahar, J. and Sportisse, B. (2002). Reduction methods for atmospheric chemistry. In Global
and Regional Atmospheric Modeling Meeting, pages 133–139. Un. of Aveiro, Portugal.

Boutahar, J. and Sportisse, B. (2005). Reduction methods and uncertainty propagation: appli-
cation to a Chemistry-Transport-Model. In Proceedings of the TAM-TAM conference.

Byrne, G. and Hindmarsh, A. (1987). Stiff ODE solvers: a review of current and coming
attractions. J. Comp. Phys., 70:1–62.

Byun, D. and Lee, S. (2002). Numerical solution of trace species advection under non-
uniform density distribution: experiment with two-dimensional linear flows. In Chock, D.
and Carmichael, G., editors, Atmospheric Modeling, IMA volumes in Mathematics and its
interactions, pages 109–151. Springer.

Capaldo, K., Pilinis, C., and Pandis, S. (2000). A computationally efficient hybrid approach for
dynamic gas/aerosol transfer in air quality models. Atmos. Env., 34:3617–3627.

Carmichael, G. and al (1997). Sensitivity analysis for atmospheric chemistry models via auto-
matic differentiation. Atmos. Env., 31(3):475–489.

Cesa-Bianchi, N., Long, P. M., and Warmuth, M. K. (1996). Worst-case quadratic loss bounds
for prediction using linear functions and gradient descent. IEEE Trans. Neural Net., 7(3):604–
619.

Chang, M., Hartley, D., Cardelino, C., and Chang, W.-L. (1996). Inverse modeling of biogenic
emissions. Geophys.Res.Lett., 23:3007.

Chang, M., Hartley, D., Cardelino, C., Haas-Laursen, D., and Chang, W.-L. (1997). On using
inverse methods for resolving emissions with large spatial inhomogeneities. J. Geophys. Res.,
102(D13):16023–16036.

Constantinescu, E., Sandu, S., Chai, T., and Carmichael, G. (2006). Investigation of ensemble
based chemical data assimilation in an idealized setting. Atmos. Env. Accepted for publication.

Coron, F., Pallegoix, J., and Sportisse, B. (1994). DSMC computations of complex test cases
and industrial configuration. In 19th International Symposium on Rarefied Gas Dynamics,
Oxford.

Dabdub, D. and Nguyen, K. (2002). Semi-lagrangian flux scheme for the solution of the aerosol
condensation/ evaporation equation. Aerosol Sci. and Technol., 36:407–418.

Daescu, D., Sandu, A., and Carmichael, G. (2003). Direct and adjoint sensitivity analysis of
chemical kinetics with KPP: II Numerical validation and applications. Atmos. Env.

Debry, E. (2004). Numerical simulation of an atmospheric aerosol distribution. PhD thesis,


Ecole Nationale des Ponts et Chaussées. Sous la direction de Bruno Sportisse.

Debry, E., Fahey, K., Sartelet, K., Sportisse, B., and Tombette, M. (2007). A new SIze REsolved
Aerosol Model: SIREAM. Atmos. Chem. Phys., 7(6):1537–1547.

Debry, E., Jourdain, B., and Sportisse, B. (2001). Modeling of aerosol dynamics: a stochastic
algorithm. In Sportisse, B., editor, APMS 2001, Geosciences, pages 308–319. Springer.
BIBLIOGRAPHY 65

Debry, E. and Sportisse, B. (2006). Reduction of the condensation/evaporation dynamics for


atmospheric aerosols: theoretical and numerical investigation of hybrid methods. J. Aerosol
Sci., 37(8):950–966.

Debry, E. and Sportisse, B. (2007a). Numerical simulation of the General Dynamics Equation
(GDE) for aerosols with two collocation methods. Appl. Numer. Math., 57(8):885–898.

Debry, E. and Sportisse, B. (2007b). Solving aerosol coagulation with size-binning methods.
Appl. Numer. Math., 47:1008–1020.

Debry, E., Sportisse, B., and Jourdain, B. (2003). A stochastic approach for the numerical
simulation of the General Dynamics Equation for aerosols. J. Comp. Phys., 184:649:689.

Dhaniyala, S. and Wexler, A. (1995). Numerical schemes to model condensation and evaporation
of aerosols. Atmos. Env., 30(6):919–927.

Djouad, R., Audiffren, N., and Sportisse, B. (2003a). Sensitivity analysis using automatic
differentiation applied to a multiphase chemical mechanism. Atmos. Env., 37(22):3029:3038.

Djouad, R. and Sportisse, B. (2002a). Partitioning techniques for reduction in chemical ki-
netics. APLA: an Automatic Partitioning and Lumping Algorithm. Appl. Numer. Math.,
43(4):383:398.

Djouad, R. and Sportisse, B. (2002b). Some reduction techniques for simplifying air pollution
models. In Sportisse, B., editor, APMS 2001, Geosciences, pages 235–245. Springer.

Djouad, R. and Sportisse, B. (2003). Solving reduced models in air pollution modelling. Appl.
Numer. Math., 44(1):49:61.

Djouad, R., Sportisse, B., and Audiffren, N. (2002). Numerical simulation of aqueous-phase
atmospheric models : use of a non-autonomous Rosenbrock method. Atmos. Env., 36:873–
879.

Djouad, R., Sportisse, B., and Audiffren, N. (2003b). Reduction of multiphase atmospheric
chemistry. Journal of Atmospheric Chemistry, 46:131–157.

Domilovskii, E., Lushnikov, A., and Piskunov, V. (1978). Monte Carlo simulation of coagulation
processes. Dockl. Akad. Nauk SSSR, Ser. Phys. Chem., 240(1).

Dumas, L., Coron, F., Dupuis, D., Pallegoix, J., and Sportisse, B. (1994). Applications of DSMC
methods to aerodynamic problems in an engineering context. In Méthodes de Monte-Carlo et
applications aux gaz raréfiés. Ecole CEA-INRIA-EDF, 1994.

Eibeck, A. and Wagner, W. (2000). Stochastic particle approximations for Smoluchowski’s coag-
ulation equation. Technical report, Weierstrass-Institut for Applied Analysis and Stochastics.
Preprint No. 585.

Elbern, H. and Schmidt (1999). A four dimensional variational chemistry data assimilation
scheme for Eulerian chemistry transport modeling. J. Geophys. Res., 104:18 583–18598.

Elbern, H. and Schmidt, H. (2001). Ozone episode analysis by four dimensional variational
chemistry data assimilation. J. Geophys. Res., 106:3569–3590.

Elbern, H., Schmidt, H., and Ebel, A. (1997). Variational data assimilation for tropospheric
chemistry modeling. J. Geophys. Res., 102:15967:15985.
66 BIBLIOGRAPHY

Elbern, H., Schmidt, H., Talagrand, O., and Ebel, E. (2000). 4D variational data assimilation
with an adjoint air quality model for emission analysis. Env. Mod. Soft., 15:539–548.

Fahey, K., Debry, E., Foudhil, H., and Sportisse, B. (2005). Formulation, development and pre-
liminary validation of the SIze REsolved Aerosol Model, SIREAM. In Proceedings GLOREAM
2004, pages 7–17.

Fife, P. (1979). Mathematical aspects of reacting and diffusing systems, volume 28. Springer
Verlag.

Fiore, A. and Jacob, D. (2003). Application of empirical orthogonal functions to evaluate ozone
simulations with regional and global models. J. Geophys. Res., 108(D14).

Frohn, L., Christensen, J., and Brandt, J. (2002). Development of a high-resolution nested air
pollution model. J. Comp. Phys., 179:68–94.

Galmarini, S., Bianconi, R., Klug, W., Mikkelsen, T., Addis, R., Andronopoulos, S., Astrup,
P., Baklanov, A., Bartniki, J., Bartzis, J. C., et al. (2004). Ensemble dispersion forecasting –
part I: concept, approach and indicators. Atmos. Env., 38(28):4,607–4,617.

Gelbard, F. and Seinfeld, J. (1977). Numerical solution of the dynamic equation for particulate
systems. J. Comp. Phys., 28:357–375.

Gelbard, F. and Seinfeld, J. (1979). The General Dynamics Equation for aerosols. J. Colloid
Interface Sci., 68(2):363–382.

Giering, R. and Kaminski, T. (1998). Recipes for adjoint code construction. ACM Trans. Math.
Software, 24(4):437–474.

Gong and Cho (1993). A numerical scheme for the integration of the gas-phase chemical rate
equations in 3D atmospheric models. Atmos. Env., 27A:2591–2611.

Graf, J. and Moussiopoulos, N. (1991). Intercomparison of two models for the dispersion of
chemically reacting pollutants. Beitr.Phys.Atmosph., 64(1):13–25.

Griewank, A. and Walther, A. (1999). Revolve: An implementation of checkpoint for the reverse
or adjoint mode of computational differentiation. ACM Trans. Math. Software.

Griffin, R., Dabdub, D., and Seinfeld, J. (2002). Secondary organic aerosol 1. Atmospheric
chemical mechanism for production of molecular constituents. J. Geophys. Res., 107(D17).

Guerova, G., Bey, I., Attié, J., Cui, J., and Sprenger, M. (2006). Impact of transatlatic transport
episodes on summertime ozone in Europe. Atmos. Chem. Phys., 6(2057-2072).

Hanna, S. R., Lu, Z., Frey, H. C., Wheeler, N., Vukovich, J., Arunachalam, S., Fernau, M., and
Hansen, D. A. (2001). Uncertainties in predicted ozone concentrations due to input uncer-
tainties for the UAM-V photochemical grid model applied to the July 1995 OTAG domain.
Atmos. Env., 35(5):891–903.

He, S. and al (2000). Application of ADIFOR for air pollution model sensitivity studies. Envi-
ronmental Modeling and Software, (15):475–489.

Heard, A., Pilling, M., and Tomlin, A. (1998). Mechanism reduction techniques applied to
tropospheric chemistry. Atmos. Env., 32(6):1059–1073.
BIBLIOGRAPHY 67

Hertel, O., Berkowicz, R., and Christensen, J. (1993). Test of two numerical schemes for use in
atmospheric transport-chemistry models. Atmos. Env., 27A(16):2591–2611.

Higdon, R. and Bennett, A. (1996). Stability analysis of operator splitting for large scale ocean
modeling. J. Comp. Phys., 123:311–329.

Hourdin, F. and Issartel, J.-P. (2000). Sub-surface nuclear tests monitoring through the CTBT
133 Xe network. Geophys. Res. Lett., 27:2245–2248.

Hu, Y., Odman, M., and Russel, A. (2006). Mass conservation in the CMAQ model. Atmos.
Env., 40:1199–1204.

Hudman, R. and al (2004). Ozone production in transpacific Asian pollution plumes and impli-
cations for ozone in air quality in California. J. Geophys. Res., 109(D23S10).

Hundsdorfer, W. and Verwer, J. (2003). Numerical Solution of Time-Dependent Advection-


Diffusion-Reaction Equations, volume 33 of Springer Series in Computational Mathematics.
Springer Verlag.

Isnard, O., Krysta, M., Bocquet, M., Dubiau, P., and Sportisse, B. (2005). Data assimilation
of radionuclides atmospheric dispersion at small scale: a tool to assess the consequences of
radiological emergencies. In Proceedings of the IAEA Conference. Rio Conference.

Jacob, D. J., Logan, J., and Murti, P. (1999). Effect of rising Asian emissions on surface ozone
in the United States. Geophys. Res. Lett., 26:2175–2178.

Jacobson, M. and Turco, R. (1995). Simulating condensational growth, evaporation and coagu-
lation of aerosols using a combined moving and stationary size grid. Aerosol Sci. and Technol.,
22:73:92.

Jacobson, M., Turco, R., Jensen, E., and Toon, O. (1994). Modeling coagulation among particles
of different composition and size. Atmos. Env., 28(7):1327–1338.

Jaubertie, A., Plion, P., Sportisse, B., Aumont, B., and Toupance, G. (1998). Reduction of
kinetic schemes for atmospheric chemistry. In Proceedings of the Conference APMS98. INRIA,
France.

Jay, L., Sandu, A., Potra, F., and Carmichael, G. (1997). Improved QSSA methods for atmo-
spheric chemistry integration. SIAM J. Sci. Comp.

Kalachev, L. and Field, R. (2001). Reduction of a model describing ozone oscillations in the tro-
posphere: example of an algorithmic approach to model reduction in atmospheric chemistry.
J. Atmos. Chem., 39:65–93.

Kalnay, E. (2003). Atmospheric Modeling, Data Assimilation and Predictability. Cambridge


University Press.

Karamchandani, P., Santos, L., Sykes, I., Zhang, Y., Tonne, C., and Seigneur, C. (2000). Devel-
opment and application of a state-of-the-science reactive plume model. Environ. Sci. Tech.,
34:870–880.

Kim, Y. and Seinfeld, J. (1990). Simulation of multicomponent aerosol condensation by the


moving sectional method. J. Colloid Interface Sci., 135(1):185–199.
68 BIBLIOGRAPHY

Korsakissok, I. and Sportisse, B. (2006). Simulations et évaluation de réseaux pour la détection


de traceurs à l’échelle locale. Technical Report 2006-46, CEREA. Rapport de contrat DGA.

Krishnamurti, T., Kishtawal, T., LaRow, T., Bachiocchi, D., Zhang, Z., Williford, C., Gadhil, S.,
and Surendran, S. (1999). Improved weather and seasonal climate forecasts from multi-model
superensemble. Science, 285:1548–1550.

Krysta, M., Bocquet, M., Isnard, O., Issartel, J.-P., and Sportisse, B. (2004). Data assimilation
of radionuclides at short and regional scales. In Farago, I., Georgiev, K., and Havasi, A.,
editors, Proceedings of the NATO Advanced Research Workshop on Advances in Air Pollution
Modeling for Environmental security, pages 275–285. Kluiwer.

Krysta, M., Bocquet, M., Sportisse, B., and Isnard, O. (2006). Data assimilation for short-range
dispersion of radionuclides: an application to wind tunnel data. Atmos. Env., 40(38):7267–
7279.

Lacour, S., Megueulle, C., Carlotti, P., Soulhac, L., and Sportisse, B. (2003). Unsteady effects on
pollutant dispersion around a tunnel portal. In Proceedings of Urban Air Quality Conference
UAQ4, Prague.

Lacour, S. and Sportisse, B. (2007). Estimation of indoor deposition velocity for ozone with a
simplified reactive box model. Atmos. Env. In revision.

Lagache, R., Declercq, C., Quélo, D., Sportisse, B., Palmier, P., Quetelard, B., and Haziak, F.
(2006). Evaluation du PDU de Lille-Métropole sur le trafic, les concentrations de polluants
atmosphériques et la mortalité. In Actes de la 15ième conférence sur les transports et la
pollution de l’air.

Lamb, R. (1973). Note on the application of K-Theory to diffusion problems involving nonlinear
chemical reactions. Atmos. Env., 7:257–263.

Lanser, L. and Verwer, J. (1998). Analysis of operator splitting for advection-diffusion-reaction


problems from air pollution modelling. In Proceedings 2nd Meeting on Numerical methods for
differential equations. Coimbra, Portugal.

Larrouturou, B. and Sportisse, B. (1997). Some mathematical and numerical aspects of reduction
in chemical kinetics. In Bristeau, M. and al, editors, Computational Science for the 21st
century, conference in honor of Professor R.Glowinski. John Wiley and Sons.

Larson, V., Golaz, J.-C., Jiang, H., and Cotton, W. (2005). Supplying local microphyisics
parameterizations with information about subgrid variability: Latin Hypercube sampling. J.
Atmos. Sci., 62:4010–4026.

Le Dimet, F. and Sportisse, B., editors (2002). Data Assimilation for Geophysical flows, INRIA-
CEA-EDF School on Applied Non-linear Problems. INRIA.

Lee, S., Yoon, S., and Byun, D. (2004). The effect of mass inconsistency of the meteorological
field generated by a common meteorological model on air quality modeling. Atmos. Env.,
38:2917–2926.

Li, Q., Jacob, D., Bey, I., Palmer, P., Duncan, B., Field, B., Martin, R., Fiore, A., Yantosca, R.,
Simmonds, P., and Oltmans, S. (2002). Transatlantic transport of pollution and its effects on
surface ozone in Europe and North America. J. Geophys. Res., 107.
BIBLIOGRAPHY 69

Louis, J. (1979). A parametric model of vertical eddy fluxes in the atmosphere. Boundary-Layer
Meteor., 17:197:202.

Lovas, T., Mastorakos, E., and Goussis, D. (2006). Reduction of the RACM scheme using
Computational Singular Perturbation analysis. J. Geophys. Res., 111(D13302).

Lowe, R. and Tomlin, A. (2000). Low-dimensional manifolds and reduced chemical models for
tropospheric chemistry simulations. Atmos. Env., 34:2425–2436.

Mallet, V. (2007). Pitfalls arising in the numerical integration of the aerosol condensational
growth. Technical report, CWI.

Mallet, V., Korssakissok, I., Quélo, D., and Sportisse, B. (2007a). Polyphemus : un système
modulaire de modélisation pour la dispersion atmosphérique et l’évaluation des risques. Pol-
lution Atmosphérique.

Mallet, V., Pourchet, A., Quélo, D., and Sportisse, B. (2007b). Investigation of some nu-
merical issues in a Chemistry-Transport Model: gas-phase simulations. J. Geophys. Res.,
112(D15301). doi:10.1029/2006JD008373.

Mallet, V., Quélo, D., Sportisse, B., Ahmed de Biasi, M., Debry, É., Korsakissok, I., Wu, L.,
Roustan, Y., Sartelet, K., Tombette, M., and Foudhil, H. (2007c). Technical Note: The air
quality modeling system Polyphemus. Atmos. Chem. Phys. Discuss., 7(3):6,459–6,486.

Mallet, V. and Sportisse, B. (2004). 3D Chemistry-Transport Model Polair3D: numerical


issues, validation and automatic-differentiation strategy. Atmos. Chem. Phys. Discuss.,
(1):1371:1392.

Mallet, V. and Sportisse, B. (2005). A comprehensive study of ozone sensitivity with respect to
emissions over Europe with a chemistry-transport model. J. Geophys. Res., 110(D22).

Mallet, V. and Sportisse, B. (2006a). Air quality modeling: from deterministic to stochastic
modeling. Computers and Mathematics with Application. In press.

Mallet, V. and Sportisse, B. (2006b). Ensemble-based air quality forecasts: a multi-model


approach applied to ozone. J. Geophys. Res., 111(D18):18302.

Mallet, V. and Sportisse, B. (2006c). Peut-on modéliser la qualité de l’air de manière déterministe
? In Actes 38ièmes Journées de Statistiques. Société Française de Statistique.

Mallet, V. and Sportisse, B. (2006d). Uncertainty in a chemistry-transport model due to physical


parameterizations and numerical approximations: an ensemble approach applied to ozone
modeling. J. Geophys. Res., 111(D01302).

Marchuk, G. (1986). Mathematical models in environmental problems, volume 16. North Holland.

Molemaker, M. and Vila-Guerau de Arellano, J. (1998). Control of chemical reactions by con-


vective turbulence in the boundary layer. J. Atmos. Sci., 55:568–579.

Mott, D., Oran, E., and Van Leer, B. (2000). A Quasi-Steady-State solver for the stiff Ordinary
Differential Equations of reaction kinetics. J. Comp. Phys., 164:407–428.

Müller, F. (2001). Splitting error estimation for microphysical-multiphase chemical systems in


meso-scale air quality models. Atmos. Env., pages 5749–5764.
70 BIBLIOGRAPHY

Naumann, K. (2003). COSIMA: a computer program simulating the dynamics of fractal aerosols.
J. Aerosol Sci., 34:1371–1397.

Neophytou, M., Goussis, D., Van Loon, M., and Mastorakos, E. (2004). Reduced chemical
mechanisms for atmospheric pollution using Computational Singular Perturbation analysis.
Atmos. Env., 38:3661–3673.

Nguyen, K. and Dabdub, D. (2003). Development and analysis of a non-splitting solution for
three-dimensional air quality models. Atmos. Env., 37:3741–3748.

Odman, M. and Russel, A. (2000). Mass conservative coupling of non-hydrostatic meteorological


models with air quality models. In Gryning, S. and Batchvarova, E., editors, Air pollution
modeling and its applications, pages 651–660. Kluwer.

Pagowski, M., Grell, G. A., McKeen, S. A., Dévényi, D., Wilczak, J. M., Bouchet, V., Gong, W.,
McHenry, J., Peckham, S., McQueen, J., Moffet, R., and Tang, Y. (2005). A simple method
to improve ensemble-based ozone forecasts. Geophys. Res. Lett., 32.

Pandis, S. and Seinfeld, J. (1989). Sensitivity analysis of a chemical mechanism for aqueous-
phase atmospheric chemistry. J. Geophys. Res., 94(D1):1105–1126.

Petersen, A., Spee, E., Van Dop, H., and Hundsdorfer, W. (1998). An evaluation and intercom-
parison of four new advection schemes for use in global chemistry models. J. Geophys. Res.,
103.

Pham, D. (1996). A singular evolutive interpolated Kalman filter for data assimilation in
oceanography. Technical Report RT 163, LMC/IMAG.

Pham, D., Verron, J., and Roubaud, M. C. (1998). A singular evolutive extended Kalman filter
for data assimilation in oceanography. J. Marine Systems, 16:323–340.

Pilinis, C. (1990). Derivation and numerical solution of the species mass distribution equations
for multicomponent particulate systems. Atmos. Env., 24A(7):1923–1928.

Pilinis, C., Capaldo, K., Nenes, A., and Pandis, S. (2000). MADM - a new Multi-component
Aerosol Dynamic Model. Aerosol Sci. and Technol., 32:482–502.

Preussner, P. and Brand, K. (1981). Application of a semi-implicit Euler method to mass action
kinetics. Chem. Eng. Sci., 36(10):1633–1641.

Pun, B., Griffin, R., Seigneur, C., and Seinfeld, J. (2002). Secondary Organic Aerosol 2. Ther-
modynamic model for gas/particle partitioning of molecular constituents. J. Geophys. Res.,
107(D17).

Putaud, J., Dingenen, R., Baltensperger, U., Brüggemann, E., Charron, A., Facchini, M., Dece-
sari, S., Fuzzi, S., Gehrig, R., Hansson, H., Harrison, R., Jones, A., Laj, P., Lorbeer, G.,
Maenhaut, W., Mihalopoulos, N., Müller, K., Palmgren, F., Querol, X., Rodriguez, S., Schnei-
der, J., Spindler, G., Brink, H., Tunved, P., Torseth, K., Weingartner, E., Wiedensohler, A.,
Wahlin, P., and Raes, F. (2003). A european aerosol phenomenology. Technical report, Joint
Research Centre, Institute for Environment and Sustainability.

Quélo, D., Krysta, M., Bocquet, M., Isnard, O., Minier, Y., and Sportisse, B. (2007). Validation
of the polyphemus system: the ETEX, Chernobyl and Algeciras cases. Atmos. Env., 41:5300–
5315.
BIBLIOGRAPHY 71

Quélo, D., Mallet, V., and Sportisse, B. (2005a). Inverse modeling of NO x emissions at re-
gional scale over Northern France. Preliminary investigation of the second-order sensitivity.
J. Geophys. Res., 110(D24310).

Quélo, D., Sportisse, B., Berroir, J., and Charpentier, I. (2002). Some remarks concerning inverse
modeling and data assimilation for slow-fast atmospheric chemical kinetics. In Sportisse, B.,
editor, APMS 2001, Geosciences, pages 499–513. Springer.

Quélo, D., Sportisse, B., and Isnard, O. (2005b). Data assimilation for short-range dispersion of
radionuclides: a case study for second-order sensitivity. J. Environ. Radioactivity, 84:393–408.

Quélo, D., Lagache, R., and Sportisse, B. (2004). Etude de l’impact qualité de l’air des scénarios
PDUs sur Lille à l’aide du modèle de Chimie-Transport Polair3D. Technical Report 2004-28,
CEREA. Rapport PREDIT.

Rabitz, H. and Alis, O. (1999). General foundations of high-dimensional model representations.


Journal of Mathematical Chemistry, 25:197–223.

Ramabhadran, T., Peterson, T., and Seinfeld, J. (1976). Dynamics of aerosol coagulation and
condensation. AIChe Journal, 22(5):840–851.

Rodriguez, M. and Dabdub, D. (2003). Monte Carlo uncertainty and sensitivity analysis of the
CACM chemical mechanism. J. Geophys. Res., 108(D15).

Rouchon, P. and Sportisse, B. (September 1997). Slow and fast kinetic scheme with slow diffu-
sion. In Proceedings of the Workshop Numerical Aspects of Reduction in Chemical Kinetics,
ENPC-CERMICS.

Roustan, Y. and Bocquet, M. (2006). Sensitivity analysis for mercury over Europe. J. Geophys.
Res., 111(D14304).

Roustan, Y., Bocquet, M., Musson-Genon, L., and Sportisse, B. (2005). Modeling of mercury at
the European scale with the chemistry-transport model Polair3D. In Proceedings GLOREAM
2004, pages 121–130. Also as Technical Report CEREA 2005-04.

Roustan, Y., Bocquet, M., Musson-Genon, L., and Sportisse, B. (2006). Modélisation du mer-
cure, du plomb et du cadmium à l’échelle du continent européen. Pollution Atmosphérique,
(191):317–327.

Roustan, Y., Lasry, F., and Sportisse, B. (2007). Modélisation de l’impact des émissions des cen-
trales EDF et SNET en France sur le transport atmosphérique des PM 10 et PM2.5 . Technical
Report 2007-9, CEREA.

Russell, A. and Dennis, R. (2000). NARSTO critical review of photochemical models and
modeling. Atmos. Env., 34:2,283–2,234.

Sandu, A. (2001a). Discretizing aerosol dynamics with B-splines. Technical report, Michigan
Technological University, Houghton, MI 49931.

Sandu, A. (2001b). Positive numerical integration methods for chemical kinetic systems. J.
Comp. Phys., 170:1–14.

Sandu, A. (2006). Targeted observations for atmospheric chemistry and transport models. In
International Conference for Computational Science (ICCS-2006).
72 BIBLIOGRAPHY

Sandu, A. and Borden, C. (2003). A framework for the numerical treatment of aerosol dynamics.
Appl. Numer. Math., 45:475:497.

Sandu, A., Daescu, D., and Carmichael, G. (2003). Direct and adjoint sensitivity analysis of
chemical kinetics with KPP: I Theory and software tools. Atmos. Env.

Sandu, A., Daescu, D., Carmichael, G., and Chai, T. (2005). Adjoint sensitivity analysis of
regional air quality models. J. Comp. Phys., 204:222–252.

Sandu, A., Potra, F., Carmichael, G., and Damian, V. (1996). Efficient implementation of fully
implicit methods for atmospheric chemical kinetics. J. Comp. Phys., 129:101–110.

Sandu, A., Verwer, J., Blom, J., Spee, E., and Carmichael, G. (1997a). Benchmarking stiff ODEs
solvers for atmospheric chemistry problems II: Rosenbrock solvers. Atmos. Env., 31:3459–3472.

Sandu, A., Verwer, J., Van Loon, M., Carmichael, G., Potra, F., Dabdub, D., and Seinfeld,
J. (1997b). Benchmarking stiff ODEs solvers for atmospheric chemistry problems I: implicit
versus explicit. Atmos. Env., 31:3151–3166.

Sartelet, K., Boutahar, J., Quélo, D., Coll, I., Plion, P., and Sportisse, B. (2002). Develop-
ment and validation of a 3d chemistry-transport model, Polair3D, by comparison with data
from ESQUIF campaign. In Proceedings of the 6th Gloream workshop: Global and regional
atmospheric modelling, page 140:146. Un. Aveiro, Portugal.

Sartelet, K., Debry, E., Fahey, K., Tombette, M., Roustan, Y., and Sportisse, B. (2007a).
Simulation of aerosols and gas phase species over Europe with the Polyphemus sys-
tem. Part I: model-to-data comparison for year 2001. Atmos. Env., 41:6116–6131.
doi:10.1016/j.atmosenv.2007.04.024.

Sartelet, K., Hayami, H., and Sportisse, B. (2007b). Dominant aerosol processes during high-
pollution episodes over Greater Tokyo. J. Geophys. Res. Accepted for publication.

Sartelet, K., Hayami, H., and Sportisse, B. (2007c). MICS-Asia Phase II: sensitivity to the
aerosol module. Atmos. Env. doi:10.1016/j.atmosenv.2007.03.005.

Sartelet, K. N., Hayami, H., Albriet, B., and Sportisse, B. (2005). Development and preliminary
validation of a Modal Aerosol Model for tropospheric chemistry: MAM. Aerosol Sci. and
Technol., 40(2):118–127.

Schmidt, H., Derognat, C., Vautard, R., and Beekmann, M. (2001). A comparison of simulated
and observed ozone mixing ratios for the summer of 1998 in Western Europe. Atmos. Env.,
35:6277:6297.

Segers, A. (2002). Data assimilation in atmospheric chemistry models using Kalman filtering.
PhD thesis, Delft University.

Seigneur, C., Stephanopoulos, G., and Carr Jr, R. (1982). Dynamic sensitivity analysis of
chemical reaction systems. Chem. Eng. Sci., 37(6):845–853.

Seinfeld, J. and Pandis, S. (1998). Atmospheric Chemistry and Physics. Wiley-interscience.

Shorter, J. and Precila, C. (1999). An efficient chemical kinetics solver using high dimensional
model representation. Journal of Physical Chemistry (1999), 103:7192–7198.
BIBLIOGRAPHY 73

Sirkes, Z. and Tziperman, E. (1997). Finite difference of adjoint or adjoint of finite difference?
Mon. Wea. Rev., 125:3373–3378.

Sommeijer, B., Van der Houwen, P., and Verwer, J. (1981). On the treatment of time-dependent
boundary conditions in splitting methods for parabolic differential equations. Int.J.for Num.
Met. in Eng., 17:335–346.

Sportisse, B., editor (1998). Air Pollution Modeling and Simulation APMS98, INRIA/ENPC
APMS Conference. INRIA.

Sportisse, B. (1999). Contribution à la modélisation des écoulements réactifs: réduction des


modèles de cinétique chimique et simulation de la pollution atmosphérique. PhD thesis, Ecole
Polytechnique.

Sportisse, B. (2000). An analysis of operator splitting techniques in the stiff case. J. Comp.
Phys., 161:140–168.

Sportisse, B., editor (2001a). Air Pollution Modeling and Simulation APMS01. Springer Verlag.

Sportisse, B. (2001b). Box models versus Eulerian models in air pollution modeling. Atmos.
Env., 35:173–178.

Sportisse, B. (2003a). Ozone des villes, ozone des champs: qualité de l’air et pollution pho-
tochimique. PCM Le Pont.

Sportisse, B. (2003b). Quelques aspects de la modélisation de la pollution atmosphérique.


Annales des Ponts et Chaussées, 107-108:44–55.

Sportisse, B. (2004). Assimilation de données. I Eléments théoriques. Technical Report 2004-24,


CEREA. Notes de cours ENSTA et DEA M2SAP.

Sportisse, B. (2007a). Management de la recherche. Enjeux et perspectives, chapter Partenariat


recherche publique/entreprise : l’exemple du CEREA, Laboratoire Commun ENPC/EDF
R&D. Chapitre 6. De Boeck.

Sportisse, B. (2007b). Pollutions atmosphériques. Enjeux, processus, modélisation. Springer.


350 pages.

Sportisse, B. (2007c). A review of current issues in air pollution modeling and simulation.
Computational Geosciences, 11(2):159–181.

Sportisse, B. (2007d). A review of parameterizations for modeling dry deposition and scavenging
of radionuclides. Atmos. Env., (41):2683–2698.

Sportisse, B. (September 1997). Reducing chemical kinetics: a mathematical point of view. In


Proceedings of the Workshop Numerical Aspects of Reduction in Chemical Kinetics, ENPC-
CERMICS.

Sportisse, B., Bencteux, G., and Plion, P. (2000). Method of Lines versus Operator Splitting
methods for reaction-diffusion systems in air pollution modelling. Env. Mod. Soft., 15(6-
7):673–679.

Sportisse, B., Boutahar, J., Debry, E., Quélo, D., and Sartelet, K. (2002). Some tracks in air
pollution modeling. RACSAM Journal of the Spanish Science Academy, Real Academia de
Ciencias de Espana, 96(3).
74 BIBLIOGRAPHY

Sportisse, B., Debry, E., Fahey, K., Roustan, Y., Sartelet, K., and Tombette, M. (2006a). PAM
project (Multiphase Air Pollution): description of the aerosol models SIREAM and MAM.
Technical Report 2006-08, CEREA. Available at http://www.enpc.fr/cerea/polyphemus/.

Sportisse, B., Debry, E., Fahey, K., Roustan, Y., Sartelet, K., and Tombette, M.
(2006b). Rapport du projet Primequal-Predit PAM (Pollution Atmosphérique Multiphasique).
Modélisation. Technical Report 2006-1, CEREA. 171 pages.

Sportisse, B. and Djouad, R. (2000). Reduction of chemical kinetics in air pollution modelling.
J.Comp.Phys., 164:354–376.

Sportisse, B. and Djouad, R. (2002). Some aspects of multi-timescales issues for the numerical
modeling of atmospheric chemistry, volume 130, chapter 4, page 39:61. IMA Volume in
Mathematics and its applications. Springer New York.

Sportisse, B. and Djouad, R. (2003). Mathematical investigation of mass transfer for atmospheric
pollutants into a fixed droplet with aqueous chemistry. J. Geophys. Res., 108(D2):4073.

Sportisse, B. and Djouad, R. (2007). Use of Proper Orthogonal Decompositions for the reduction
of atmospheric chemistry. J. Geophys. Res., 112(D06303).

Sportisse, B. and Du Bois, L. (2002). Numerical and theoretical investigation of a simplified


model for the parameterization of below-cloud scavenging by falling raindrops. Atmos. Env.,
36:5719–5727.

Sportisse, B., Jaubertie, A., and Plion, P. (1997). Reducing mechanisms in chemical kinetics for
the simulation of reactive transport; an application to air pollution modelling. In Proceedings
of the St Venant Symposium.

Sportisse, B. and Mallet, V. (2005). Calcul Scientifique pour l’Environnement. Cours ENSTA.
86 pages.

Sportisse, B. and Mallet, V. (2007). Data assimilation and inverse modelling of slow/fast atmo-
spheric chemistry. Journal of Computational Geosciences. Submitted.

Sportisse, B. and Quélo, D. (2002). Observational network design, adaptive observations and
targeting strategies in atmospheric data assimilation: sketch of a methodological review. In
INRIA-CEA-EDF Summer School on Data Assimilation for Geophysical flows.

Sportisse, B. and Quélo, D. (2003). Data assimilation and inverse modeling of atmospheric
chemistry. Proc. of Indian National Science Academy. Part A Physical Sciences, 69.

Sportisse, B. and Quélo, D. (2004). Assimilation de données. II Implémentation des approches


variationnelles: le cas d’un modèle de Chimie-Transport. Technical Report 2004-25, CEREA.
Notes de cours ENSTA et DEA M2SAP.

Sportisse, B., Quélo, D., and Mallet, V. (2007). Impact of mass consistency errors for atmo-
spheric dispersion. Atmos. Env., 41:6132–6142.

Sportisse, B., Sartelet, K., Debry, E., Fahey, K., Roustan, Y., and Tombette, M. (2006c). The
SIze REsolved Aerosol Model (SIREAM) and the Modal Aerosol Model (MAM). Technical
documentation. Technical Report 2006-8, CEREA. 115 pages.

Srivastava, R., Mc Rae, D., and Odman, M. (2001). Simulation of dispersion of a power plant
plume using an adaptive grid algorithm. Atmos. Env., 35:4801–4818.
BIBLIOGRAPHY 75

Stockwell, W., Kirchner, F., Kuhn, M., and Seefeld, S. (1997). A new mechanism for Regional
Atmospheric Chemistry Modeling. J. Geophys. Res., 95(D10):16343:16367.

Strang, G. (1968). On the construction and comparison of difference schemes. SIAM J. Numer.
Anal., 5:506–517.

Stratton, D., Gans, J., and Williams, E. (1994). Coagulation algorithms with size binning. J.
Comp. Phys., 112:364–369.

Sun, P. (1996). A pseudo Non-Time Splitting Method in Air Quality modeling. J. Comp. Phys.,
(127):152–157.

Taghavi, M., Musson-Genon, L., and Sportisse, B. (2004). Impact des rejets de la centrale
thermique de Martigues sur la qualité de l’air. Modélisation avec Polair3D. Technical Report
2004-27, CEREA. Rapport pour la Mission Thermique EDF.

Tarantola, A. (1987). Inverse Problem Theory. Elsevier Amsterdam.

Tatang, M., Pan, W., Prinn, R., and McRae, G. (1997). An efficient implementation for paramet-
ric uncertainty analysis of numerical geophysical models. J. Geophys. Res., 102:21,925–21,932.

Tikhonov, A., Vasileva, A., and Sveshnikov, A. (1985). Differential equations. Springer Verlag.

Tombette, M., Chazette, P., and Sportisse, B. (2007). Simulation of the Aerosol Optical Thick-
ness over Europe with a 3D size-resolved aerosol model: Comparisons with AERONET data.
Atmos. Chem. Phys. Submitted.

Tombette, M., Fahey, K., Sartelet, K., and Sportisse, B. (2005). Aerosol modelling at regional
scale: a sensitivity study with the Polyphemus platform. In Proceedings of GLOREAM 2005.
Also as Report CEREA 2005-50.

Tombette, M. and Sportisse, B. (2006). Rapport de contrat préliminaire: évaluation des impacts
du CPT Porcheville. Technical Report 2006-55, CEREA. Rapport pour la Mission Thermique
EDF.

Tombette, M. and Sportisse, B. (2007). Aerosol modeling at regional scale: Model-to-data


comparison and sensitivity analysis over Greater Paris. Atmos. Env. In press.

Tomlin, A., Berzins, M., Ware, J., Smith, J., and Pilling, M. (1997). On the use of adaptative
gridding methods for modelling chemical transport from multi-scales sources. Atmos. Env.,
31(18):2945–2959.

Troen, I. B. and Mahrt, L. (1986). A simple model of the atmospheric boundary layer; sensitivity
to surface evaporation. Boundary-Layer Meteor., 37:129–148.

Tzivion, S., Reisin, T., and Levin, Z. (1999). A numerical solution of the kinetic collection
equation using high spectral grid resolution : A proposed reference. J. Comp. Phys., 148:527–
544.

Van Loon, M. (1996). Numerical methods in smog prediction. PhD thesis, Univ. Amsterdam.

Vasileva, A., Butuzov, V., and Kalachev, L. (1994). The boundary function method for singular
perturbation problems. SIAMs monographs.
76 BIBLIOGRAPHY

Vercauteren, N. and Bocquet, M. (2006). Projet de développement d’un réseau automatisé de


télésurveillance des particules radioactives dans l’air - Etude pour l’optimisation du réseau.
Technical Report 2006-41, CEREA. Rapport de contrat IRSN.

Verlaan, M. (1998). Efficient Kalman filtering algorithms for hydrodynamic models. PhD thesis,
TU Delft.

Verwer, J. (1994). Gauss-Seidel iteration for stiff ODEs from chemical kinetics. SIAM J. Sci.
Comp., 15:1243–1250.

Verwer, J., Blom, J., and Hundsdorfer, W. (1996a). An implicit-explicit approach for atmo-
spheric transport-chemistry problems. Appl. Numer. Math., 20:191–209.

Verwer, J., Blom, J., Van Loon, M., and Spee, E. (1996b). A comparison of stiff ODEs solvers
for atmospheric chemistry problems. Atmos. Env., 30:49–58.

Verwer, J., Hundsdorfer, W., and Blom, J. (1998). Numerical time integration for air pollution
models. In Proceedings of the Conference APMS’98. ENPC-INRIA.

Verwer, J. and Simpson, D. (1995). Explicit methods for stiff ODEs from atmospheric chemistry.
Appl. Num. Math., 18:413–430.

Verwer, J., Spee, E., Blom, J., and Hundsdorfer, W. (1999). A second order Rosenbrock method
applied to photochemical dispersion problem. SIAM J. Sci. Comp., 20(4):1456–1480.

Verwer, J. and Sportisse, B. (1998). A note on operator splitting analysis in the stiff linear case.
Technical Report MAS-R9830, CWI.

Verwer, J. and Van Loon, M. (1994). An evaluation of explicit pseudo-steady state approximation
schemes for stiff ODE systems from chemical kinetics. J. Comp. Phys., 113:347–352.

Vila-Guerau de Arellano, J., Dosio, A., Vinuesa, J., Holtslag, A., and Galmarini, S. (2004).
The dispersion of chemically reactive species in the atmospheric boundary layer. Me-
teor.Atmos.Phys., 87(23-38).

Vila-Guerau de Arellano, J., Duynkerke, P., Jonker, P., and Builtjes, P. (1993). An observational
study on the effects of time and space averaging in photochemical models. Atmos. Env.,
27A(3):353–362.

Vinuesa, J.-F. and Vila-Guerau de Arellano, J. (2003). Fluxes and covariances of reacting scalars
in the convective boundary layer. Tellus, B 55(4):935–949.

Vinuesa, J.-F. and Vila-Guerau de Arellano, J. (2005). Introducing effective reaction rates to
account for the inefficient mixing of the convective boundary layer. Atmos. Env., 39:445–461.

Wang, S., Levy II, H., Li, G., and Rabitz, H. (1999). Fully equivalent operational models for
atmospheric chemical kinetics within global chemistry-transport models. J. Geophys. Res.,
104(D23):30417–30426.

Wang, Z., Navon, I., Le Dimet, F., and Zou, X. (1992). The second order adjoint analysis:
theory and applications. Meteor.Atmos.Phys., (50):3:20.

Whitby, E. and McMurry, P. (1997). Modal aerosol dynamics modeling. Aerosol Sci. and
Technol., 27:673–688.
BIBLIOGRAPHY 77

Whitby, E., Stratmann, F., and Wilck, M. (2002). Merging and remapping modes in a modal
aerosol dynamics models: a ”dynamic mode manager”. J. Aerosol Sci., 33:623–645.

Williamson, D. (2001). Time-split versus process-split coupling of parameterizations and dy-


namical core. J. Climate.

Wolke, R. and Knoth, O. (1998). An explicit-implicit numerical approach for atmospheric


chemistry-transport modeling. Atmos. Env., 32:1785–1797.

Yanenko, N. (1971). The method of fractional steps. New-York.

Yang, B. and Pope, S. (1998). An investigation of the accuracy of manifold methods and
splitting schemes in the computational implementation of combustion chemistry. Combustion
and Flame, 112:16–32.

Young, T. and Boris, J. (1977). A numerical technique for solving stiff ODE associated with the
chemical kinetics of reaction flow problems. J.Phys.Chem., 81:2424.

Zhang, Y., Bischof, C., Easter, R., and Wu, P. (1998a). Sensitivity analysis of a mixed-phase
chemical mechanism using automatic differentiation. J. Geophys. Res., 103(D15):18,953–
18,979.

Zhang, Y., Bischof, C., Easter, R., and Wu, P. (2005). Sensitivity analysis of photochemical
indicators for O3 chemistry using automatic differentiation. J. Atmos. Chem., 51:1–41.

Zhang, Y., Seigneur, C., Seinfeld, J., Jacobson, M., and Binkowski, F. (1999). Simulation of
Aerosol dynamic: A comparative Review of Algorithms used in Air Quality Models. Aerosol
Sci. and Technol., 31:487–514.

Zhang, Y., Seinfeld, J., Jacobson, M., Clegg, S., and Binkowski, F. (1998b). A comparative
review of inorganic aerosol thermodynamic equilibrium modules: similarities, differences and
their likely causes. Atmos. Env.

You might also like