First- and Second-Order Sensitivity Analysis of a Photochemically

Specifically, the variation of photolysis rate coefficients for NO2 (a) (reaction 1) ..... The range of concentration variations is bounded by the nat...
0 downloads 0 Views 374KB Size
Download PDF
Environ. Sci. Technol. 1997, 31, 1206-1217

First- and Second-Order Sensitivity Analysis of a Photochemically Reactive System (a Green’s Function Approach) †

LAURENT VUILLEUMIER, ROBERT A. HARLEY,‡ AND N A N C Y J . B R O W N * ,† Energy and Environment Division, Lawrence Berkeley National Laboratory, and Civil and Environmental Engineering Department, University of California, Berkeley, California 94720

Local sensitivity analysis of air quality models permits an estimation of the effects of parameter variations, can help assess uncertainties, and distinguishes important from unimportant model features. While local sensitivity analysis in atmospheric chemistry is usually conducted to the firstorder, this work includes second-order sensitivity analysis. A Green’s function method is used to compute the sensitivities of a photochemically reactive system. Timedependent photolytic reaction rate coefficients are specified using two-parameter functions of time. The Green’s function is used directly to study the temporal dependence of ozone concentrations on the NOx concentrations. ROGlimited (ROG/NOx ) 6) and NOx-limited (ROG/NOx ) 24) conditions are both evaluated. The importance of NOx, formaldehyde, and aromatic chemistry is emphasized by the high ozone sensitivity to the initial concentrations of these species and to the reaction rate parameters of NO2 photolysis, NO recombination with ozone (NO + O3 f NO2 + O2), formaldehyde photolysis (HCHO + hν f 2HO2 + CO), and NO2 + OH f HNO3. The parameters that ozone concentration is most sensitive to were treated with second-order sensitivity analysis. A comparison of ozone concentration changes predicted by second-order sensitivity analysis with concentration changes computed directly shows that second-order sensitivity predictions are valid for parameter variations as large as (50% of the parameter’s nominal value, whereas first-order predictions are valid for parameter variations of at most (25%. Even though all of the NOx in these simulations is introduced via initial concentrations, the Green’s function reveals that the temporal pattern of NOx emission is likely to influence ozone concentrations.

1. Introduction The difficulty of implementing successful regulations to reduce ozone concentrations in urban areas throughout the United States highlights the importance of understanding the mechanism of ozone formation. Sensitivity analysis is a powerful tool for improving our understanding of the relevant chemistry by focusing our attention on the most significant reactions and guiding future research to reduce critical uncertainties. † ‡

Lawrence Berkeley National Laboratory. University of California.

1206

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

Sensitivity analysis quantifies the effect that single parameter (or independent variables, Rj) variations have on model outputs (or dependent variables, y). This can be achieved in several ways. First, the model can be re-run with different parameter values to compute sensitivity coefficients with a finite difference approximation:

yRj )

∂y y(Rj + ∆Rj) - y(Rj) ≈ ∂Rj ∆Rj

Second, the sensitivity coefficients can be computed in a more analytical way by solving differential equations for the sensitivity coefficients as functions of time. The sensitivity equations are solved using the direct method either with the sensitivity equations solved together with the model equations [coupled method (1)] or separately [decoupled method (2)] or with the Green’s function method (3). Solving the sensitivity equations minimizes the computational burden when a large number of parameters are involved because the sensitivities for all dependent-independent variable combinations are computed in a single run. Some other techniques used in sensitivity analysis do not compute local sensitivity coefficients but rather use global sensitivity analysis (for a review of sensitivity analysis methods see ref 4). Sensitivity analysis can be used to help assess model uncertainties, to distinguish important from unimportant features of a model, and to identify input parameters that most significantly influence model predictions. Sensitivity analysis techniques have been used to study combustion chemistry (5-7) as well as atmospheric chemistry (1, 8-15). In this work, the Green’s function method is applied to perform sensitivity analysis on a photochemically reactive system. This approach allows the computation of first- and second-order sensitivities to both initial conditions and chemical reaction rate coefficients. While first-order sensitivities are used to estimate the effect of parameter variations on model outputs, second-order sensitivities indicate the effect of parameter variations on first-order sensitivity coefficients. Their analysis enables the study of the interactions, within the frame of the model, among all kind of parameters used for first-order sensitivity analysis (e.g., initial concentrations or reaction rate parameters). Second-order sensitivity coefficients are also needed to improve the accuracy of Taylor expansion predictions when local sensitivities are used to predict the effect of large parameter variations. Time-varying photolytic reaction rate coefficients are specified using twoparameter functions of time, and sensitivity analysis is performed for the parameters specifying the photolytic reaction rates. The limits of applicability of first- and secondorder sensitivity analysis are also determined for this model. Green’s function values are used to explore the dependence of one chemical species concentration at any time to another species concentration at an earlier time. Another goal of this study is to determine whether findings from previous sensitivity analysis of combustion chemistry (5) concerning second-order sensitivity also apply to atmospheric chemistry. Specifically, Yetter et al. determined that large second-order sensitivity coefficients are more likely to occur for parameters to which the model outputs already show large first-order sensitivity and that large second-order diagonal coefficients reflect the bounded nature of the system.

2. AIM Method The Analytically Integrated Magnus method (AIM) was developed by Kramer et al. (16-19) to compute sensitivity coefficients in combustion modeling using a Green’s function approach. This method allows the computation of first- and

S0013-936X(96)00727-4 CCC: $14.00

 1997 American Chemical Society

second-order sensitivities to parameters like rate coefficients and initial conditions. This method is also designed to minimize the computational requirements for problems where the number of parameters exceeds the number of ordinary differential equations (ODE) involved in the model description and, as such, is well suited for our purpose. AIM solves a system of ODE describing the sensitivity problem. The rate equation for the time evolution of species concentrations, y, is given by

y3 (t) ) f(y, r, t)

and

y(t0) ) y0

(1)

where y is a vector of n species concentrations of interest, r is a vector of m parameters, and t is time. The corresponding first-order sensitivity equation is

y3 Rj(t) ) J(t)yRj(t) + fRj(t)

and

yRj(t0) ) δRjy0i

(2)

where yRj ) ∂y/∂Rj, J(t) (an n × n matrix) is the Jacobian, Jij ) ∂fi/∂yj, fRj ) ∂f/∂Rj, and δRjy0i is an n vector with all components being 0, except the ith component being 1 if Rj ) y0i. The solution of such a problem is given by

yRj(t) ) K(t,t0)yRj(t0) +

∫ K(t,τ)f t

Rj(τ)

t0

dτ

(3)

