KPP A Software Environment for the Simulation of Chemical Kinetics - - PowerPoint PPT Presentation

KPP – A Software Environment for the Simulation of Chemical Kinetics

Adrian Sandu Computer Science Department Virginia Tech NIST, May 25, 2004

May 25, 2004

Acknowledgements

• Dr. Valeriu Damian
• Dr. Florian Potra
• Dr. Gregory Carmichael
• Dr. Dacian Daescu
• Dr. Tianfeng Chai
• Dr. Mirela Damian

Part of this work is supported by

NSF ITR AP&IM-0205198

May 25, 2004

Related Work

CHEMKIN (http://www.reactiondesign.com/lobby/open) Translates symbolic chemical system in a data file that is then used by internal libraries for

• simulation. Gas-phase kinetics, surface kinetics, reversible equations, transport, mixing,

deposition for different types of reactors, direct sensitivity analysis (Senkin). Database of reaction data, graphical postprocessor for results. KINALC (http://www.chem.leeds.ac.uk/Combustion/kinalc.htm) Postprocessor to CHEMKIN for sensitivity and uncertainty analysis, parameter estimation, and mechanism reduction; etc. KINTECUS (http://www.kintecus.com) Compiler/ chemical modeling software. Can run heterogeneous and equilibrium chemistry, generates analytical Jacobians, fit/optimize rate constants, initial concentrations, etc. from data; sensitivity analysis; Excel interface. Can use Chemkin models and databases. CANTERA (http://rayleigh.cds.caltech.edu/~goodwin/cantera) Object-oriented package for chemically-reacting flows. C++ kernel, STL, standard numerical libraries, Interfaces for MATLAB, Python, C++, and Fortran. Capabilities: homogeneous and heterogeneous kinetics, equilibria, reactor modeling, multicomponent transport. LARKIN/LIMKIN (http://www.zib.de/nowak/codes/limkin_1.0/full) Simulation of LARge systems of chemical reaktion KINetics and parameter identification. Parser generates Fortran code for function and Jacobian, or internal data arrays. DYNAFIT (http://www.biokin.com/dynafit) Performs nonlinear least-squares regression of chemical kinetic, enzyme kinetic, or ligand- receptor binding data using experimental data. Parses symbolic equations.

May 25, 2004

KPP in a Nutshell

• The Kinetic PreProcessor
• Purpose: automatically implement building blocks for large-scale

simulations

• Parses chemical mechanisms
• Generates Fortran and C code for simulation, and direct and adjoint

sensitivity analysis

• Function, Jacobian, Hessian, Stoichiometric matrix, derivatives of

function&Jacobian w.r.t. rate coefficients

• Treatment of sparsity
• Comprehensive library of numerical integrators
• Used in several countries by academia/research/industry
• Free! http://www.cs.vt.edu/~asandu/Software/KPP

May 25, 2004

KPP Architecture

Fortran C Fortran C Fortran C Fortran C Scanner / Parser Error checking Compute best sparsity structure Reorder species for best sparsity Compute expression trees for assignments Substitution pre-processor prod/destr function jacobian code sparsity structure integrator function driver function(main) utility functions I/O functions Fortran C Driver file Integrator file Utility file I/O file kinetic description files atoms species equations

• ptions

inline code

KPP

Fortran C Code generation

May 25, 2004

KPP Example

SUBROUTINE FunVar ( V, R, F, RCT, A_VAR ) INCLUDE 'small.h' REAL*8 V(NVAR), R(NRAD), F(NFIX) REAL*8 RCT(NREACT), A_VAR(NVAR) C A - rate for each equation REAL*8 A(NREACT) C Computation of equation rates A(1) = RCT(1)*F(2) A(2) = RCT(2)*V(2)*F(2) A(3) = RCT(3)*V(3) A(4) = RCT(4)*V(2)*V(3) A(5) = RCT(5)*V(3) A(6) = RCT(6)*V(1)*F(1) A(7) = RCT(7)*V(1)*V(3) A(8) = RCT(8)*V(3)*V(4) A(9) = RCT(9)*V(2)*V(5) A(10) = RCT(10)*V(5) C Aggregate function A_VAR(1) = A(5)-A(6)-A(7) A_VAR(2) = 2*A(1)-A(2)+A(3)-A(4)+A(6)- &A(9)+A(10) A_VAR(3) = A(2)-A(3)-A(4)-A(5)-A(7)-A(8) A_VAR(4) = -A(8)+A(9)+A(10) A_VAR(5) = A(8)-A(9)-A(10) RETURN END #INCLUDE atoms #DEFVAR O = O; O1D = O; O3 = O + O + O; NO = N + O; NO2 = N + O + O; #DEFFIX O2 = O + O; M = ignore; #EQUATIONS { Small Stratospheric } O2 + hv = 2O : 2.6E-10*SUN**3; O + O2 = O3 : 8.0E-17; O3 + hv = O + O2 : 6.1E-04*SUN; O + O3 = 2O2 : 1.5E-15; O3 + hv = O1D + O2 : 1.0E-03*SUN**2; O1D + M = O + M : 7.1E-11; O1D + O3 = 2O2 : 1.2E-10; NO + O3 = NO2 + O2 : 6.0E-15; NO2 + O = NO + O2 : 1.0E-11; NO2 + hv = NO + O : 1.2E-02*SUN;

May 25, 2004

Sparse Jacobians

1 10 20 30 40 50 60 70 74 1 10 20 30 40 50 60 70 74

#JACOBIAN [ ON | OFF | SPARSE ]

IDEAS:

• Chem. interactions: sparsity pattern (off-line)
• Min. fill-in reordering
• Expand sparsity structure
• Row compressed form
• Doolitle LU factorization
• Loop-free substitution

JacVar(…), JacVar_SP(…), JacVar_SP_Vec(…), JacVarTR_SP_Vec(…) KppDecomp(…) KppSolve(…), KppSolveTR(…) E.g. SAPRC-99 74+5 spc./211 react. NZ=839, NZLU=920

May 25, 2004

Computational Efficiency

8 9 10 20 30 50 80 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 CPU time [seconds] Number of Accurate Digits qssa chemeq vode sdirk4 rodas S−rodas S−sdirk4 S−vode Lsodes

0.12 0.06 0.23 KPP 0.35 0.21 0.61 Harwell Dec+7Sol Sol Dec (1/Lapack) Linear Algebra Stiff I ntegrators

May 25, 2004

Sparse Hessians

#HESSIAN [ ON | OFF ]

k j i k j i

y y f H ∂ ∂ ∂ =

2 , ,

• 3-tensor
• sparse coordinate format
• account for symmetry

HessVar(…) HessVar_Vec(…) HessVarTR_Vec(…) E.g. SAPRC99. NZ = 848x2 (0.2%)

May 25, 2004

Stoichiometric Form

1 20 40 60 80 100 120 140 160 180 200 211 1 20 40 60 79

Reaction SAPRC−99. Stoichiometric Matrix. NZ = 1024 ( 6 % ) Species

STOICM (column compressed) ReactantProd(…) JacVarReactantProd(…) dFunVar_dRcoeff(…) dJacVar_dRcoeff(…)

#STOICMAT [ ON | OFF ]

May 25, 2004

Requirements for Numerics

Numerical stability (stiff chemistry) Accuracy: medium-low (relerr~10-6-10-2) Low Computational Time Mass Balance Positivity

May 25, 2004

Stiff Integration Methods

( )

∑ ∑

= − = − = k i i n n i n k i i n n i

y f h y

] [ ] [

β α

BDF

( ) ( )

j s j j i n i s j j j n n

Y f a y Y Y f b y y

∑ ∑

= = +

+ = + =

1 , 1 1

,

Implicit Runge-Kutta

May 25, 2004

Rosenbrock Methods

( )

( )

j i j j i h i i n h i j j j i n i s j j j n n

k c Y f k J I k a y Y k m y y

∑ ∑ ∑

− = − = = +

+ = ⋅ − + = + =

1 1 , 1 1 1 1 , 1 1 γ

No Newton Iterations Suitable for Stiff Systems Mass Conservative Efficient: Low/Med Accuracy

0.6 0.8 9 1 2 3 4 5 6 7 8 10 0.5 1 1.5 2 2.5 3 3.5 4 SEULEX VODE Rodas3 Ros3 RODAS TWOSTEP CPU time [seconds] SDA

May 25, 2004

Direct Decoupled Sensitivity

( ) ( ) ( )

⎩ ⎨ ⎧ ≤ ≤ + ⋅ = ′ ∂ ∂ = = ′ m p y t f S p y t J S p y S p y t f y

p

l

l

l l l l

1 , , , , , ,

( )

[ ]

( )

[ ]

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⋅ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + − + − = ℑ − ⇒ = −

− −

U U L P U J JS h L P U J JS h L P h I LU P J h I

T p y m T p y T T

m

L M O M L L M O M L L

1 1 1

1

γ γ γ γ

Problem

May 25, 2004

DDM with KPP

[ ]

∑ ∑

− = + − + + − = + − +

≤ ≤ + = ⋅ − + =

1 1 1 1 1 1 1

1 ,

k i n p i n n n k i n i n n

m f h S S J h I f h y y l

l

l l

β β β

[ ]

( )

[ ]

( ) ( ) ( )

∑ ∑ ∑ ∑ ∑ ∑

− = − = − = = − = + = +

+ + ⋅ × + ⋅ + + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⋅ = ⋅ − + + = ⋅ − + = + = + =

1 1 , 1 1 1 1 1 1 1 1 1 1 1 1 1 1

, , , , , , , ,

i j n t p i n n t i i n n i n p j ij h i i p i j j ij n i i i n h i j n t i j ij h i i i n h s i i j j ij n i i i n n s i i i n n

f h S J h k S H k J k c p Y T f k a S p Y T J k J I f h k c p Y T f k J I k a y Y k m S S k m y y

l l l

l l l l l l l l l

γ γ γ

γ γ

BDF

(Dunker, ’84)

Rosenbrock

(Sandu et. al. ‘02)

May 25, 2004

Adjoint Sensitivity Analysis

( ) ( )

( )

) ( ) ( , , , '

f y f

t y t y t t t y t f y ψ ψ ∇ ⇒ ≤ ≤ =

( )

( ) ( )( ) ( )

) ( , , , ' t t y t t t t y t J

y f y f f T

λ ψ ψ λ λ λ = ∇ ⇒ ∇ = ≥ ≥ ⋅ − =

Problem: Stiff ODE, scalar functional. Continuous: Take adjoint of problem, then discretize. Discrete: Discretize the problem, then take adjoint: Note: In both approaches the forward solution needs to be precomputed and stored.

( ) ( ) ( ) ( )( ) ( )

1

, 1 , , 1 , , λ ψ ψ λ λ = ∇ ⇒ ∇ = ∇ = − = =

+ N y N y N k T k y k k k k

y y y F N k y F y L

May 25, 2004

Linear Multistep Methods

( )

i m k i i m i i m m T k i i m i m i

h y J

+ = + + = + +

∑ ∑

⋅ = λ β λ α

] [ ] [

( )

∑ ∑

= − = − = k i i n n i n k i i n n i

y f h y

] [ ] [

β α

Method Continuous Adjoint Discrete Adjoint

( )

i m k i i m T m i m k i i m m i

y J h

+ = + = +

⋅ =

∑ ∑

λ β λ α

] [ ] [

Consistency: ~one-leg method, in general not consistent with continuous adj. eqn. I mplementation: with KPP

May 25, 2004

Runge Kutta Methods

( )

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⋅ = + =

∑ ∑

= + = + s j j i j n i i T i s i i n n

a b Y hJ

1 , 1 1 1

, θ λ θ θ λ λ

( ) ( )

∑ ∑

= = +

+ = + =

s i j j i n i s i i i n n

Y f a h y Y Y f b h y y

1 , 1 1

,

Method Continuous Adjoint Discrete Adjoint (Hager, 2000)

( ) ( )

j i n T s j j i n i i i n T s i i n n

h c t y J a h h c t y J b h Λ ⋅ − + = Λ Λ ⋅ − + =

+ = + + = +

∑ ∑

) ( , ) (

1 1 , 1 1 1 1

λ λ λ

Consistency (Sandu, 2003) Consider a Runge-Kutta method of order p. Its discrete adjoint is an order p numerical discretization of the continuous adjoint equation. (Proof: use elementary differentials of transfer fcns). I mplementation: with KPP

May 25, 2004

Singular Perturbation

(Sandu, 2003) I f RK with invertible coefficient matrix A and R(∞) = 0; and the cost function depends only on the non-stiff variable y Then µ= 0 and λ are solved with the same accuracy as the original method, within O(ε). A similar conclusion holds for continuous RK adjoints.

f

t t t z y g z z y f y ≤ ≤ = = ), , ( ' ), , ( ' ε

( ) ( )

) ( ), ( ) ( , ) ( ), ( ) (

) ( ) ( f f t z f f t y

t z t y t t z t y t ψ µ ψ λ ∇ = ∇ =

Test problem relevant for stiff systems: Distinguish between derivatives w.r.t. stiff/non-stiff variables

May 25, 2004

Rosenbrock Methods

( )

[ ]

( )

( )

( )

∑ ∑ ∑

= = + + = +

+ ⋅ × + = ≥ ≥ ⋅ = + + = ⋅ −

s i i i s i T i n n n i i T i s i j j i j h j i j n i i T n h

v u k H i s u Y J v u c v a m u J I

1 1 1 1 , 1 , 1 1

1 λ λ λ

γ

( )

[ ]

( )

s i z c Y J z J I z m h t y Y z a

j j i i j h i i T i T n h j s j j n n i j i n i j j i n i

≤ ≤ + Λ ⋅ = ⋅ − + = − = + = Λ

∑ ∑ ∑

− = + = + − = + +

1 , , ) ( ,

, 1 1 1 1 1 1 1 1 1 1 , 1 γ

λ λ α λ

Discrete Adjoint Continuous Adjoint Consistency (Sandu, 2003) Similar to Runge-Kutta I mplementation: with KPP

May 25, 2004

Computational Efficiency

Tdiscrete adjoint≤ 5 Tforward (Griewank, 2000)

KPP/Rodas-3 on SAPRC-99:

4.4 2.3 3.3 1.2 Tdiscr-grad /Tfwd Tdiscr-adj /Tfwd Tcont-grad /Tfwd Tcont-adj /Tfwd

May 25, 2004

KPP Numerical Library

Simulation

Sparse: BDF (VODE, LSODES), Runge-Kutta (Radau5), Rosenbrock (Ros-{1,2,3,4},Rodas- {3,5}. QSSA. Drivers

Direct Decoupled Sensitivity

Sparse: ODESSA, Rosenbrock, I.C. and R.C.

Adjoint Sensitivity

Continuous Adj. with any simulation method Discrete Adj. Rosenbrock. Drivers.

May 25, 2004

Swedish Meteo & Hydro Inst.

joakim.langner@smhi.se, www.smhi.se

May 25, 2004

Max Planck Institute

May 25, 2004

MISTRA-MPIC, U. Heidelberg

“Chemical reactions in the gas phase are considered in all model layers, aerosol chemistry only in layers where the relative humidity is greater than the crystallization humidity. […] The set of chemical reactions is solved using KPP.”

(http://www.iup.uni-heidelberg.de/institut/forschung/groups/atmosphere/modell/glasow)

May 25, 2004

User Contributions to KPP

Artenum SARL ParisCyberVillage 101-103 Bld Mac Donald 75019 Paris 19ème France www.artenum.com (Consulting company)

May 25, 2004

Example: Modeling Air Pollution

May 25, 2004

The Forward Model

( ) ( )

( )

( ) ( ) ( )

GROUND i i DEP i i OUT i IN IN i i F i i i i i i i

• n

Q C V n C K

• n

n C K

• n

x t C x t C t t t x C x t C E C f C K C u dt dC Γ − = ∂ ∂ Γ = ∂ ∂ Γ = ≤ ≤ = + + ∇ ⋅ ∇ + ∇ ⋅ − = , , , , 1 1 ρ ρ ρ ρ r

Mass Balance Equations. C = mole fraction. ρ = air density.

May 25, 2004

Discrete Forward Model

t HOR t VERT t CHEM t VERT t HOR t t t k t t t k

T T R T T N C N C

∆ ∆ ∆ ∆ ∆ ∆ + ∆ + +

= =

• ]

, [ ] , [ 1

Operator Splitting:

Conservative Methods for Transport Stiff Methods for Chemistry (KPP), Specific Methods for Aerosols, Different Time Steps.

May 25, 2004

Chemical Data Assimilation

• e.g. clouds, winds.
• fthe atmosphere

Physical aspects and optical properties chemical components destruction of the evolution, including (CTMs) provide the production and provide

Chemical Models Meteorological Models

Atmospheric Models

New Approaches

Assimilation Techniques:

Improved repre− transport features Chemical & Optical Observations: Sattelite, ground, aircraft Column (integral) quantities Algorithms Better Retrieval Meteorological Observations Winds, Temperature, Clouds sentation of cloud fluxes

Emission Estimates

Adjoint Modeling Variational (4D−Var) Fluxes Budgets Concentrations

Model Outputs:

Optimal Analysis State

analysis fields assimilated forecasts; Meteorological Improved for Improved Inverse Modeling Emission Estimates measurement strategies Model Performance; Feedbacks Applications Better design

• f observation networks/

weather forecasts and Chemical Improved Air Quality Evaluation;

May 25, 2004

4D-Var Data Assimilation

( )

( ) ( ) ( ) ( )

b T b k k k k T N k k k k

p p B p p

• y

H R

• y

H p − − + − − =

− − =

∑

1 2 1 1 1 2 1

min ψ

y0 yN

FORWARD ADJOINT

λN y1 y2 y3 λ1 λ2 λ3

(Note: Need the gradient of J)

λ ψ = ∇ p

May 25, 2004

Continuous Adjoint Model

( )

( ) ( )

( ) ( ) ( ) ( )

GROUND i DEP i i OUT i i IN i

• F

F i F i i i T i i i

• n

V n K

• n

n K u

• n

x t t t t x x t C F K u dt d Γ = ∂ ∂ Γ = ∂ ∂ + Γ = ≥ ≥ = − ⋅ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ ⋅ ⋅ ∇ − ⋅ ∇ − = λ ρ λ ρ ρ λ ρ λ λ λ λ φ λ ρ ρ λ ρ λ λ , , , ) ( r r

Note: Linearized chemistry generated by KPP.

May 25, 2004

Discrete Adjoint Model

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

'* '* '* '* '* ] , [ '* 1 ] , [ '* ' ' ' ' ' ] , [ ] , [ 1 t HOR t VERT t CHEM t VERT t HOR t t t k t t t k t HOR t VERT t CHEM t VERT t HOR t t t k t t t k

T T R T T N N T T R T T N C N C

∆ ∆ ∆ ∆ ∆ ∆ + + ∆ + ∆ ∆ ∆ ∆ ∆ ∆ + ∆ + +

= = = ′ ′ =

• λ

λ δ δ

Operator Split Tangent Linear and Adjoint Discrete Models Note: Chemical Model Discrete/Continuous Adjoints automatically generated by KPP

May 25, 2004

Adjoint STEM-III

Measurement info used to

adjust initial fields and improve predictions;

East Asia. Grid: 80 x 80Km.

Time: [0,6] GMT, 03/01/01.

SAPRC 99 (Ros-2); 3rd order

upwind FD; LBFGS;

Parallelization with PAQMSG Distributed Level-2

Checkpointing Scheme

I/O Data Local Chkpt Local Chkpt Local Chkpt Local Chkpt Node Node Node Node Master

5 10 15 20 5 10 15 20

• No. Workers

Relative Speedup Ideal Speedup Relative Speedup

May 25, 2004

NASA Trace-P Experiment

Transport and Chemical Evolution over the Pacific March-April 2001 Quantify Asian transport Understand chemistry over West Pacific

May 25, 2004

Computational domain:

• Area: 7.200 x 4.800 km
• Horizontal grid: 80 x 80 Km
• Vertical grid: 18 layers,10 km.

Computational Setting

Model Size ~ 8,000,000 variables

May 25, 2004

Trace-P Simulation

March 01-04, 2001 O3 NO2 SO2

May 25, 2004

Model vs. Observations

10 20 30 40 50 60 70 80 90 24 25 26 27 28 29 30 31 32 33 2000 4000 6000 8000 10000 12000 O3 (ppbv) Altitude (m) TIME (GMT) Observed O3 Modeling O3 Flight Height (m)

Modeled O3 vs. Measured O3

• Cost functional =

model-observation gap.

• The analysis produces

an optimal state of the atmosphere using: Model information consistent with physics/chemistry Measurement information consistent with reality

May 25, 2004

More Observations

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 24 25 26 27 28 29 30 31 32 33 2000 4000 6000 8000 10000 12000 NO (ppbv) Altitude (m) TIME (GMT) DC-8 Flight #6 NO on 3/3/2001 Observed NO Modeling NO Flight Height (m) 0.5 1 1.5 2 2.5 3 3.5 4 24 25 26 27 28 29 30 31 32 33 2000 4000 6000 8000 10000 12000 SO2 (ppbv) Altitude (m) TIME (GMT) DC-8 Flight #6 SO2 on 3/3/2001 Observed SO2 Modeling SO2 Flight Height (m)

NO: Observation vs. Model SO2: Observation vs. Model

May 25, 2004

Target: O3 at Cheju Island

March 4-6, 2001: Strong NE flow, Beijing-Cheju March 22-25, 2001: Weaker Flow

May 25, 2004

Cones of Influence

March 4-6, 2001 March 22-25, 2001 O3 NO2 HCHO Isosurfaces of time integrals of adjoint vars. (ψ = O3 at Cheju).

May 25, 2004

Areas of Influence

Ψ = Ozone at Cheju, 0GMT, 03/04/2001

Influence areas (adjoint isosurfaces) depend on meteo, but cannot be determined solely by them (nonlinear chemistry). Boundary condition uncertainties 3 days before.

d ψ / d O3 d ψ / d HCHO d ψ / d NO2

May 25, 2004

Flights on March 07, 2001

DC-8 P3-B

Taiwan (120E-122E, 22N-25N)

May 25, 2004

Assimilation of DC-8 O3

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x time (GMT hr)

• bservation (ppbv), model result (ppbv)

2 4 6 8 10 20 30 40 50 60 70 80 90

• bservation (ppbv)

model result (ppbv) (before DA) model result (ppbv) (After DA)

+ x

Control = Initial O3

• Assim. window [0,12] GMT

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x time (GMT hr)

• bservation (ppbv), model result (ppbv)

2 4 6 8 10 20 40 60 80 100 120

• bservation (ppbv)

model result (ppbv) (before DA) model result (ppbv) (After DA)

+ x

O3 along P3-B flight O3 along DC-8 flight

May 25, 2004

Assimilation of P3-B O3

Control = Initial O3

• Assim. Window [0,12] GMT

O3 along P3-B flight O3 along DC-8 flight

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x time (GMT hr)

• bservation (ppbv), model result (ppbv)

2 4 6 8 10 20 40 60 80 100 120

• bservation (ppbv)

model result (ppbv) (before DA) model result (ppbv) (After DA)

+ x + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x time (GMT hr)

• bservation (ppbv), model result (ppbv)

2 4 6 8 10 20 40 60 80 100

• bservation (ppbv)

model result (ppbv) (before DA) model result (ppbv) (After DA)

+ x

May 25, 2004

Assimilation of Multiple Species

P3-B Obs: O3(8%); NO, NO2(20%); HNO3, PAN, RNO3(100%)

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) O3

4 6 8 10 20 40 60 80 100 120

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + ++ + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) RNO3

4 6 8 10 0.02 0.04 0.06 0.08 0.1 0.12

+ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) PAN

4 6 8 10 0.2 0.4 0.6 0.8 1 1.2

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) HNO3

