Fast, Direct Sensitivity Analysis of Multidimensional ... - ACS Publications

Environ. Sci. Technol. 1997, 31, 2859-2868

Fast, Direct Sensitivity Analysis of Multidimensional Photochemical Models Y U E H - J I U N Y A N G , * ,† J A M E S G . W I L K I N S O N , †,‡ A N D ARMISTEAD G. RUSSELL† School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0512, and Department of Engineering and Public Policy, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213

Sensitivity analysis is a powerful tool for understanding the response of an environmental system (e.g., air quality) or model of that system to both inputs and system parameters. A fast, formal sensitivity analysis method was developed for application to multidimensional, chemicallyactive environmental models. The technique was then implemented in a three-dimensional air quality model and applied to the South Coast Air Basin (SoCAB) of California. Using direct derivatives of the equations governing the evolution of species concentrations, the local sensitivities to a variety of model parameters (e.g., rate constants, dry deposition velocities, wind speed) and inputs (e.g., initial concentrations, ground-level emissions, wind speed) are computed simultaneously. Since the equations governing sensitivity coefficients have a structure similar to that of the pollutant concentrations, the implementation using this technique is straightforward and computationally efficient.

Introduction Photochemical air quality models increasingly are being used to understand the atmospheric dynamics of air pollutants and as the basis of emission control regulations. The response of the these model predictions to system parameters or emission controls provides valuable information for control strategy design to improve air quality. Such information can be pursued via sensitivity analysis, the systematic calculation of sensitivity coefficients, to quantitatively measure these dependencies. However, sensitivity analysis has not seen as wide of use as desired, in part because of the implementation complexity as well as computational limitations. For these reasons, sensitivity analysis has been applied primarily to subsystems of air quality models [e.g., Koda et al. (1), Rabitz et al. (2), and Milford et al. (3)] or to limited aspects in air quality models (4). The “brute-force” method has been the most frequently used to determine model sensitivities. The brute-force method involves one at-a-time parameter perturbation and repeated solution of the model [e.g., Seigneur et al. (5) and Dechaux et al. (6)]. Regardless of a model formulation and complexity, the brute-force method provides an easy way to calculate the sensitivity coefficients, but it * Author to whom correspondence should be addressed; e-mail: [email protected]. † Georgia Institute of Technology. ‡ Carnegie Mellon University.

S0013-936X(97)00117-X CCC: $14.00

 1997 American Chemical Society

rapidly becomes less viable and prohibitively inefficient for a model when a large number of sensitivity coefficients need to be computed. Also, there are significant numerical limitations. These limitations, however, can be overcome by directly solving sensitivity equations using a formal sensitivity analysis. A number of approaches have been developed to calculate sensitivity coefficients. One class of techniques is built on the coupled direct method [e.g., Dickerson et al. (7)]. In this method, the sensitivity equations are derived from the model equations and solved simultaneously with the model equations. The total number of equations for the solution is N × (M + 1), where N is the species number and M is the number of parameters for the sensitivity equations. This method was found to be unstable and inefficient when applied to stiff equations found in many problems (8). A second set of techniques relies on the Green’s function (2, 4, 9) or adjoint method in which the sensitivity coefficients are computed from integrals of the Green’s function of sensitivity equations derived from the model equations. The implementation of this method is unwieldy for most current, comprehensive air quality models, especially for those having a modular algorithm structure. Another approach for computing sensitivity coefficients is the decoupled direct method (DDM) (8) in which the sensitivity equations are derived from the model equations but solved separately from the model equations. DDM does not share the instability problem found with the direct methods or adjoint methods (8). Furthermore, the implementation of this method is more straightforward than the coupled direct or adjoint methods since the sensitivity equations are similar to the model equations. Therefore, only minimal modifications in the model algorithms are required to simultaneously calculate the sensitivity coefficients. DDM has been applied to zero-dimensional models (3, 8, 10, 11), and a somewhat different approach in a threedimensional model (12). Another technique is ADIFOR (Automatic DIfferentiation in FORtran) (13), which automatically translates large FORTRAN codes to a subprogram that includes the original functions as well as those for the desired sensitivity coefficients. This method has been used in past studies, notably for sensitivity analysis of the advection equation as used for atmospheric modeling (14) or for atmospheric chemistry alone (15). Because ADIFOR is designed for general purpose sensitivity analysis, the expanded codes do not take advantage of the program structure and re-use of calculations. Also, to compute some sensitivity coefficients, such as those with respect to the subdomain emissions or the boundary conditions, requires additional modifications that can be cumbersome. In this paper, a fast, direct, multidimensional sensitivity analysis similar to the DDM was developed and applied to the CIT (California/Carnegie Institute of Technology) airshed model, which is a three-dimensional, photochemical air quality model. The method is referred to as DDM-3D. The technique, as applied here, computes sensitivity coefficients, including those with respect to initial conditions, dry deposition velocities, rate constants, and ground-level emissions in the SoCAB during the August 27-29, 1987, ozone episode. Furthermore, the approach can also be applied to other multidimensional air quality models. Below, the derivation of the method is given, followed by results of this analysis.

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