K(t,τ) is an n × n matrix called the kernel or Green’s function of the problem. The AIM method is based upon the following approach to solve for K(t,τ):

K(t,τ) ) exp(Ω(t,τ)) Ω(t,τ) )

∫ J(t ) dt + 21∫ [J(t ), ∫ J(t ) dt ] dt + 1 [J(t ), ∫ [J(t ), ∫ J(t ) dt ] dt ] dt + ... 2∫ t

t

1

τ

1

t

t1

1

τ

t1

1

τ

Ω(t,τ) ≈

(4)

τ

2

τ

2

1

t2

2

τ

3

3

2

1

∫ J(t ) dt t

τ

1

1

(5) (6)

Once Ω is obtained by numerical integration of the Jacobian, it is decomposed into its eigenvalues and eigenvectors to facilitate its exponentiation:

exp(Ω(t,τ)) ) exp(ZD∆tZ-1) ) Z exp(D∆t)Z-1

(7)

where Z is an n × n matrix containing the right-hand eigenvectors of Ω as its columns, and D is a diagonal n × n matrix containing the eigenvalues of Ω/∆t. The main approximation underlying the AIM approach occurs when eq 5 is approximated as eq 6 to obtain Ω(t,τ). Errors resulting from this approximation are limited by controlling the step size (∆t). The criterion used to determine an acceptable step size is

|(Ω(t + ∆t,t)/∆t - J(t + ∆t))K(t + ∆t,t)| e |J(t + ∆t)K(t + ∆t,t)|

(8)

where  is a user-defined tolerance. Numerical errors may also be introduced via eq 7 because the decomposition of a matrix into its eigenvalues and eigenvectors is not exact. Since the eigenvalue diagonal matrix is exponentiated and multiplied by the eigenvector matrix and its inverse, a small error in the decomposition process can result in a larger error in K(t + ∆t,t). Furthermore, since the solution is advanced in time by multiplying K(t + ∆t,t) by its predecessors, this error may be magnified as the calculations proceed. It is unlikely that an error of this type will be prevented using the step size selection criterion (eq 8), because K(t + ∆t,t) is present in both the numerator and the denominator. Thus, special care should be taken that the eigenvalue-eigenvector decompo-

TABLE 1. Chemical Species Defined in LCC Mechanism Condensed Version no.

code

description

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

Normal Species NO nitric oxide NO2 nitrogen dioxide O3 ozone HONO nitrous acid HNO3 nitric acid HNO4 pernitric acid N2O5 dinitrogen pentoxide NO3 nitrate radical HO2 hydroperoxyl radical CO carbon monoxide HCHO formaldehyde ALD2 acetaldehyde MEK methyl ethyl ketone MGLY methyl glyoxal PAN peroxyacetyl nitrate RO2 total peroxy radicals MCO3 peroxyacetyl radical ALKN alkyl nitrate ALKA C4+ alkanes ETHE ethene ALKE C3+ alkenes TOLU monoalkyl benzenes AROM di/trialkyl benzenes DIAL unknown dicarbonyls CRES o-cresol NPHE nitrophenols O(1D) oxygen (1D) O(3P) oxygen (3P) OH hydroxyl radical RO2R general RO2 radical 1 R2O2 general RO2 radical 2 RO2N alkyl nitrate RO2 RO2P phenol RO2 BZN2 benzaldehyde N-RO2 BZO phenoxy radical

36 37 38

BuildUp Product Species H2O2 hydrogen peroxide H2 hydrogen LN lost nitrogen H2O CO2 O2

Constant Species water vapor carbon dioxide oxygen

sition is accomplished as accurately as possible. Therefore, we replaced the original matrix decomposition code provided with AIM by more accurate algorithms [from the LAPACK library (20)]. With these changes, sensitivities calculated using AIM agreed with those calculated by re-running the model with different parameter values as described further in section 5. To employ the AIM method, Fortran code for the rate equations that describe the time evolution of the species concentrations (eq 1) and the Jacobian J used in eq 2 must be provided. Additionally for second-order sensitivities, the partial derivatives of the Jacobian with respect to the vector y must be included as well. The user may also provide the derivative of the vector function f with respect to the parameters R (fRj used in eq 2). In order to facilitate changing the photochemical model to accommodate different chemistry or a complete change of mechanism (on the order of 50 species and 100 reactions, see section 3), a symbolic chemical interpreter is desirable. Therefore, a new symbolic interpreter called CHESS (CHemical mEchaniSm Symbolic interpreter) was developed. CHESS is able to read a text file describing the chemical mechanism and to generate the Fortran subroutines that are required for AIM.

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1207

TABLE 2. Listing of Condensed Version of LCC Mechanism no.

reactants

productsa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

NO2 + hν O(3P) + O2 + M O(3P) + NO2 O(3P) + NO2 NO + O3 NO2 + O3 NO + NO3 NO + NO + O2 NO2 + NO3 + M N2O5 N2O5 + H2O NO2 + NO3 NO3 + hν NO3 + hν O3 + hν O3 + hν O(1D) + H2O O(1D) + M NO + OH + M HONO + hν NO2 + H2O NO2 + OH + M HNO3 + OH CO + OH O3 + OH NO + HO2 NO2 + HO2 HNO4 HNO4 + OH O3 + HO2 HO2 + HO2 HO2 + HO2 + H2O NO3 + HO2 NO3 + HO2 + H2O RO2 + NO RO2 + HO2 RO2 + RO2 RO2 + MCO3 HCHO + hν HCHO + hν HCHO + OH HCHO + NO3 HCHO + HO2 ALD2 + OH ALD2 + hν ALD2 + NO3 MCO3 + NO MCO3 + NO2 MCO3 + HO2 MCO3 + MCO3 PAN MEK + hν MEK + OH MGLY + hν MGLY + OH MGLY + NO3 ALKA + OH

58

ALKN + OH

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

RO2N + NO RO2N + HO2 RO2N + RO2 RO2N + MCO3 R2O2 + NO R2O2 + HO2 R2O2 + RO2 R2O2 + MCO3 RO2R + NO RO2R + HO2 RO2R + RO2 RO2R + MCO3 ETHE + OH ETHE + O3 ETHE + O(3P)