4 6 8 10 0.5 1 1.5 2

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) NO2

4 6 8 10 0.1 0.2 0.3 0.4

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) NO

4 6 8 10 0.02 0.04 0.06 0.08 0.1 0.12 0.14

O3 NO2 PAN NO HNO3 RNO3

May 25, 2004

Sensitivity Analysis

Vertical level

5 10 15 20 25 30 35 40 45 50 55 60 65 5 10 15

1E-09 1E-08 1E-07 1E-06 1E-05 0.0001 0.001 0.01

Species Index

1. O3 2. H2O2 3. NO 4. NO2 5. NO3 6. N2O5 7. HONO 8. HNO3 9. HNO4

• 10. SO2

11. H2SO4 12. CO 13. HCHO 14. CCHO 15. RCHO 16. ACET 17. MEK 18. HCOOH 19. MEOH 20. ETOH 51. PHEN 52. PAN 53. PAN2 54. PBZN 55. MA_PAN 56. BC 57. OC 58. SSF 59. SSC 60. OPM10 41. ALK3 42. ALK4 43. ALK5 44. ARO1 45. ARO2 46. OLE1 47. OLE2 48. TERP 49. RNO3 50. NPHE 31. PROD2 32. DCB1 33. DCB2 34. DCB3 35. ETHENE 36. ISOPRENE 37. C2H6 38. C3H8 39. C2H2 40. C3H6 21. CCO_OH 22. RCO_OH 23. GLY 24. MGLY 25. BACL 26. CRES 27. BALD 28. ISOPROD 29. METHACRO 30. MVK 61. OPM25 62. DUST1 63. DUST2 64. DUST3 65. DMS 66. CO2 Sensitivity magnitude