9

2859

Methodology Formulation of the Multidimensional Direct Differentiation Method. The formation and transport of chemically reacting atmospheric contaminants is described by the atmospheric diffusion equation (ADE) (16):

∂ci ) -∇(uci) + ∇(K∇ci) + Ri[c1, c2, ..., cn; T, t] + ∂t i ) 1, ..., N (1) Si where ci(x, t) is the ensemble average concentration of species i, u(x, t) is the wind field, K(x, t) is a second-order turbulent diffusivity tensor, Ri ) Ri[c1, c2, ..., cn, t] is the net rate of chemical production, Si(x, t) is the elevated source rate of emissions of specie i, and N is the number of total species. This equation is solved numerically subject to initial conditions (IC) and boundary conditions (BC):

IC

nominal values

parameters O3 initial condition, ppb O3 dry deposition velocity, cm/s HO + NO2 f HNO3 rate constant, ppm-1 min-1 b ground-level emissions,c t/day NOx(M) NOx(A) VOC(W) VOC(E) wind speed, m/s

60 0.02-2.0a 1.7 × 10-4 668 365 1447 297 0.0-7.0

a Depends on terrain difference. b Evaluated at 300 K but varies at locations in the model region. c M and A in the parentheses indicate mobile and area source; W and E specify the western and eastern modeling domain.

and initial and boundary conditions are obtained as follows:

∂Ri ∂Si ∂s*ij + ) -∇(us*ij) + ∇(K∇s*ij) + Jiks*kj + ∂t ∂j ∂j ∇(uci)δij + ∇(K∇ci)δij (3)

ci(to) ) coi BCs (1) uci - K∇ci ) ucbi

(3) vigci - Kzz

(4) -

IC

horizontal inflow

s*ij(to) ) coj δij

(2) - ∇ci ) 0

horizontal outflow BCs ∂ci ) Ei ∂z

∂ci )0 ∂z

z)0

(1) us*ij - K∇s*ij ) ucbi δij - uciδij + K∇ciδij horizontal inflow

z ) H (top of model domain)

(2) -∇s*ij ) 0

where cbi is the concentration of compound i at the boundary, vig is the dry deposition velocity, Ei is the ground-level emission rate, and H is the height of the model domain. The local sensitivity of a model output to a parameter can be calculated as the partial derivative of the output with respect to the parameter:

∂ci(t) sij(t) ) ∂pj where pj is a model parameter (e.g., rate constant, dry deposition velocity) or input (e.g., initial condition, emissions). The magnitudes of model parameters and inputs can range dramatically in space and time, and, likewise, sensitivity coefficients to different parameters may differ by several orders of magnitude. Thus, it is useful to define and use semi-normalized sensitivity coefficients. Given a model parameter, pj, its variation is defined here as pj(x, t) ) jPj(x, t), where Pj(x, t) is the unperturbed field, which can vary in time and space, and j is a scaling variable with a nominal value of 1. Here, what is termed the semi-normalized sensitivity coefficient, s*ij(t), is calculated by the partial derivative of a species concentration, ci, to the scaling variable of parameter j, j:

s*ij(t) ) Pj

∂ci(t) ∂ci(t) ∂ci(t) ) Pj ) ∂pj ∂j ∂(jPj)

(2)

for non-zero pj. For pj ) 0, s*ij(t) is equal to zero. As formulated, the parameter perturbation is uniform relative to its nominal value over the computational domain. By substituting eq 2 into eq 1, the auxiliary equations of sensitivity coefficients