f NO + O(3P) f O3 + M f NO + O2 f NO3 f NO2 + O2 f NO3 + O2 f 2 NO2 f 2 NO2 f N2O5 + M f NO2 + NO3 f 2 HNO3 f NO + NO2 + O2 f NO + O2 f NO2 + O(3P) f O(3P) + O2 f O(1D) + O2 f 2 OH f O(3P) + M f HONO + M f NO + OH f HONO + HNO3 - NO2 f HNO3 + M f NO3 + H2O f HO2 + CO2 f HO2 + O2 f NO2 + OH f HNO4 f NO2 + HO2 f NO2 + O2 + H2O f OH + 2 O2 f H2O2 + O2 f H2O2 + O2 + H2O f HNO3 + O2 f HNO3 + H2O f NO f HO2 f ... f MCO3 f 2 HO2 + CO f CO + H2 f HO2 + CO + H2O f HNO3 + HO2 + CO f RO2R + RO2 f MCO3 + H2O f CO + HCHO + HO2 + RO2R + RO2 f HNO3 + MCO3 f NO2 + HCHO + CO2 + RO2R + RO2 f PAN f HCHO f 2 HO2 + 2 HCHO f MCO3 + NO2 f ALD2 + MCO3 + RO2R + RO2 f 3/2R2O2 + 3/2RO2 + MCO3 + 1/2ALD2 + 1/2HCHO f MCO3 + HO2 + CO f MCO3 + CO f HNO3 + MCO3 + CO f 0.1533HCHO + 0.4268ALD2 + 0.5378MEK + 0.102RO2N + 0.9RO2R + 0.64R2O2 + 1.64RO2 f NO2 + 0.15MEK + 1.53ALD2 + 0.16HCHO + 1.39R2O2 + 1.39RO2 f ALKN f MEK f RO2 + 1/2 HO2 + MEK f HCHO + HO2 + MEK f NO2 f ... f RO2 f HCHO + HO2 f NO2 + HO2 f ... f 1/2HO2 + RO2 f HCHO + HO2 f RO2R + RO2 + 1.56HCHO + 0.22ALD2 f 1.12HCHO + 0.42CO f HCHO + HO2 + CO + RO2R + RO2

1208

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

A

Bb

5.05E-01 and 64.8c 1.05E+03 -1282 9.30E-12 1.11E-13 -894 1.80E-12 1370 1.20E-13 2450 8.00E-12 -252 1.64E-20 -529 4.62E-13 -273 1.33E+15 11379 1.00E-21 2.50E-14 1229 1.11E+00 and 73.0 1.01E+01 and 72.1 2.73E-02 and 70.9 2.27E-03 and 41.9 2.20E-10 7.20E+08 4.03E-13 -833 9.95E-02 and 63.6 4.00E-24 9.58E-13 -737 9.40E-15 -778 2.18E-13 1.60E-12 942 3.70E-12 -240 1.02E-13 -773 4.35E+13 10103 4.00E-12 1.40E-14 579 2.27E-13 -771 3.26E-34 -2971 2.27E-13 -771 3.26E-34 -2971 4.20E-12 -180 3.00E-12 1.00E-15 3.00E-12 1.86E-03 and 51.7 2.85E-03 and 57.7 9.00E-12 6.00E-13 2060 1.00E-14 6.90E-12 -250 4.52E-04 and 49.0 3.00E-13 1427 4.20E-12 -180 2.80E-12 -180 3.00E-12 2.50E-12 2.00E+16 13542 9.71E-05 and 48.3 1.20E-11 745 8.58E-03 and 69.9 1.70E-11 3.00E-13 1427 3.84E-12 2.19E-11 4.20E-12 3.00E-12 1.00E-15 3.00E-12 4.20E-12 3.00E-12 1.00E-15 3.00E-12 4.20E-12 3.00E-12 1.00E-15 3.00E-12 2.15E-12 1.20E-14 1.04E-11

709 -180

-180

-180

-411 2634 792

TABLE 2 (Continued) no.

productsa

A

f NO2 + 2HCHO + R2O2 + RO2 f HCHO + ALD2 + RO2R + RO2 f 0.64HCHO + 0.50ALD2 + 0.13RO2R + 0.13RO2 + 0.17HO2 + 0.06OH + 0.28CO f 0.40CO + 0.40MEK + 0.40HCHO + 0.20ALD2 + 0.20HO2 + 0.60RO2R + 0.60RO2 f NO2 + HCHO + ALD2 + R2O2 + RO2 f 0.16CRES + 0.16HO2 + 0.84RO2R + 0.84RO2 + 0.40DIAL + 0.144MGLY + 0.114HCHO + 0.114CO f 0.17CRES + 0.17HO2 + 0.095HCHO + 0.095CO + 0.65DIAL + 0.316MGLY + 0.83RO2R + 0.83RO2 f MCO3 f HO2 + CO + MCO3 f 0.2MGLY + 0.15RO2P + 0.85RO2R + RO2 f HNO3 + BZO f NPHE f ... f 1/2HO2 + RO2 f HCHO + HO2 f NPHE f ... f ... f HNO3 + BZN2 f 2 LN f NPHE f NPHE

2.00E-12 2.63E-11 1.13E-17

reactants

74 75 76

ETHE + NO3 ALKE + OH ALKE + O3

77

ALKE + O(3P)

78 79

ALKE + NO3 TOLU + OH

80

AROM + OH

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

DIAL + OH DIAL + hν CRES + OH CRES + NO3 RO2P + NO RO2P + HO2 RO2P + RO2 RO2P + MCO3 BZO + NO2 BZO + HO2 BZO NPHE + NO3 BZN2 + NO2 BZN2 + HO2 BZN2

Bb 2925

3.98E-12 7.57E-15 2.10E-12

-322

2.45E-11 3.00E-11 3.25E-02 and 58.3 4.00E-11 2.20E-11 4.20E-12 -180 3.00E-12 1.00E-15 3.00E-12 1.50E-11 3.00E-12 1.00E-03 3.80E-12 1.50E-11 3.00E-12 1.00E-03

a “...” means the reaction produces species considered as inert and not accounted for in the LCC mechanism. b Arrhenius constants, k(T) ) A exp(-B/T)(T n). Units are [cm3 molecule-1 s-1] (A) and [K] (T). Values not given are assumed to be zero, n (always zero) is not listed. c For photolytic reactions, the two parameters described in section 3A are given.