ψ = NOY (P3-B) Averaged gradients help with choice of control variables

May 25, 2004

Assimilation of DC-8 CO

Control: Initial Conc. of 50 spc.

+ + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

Time (GMT hr) CO concetrations (ppbv)

2 4 6 8 10 12 50 100 150 200 250 300 350 400

Observations Model (before DA) Model (after 10 iterations)

+

X

+ + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X XX X X X X X X X X X X X X X X XX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

Time (GMT hr) CO Observation (ppbv), CO model prediction (ppbv)

2 4 6 8 10 50 100 150 200 250 300 350 400

CO Observation (ppbv) CO model prediction (ppbv) (before DA) CO model prediction (ppbv) (after 4 iterations)

+

X

Control: Initial CO Conc.

May 25, 2004

Assimilation of DC-8 CO

Control: Initial Conc. 50 spc. CO along P3-B flight CO along DC-8 flight

+ + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

Time (GMT hr) CO concetrations (ppbv)

2 4 6 8 10 12 50 100 150 200 250 300 350 400

Observations Model (before DA) Model (after 10 iterations)

+

X

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

Time (GMT hr) CO concetrations (ppbv)

2 4 6 8 10 50 100 150 200 250 300 350 400 450

p3 Observations Model (before DA) Model (after 10 iterations)