2860

TABLE 1. Model Parameters Used for Ozone Sensitivity Analysis

9

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

(3) vigs*ij - Kzz

(4) -

horizontal outflow ∂s*ij ∂ci ) - vigciδij + Kzz δ + Eiδij ∂z ∂z ij z)0

∂s*ij )0 ∂z

z)H

where J is the Jacobian matrix defined by Jik ) ∂Ri/∂ck. δij is a binary variable, and depending on whether parameter j is a constituent of the term, the δij on the right-hand side will be either one or zero. If parameter j exists in the term, then δij is one; otherwise δij will be zero. For example, for the inflow boundary condition, if parameter j refers to either the wind speed or the boundary condition of ci, the term ucbi δij will be ucbi with δij ) 1 since both u and cbi are in that term. As can be seen in eq 3, the sensitivity coefficients are functions of the concentrations and are linear in the sensitivity coefficients although coupled. Also, while the sensitivity coefficients do depend on the concentrations, the converse is not the case; so the sensitivities can be integrated separately from the concentrations over the time step. Furthermore, the temporal evolution of the concentrations, which are found before the sensitivity equations are integrated, can be used to stabilize the calculation of sensitivity coefficients and to use larger time steps when integrating eq 3. Numerical Implementation. Like most three-dimensional air quality models, the CIT model uses an operatorsplitting technique to reduce the multidimensional problem to a sequence of one-dimensional equations. For the CIT model, the operator splitting for concentrations ci at t ) t + ∆t is

(∆t2 )L (∆t)L

ci(t + ∆t) ) LH

V

RS(∆t)LH

(∆t2 )c (t) i

(4)

FIGURE 1. NOx emissions (t/day) from (a, top) mobile sources and (b, bottom) area sources. where

LH(ci) ) -∇H(uci) + ∇H(K∇Hci)

port and a separate technique for chemical production with diffusion-dominated vertical transport (17). For the sensitivity coefficients of eq 3, a similar operatorsplitting sequence is used:

LV(ci) + -∇V(uzci) + ∇V(K∇Vci)

(∆t2 )L (∆t)L

s*ij(t + ∆t) ) LH

LRS(ci) ) Ri + Si LH and LV are the operators applied to horizontal and vertical transport, respectively, and LRS is an operator that combines the chemistry and source processes. The horizontal transport operator is further decomposed into operators in the x and y directions. However, the time scale for the vertical diffusion transport is very short during the day in the mixed layer and is comparable to the chemical reaction time scales of many compounds. Moreover, the solution of a diffusion-dominated process has an exponential structure similar to chemical decay. Thus, the vertical diffusion is implemented together with the chemistry and source process in the LRS operator. This structure allows the application of very fast and accurate solution techniques to solve the advection-dominated trans-

V

p RS

(∆t)LH

(∆t2 )s*(t)

(5)

ij

where

LH(s*ij) ) -∇H(us*ij) + ∇H(K∇Hs*ij) LV(s*ij) ) -∇V(uzs*ij) + ∇V(K∇Vs*ij)

(

p LRS (s*ij) ) Jiks*kj +

)

∂Ri ∂Si + - ∇(uci)δij + ∇(K∇ci)δij ∂j ∂j

∇H and ∇V denote the del operator applied to horizontal and p vertical transport, respectively; LRS is the production term arising from the sensitivity to the chemistry, emission sources, and transport. As seen, the difference in solution procedure


9

2861

FIGURE 2. Distributed emissions of volatile organic compounds (t/day). Dashed line indicates the demarcation between the western and eastern domains.

FIGURE 3. Spatial ozone distribution (ppb) at 1500 h, August 29, 1987. between the concentrations and sensitivity coefficients is the calculation of the sensitivity production terms in the operator p . Thus, the same algorithm routines can be used for the LRS transport process of the sensitivities, and the same general model structure is maintained. The last three terms in the p operator are independent of the sensitivity coefficients LRS and can be treated as source/sink terms. In the solution procedure, the sensitivity coefficients are calculated by alternating the operator solution of eq 5 with the operator solution of eq 4 in sequence, i.e., in a given time step, each concentration operator is solved and followed by the corresponding sensitivity operator. This procedure allows the use of the same time step for both concentrations and sensitivity coefficients. In this manner, calculating the production term -∇H(uci)δij + ∇H(K∇Hci)δij of sensitivity coefficients is equivalent to computing the concentration differences resulting from applying the operator LH(ci):