3. Chemistry and Initial Conditions The model used here to perform sensitivity analysis of a photochemically reactive system is a single-cell model appropriate for simulating well-mixed reactor experiments. All reactants enter this system at time zero as initial concentrations, temperature is constant, and deposition is ignored. Temporally dependent variables such as photolytic rates can be parametrized as functions of time for sensitivity analysis, instead of treating them as constant over the course of the simulation. The condensed version of the LCC [Lurmann-CarterCoyner (21)] chemical mechanism was used in this work. The mechanism includes 35 normal species, 3 constant species, and 3 buildup product species that participate in 95 chemical reactions (see Tables 1 and 2). Concentrations of species 27-35 in the LCC mechanism are usually computed using steady-state approximations; however, in this work, rate equations for these species are integrated in the same manner as other active chemical species (1-26). Concentrations of buildup product species (36-38) are also integrated in this same manner, but these species do not reactsthey are only formed as product species. While constant species concentrations do not vary in time, their concentrations do affect the rates of chemical reactions. The LCC mechanism has been used in previous photochemical modeling (e.g., ref 22) and sensitivity analysis studies (e.g., ref 9). The mechanism is of the same vintage and comparable detail to the carbon bond IV mechanism. The CB-IV and LCC chemistry have been compared by Milford et al. (9) and other studies (23) and found to be similar. A. Photolytic Reactions. Twelve of the reactions in the LCC mechanism are photolytic with rate coefficients that depend on the solar actinic flux and are thus variable with time. These reactions are included in Table 2 and numbered 1, 13-16, 20, 39, 40, 45, 52, 54, and 82. To compute photolytic rate coefficients as functions of the solar zenith angle, most atmospheric chemistry studies use results from detailed radiative transfer models. These models give sets of data points that are interpolated or sometimes fitted to give a

continuous prediction. Lurmann et al. (21) used model predictions of the solar actinic flux by Peterson (24) to compute the photolytic reaction rate coefficients for different solar zenith angles using the expression:

kiph(θz) )

∑F (θ )σ λ

z

(9)

λ,iφλ,i

λ

where F is the actinic flux for wavelength λ and solar zenith angle θz, σ is the absorption cross-section, and φ is the quantum yield. This calculation method is not suitable here since constant parameters are required in the AIM method for parametric sensitivity analysis. For local sensitivity analysis, prior studies (9, 12) have treated photolysis rate coefficients as constant over the entire day. Uncertainty studies (13-15, 25) treated the sensitivity parameter as an overall multiplicative factor to be applied uniformly to the photolysis rate coefficient. Here, this constraint is overcome, while retaining the temporal dependence of the photolytic reaction rate coefficients, by fitting the estimates by Lurmann et al. with a simple function of time using constant parameters. Under clear-sky conditions, the light input depends mainly on the sun’s position in the sky and atmospheric optical depth in near-UV wavelengths. Taking light input as the dominant factor for photolytic processes, we only need to consider the solar zenith angle (θz) dependence. The optical depth increases linearly with the path length of the solar beam through the atmosphere. However, increasing optical depth has greater effect at shorter wavelengths in the near-UV portion of the solar spectrum. Therefore, it is not sufficient to use only a linear dependence in the solar beam path length through the atmosphere for predicting the dependence of kiph on θz. A good approximation is

kiph ) ai cos (θz) + bi cos2 (θz)

(10)

where θz can be readily computed as a function of date, time, and latitude (26). The two free-fitting parameters ai and bi

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1209

FIGURE 1. Parametrization of the solar zenith angle dependence of photolytic reaction rate coefficients. Specifically, the variation of photolysis rate coefficients for NO2 (a) (reaction 1) and O3 to O(1D) (b) (reaction 16) with θz are given as examples. can be related to the geometric parameters βi and γi shown in Figure 1:

kiph )

(

)

βi 1 - 2 cos (γi) cos (θz) + 2 cos2 (γ ) - cos (γ ) i i

} } ai

(

)

βi 2 cos2 (γi) - 1 cos2 (θz) 2 cos2 (γ ) - cos (γ ) i i

used by Milford et al. (9). When variable rate coefficients were used, the sensitivities to the parameters characterizing the time dependence (βi and γi in eq 11) were computed. Parameters to which ozone showed large first-order sensitivity were selected, and second-order sensitivities were then computed for all pair combinations of the selected parameters. First- and second-order sensitivities were normalized as

(11)

∂[A] R ∂R [A]

and

∂2[A] R1R2 ∂R1∂R2 [A]

(12)

bi This way βi is the height (value of kiph at θz ) 0°) and γi is the width at half maximum. In an attempt to explore the effect of diurnal variations in the rate coefficients on the results of sensitivity analysis, we compared calculations performed with the fitted rate coefficients with ones performed with constant photolytic rate coefficients that have been used in past studies (9). The data used to compute the solar zenith angle are appropriate for Los Angeles, CA, in late summer. B. Initial Conditions. While most of the species shown in Table 1 are absent initially, nitrogen oxides (NOx) and seven reactive organic gases (ROG) concentrations were set to nonzero initial values given in Table 3. Two sets of initial concentrations were chosen: one with a low ROG to NOx ratio (6 ppm of C/ppm of NOx) and one with a high ratio (24 ppm of C/ppm of NOx). These initial conditions represent cases where O3 formation is ROG-limited and NOx-limited, respectively. Temperature (298 K), pressure (101 325 Pa), and water vapor concentration (9824 ppm) were held constant throughout the simulations. The simulation time is 8:00 a.m. to 4:00 p.m. (PST) on August 28, 1995, and the location is Los Angeles, CA, (latitude 34.058° N, longitude 118.250° W). With these data, the sun rises at 5:27 a.m. and sets at 6:20 p.m.

4. Results and Discussion Comparisons were made between results of calculations with low (6) and high (24) initial ROG to NOx ratio. Calculations with variable and constant photolytic reaction rate coefficients were also compared. First-order sensitivities to both initial concentrations and reaction rate coefficients were examined, and sensitivities were ranked to identify the most important parameters for further study. Time-varying photolytic reaction rate coefficients were parametrized as described previously. Three sets of constant photolytic reaction rate coefficients were also used: 9:00 a.m. values, 12:00 p.m. values, and the daytime average values

1210

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

First-order sensitivities to initial concentrations and reaction rate coefficients were rank-ordered separately. Sensitivities to the nine non-zero initial concentrations and the full set of reaction rate parameters were considered. Normalized sensitivities to reaction rate parameters were ranked according to their maximum absolute values. In addition to sensitivities of O3, sensitivities to all input variables were also computed for NO, NO2, HONO, HO2, HCHO, RO2, MCO3, and OH. These species were selected because of their importance to ozone formation. Considering sensitivities of the nine species listed above (including ozone) or considering ozone sensitivities alone gave similar results, with seven reactions (see Figure 2) consistently emerging as the most important ones. The set of selected reactions is consistent with findings of previous studies (9, 12, 13). The nine nonzero initial concentrations and the seven key reaction rate parameters constituted the set of important parameters selected for more detailed study. Such a selection scheme is useful for choosing parameter combinations for secondorder sensitivities that can be use to extend the range of sensitivity analysis predictions. A. Comparison of Model Scenarios. Sensitivities computed with time-dependent photolytic reaction rate coefficients were similar to those computed with constant photolytic reaction rate coefficients (see Figure 3). However, at ROG/NOx ) 24, the shape of the sensitivities is slightly different when using constant photolytic reaction rate coefficients that are averaged over the day (9). The high ROG/ NOx ratio simulation is divided in two parts: initially the sensitivities undergo large changes, then they remain approximately constant. The initial period lasts longer when using averaged photolytic reaction rate coefficients than in the two other cases. The parametrization of photolysis rate using eq 11 allowed direct comparison of sensitivity to the parameters βi (the maximum values of the photolytic reaction rate coefficients) with sensitivity to the reaction rate coef-