+

X

May 25, 2004

Different Control Sets

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

TIME (GMT hr) NOy (ppbv)

2 4 6 8 10 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Observation Model (Without 4D-Var) Model (Ctrls: NO + NO2 + HNO3 + PAN +PAN2) Model (50 Ctrls)

+

X

NOY= NO+NO2+NO3+2*N2O5+HONO+HNO3+HNO4+RNO3+PAN+PAN2+PBZN+MA_PAN

Observed: NOY (P3-B)

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

TIME (GMT hr) NOy (ppbv)

2 4 6 8 10 0.5 1 1.5 2 2.5 3 3.5

Observation Model (Without 4D-Var) Model (Ctrls: NO + NO2 +PAN) Model (Ctrls: NO + NO2 + O3)

+

X

(NO,NO2,PAN) > (NO,NO2,O3) 50 spc > (NO,NO2,HNO2,PAN,PAN2)

May 25, 2004

Effect on Unobserved Species

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) NO2

4 6 8 10 0.1 0.2 0.3 0.4

+ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Time GMT (hr) PAN

4 6 8 10 0.2 0.4 0.6 0.8 1 1.2

Observed: NOY (P3-B) Control: Initial NO, NO2, HNO3, PAN, PAN2 NO2 (unobserved) PAN (unobserved)

May 25, 2004

Conclusions

KPP software tool for the simulation of

chemical kinetics

Code generation Useful (and widely used!) to build

blocks for large-scale simulations

Examples of chemical data assimilation

which allows enhanced chemical weather forecasts

May 25, 2004

Quote of the Day

“Persons pretending to forecast the future shall be considered disorderly under subdivision 3, section 901 of the criminal code and liable to a fine of $250 and/or 6 months in prison.”

Section 889, New York State Code of Criminal Procedure (after M.D. Webster)

May 25, 2004

Thank you!