2862

9


-∇H(uci)δij +∇H(K∇Hci)δij ) (cn+1 - cni )/dt i

(6a)

represent the concentration of species i where cni and cn+1 i before and after solving the horizontal transport operator. However, the effect of the diffusion process on model predictions is relatively small, and hence the production resulting from the sensitivity due to diffusivity is negligible in magnitude, especially in the horizontal transport. Thus, eq 6a can be approximated as

-∇H(uci)δij ≈ (cn+1 - cni )/dt i

(6b)

i.e., the concentration change of applying LH(ci). For assessing the sensitivities of the concentrations to the wind velocities, the production term in eq 6b can be added in a computationally effective manner by recognizing that the corresponding production term of sensitivity (s*ijn+1 - s*ijn), is equivalent to the concentration change or

FIGURE 4. Ozone sensitivity (ppb of O3) at 1500 h, August 29, 1987, to the (a, top) O3 initial condition, (b, second from the top) O3 dry deposition velocity, (c, second from the bottom) HO + NO2 rate constant, and (d, bottom) HO + NO2 rate constant by the brute force method.


9

2863

FIGURE 5. Time series plots of ozone sensitivity (ppb O3) to the O3 initial condition (solid line) and HO + NO2 rate constant (dashed line) shown at Claremont.

s*ijn+1 ) s*ijn + (cn+1 - cni ) i

(7)

Similarly, the production resulting from the sensitivity due to vertical wind speed can be implemented in the sensitivity p equations using eq 7 before solving for the LRS (s*ij). The most intensive computational aspect of atmospheric transport chemistry problems is the task of solving the differential equations of chemical transformation. Unlike the equations governing the evolution of species concentrations, in which the integration is dominated by the nonlinear chemistry, the linear system of sensitivity equations with respect to chemistry can be solved very efficiently by LU decomposition using larger time steps than used during chemical integration. In particular, the time step used for integration of the LRS step is twice the advection time step (which is determined by the Courant-Friedrichs-Lewy limit (15)). The linear algebraic system of sensitivity equations associated with the chemistry is represented by

(

)

n+1

∂R(c s* n+1 - s* n s* n+1 + s* n + ) ˜Jn+1 ∆t 2 ∂

)

(8)

where

˜Jn+1 ) J(c˜)

and

c˜ )

cn+1 + cn 2

The solution of sensitivity coefficients in eq 8 at t ) t + ∆t is

(

s* n+1 ) I -

∆t n+1 ˜J 2

) {( -1

I+

)

n+1

∆t n+1 *n ∂R(c s + ˜J 2 ∂

}

)

(9)

Of note, the LU factorization of (I - (∆t/2)J˜n+1) need only be carried out once for all sensitivity coefficients calculated in each computational cell because both the Jacobian and the local parametric derivation matrix (i.e., ∂R/∂) are independent of sensitivity coefficients, only being functions of the concentrations and rate constants (see eq 3) (Note: since the hybrid solver is used to find cn+1 , the LU decomposition is i required here for the sensitivity analysis. If a method is chosen to find cn+1 that requires finding Jn+1, it may not be necessary i to do the LU decomposition again). Because a large number of compounds are involved in the chemical mechanism used in this study, the sparsity of the Jacobian matrix has been

2864

9


considered for more efficient factorization. Here, a numerical package for sparse matrices from the Harwell Subroutine Library (HSL) (18) was adopted for the LU decomposition in eq 9. Another feature of solving sensitivities in this study is the availability of an automatic code-generation program which can write the FORTRAN code for calculating the reaction rates, the Jacobian matrix, and ∂R/∂ matrix required by the integration. The program has been enhanced from its original design (19) for the DDM to provide a suitable format of ∂R/∂ for the absolute and semi-normalized sensitivity calculations, and the sensitivities to the stoichiometric coefficients of reactions. This enables new chemical mechanisms or updated reactions to be implemented without the need for manual coding.