TABLE 3. Initial Concentrations for Low and High ROG to NOx Ratio ROG/NOx ) 6

ROG/NOx ) 24

species

initial concn (ppmv)

species

initial concn (ppmv)

NO NO2 HONO HCHO ALKA ETHE ALKE TOLU AROM H2O

0.12895 0.03224 0.00163 0.0413 0.1059 0.0275 0.0184 0.0197 0.0172 9824.0

NO NO2 HONO HCHO ALKA ETHE ALKE TOLU AROM H2O

0.03224 0.00806 0.00041 0.0413 0.1059 0.0275 0.0184 0.0197 0.0172 9824.0

FIGURE 3. Comparison of largest sensitivities of O3 to reaction rate parameters computed with variable and constant photolytic reaction rate coefficients (ROG/NOx ) 24).

FIGURE 2. Largest normalized sensitivities of O3 to reaction rate coefficients vs time. Results are given for low and high ROG to NOx ratio with the photolytic reaction rate coefficients parametrized with respect to time. ficients for non-photolytic reactions or for the case where photolytic rate coefficients were held constant. Large changes in ozone sensitivities occurred when the ROG to NOx ratio was increased from 6 to 24. At low ROG to NOx ratio, the ozone sensitivities to the set of important parameters have a smooth shape for the entire 8-h simulation (see Figures 2 and 4). Accordingly, the rate of ozone formation does not vary much during the simulation (see Figure 5). Radical concentrations are plotted in Figure 6: [OH] builds up very quickly to reach a maximum 1 h after the beginning of the simulation and then slowly decreases; whereas [HO2] and [RO2] increase at a rate nearly constant from the beginning to the end of the simulation. In contrast, at a high ROG to NOx ratio, the sensitivities undergo rapid changes in the first 200 min and then remain nearly constant until the end of the simulation. This can be understood by noticing (see Figure 5) that during the first 200 min, almost all the NO and NO2 are consumed and most of the ozone is formed. The corresponding OH concentration reaches a

FIGURE 4. Most important sensitivity of O3 to initial conditions vs time. Results are given for low and high ROG to NOx ratio with the photolytic reaction rate coefficients parametrized with respect to time. maximum during the period of strong ozone production and decreases rapidly afterward, whereas [HO2] and [RO2] peak at the end of this initial period and slowly decrease afterward (see Figure 6).

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1211

FIGURE 5. Concentration of species important for O3 formation vs time. Results are given for low and high ROG to NOx ratio with the photolytic reaction rate coefficients parametrized with respect to time.

FIGURE 6. Radical concentrations vs time. Results are given for low and high ROG to NOx ratio with the photolytic reaction rate coefficients parametrized with respect to time. OH concentration is scaled. B. Sensitivities to Initial Conditions. Normalized sensitivities of ozone to initial conditions shown in Figure 4 indicate that NO, NO2, HCHO, and lumped aromatics (AROM)

1212

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

are the species whose initial concentrations are the most important for both ROG/NOx ) 6 and ROG/NOx ) 24. In both cases, the largest [O3] sensitivity to initial conditions is the sensitivity to [NO]0. The distinction between ROG-limited and NOx-limited cases is clearly demonstrated in Figure 4 by these sensitivities since their shape differs between the two cases; some sensitivities even change sign (after 200 min) when switching from the ROG to the NOx-limited case. Similar results have been obtained and explained in previous studies (see refs 8 and 9). When ozone production is ROG-limited, the shape of the [O3] sensitivities is smooth. When the NOx concentration is limiting, the chemistry proceeds in two phases. During the first 3 h, when ozone formation is most rapid, ozone sensitivity to initial concentrations exhibits a pattern similar to the one found in the ROG-limited case, but with a shorter time scale. For this initial period, the shape of the ozone sensitivities can be explained in the same way as in the ROG-limited case. Then the situation reverses, the sensitivities to [NO]0 and [NO2]0 become positive and the sensitivities to [HCHO]0 and [AROM]0 become very small or negative. The consumption of nearly all NOx during the first 3 h explains the ozone sensitivity sign change. The ozone concentration growth is largest between t ) 20 min and t ) 150 min when most of the ozone is produced. During this time, the sensitivities of [O3] to [NO]0 and [NO2]0 are negative because NO scavenges ozone and NO2 reacts with hydroxyl radicals whose concentration is, at this time, relatively high. Hence the presence of NO and NO2 still diminishes ozone production; but the smaller the NOx concentrations, the more rapidly the ozone production will cease. Therefore, an increase in NOx initial concentrations will result in a smaller ozone production rate during the period between 20 and 150 min (negative ozone sensitivity), but will result in a longer ozone production period that will balance the initial smaller production rate and result in greater final ozone concentrations (positive final ozone sensitivity). This suggests that the timing, not just the total mass of NOx emissions, may be important when studying ozone production because the picture described above is likely to change if the initial NOx concentration was lower and continuous NOx emissions were maintained throughout the run. The Green’s function can reveal the dependence of species concentration at any time to a variation of another species concentration at an earlier time. It is apparent from eq 3 that Kij(t,t′) is the first-order sensitivity of the concentration of species i at time t to the concentration of species j at time t′. The dependence of [O3](t) to [NO2](t′) is shown in Figure 7. The sensitivity of [O3] at various times to a change of [NO2] at a given time is given by slices parallel to the t-axis (the leftmost slice being the sensitivity to initial condition). Slices parallel to the t′-axis give the sensitivity of [O3] at a given time to [NO2] change at various times. The rearmost slice gives the dependence of the final O3 concentration to changes of [NO2] at different times. The Green’s function is plotted unnormalized: it gives the variation of [O3] (according to the first-order sensitivity) (in ppmv) for a fixed (1 ppmv) variation of [NO2]. The maximal sensitivity of the final [O3] occurs for a change in [NO2] at t ) 145 min. For changes at later t′, the variation of the final ozone concentration decreases to zero (a change of [NO2] at t′ cannot influence [O3] at t ) t′ because the chemistry is not instantaneous). Thus, the largest effect of adding NO2 at different times occurs when it is added as the NO concentration approaches zero. Thus, model calculations will yield different results if the temporal nature of the NOx emissions is changed. C. Sensitivities to Reaction Rate Coefficients. On the basis of the sensitivities to the reaction rate coefficients, seven reactions appeared consistently in the top 10 most important reactions for all the simulations that were examined. The [O3] sensitivities to the rate parameters of these reactions are shown in Figure 2. When the reactions were photolytic, the