Application to Southern California Ozone formation in an airshed usually is a multiday event associated with transport and chemical transformation. In this study, the sensitivity analysis technique developed in the previous section is applied to the August 27-29, 1987, episode, an intensive monitoring period of the Southern California Air Quality Study (SCAQS). The Los Angeles basin was selected because of the severity of its ozone problem; therefore, finding an effective ozone control strategy is critical for this area. While the results presented in the next section are only for surface level ozone sensitivities to a set of selected parameters, the sensitivity coefficients of other compounds, including aloft concentrations, are computed at the same time. The parameters for which sensitivity coefficients were computed are summarized in Table 1. Of these parameters, the rate constant of HO + NO2 f HNO3 was chosen because predicted ozone concentrations were found to be significantly sensitive to its variation (10). In addition, recent work suggests that the rate constant which has been used in past modeling studies may be overestimated by up to 40% (20). NOx emissions were split into two classes: mobile (M) and area (A) sources, which together account for 92% of the total NOx in the emission inventory during the episode. The spatial distribution of the two NOx emission classes are shown in Figure 1. For the ground level VOC emissions, the modeling region was divided into two subdomains: the western region (W), which includes most of urban Los Angeles, and the eastern

FIGURE 6. Ozone sensitivity (ppb of O3) at 1500 h, August 29, 1987, to the (a, top) NOx(M) emissions and (b, bottom) NOx(A) emissions. region (E), which covers a large part of the downwind area. Figure 2 shows the distribution of the daily total VOC emissions over the model region with the two subdomains demarcated by a dashed line. In order to calculate sensitivity coefficients to these emissions, four emission classessNOx(M), NOx(A), VOC(W), and VOC(E)swere defined in the SAPRC90 chemical mechanism (21). Because each class contains two (for NOx) or more (for VOC) detailed compounds, a transformation reaction was defined to account for the detailed emission species in the mechanism, i.e. ki

NOx (or VOC) 98

∑γ (detailed species i) i

i

where γi is the fraction of detailed species i in the emission class, e.g., NOx f 0.95NO + 0.05NO2. These reactions are defined fast enough so that the resulting dependency on these rate constants is negligible.

Results and Discussion Figure 3 shows the surface ozone field at the time of peak predicted ozone (1500 h) on August 29, 1987, during the multiday episode. As seen in the figure, the high ozone levels occur outside the urban core of Los Angeles, with the predicted peak value of 175 ppb of O3 downwind of the city. Sensitivity analysis results are given as the change in ozone levels (ppb) per 1% increase (domain-wide) in the nominal value of the parameters, and results are shown at the time of predicted peak ozone occurrence on August 29. Figure 4 shows the ozone sensitivities to the ozone initial condition (Figure 4a), the ozone dry deposition velocity (Figure 4b), and the rate constant of HO + NO2 f HNO3 (Figure 4c), respectively. As expected, an increase in initial ozone concentrations over the domain increases predicted ozone throughout the modeling region with higher sensitivities in the downwind area where the maximum ozone sensitivity is about 0.01 ppb of O3 (Figure 4a). By the third day of the simulation, only slight effects over a small portion of the


9

2865

FIGURE 7. Ozone sensitivity (ppb of O3) at 1500 h, August 29, 1987, to the (a, top) VOC(W) emissions and (b, bottom) VOC(E) emissions. domain are predicted. The time series plot of ozone peak sensitivity to ozone initial condition is presented in Figure 5. The plot shows the temporal peak sensitivity evolution at Claremont (location shown in Figure 3), one of the intensive monitoring sites for detailed VOCs during SCAQS. As shown in Figure 5, starting from its maxima, the peak sensitivity decreases dramatically before photochemical reactions can be driven by the sunlight during the first day simulation. Aftersun rise, the peak sensitivity increases to its diurnal maxima at about noon then decreases. A similar diurnal pattern also is shown on the second day, but with a much lower amplitude than on the first day. Finally, the magnitude of the peak sensitivity on the third day becomes very small. In Figure 4b, increasing the ozone deposition velocity decreases predicted ozone everywhere because ozone removal is enhanced on the surface boundary. In this case, ozone concentrations downwind have the greatest decrease. Figure 4c presents the spatial distributions of ozone sensitivities due to the reaction rate constant of HO + NO2 f HNO3. This reaction rate is one of the parameters in the chemical

2866

9