FIGURE 7. K3,2 (t,t′): Sensitivity of [O3](t) to [NO2](t′). ROG/NOx ) 24, and the photolytic reaction rate coefficients are parametrized with respect to time. parameter βi, specifying the maximum value (i.e., for a solar zenith angle of 0°) of the time-dependent reaction rate coefficient was considered. Inspection of Figure 2 reveals cases where sensitivities to some pairs of rate coefficients are symmetric with respect to the time axis. This phenomenon has been identified previously in sensitivity analysis applied to combustion kinetics (5) and can also be deduced from data presented in some previous atmospheric chemistry sensitivity analyses (9, 15, 25). Such effects occur when pairs of reactions have equal but opposite effects on the dependent variable (here ozone concentration) and are balanced with nearly equal rates. For example, reactions 1 (NO2 photolysis) and 5 (NO + O3 f NO2 + Ο2) have opposing effects on the distribution of NOx between NO and NO2 that affects ozone concentration directly. Reactions 19 (NO + OH f HONO) and 20 (HONO photolysis) have exactly opposite effects, and the sensitivity of [O3] to these two reactions shows opposite values at all times. Nitrous acid reactions only play an important role very early in the simulation because HONO photolysis is an early morning source of hydroxyl radical. Ozone sensitivity to reactions 22 (NO2 + OH f HNO3) and 39 (HCHO photolysis) are also nearly symmetric with respect to the time axis, although less obviously than in the two examples above. Reaction 22 is a termination step since it recombines both NO2 and OH into HNO3. In contrast, reaction 39 is an initiation step that generates radicals. However, reactions 22 and 39 do not influence ozone formation in a manner that is as obviously opposite as the reaction pairs 1 and 5 and 19 and 20. The strong sensitivity of ozone concentration to formaldehyde photolysis and to formaldehyde initial concentration illustrates the important role formaldehyde plays in ozone formation. The importance

of aromatic hydrocarbon chemistry is also emphasized by sensitivity analysis: both the sensitivities to initial aromatic concentration and to its oxidation reaction rate coefficient were among the most important sensitivities. While further study of formaldehyde photolysis and aromatic chemistry has occurred in recent years and has led to an update of chemical mechanisms, large uncertainties still remain (12, 13), and more study of HCHO photolysis and aromatic chemistry is needed. D. Second-Order Sensitivities. Second-order sensitivities are useful for extending the domain of applicability of local sensitivity analysis. First-order sensitivities allow the prediction of the effect of a parameter variation over the limited range where the system response can be linearized. With second-order sensitivity, the prediction can be extended beyond this limit to include non-linear effects. Second-order sensitivities also allow formal consideration of the combined effect of two-parameter variations. Yetter et al. (5) noticed that when concentrations exhibit large first-order sensitivities to one parameter, it is more likely that the second-order sensitivities to that parameter will be large. We confirmed this to be also true in our case for the diagonal second-order sensitivities. Mixed second-order sensitivities are too numerous to be checked exhaustively. Therefore, only the second-order sensitivities to the nine nonzero initial concentrations and seven key reaction rate parameters identified previously are studied here. When combined in pairs for second-order sensitivity analysis, these 16 parameters yield 136 [n(n + 1)/2] potentially large secondorder sensitivities for each species. Normalized second-order sensitivities of [O3] and also [NO], [NO2], [HONO], [HO2], [HCHO], [RO2], [MCO3], and [OH] were ranked by their

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1213

FIGURE 9. Variation of ∆[O3] ) [O3](r + ∆r) - [O3](r) with respect to NO initial concentration variation. The dashed curve is obtained by varying the parameter and re-running the model, while the dashdotted and dotted curves are obtained using the first- and first- plus second-order sensitivities, respectively, at the parameter’s nominal value. FIGURE 8. Largest second-order [O3] sensitivities vs time. Results are given for low and high ROG to NOx ratio with the photolytic reaction rate coefficients parametrized with respect to time. absolute maximum values. Seven second-order sensitivities consistently appeared in the 10 largest sensitivities of ozone and several other species considered for both low and high ROG to NOx ratios. In addition, the diagonal second-order sensitivity of ozone to NO recombination with ozone (reaction 5) was large for ozone but less important for the other species considered. The eight largest second-order sensitivities for ozone are plotted in Figure 8. Ozone sensitivities to [NO]0, [HCHO]0, and [AROM]0 were the largest normalized first-order sensitivities to initial conditions. Similarly, first-order sensitivities to reaction of NO2 with OH (reaction 22) and to formaldehyde photolysis (reaction 39) showed the largest peak value among sensitivities to reaction rate parameters. Six of the eight second-order sensitivities displayed in Figure 8 are sensitivities to a combination of the five parameters enumerated here. This also is consistent with the largest second-order sensitivities being sensitivities to the parameters already inducing the largest first-order sensitivities. Inspection of the second-order sensitivities plotted in Figure 8 shows a trend already indicated by Yetter et al. (5). Diagonal second-order sensitivities balance the effects of firstorder sensitivities to prevent the prediction of unphysical (e.g., negative) concentrations. The range of concentration variations is bounded by the nature of the model. Concentrations cannot be negative, and there is an upper limit set by the total mass availability. Physical and chemical considerations also impose limits on parameters (initial concentrations and reaction rate coefficients). However some parameters, such as initial concentrations, can vary over an order of magnitude or more, and predictions of the effect of large parameter variations with first-order sensitivity analysis can become meaningless. For example, first-order sensitivity analysis could predict ozone concentrations negative or in excess of the total amount of ozone produced if all available ROG and NOx had been consumed solely to produce ozone. The sign of large diagonal second-order sensitivities can be guessed by assuming that contributions by higher-order

1214

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