mechanism to which ozone is most sensitive. Increasing the rate decreases ozone concentrations in the downwind area because more hydroxyl radical and NO2 in system are removed through the reaction. The predicted maximum seminormalized sensitivity is -0.78 ppb of O3. If, as suggested by Donahue et al. (20), the rate constant is 40% lower than previously assumed, the ozone levels could increase by 30 ppb or more. The time series evolution of the ozone sensitivity to this rate constant at Claremont is presented in Figure 5. As can be seen, there is a strong diurnal cycle showing the peak sensitivity occurring late afternoon when NO2 concentrations becoming relatively low. In contrast to the sensitivity to initial conditions, the sensitivity to the HO + NO2 rate constant does not reduce greatly from day to day. To compare DDM-3D with other methods, a brute-force sensitivity to the HO + NO2 rate constant for the same episode is computed. As shown in Figure 4d, the peak ozone sensitivity is predicted as -1.2 ppb of O3, as compared to -1.4 ppb of O3 from the previous study of box model (10) using a 10-h average scenario of 39 city conditions (22).

FIGURE 8. Ozone sensitivity (ppb of O3) at 1500 h, August, 1987, to the wind speed. Figure 6 shows ozone sensitivities to the NOx emissions from mobile and area sources. In these cases, increasing NOx emissions has the effect of decreasing predicted ozone concentrations immediately downwind of the urban core (0100 km). However, ozone concentrations tend to increase much further downwind of the urban core (g100 km). The similar spatial distributions of ozone sensitivity to the two sources are attributed to the local abundance of VOC and NOx in the two regions (immediately downwind and much further downwind). In urban Los Angeles, the NOx is relatively abundant and ozone production is VOC-limited. On the other hand, downwind of the urban core, the NOx is less abundant. The maximum ozone sensitivities found for a 1% increase in mobile and area source NOx emissions are about -0.9 and -0.4 ppb of O3, respectively. Figure 7 shows the ozone sensitivities to increasing VOCs in the western and eastern subdomains. In both cases, ozone levels tend to increase downwind of the emissions. An increase in VOC emissions results in more peroxy radical generation through the VOC oxidation to convert NO to NO2 and, consequently, increases ozone formation. The area of ozone increase due to VOC(W) is much larger than that due to VOC(E) because the total VOC emissions are dominated by the western subdomain (see Table 1 and Figure 2), approximately accounting for 83% of total VOCs. The maximum ozone sensitivities for both cases are about 1.4 and 0.3 ppb O3, respectively. Figure 8 shows the sensitivity of ozone concentrations to the wind speed over the domain. In this case, increasing the wind speed potentially increases dilution, though this leads to both positive and negative sensitivities. The effect on the ozone levels at the inflow boundary is to lower concentrations (a negative sensitivity) because the air mass carried over the boundary is relatively clean. In the urban core, however, the dilution decreases the ambient NOx and ozone levels tend to increase as a result of decreased NOx scavenging. Somewhat downwind, the negative sensitivity is caused by the net effect of dilution of VOC and NOx (see Figures 6a and 7a). Diluting the VOC has a stronger impact than diluting the NOx emissions, which would increase ozone. On the other hand, increasing the wind speed will decrease the ozone removal by the dry deposition process; consequently, there is a positive sensitivity in the far downwind area. One of the attractive features of this method is the computational efficiency. Table 2 provides a comparison of

TABLE 2. Relative Execution Time for Sensitivity Analysisa concentrations alone (base case simulation) sensitivity coefficients to one parametersb,c sensitivity coefficients to 10 parametersb,d sensitivity coefficients to 20 parametersb,e

1.0 1.30 1.52 1.81

a A set of sensitivity coefficients represent all compound sensitivities to a given parameter or input. b Includes time needed to calculate concentrations. c Ozone initial concentration. d Five initial conditions and five rate constants. e Five initial conditions and 15 rate constants.