sensitivity will assure the prediction of meaningful concentration changes. Figure 8 contains three examples of [O3] diagonal second-order sensitivity to (1) [NO]0; (2) the rate coefficient of NO recombination with O3 (reaction 5); and (3) the rate coefficient of NO2 reaction with OH (reaction 22). At a low ROG to NOx ratio and during the first 3 h at a high ROG to NOx ratio, the first-order sensitivity of [O3] to [NO]0 is negative. After the first 3 h at a high ROG to NOx ratio, it is positive. Thus, an arbitrarily large positive change of the NO initial concentration could lead to a first-order sensitivity prediction of a negative ozone concentration at low ROG to NOx ratio and at high ratio during the first 3 h. After the first 3 h at a high ROG to NOx ratio, an arbitrarily high ozone concentration could be predicted. The second-order sensitivity of [O3] to [NO]0 is always of the opposite sign to the first-order sensitivity to prevent this. In Figure 9, the firstand second-order prediction of [O3] variation resulting from an [NO]0 variation at low ROG to NOx ratio is shown. For a variation of +45% of the parameter’s nominal value, the firstorder analysis predicts a large decrease in [O3], whereas the prediction using both first- and second-order is less negative and much closer to the ozone concentration variation determined by varying the [NO]0 and re-running the model. The same relation is found between the first-order [O3] sensitivities to k5 and k22 and the corresponding diagonal second-order sensitivities. The first-order sensitivities of [O3] to k5 and k22 are always negative, and the diagonal secondorder sensitivities are always positive. The similarity among three ozone mixed second-order sensitivities (∂2[O3]/∂[NO]0 ∂[HCHO]0, ∂2[O3]/∂[NO]0 ∂[AROM]0, and ∂2[O3]/∂[NO]0 ∂β39) is noticeable. These second-order sensitivities are nearly equal at all times for both high and low ROG to NOx ratios. Another mixed second-order sensitivity (∂2[O3]/∂[NO]0 ∂k22) has a similar shape but is the opposite sign. These are cases of self-similarity among second-order sensitivities. Self-similarity implies that a single dependent variable dominates the behavior of the selfsimilar sensitivities (27). In this case, the dominant factor for the four self-similar second-order sensitivities is the firstorder ozone sensitivity to the NO initial concentration. This implies that a change in NO initial concentration will affect

FIGURE 10. Variation of ∆[O3] ) [O3](r1 + ∆r1, r2 + ∆r2) - [O3](r1, r2) with respect to two parameter variations. The dashed, dash-dotted and dotted curves are obtained in the same way as in Figure 9. equally the ozone sensitivity to [HCHO]0, [AROM]0, and β39 and, in an equal but opposite manner, the ozone sensitivity to k22. The normalized second-order sensitivities shown in Figure 8 often have peak absolute value exceeding 1. On the other hand, the absolute peak value of the corresponding firstorder sensitivity rarely exceeds 1. Thus, when predicting the influence of a parameter variation with a Taylor expansion, the contribution of second-order terms will dominate the first-order contribution if the parameter variation is of the order of (or greater than) the parameter nominal value. While the sign of diagonal second-order terms can be predicted, there is no rule of thumb to guess the sign of the off-diagonal (or mixed) second-order sensitivities, and it is not possible to predict their influence. Computing second-order sensitivity is necessary if local sensitivity analysis results are used to predict responses to large parameter variations. 5. Validation of the Taylor Expansion. Sensitivity analysis using first-order derivatives is called local because it allows the prediction of the effect of parameter variations only over a limited range where this effect can be linearized. The prediction of the effect of parameter variations can be extended if higher order derivatives are used, and the response of the model could be predicted exactly for any value of

parameter variation if an infinite-order Taylor expansion were used. The reliability of the sensitivity analysis in first- and second-order was established by computing the species concentration changes induced by parameter variation with a Taylor expansion using the first- and second-order sensitivities. Taylor series estimates were compared with results of direct calculations with perturbed parameter values. The calculation with perturbed parameter values would in theory give the same results as an infinite-order Taylor expansion. The direct calculations are based on numerical methods rather than analytical evaluation and, as such, are subject to inaccuracies that can be large relative to the concentration changes predicted. Fortunately, such errors are only of consequence for cases of small sensitivities and small parameter variations. To assess the maximum range for which first- or second-order sensitivity analysis is valid, the parameters were perturbed one at a time and the following expressions were compared:

∆[A] ) [A](R + ∆R) - [A](R)

(13)

∂[A] (R) ∂R

(14)

with

∆[A]T1 ) ∆R

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1215

and

∂ 2[A] ∂[A] 1 ∆[A]T2 ) ∆R (R) (R) + (∆R)2 ∂R 2 ∂ R2

(15)

where [A] is a chemical species concentration, R is a parameter, ∆[A]T1 and ∆[A]T2 are the estimated changes in [A] using first- and second-order Taylor series, respectively. These comparisons were performed for every species and for perturbation of every parameter. A similar procedure was repeated with the parameters selected for detailed sensitivity analysis perturbed two by two in all possible combinations. Every coefficient, including the diagonal second-order terms as well as the mixed second-order terms, was considered to compare the following:

∆[A] ) [A](R + ∆R) - [A](R)

(16)

with

∂[A] ∂[A] ∆[A]T1 ) ∆R1 (R ,R ) + ∆R2 (R ,R ) ∂ R1 1 2 ∂ R2 1 2

(17)

and

∂[A] ∂[A] ∆[A]T2 ) ∆R1 (R ,R ) + ∆R2 (R ,R ) + ∂ R1 1 2 ∂ R2 1 2

[

]

∂ 2[A] ∂ 2[A] 1 (∆R1)2 (R1,R2) + (∆R2)2 (R1,R2) + 2 2 ∂R ∂R2 1

2

2