execution time for four different calculations where the number of parameters chosen for sensitivity analysis were increased from 0 to 20. All comparisons were conducted on a SUN Ultra 1-170 workstation with 64 Mbyte of RAM. A one parameter sensitivity analysis increases execution time 30% relative to concentration calculations alone. For the cases where 10 and 20 parameters are calculated, the execution time increases about 50% and 80%, respectively. Recall that sensitivity coefficients are independent of the nonlinear coupling in the reactions and the factorization of the matrix (I - (∆t/2)J˜n+1) needs to be done once when all the sensitivity coefficients are advanced from tn to tn+1. Therefore, solving additional sensitivity coefficients requires relatively less CPU time as compared to solving additional species concentrations. As shown in Table 2, the average execution time for each additional set of sensitivity coefficients decreases as the number of coefficients increases. As described, the sensitivity equations directly derived from the model equations have a similar structure to that of the concentration equations, and the general structure of model is maintained while calculating the sensitivity coefficients. Therefore, the effort to solve the sensitivity equations separately only requires limited programing modifications from the original model code. The results of an application presented here indicate that the new method is viable and can efficiently compute the sensitivity coefficients of all species simultaneously. Moreover, the implementation of DDM-3D presented here can also be conceptually applied to other multidimensional air quality models.

Acknowledgments This work was supported, in part, by the U.S. Environmental Protection Agency, Office of Exploratory Research under Grant


9

2867

R819356, the National Science Foundation under Grant ASC9217365, and the Georgia Power Company.

Literature Cited (1) Koda, M.; Dogru, A. H.; Seinfeld, J. H. J. Comput. Phys. 1979, 30, 259. (2) Rabitz, H.; Kraman, M. A.; Dacol, D. Annu. Rev. Phys. Chem. 1983, 34, 419. (3) Milford, J. B.; Gao, D.; Russell, A. G.; McRae, G. J. Environ. Sci. Technol. 1992, 26, 1179. (4) Cho, S.-Y.; Carmichael, G. R.; Rabitz, H. Atmos. Environ. 1987, 12, 2589. (5) Seigneur, S.; Teche, T. W.; Roth, P. M.; Reid, L. E. J. Appl. Meteorol. 1981, 20, 1020. (6) Dechaux, J. C.; Zimmerman, V.; Nollet, V. Atmos. Environ. 1994, 28, 195. (7) Dickerson, R. R.; Stedman, D. H.; Delany, A. C. J. Geophys. Res. 1982, 87, 4933. (8) Dunker, A. M. J. Chem. Phys. 1984, 81, 2385. (9) Dougherty, E. P.; Hwang, J. T.; Rabitz, H. J. Chem. Phys. 1979, 71, 1794. (10) Yang, Y.-J.; Stockwell, W. R.; Milford, J. B. Environ. Sci. Technol. 1995, 29, 1336. (11) Gao, D.; Stockwell, W. R.; Milford, J. B. J. Geophys. Res. 1995, 100, 23153. (12) Dunker, A. M. Atmos. Environ. 1981, 15, 1155. (13) Bischof, C. A.; Carle, A.; Corliss, G.; Griewank, A.; Hovland, P. Sci. Prog. 1992, 1, 11.

2868

9


(14) Hwang, D.; Byun, D. W. Ninth Joint Conference on Applications of Air Pollution Meteorology with A&WMA, Atlanta, GA; American Meterological Society: Boston, MA, 1996; Paper 9.7. (15) Carmichael, G. R.; Sandu, A.; Potra, F. A. Atmos. Environ. 1997, 31, 475. (16) McRae, G. J.; Goodin, W. R.; Seinfeld, J. H. Atmos. Environ. 1982, 16, 679. (17) Russell, A. G.; McCue, K. F.; Cass, G. R. Environ. Sci. Technol. 1988, 22, 263. (18) Harwell Subroutine Library. http://www.rl.ac.uk/departments/ ccd/numerical/hsl/hsl.html, 1996. (19) McCrosky, P. S.; McRae, G. J. Documentation for the Direct Decoupled Sensitivity Analysis MethodsDDM; Technical Report; Department of Chemical Engineering, Carnegie Mellon University: Pittsburgh, PA, 1987. (20) Donahue, N. M.; Dubey, M. K.; Mohrschladt, R.; Demerjian, K. L.; Anderson, J. G. J. Geophys. Res. 1997, 102, 6159. (21) Carter, W. P. L. Atmos. Environ. 1990, 24A, 481. (22) Carter, W. P. L. J. Air Waste Manage. Assoc. 1994, 44, 881.

Received for review February 11, 1997. Revised manuscript received May 15, 1997. Accepted May 20, 1997.X ES970117W X

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

Fast, Direct Sensitivity Analysis of Multidimensional ... - ACS Publications

Recommend Documents