∂ [A] (R ,R ) (18) ∆R1 ∆R2 ∂ R1∂ R2 1 2 Plots of ∆[A], ∆[A]T1, and ∆[A]T2 were compared for different values of ∆R. The first-order expansion is adequate if the prediction made using it (∆[A]T1) is nearly equal to the real model prediction (∆[A]). The second-order expansion is required if the second-order prediction (∆[A]T2) is significantly closer to ∆[A] than ∆[A]T1. When only the first-order and the diagonal elements of second-order sensitivities were tested, the parameter R was changed by -45%, -30%, ..., 30% and 45% as shown in Figure 9. When first-order and all elements of second-order sensitivities were tested, two parameters had to be varied, and all combinations of both parameters changed by -50%, -25%, 0%, 25%, and 50% were considered as shown in Figure 10. The second-order Taylor expansions almost always remained accurate for parameter variations as large as (50% of the nominal value. First-order sensitivities alone were sometimes sufficient to predict changes induced by (50% variations of the parameters, but remained accurate, in general, only when parameter variations were limited to at most (25%. Therefore the ability to compute second-order sensitivities is important when the effects of larger variations in the parameters are being studied. Not surprisingly, the largest differences were found for the parameters to which [O3] is the most sensitive. Variations of the parameters γi of eq 11 resulted in some of the largest differences. The variation of γi strongly influences the time dependence of the photolytic reactions, for example, the rate of photolytic reactions at large solar zenith angle is highly sensitive to the width parameter. Fortunately, the γi parameters are not uncertain by large amounts, and the range of parameter variation used in this test (up to (50%) is not suitable for such parameters. The next largest differences were found when the NO initial concentration was varied (see Figure 9). The prediction of the variation made with the first-order sensitivity computed at the nominal value of the parameter is of the same sign and order of magnitude as the calculation made by perturbing the parameter. Sometimes the results differ by a factor of 2, and the second-order correction brings the sensitivity prediction in much better agreement with the perturbed calculation. This conclusion is valid for parameter variations over the whole range tested

1216

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 4, 1997

((45% of the nominal value). The parameters varied in the example of two-parameter variation presented in Figure 10 also have a large influence on ozone concentration. Parameters varied are the rate coefficient of the NO2 reaction with hydroxyl radical (reaction 22) and one of the parameters (β39, the maximum height) used to describe the rate coefficient for formaldehyde photolysis. When two parameters are varied, the secondorder contribution consists of two diagonal terms (∂2[O3]/∂ k222 and ∂2[O3]/∂ β239) and a mixed one (∂2[O3]/∂ k22 ∂β39). When both parameters are varied by (50% (lower-left and upper-right plots), the contributions of the first-order sensitivities for both parameters almost cancel, as do the secondorder contributions. When both parameters are varied in opposite directions (upper-left and lower-right plots), the firstorder contributions add together. In the upper-left plot, the first-order prediction is less than the concentration change predicted by direct calculation, but the second-order contribution added on top of the first-order prediction brings the second-order prediction in much better agreement. In the lower-right plot, the second-order contribution again brings the Taylor expansion prediction in better agreement with the concentration change predicted by direct calculation, but this time by canceling part of the excessive first-order prediction. Second-order contributions were usually small or negligible are compared to the first-order contribution. However, there are cases when the second-order contribution is as important as the first-order one. This is more likely if the parameters are especially important for predicting species under consideration (see section 4D). A good example is the one shown in Figure 10.

Acknowledgments The authors wish to thank Thomas Kirchstetter of UC Berkeley. Financial support was provided by the Swiss National Fund for Scientific Research, and by the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences Division of the U.S. Department of Energy, under Contract DE-AC03-76SF00098.

Literature Cited (1) Dickinson, R. P.; Gelinas, R. J. J. Comp. Phys. 1976, 21, 123143. (2) Dunker, A. M. J. Chem. Phys. 1984, 81, 2385-2393. (3) Hwang, J.-T.; Dougherty, E. P.; Rabitz, S.; Rabitz, H. J. Chem. Phys. 1978, 69, 5180-5191. (4) Hamby, D. M. Environ. Monit. Assess. 1994, 32, 135-154. (5) Yetter, R. A.; Dreyer, F. L.; Rabitz, H. Combust. Flame 1985, 59, 107-133. (6) Vajda, S.; Rabitz, H. Chem. Eng. Sci. 1991, 47, 1063-1078. (7) Brown, N. J.; Koszykowski, M. L. Int. J. Chem. Kinet. In press. (8) Falls, A. H.; McRae, G. J.; Seinfeld, J. H. Int. J. Chem. Kinet. 1979, 11, 1137-1162. (9) Milford, J. B.; Gao, D.; Russell, A. G.; McRae, G. J. Environ. Sci. Technol. 1992, 26, 1179-1189. (10) Carter, W. P. L. J. Air Waste Manage. Assoc. 1994, 44, 881-899. (11) Chock, D. P.; Yarwood G.; Dunker A. M.; Morris R. E.; Pollack A. K.; Schleyer C. H. Atmos. Environ. 1995, 29, 3067-3084. (12) Gao, D.; Stockwell, W. R.; Milford, J. B. J. Geophys. Res. 1995, 100, 23153-23166. (13) Gao, D.; Stockwell, W. R.; Milford, J. B. J. Geophys. Res. 1996, 101, 9107-9119. (14) Yang, Y. J.; Milford, J. B. Environ. Sci. Technol. 1996, 30, 1591-1599. (15) Yang, Y. J.; Stockwell, W. R.; Milford, J. B. Environ. Sci. Technol. 1996, 30, 1392-1397. (16) Kramer, M. A.; Calo, J. M.; Rabitz, H. Appl. Math. Modell. 1981, 5, 432-441. (17) Kramer, M. A.; Calo, J. M.; Rabitz, H.; Kee, R. J. AIM: The Analytically Integrated Magnus Method for Linear and Second Order Sensitivity Coefficients; Sandia National Laboratories: Albuquerque, NM, and Livermore, CA, 1982; SAND82-8231 for the United States Department of Energy under Contract DEAC04-76DP00789.

(18) Kramer, M. A.; Rabitz, H.; Calo, J. M.; Kee, R. J. Int. J. Chem. Kinet. 1984, 16, 559-578. (19) Kramer, M. A.; Rabitz, H.; Calo, J. M. Appl. Math. Modell. 1984, 8, 341-350. (20) Anderson, E.; et al. LAPACK Users’ Guide, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1995. (21) Lurmann, F. W.; Carter, W. P. L.; Coyner, L. A. A Surrogate Species Chemical Reaction Mechanism for Urban-Scale Air Quality Simulation Models; ERT Inc., Newbury Park, CA, and Statewide Air Pollution Research Center, University of California: Riverside, CA, 1987; Vols. I and II; Report to the U.S. Environmental Protection Agency under Contract 68-02-4104. (22) Harley R. A.; Russell A. G.; Mcrae G. J.; Cass G. R.; Seinfeld J. H. Environ. Sci. Technol. 1993, 27, 378-388. (23) Dodge, M. C. J. Geophys. Res. 1990, 95, 3635-3648. (24) Peterson, J. T. Calculated Actinic Fluxes (290-700 nm) for Air Pollution Photochemistry Applications; U.S. Environmental

Protection Agency, Research Triangle Park, NC, 1976; EPA-600/ 4-76-025. (25) Yang, Y. J.; Stockwell, W. R.; Milford, J. B. Environ. Sci. Technol. 1995, 29, 1336-1345. (26) McRae, G. J.; Goodin, W. R.; Seinfeld, J. H. Mathematical Modeling of Photochemical Air Pollution; Environmental Quality Lab Report 18; California Institute of Technology: Pasadena, CA, 1982. (27) Rabitz, H.; Smooke, M. D. J. Phys. Chem. 1988, 92, 1110-1119.

Received for review August 22, 1996. Revised manuscript received December 5, 1996. Accepted December 6, 1996.X ES960727G X

Abstract published in Advance ACS Abstracts, February 15, 1997.

VOL. 31, NO. 4, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1217