Implementing the Decoupled Direct Method for Sensitivity Analysis in

quality Model with extensions (CAMx) to calculate first-order sensitivities for gas-phase species (e.g., ozone) with respect to initial concentrations...
9 downloads 0 Views 183KB Size
Download PDF
Environ. Sci. Technol. 2007, 41, 2847-2854

Implementing the Decoupled Direct Method for Sensitivity Analysis in a Particulate Matter Air Quality Model B O N Y O U N G K O O , * ,† A L A N M . D U N K E R , ‡ AND GREG YARWOOD† ENVIRON International Corporation, 101 Rowland Way, Novato, California 94945, and General Motors R&D Center, 30500 Mound Road, Warren, Michigan 48090

The decoupled direct method (DDM) is an efficient and accurate way of performing sensitivity analysis to model inputs. As the impact of atmospheric particulate matter (PM) on human health and visibility became evident, the need to extend the DDM to PM sensitivity has grown. In this work, the DDM is implemented in the three PM modules employed in the Comprehensive Air-quality Model with extensions (CAMx): ISORROPIA for inorganic gas/aerosol partitioning, SOAP for secondary organic gas/aerosol partitioning, and RADM-AQ for aqueous-phase chemistry. The PM modules are complex and the DDM implementation is discussed in detail. Stand-alone tests are performed for each PM module in which first-order sensitivities computed by the DDM are compared to the traditional brute-force method (BFM). The DDM sensitivities are shown to be accurate and agree well with the BFM within the linear response range. The SOAP module showed nearly linear response for up to (30% changes in concentration inputs. The RADMAQ module showed moderately nonlinear response in some tests but first-order sensitivities accounted for most of the response for input changes up to (20%. ISORROPIA shows greater deviation from linear response than the other PM modules and the near-linear range can be as restricted as (10% changes in concentration inputs. Nonlinearity in ISORROPIA results both from the equations that describe thermodynamic equilibrium and the computational approaches within ISORROPIA that are employed to boost efficiency.

Introduction Air quality models are powerful tools to understand the relationship between emissions and ambient concentrations of atmospheric pollutants providing valuable information for the establishment of effective control strategies to reduce adverse effects on our environment. To date, the most frequently used technique to study the model responses to various system parameters is the brute-force method (BFM), which involves repeated model simulations with different model inputs. This approach is easy to apply and interpret. However, it is computationally demanding (a full model simulation is required for each model input of interest) and susceptible to numerical uncertainty (model response to a small perturbation is easily contaminated by numerical noise). * Corresponding author phone: (415)899-0700; fax: (415)899-0707; e-mail: [email protected]. † ENVIRON. ‡ General Motors. 10.1021/es0619962 CCC: $37.00 Published on Web 03/20/2007

 2007 American Chemical Society

The decoupled direct method (DDM) offers an alternative to the BFM by directly solving sensitivity equations derived from the governing equations of the model (1, 2). The DDM decouples the sensitivity equations from the model equations, thereby providing stability and computational efficiency. The DDM has been implemented in the Comprehensive Airquality Model with extensions (CAMx) to calculate first-order sensitivities for gas-phase species (e.g., ozone) with respect to initial concentrations, boundary concentrations, and emissions (3, 4). A variant of the DDM called DDM-3D has also been developed and applied to other three-dimensional air quality models (5-7). The main difference between the DDM and DDM-3D is that the DDM closely follows the host model’s chemistry solution algorithm and uses the same time steps as the chemistry routine, while DDM-3D uses a semi-implicit finite difference scheme with different time steps to solve the chemistry sensitivity equations. The latter can help to improve the computational efficiency, but may introduce inconsistencies between the sensitivities and concentrations (3). Generally the DDM has been used for first-order sensitivity analysis for ozone. There have been a few attempts to extend the capability of the DDM. Hakami et al. (8, 9) and Cohan et al. (10) extended DDM-3D to calculate higher-order sensitivities of ozone concentrations. Recently, Napelenok et al. (11) published the implementation of DDM-3D in the PM modules of the Community Multiscale Air Quality model (12). We will discuss their implementation later. In this paper, the DDM is implemented in the three PM chemistry modules of the CAMx. Detailed solution algorithms for each module and the derivation of the sensitivity equations are described. The DDM performance is then compared with the BFM in stand-alone tests for each module. Full 3-D tests of CAMx with DDM implemented for the PM modules also are presented.

Model Description There are three main PM chemistry modules in version 4.20 of CAMx (13). The Regional Acid Deposition Model-Aqueous Chemistry (RADM-AQ) module performs aqueous-phase SO2 oxidation reactions in cloudwater using the algorithm developed for the RADM (14). Partitioning of condensable organic gases (CGs) to secondary organic aerosols (SOAs) to form a condensed “organic solution phase” is performed using a semi-volatile, equilibrium scheme called Secondary Organic Aerosol Partitioning (SOAP) (15). Gas/aerosol partitioning for inorganic aerosol constituents (sulfate, nitrate, ammonium, sodium, and chloride) is performed by the ISORROPIA (Greek for “equilibrium”) thermodynamic model (16). In applying the DDM to the PM modules, we are primarily concerned with the first-order sensitivity coefficients (S), which represent the change in concentration (Y) with respect to some input parameter (p), evaluated relative to the base case (p ) p0):

S)

|

|

|

|

∂Y ∂Y ∂Y ∂Y ) + + (1) ∂p p0 ∂p p0;SOAP ∂p p0;RADM-AQ ∂p p0;ISORROPIA

The right side of eq 1 is shorthand to indicate that the sensitivity is propagated first through SOAP, then through RADM-AQ, and last through ISORROPIA, analogous to the operator splitting method used for the PM concentrations. In this implementation, the sensitivities with respect to meteorological inputs such as temperature are not derived. VOL. 41, NO. 8, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

2847

TABLE 1. Species Involved in the PM Modules of CAMx ISORROPIA gases: NH3, HNO3, HCl ions: H+, Na+, NH4+, SO42-, HSO4-, NO3-, Cl-, OHsolid salts: Na2SO4, NaHSO4, NaNO3, NaCl, (NH4)2SO4, NH4HSO4, NH4NO3, NH4Cl, (NH4)3H(SO4)2 aerosol water: H2O a

RADM-AQ

SOAP

gases: SO2, H2O2, O3, CO2, HNO3, NH3, HCOOH, MHPa, PAAa dissolved gases: SO2(aq), H2O2(aq), O3(aq), CO2(aq), HNO3(aq), NH3(aq), HCOOH(aq), MHP(aq), PAA(aq) ions: H+, OH-, HSO3-, SO32-, HCO3-, CO32-, NO3-, NH4+, HCOO-, HSO4-, SO42-

gases: CG1, CG2, CG3, CG4, CG5 aerosols: SOA1, SOA2, SOA3, SOA4, SOA5

MHP and PAA represent methylhydroxyperoxide and peroxyacetic acid, respectively.

We describe the DDM implementation for each PM module in detail below. Implementing the DDM in SOAP. Volatile organic compounds react with atmospheric oxidants to produce condensable gases (CGs) that can condense to secondary organic aerosols (SOAs). CAMx version 4.20 includes 5 pairs of CG/SOA species from various precursors (aromatics, monoterpenes, etc.). SOAP determines the equilibrium partitioning between the gas and aerosol phases assuming that the organic mixture forms an ideal solution (15):

Y0(tSOAi) ) Y1(tSOAi)

(2a)

* Y1(CGi) ) xSOAicCGi

(2b)

where i ) 1,...5; Y0(A) and Y1(A) are initial and equilibrium mass concentrations of species A, respectively; xA is mole fraction of A in the aerosol phase; c*A is the effective saturation concentration of A; tA represents species A in both gas and aerosol phases (e.g., tSOA1 ) CG1 + SOA1). SOAP transforms eq 2 to the following form

∑ i

[

Y0(tSOAi)MOM

* MWSOAiMOM + cCGi

]

+

Y0(POA) MWPOA

- MOM ) 0 (3)

where MOM is the molar concentration of total organic matter in the particle phase, MWA is the molecular weight of species A, and POA represents primary organic aerosol, which is assumed to be entirely in the particle phase. Equation 3 is then solved for MOM using a bi-section method. The equilibrium concentrations of CG/SOA species are then determined as follows:

Y1(SOAi) )

Y0(tSOAi)MWSOAiMOM * MWSOAiMOM + cCGi

Y1(CGi) ) Y0(tSOAi) - Y1(SOAi)

(4a) (4b)

Differentiating eqs 3 and 4 with respect to p gives the sensitivity equations for the SOAP module:

∑ i

[

S0(tSOAi)MOM + Y0(tSOAi)

[

∂p

* MWSOAiMOM + cCGi

∑ i

]

∂MOM

Y0(tSOAi)MOMMWSOAi

-

∂p

* cCGi ]2

[MWSOAiMOM + S0(POA) MWPOA

2848

9

]

∂MOM

-

+

∂MOM ∂p

) 0 (5)

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 8, 2007

S1(SOAi) )

S0(tSOAi)MWSOAiMOM * MWSOAiMOM + cCGi

+

∂MOM ∂p (6a) * [MWSOAiMOM + cCGi ]2

* Y0(tSOAi)cCGi MWSOAi

S1(CGi) ) S0(tSOAi) - S1(SOAi)

(6b)

Implementing the DDM in RADM-AQ. RADM-AQ has a total of 29 internal species (Table 1) that are defined by 29 equations (Table S1 in the Supporting Information [SI]): 18 dissolution/dissociation equations, 10 mass balance equations, and a charge balance equation. The solution procedure of RADM-AQ consists of two steps. First, it solves equilibrium concentrations of the 29 internal species (Table 1) using 29 equations (Table S1). Then, based on the equilibrium ion concentrations, it advances the aqueous-phase oxidation reactions (Table S1) in time, which adjusts the concentrations of S(IV), S(VI), and the oxidants. Multiple time steps may be taken during the simulation and each step is iterated to convergence. Note that the equilibrium equations are expressed with molarities and partial pressures which need to be converted to molar concentration unit: [A] ) Y1(A)/wL; PA ) Y1(A)RT. The liquid water mixing ratio in the cloud droplet, wL (volume water/volume air), is assumed to be independent of changes of the solutes. In RADM-AQ, the mean activity coefficients of electrolytes are calculated by the Davies equation (17). Differentiating all the equations of the equilibrium, mass balance, charge balance, and oxidation reactions with respect to p gives a set of sensitivity equations that are solved by the same iteration scheme as used to solve for the concentrations (Figure 1). We use the same time step to solve for the sensitivities and the concentrations, which keeps the solutions for the concentrations and the sensitivities consistent. A similar approach was used to implement DDM for the gas-phase chemistry (3). Implementing the DDM in ISORROPIA. ISORROPIA calculates the gas/aerosol partitioning for inorganic species by assuming that thermodynamic equilibrium exists between the gas and aerosol phases. Inputs to ISORROPIA are total (gas and aerosol) concentrations for five inorganic species (sodium, sulfate, ammonium, nitrate, and chloride) plus ambient temperature and relative humidity (RH). Internally, ISORROPIA calculates equilibrium concentrations for 21 species (Table 1) although CAMx uses only 4 of these (ammonia, nitric acid, hydrogen chloride, and aerosol water) to describe the inorganic gas/aerosol system in a massconserving way (e.g., aerosol nitrate equals the difference between total nitrate and gaseous nitric acid). To determine the sensitivities of the CAMx species, however, we must calculate the sensitivities of all the ISORROPIA species because they are coupled together via equilibrium reactions. To define the concentrations of the 21 species, there must

3 2 2 2) γH , γH-HSO ) γH +γSO 2SO4 4 4

γH γHSO4 , C ) +

FIGURE 1. Diagram for the solution procedure of the DDM and RADM-AQ. be 21 equations relating the concentrations of the species. The 14 gas-aerosol phase equilibria considered by ISORROPIA define 14 equilibrium equations. Five mass balance equations for sulfate, nitrate, chloride, ammonia, and sodium and the electro-neutrality equation make 6 more equations. Last, there is the Zdanovskii-Stokes-Robinson (ZSR) relationship for the aerosol water concentration (18). Table S2 lists all 21 equations for ISORROPIA. Differentiating them with respect to the input parameter p gives 21 equations for sensitivities. The equilibrium equations are similar to those of RADM-AQ except that [A] represents molality here: [A] ) Y1(A)/W1(H2O) where W1(H2O) is aerosol water content in kg/m3 of air. For example, differentiating the first equation in Table S2 with respect to p gives the following sensitivity equation (note that equilibrium constants do not depend on species concentrations):

Y1(SO24 ) Y1(HSO4 )W1(H2O)

ΓS1(H+) + Y1(H+)

Y1(HSO4 )W1(H2O) Y1(H+)Y1(SO24 ) 2 W1(H2O)Y1(HSO4)

ΓS1(SO24 ) ΓS1(HSO4) -

Y1(H+)Y1(SO24 )

2 Y1(HSO4 )W1(H2O)

3

2

2 γH 2SO4

J

2 γH-HSO 4 3 γH 2SO4

3 γH-HSO 4

C

∂γH2SO4

∑ ∂Y (A ) S (A ) 1

j)1 J

C

ΓS1(H2O) +

∑ j)1

1

∂γH-HSO4 ∂Y1(Aj)

S1(Aj) )

S1(H2O) ) ∂K1 ∂p

∑ i

) 0 (7)

where J is the number of ISORROPIA species that are involved in calculating the activity coefficients,

Y1(HSO4 )W1(H2O)

, and Γ )

3 γH 2SO4 2 γH-HSO 4

In ISORROPIA, the binary and multicomponent activity coefficients are determined using the Kusik-Meissner method (19) and Bromley’s formula (20), respectively. To boost efficiency, ISORROPIA uses pre-calculated tables of binary activity coefficients for each salt over ranges of temperatures and ionic strength. Similar tables could be prepared for the sensitivities of binary activity coefficients. The computational advantages of such a table would be small because the sensitivity equations are linear and can be solved noniteratively. In this work, the sensitivities of the activity coefficients are calculated by direct differentiation of the Kusik-Meissner and Bromley formulas to simplify coding and save memory (see SI for the equations of the activity coefficients and the sensitivity). This may cause a discrepancy in the sensitivities calculated by the DDM and BFM. However, it would be much less significant compared to that caused by the discontinuous response surface of ISORROPIA which will be discussed later. Nine of the equilibrium reactions represent the dissolution or precipitation of solids. If the solid involved in a particular reaction is not present at equilibrium, then we do not include the sensitivity equation derived from the equilibrium reaction. That is, the total number of sensitivity equations to solve is reduced by one for every solid that is absent at equilibrium. Eliminating the sensitivity equation also eliminates one unknown sensitivity because each solid appears in only one reaction. Thus, the number of equations remains equal to the number of unknown sensitivities regardless of how many solids are present at equilibrium. The presence or absence of solids can be diagnosed from the ISORROPIA output. If a solid species is absent at equilibrium, then the sensitivity of the species is set to zero. For a sufficiently large perturbation to an input parameter, the solid species may form. However, the sensitivity coefficient, which is identically zero, cannot predict the formation of the solid when the parameter is varied. This is a limitation of local sensitivity methods and other methods that use information from only one atmospheric condition. The consequences of this limitation may be mitigated by the fact that concentrations of the CAMx aerosol species are calculated as the difference between the total (gas and aerosol) and equilibrium gas concentrations of the species. Then, the sensitivity of the aerosol species is the difference between the sensitivities of the total concentration and equilibrium gas concentration, which likely will not be highly dependent on whether the aerosol species is in liquid phase or in a solid. The ZSR relationship is used to calculate the liquid water content of the aerosols (Table S2). ISORROPIA considers 10 electrolytes (NaCl, Na2SO4, NaNO3, (NH4)2SO4, NH4NO3, NH4Cl, NH4HSO4, NaHSO4, (NH4)3H(SO4)2, and H2SO4) whose concentrations are determined from the equilibrium concentrations of ions and solid salts in different subcases (see the SI for details). Differentiating the ZSR equation with respect to p yields the last sensitivity equation:

j

j

Y1(H+)Y1(SO24 )

-

[

J

1

m0i

∂Ei

∑ ∂Y (A ) S (A ) (a ) 1

w

j)1

1

j

j

]

(8)

In equilibrium, the water activity, aw, is equal to the ambient RH (21). If only solid salts are present, the equilibrium reactions in Table S2 (except for the two solid-to-gas reactions) do not apply and neither does the ZSR equation. Furthermore, the chemical formulas for the solids force VOL. 41, NO. 8, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

2849

TABLE 2. Input Conditions for the Test Cases SOAP

RADM-AQ

ISORROPIA

temp ) 298 K press ) 1 atm CG1 ) 0.15 ppb CG2 ) 0.55 ppb CG3 ) 0.02 ppb CG4 ) 0.50 ppb CG5 ) 0.05 ppb SOA1 ) 0 SOA2 ) 0 SOA3 ) 0 SOA4 ) 0 SOA5 ) 0 POA ) 5 µg/m3

temp ) 298 K press ) 1 atm cloudwater ) 0.2 g/m3 reaction time ) 6 min SO2 ) 10 ppb H2O2 ) 1 ppb O3 ) 50 ppb N xO y ) 0 tNO3 ) 2 ppb tSO4 ) 1 ppb tNH4 ) 5 ppb (case A) tNH4 ) 2 ppb (case B) CO2 ) 330 ppm HCOOH ) 1 ppt MHP ) 1 ppt PAA ) 1 ppt NaCla ) 0.05 µg/m3 KCla ) 0 CaCO3a ) 0 MgCO3a ) 0 Fe(III)a ) 0.01 µg/m3 Mn(II)a ) 0.005 µg/m3

temp ) 298 K RH ) 25∼95% sulfate-poor tNa ) 0.001 µmol/m3 tSO4 ) 0.1 µmol/m3 tNH4 ) 0.25 µmol/m3 tNO3 ) 0.2 µmol/m3 tCl ) 0.001 µmol/m3 sulfate-rich, no free acid tNa ) 0.001 µmol/m3 tSO4 ) 0.1 µmol/m3 tNH4 ) 0.17 µmol/m3 tNO3 ) 0.2 µmol/m3 tCl ) 0.001 µmol/m3 sulfate-rich, free acid tNa ) 0.001 µmol/m3 tSO4 ) 0.2 µmol/m3 tNH4 ) 0.17 µmol/m3 tNO3 ) 0.2 µmol/m3 tCl ) 0.001 µmol/m3

a

Species that remain constant throughout the simulation.

charge balance so that the electroneutrality equation is automatically satisfied. Consequently, ISORROPIA uses only the mass balance equations (and some chemical assumptions) to determine the concentrations of solid salts. There are seven “dry” cases considered in ISORROPIA. These cases

are solved case-by-case using mass balances of the solid species that are possible in the given system. Differentiating the mass balance equations of each dry case with respect to p gives the sensitivities of the solid and gas species. The sensitivities of ionic species are identically zero in these cases.

Results and Discussion The DDM was implemented in each PM module and tested in the stand-alone mode. The input parameter, p, is formulated as follows:

Y0(A)|p ) Y0(A)|p)0 + p Y0(A)|p)0

(9)

With this formulation, p is a dimensionless parameter with base value of zero. Then, the first-order response of a model species with respect to the parameter is predicted as follows:

Y1(A)|p ) Y1(A)|p)0 + p S1(A)

(10)

The model responses predicted by the DDM were compared with those obtained by the BFM (varying p and re-running the model). Testing the DDM in SOAP. The input conditions of the test case are listed in Table 2. The input concentrations chosen for the test case are loosely based on a heavily polluted urban area where high levels of biogenic SOA precursors exist. Figure 2 shows the model responses to the changes in CG1. The model responses are almost perfectly linear within the tested

FIGURE 2. Model responses of SOA species to the change of CG1 concentration predicted by the BF and DDM; the concentrations are in µg/m3 of air. 2850

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 8, 2007

FIGURE 3. Model responses of the test cases A and B to the change of SO2 concentration predicted by the BF and DDM; the concentrations are in ppb except for PSO4 which is in µg/m3 of air. range of the input changes (-30% to +30%) and the DDM agrees almost exactly with the BFM. The model responses to changes in POA show slightly less linear results and start to deviate from linearity at around (30% input changes (see Figure S1 of the SI). The agreement between the DDM and BFM for the other two precursors with substantial initial concentrations, CG2 and CG4, is similar to the agreement for CG1 and POA, respectively. Testing the DDM in RADM-AQ. For the RADM-AQ tests, because the aqueous-phase reactions are sensitive to pH, two test cases with different pH ranges were prepared (Table 2). The test case is again loosely based on an urban area with moderate to high level of pollutants. For non-CAMx species such as CO2, HCOOH, MHP, PAA, NaCl (if Na and Cl are not simulated in CAMx), KCl, CaCO3, MgCO3, Fe(III), and

Mn(II), CAMx default background concentrations were assigned. Figure S2 shows the hydrogen ion concentration during the simulation for the two test cases. Among the many outputs that RADM-AQ provides, only SO2, H2O2, O3, H2SO4, and particulate sulfate (PSO4) are transported between grid cells by CAMx. Moreover, all H2SO4 is assumed to exist in the particle phase, thus its sensitivity coefficient is always zero. Therefore, the test results are shown only for SO2, H2O2, O3, and PSO4. Again, up to (30% input changes were applied for initial SO2, H2O2, O3, and total sulfate (H2SO4 and PSO4) concentrations. Figure 3 gives results for changes in the initial SO2 concentration. The model responses are linear or nearly linear and are closely predicted by the DDM. The agreement between the DDM and BFM for changes in the initial H2O2 or O3 concentrations (Figures S3 and S4) is similar to that in VOL. 41, NO. 8, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

2851

FIGURE 4. Sensitivity coefficients of a sulfate-poor system to the change of total SO4 (a), total NH4 (b), and total NO3 (c) predicted by BF and DDM; the unit is µmol/m3 except for water which is in mg/m3; the solid and dotted arrows indicate RH ) 65% and 70%, respectively. Figure 3. Changes to the initial total sulfate concentration give more nonlinear model responses for test case A (Figure S5). Nevertheless, for changes in this initial concentration up to (20%, the DDM-predicted response is very close to the actual model response. Testing the DDM in ISORROPIA. Unlike SOAP and RADM-AQ where mostly the same set of equations is solved for all input conditions, ISORROPIA uses different solution algorithms for various subcases. Therefore, ISORROPIA creates a discontinuous response surface while those of SOAP and RADM-AQ are continuous, thus more comprehensive tests are required for ISORROPIA. As shown in Table 2, three types of aerosols are selected and tested with RH ranging from 25% to 95% (with increments of 5%). The input concentrations, although based on a condition similar to that used for the above test cases, were specifically chosen to represent the different states of sulfate neutralization as distinguished by ISORROPIA. Figures 4 and 5 present comparisons of the BFM and DDM sensitivities to the change of total sulfate, total ammonium, and total nitrate inputs for the sulfate-poor and sulfate-rich/no-free-acid cases. The BF first-order sensitivity, SBF is estimated as follows:

S1BF(A) )

Y1(A)|p)0.1 - Y1(A)|p)-0.1 0.2

(11)

Computing SBF with the two-sided perturbation shown in eq 11 prevents second-order, but not third-order, sensitivity from biasing SBF. In most cases, the DDM agrees well with 2852

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 8, 2007

the BFM. For the sulfate-poor system, there are two cases where the DDM and BF disagree (indicated by solid and dotted arrows in Figure 4). In both cases, the systems at p ) -0.1 and p ) 0.1 fall into different solution regimes (subcases) in ISORROPIA (Table S3). Therefore, nonlinear response to the input change is not surprising, and the sensitivity coefficient predicted by the DDM at p ) 0 may not agree well with that estimated by the BFM using eq 11. These nonlinear BF responses from ISORROPIA could result either from the governing equations (i.e., true nonlinearity) or from computational assumptions (e.g., the use of subcases in the numerical solution that have a discontinuity in the concentrations or derivatives of the concentrations at the interface between subcases). For the sulfate-rich/no-free-acid system, the sensitivity coefficients of nitric acid predicted by the DDM and BFM differ at RH ) 95% (indicated by solid arrows in Figure 5). This is due to several assumptions in ISORROPIA. In the sulfate-rich cases, ISORROPIA assumes that dissolution of nitric acid is minor. Also, ammonia in the gas phase is assumed to be a minor species, and the dissolution of hydrogen chloride is assumed to be minor. Therefore, ISORROPIA first establishes aerosol phase equilibrium without considering ammonia, nitric acid, and hydrogen chloride in the gas phase. Then, the equilibria of ammonia, nitric acid, and hydrogen chloride between gas and liquid phases are calculated assuming that these equilibria do not significantly affect the equilibrium concentrations of other species established in the first step. The DDM equations, however, were derived from the full set of equilibrium reac-

FIGURE 5. Sensitivity coefficients of a sulfate-rich, no-free-acid system to the change of total SO4 (a), total NH4 (b), and total NO3 (c) predicted by BF and DDM; the unit is µmol/m3 except for water which is in mg/m3; the solid arrows indicate RH ) 95%. tions without these assumptions. Normally, the assumptions in ISORROPIA hold and the DDM and BFM agree well, but at high relative humidity, discrepancies may be noticeable as in Figure 5. Although the same assumptions are made by ISORROPIA for the subcase of a sulfate-rich aerosol with free sulfuric acid, the aerosol phase is significantly more acidic. Thus, ammonia in the gas phase and dissolution of nitric acid and hydrogen chloride are negligible. For this subcase, the DDM and BFM indeed agree well (Figure S6). Comparisons of sensitivities computed by the DDM and the BFM found very good agreement in most cases. Differences between the two methods resulted from nonlinear response contaminating the BFM calculations. Overall, the DDM sensitivities describe the ISORROPIA responses well for (20% changes in input concentrations for most conditions, but only for (10% changes for some conditions. Testing the DDM in Full 3-D Model. The DDM implementations for the three PM modules were incorporated into CAMx v4.20. The completed model was then applied for a 2-day test case for June 2002 and the Eastern United States modeling domain (Figure S7). The sensitivities to the domainwide emissions of NOx, VOC (anthropogenic and biogenic), SO2, NH3, and POA were computed by the DDM and compared to those of the BFM using eq 11. Linear regression analysis was performed for 24-hr average sensitivities by the DDM and BFM for all the surface grid cells, and the statistics are compiled in Table 3 (a series of scatter plots for the same sensitivities is shown in Figures S8-S11). The degree of agreement between the two methods is very high with the values of coefficient of determination

TABLE 3. Linear Regression Analysis Statistics for the Agreement between the 24-hr Average Sensitivities by the DDM and BF Methodsa model speciesb statistics

perturbed emission species NOx

2

VOC

SO2

NH3

POA

PSO4

R 0.9931 slope 1.0350 intercept -0.0040

0.9789 0.9906 0.0026

PNO3

R2 0.9860 slope 0.9381 intercept -0.0015

0.9659 0.9559 0.9771 c 0.8970 0.8445 0.9052 c 0.0041 -0.0027 -0.0080 c

PNH4

R2 0.9823 slope 0.9131 intercept -0.0023

0.9727 0.9055 0.0020

SOAd

R2 0.9997 1.0000 slope 0.9910 0.9993 intercept -0.0001 -0.0005

0.9927 0.9729 0.0020

0.9866 c 1.1720 c 0.0035 c

0.9758 0.9834 c 0.9621 0.9626 c 0.0046 -0.0024 c

c c c

c c c

0.9999 1.0059 0.0

a From full 3-D CAMx simulation for June 14, 2002. b Species whose sensitivities are calculated. c Regression statistics are not shown for the species whose sensitivities are negligible. d SOA ) SOA1 + SOA2 + SOA3 + SOA4 + SOA5.

(R 2) higher than 0.95 for all the sensitivities considered. Consistent with the box model tests, sensitivities of SOA (determined by SOAP module) and PSO4 (determined by RADMAQ; ISORROPIA simply puts all sulfuric acid into the aerosol phase) show better agreement with the BFM than those of PNO3 and PNH4 (determined by ISORROPIA and RADMAQ) do. VOL. 41, NO. 8, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

2853

TABLE 4. CPU Time for CAMx Simulation with DDM and Varying Numbers of Sensitivity Parametersa no. of sensitivity parameters

total CPU time (min/simulation-day)

time per parameterb

0c 1 2 4 8

32.3 78.8 93.6 124.3 194.5

1.44 0.95 0.71 0.63

a

b

Based on a Linux PC with a 2.1 GHz processor Normalized to time required for standard CAMx: (Tn - T0)/nT0 where Tn is the total CPU time to calculate the concentrations and sensitivities for n parameters together and T0 is the time for the concentrations only (standard CAMx run). c Standard CAMx run.

It should be noted that a potential source of discrepancy between the DDM and BFM sensitivities in the 3-D test lies in the deposition modules. Current DDM implementation in the CAMx deposition modules assumes that the impact of particle size on the deposition rates is negligible. If a particle absorbs a significant amount of water, however, it will change density and diameter, which in turn affect dry deposition velocity and wet scavenging rate of the particle. Further investigation is under way to assess the significance of this issue. Table 4 presents the computational efficiency of the DDM with different numbers of sensitivities. When sensitivities are calculated together for 8 input parameters, the time per sensitivity is 37% less than computing the concentrations alone. It is difficult to directly compare the results by Napelenok et al. (11) with ours because different host models and test cases were used. That said, our results exhibit comparable to higher level of agreement between the DDM and BFM sensitivities. Of note, our correlations between DDM and BFM sensitivities are substantially higher than theirs for PM species that are indirectly related to the emitted species, e.g., sensitivity of sulfate to NH3 (Table 3). Their DDM implementation in the PM modules is similar to ours with one notable distinction. Napelenok et al. ignored activity coefficients in the equilibrium equations, assuming these are insensitive to local perturbations, while we explicitly differentiated the activity coefficient formulas. Their assumption may not always be justified as activity coefficients can be highly sensitive to the solution composition. First-order sensitivity analysis has limitations and, for example, it has been shown that higher-order sensitivities may be a useful addition to first-order sensitivities for ozone modeling (8-10). Although more complicated, implementing higher-order sensitivities for PM is essentially not much different from that for gas-phase species in that higher-order (e.g., second-order) sensitivities are calculated by differentiating the lower-order (e.g., first-order) sensitivity equations with respect to input parameters. However, the number of higher-order sensitivities becomes large as the number of parameters increases, so higher-order analysis is best applied after screening the parameters for importance. Because the first-order DDM is more computationally efficient than the BFM, it can provide an efficient way to screen emission sensitivities for purposes such as improving model performance and evaluating emission control technologies.

Acknowledgments This work was funded by the Coordinating Research Council, Alpharetta, GA, under Project A-51a.

Supporting Information Available Calculating activity coefficients and concentrations of electrolytes in ISORROPIA, equations of RADM-AQ and ISORROPIA, ISORROPIA subcases for the sulfate-poor system, model responses of SOAP species to the change of initial 2854

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 8, 2007

POA concentration, pH ranges for the RADM-AQ test case A and B, model responses of RADM-AQ species to the change of initial H2O2, O3, and total sulfate concentrations, comparison of DDM and BFM sensitivities for the sulfate-rich, free-acid system of ISORROPIA, modeling domain for the 3-D test case, scatter plots of 24-hr average DDM and BF sensitivities from the 3-D model test. This material is available free of charge via the Internet at http://pubs.acs.org.

Literature Cited (1) Dunker, A. M. Efficient calculation of sensitivity coefficients for complex atmospheric models. Atmos. Environ. 1981, 15, 11551161. (2) Dunker, A. M. The decoupled direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 1984, 81, 2385-2393. (3) Dunker, A. M.; Yarwood, G.; Ortmann, J. P.; Wilson, G. M. The decoupled direct method for sensitivity analysis in a threedimensional air quality model: Implementation, accuracy, and efficiency. Environ. Sci. Technol. 2002, 36, 2965-2976. (4) Dunker, A. M.; Yarwood, G.; Ortmann, J. P.; Wilson, G. M. Comparison of source apportionment and source sensitivity of ozone in a three-dimensional air quality model. Environ. Sci. Technol. 2002, 36, 2953-2964. (5) Yang, Y.-J.; Wilkinson, J. G.; Russell, A. G. Fast, direct sensitivity analysis of multidimensional photochemical models. Environ. Sci. Technol. 1997, 31, 2859-2868. (6) Mendoza-Dominguez, A.; Wilkinson, J. G.; Yang, Y.-J.; Russell, A. G. Modeling and direct sensitivity analysis of biogenic emissions impacts on regional ozone formation in the MexicoU.S. border area. J. Air Waste Manage. Assoc. 2000, 50, 21-31. (7) Hakami, A.; Harley, R. A.; Milford, J. B.; Odman, M. T.; Russell, A. G. Regional, three-dimensional assessment of the ozone formation potential of organic compounds. Atmos. Environ. 2004, 38, 121-134. (8) Hakami, A. M.; Odman, T.; Russell, A. G. High-order, direct sensitivity analysis of multidimensional air quality models. Environ. Sci. Technol. 2003, 37, 2442-2452. (9) Hakami, A. M.; Odman, T.; Russell, A. G. Nonlinearity in atmospheric response: A direct sensitivity analysis approach. J. Geophys. Res. 2004, 109, D15303. (10) Cohan, D. S.; Hakami, A.; Hu, Y.; Russell, A. G. Nonlinear response of ozone to emissions: Source apportionment and sensitivity analysis. Environ. Sci. Technol. 2005, 39, 6739-6748. (11) Napelenok, S. L.; Cohan, D. S.; Hu, Y.; Russell, A. G. Decoupled direct 3D sensitivity analysis for particulate matter (DDM-3D/ PM). Atmos. Environ. 2006, 40, 6112-6121. (12) Byun, D. W.; Ching, J. K. S. Science algorithms of the EPA Models-3 Community Multiscale Air Quality (CMAQ) Modeling System; U.S. EPA/600/R-99/030; U.S. EPA: Research Triangle Park, NC, 1999. (13) ENVIRON. User’s Guide, Comprehensive Air Quality Model with Extensions (CAMx), Version 4.20; ENVIRON International Corporation: Novato, CA, June, 2005; available at http://www. camx.com. (14) Chang, J. S.; Brost, R. A.; Isaksen, I. S. A.; Madronich, S.; Middleton, P.; Stockwell, W. R.; Walcek, C. J. A three-dimensional Eulerian acid deposition model: Physical concepts and formulation. J. Geophys. Res. 1987, 92, 14681-14700. (15) Strader, R. F.; Lurmann, F. W.; Pandis, S. N. Evaluation of secondary organic aerosol formation in winter. Atmos. Environ. 1999, 33, 4849-4863. (16) Nenes, A.; Pandis, S. N.; Pilinis, C. ISORROPIA: A new thermodynamic equilibrium model for multiphase multicomponent inorganic aerosols. Aquat. Geochem. 1998, 4, 123-152. (17) Davies, C. W. Ion Association; Butterworths: Washington, DC, 1962. (18) Stokes, R. H.; Robinson, R. A. Interactions in aqueous nonelectrolyte solutions. I. Solute-solvent equilibria. J. Phys. Chem. 1966, 70, 2126-2130. (19) Kusik, C. L.; Meissner, H. P. Electrolyte activity coefficients in inorganic processing. AIChE Symp. Ser. 1978, 173, 14-20. (20) Bromley, L. A. Thermodynamic properties of strong electrolytes in aqueous solutions. AIChE J. 1973, 19, 313-320. (21) Seinfeld, J. H.; Pandis, S. N. Atmospheric Chemistry and Physics; J. Wiley & Sons: New York, 1998.

Received for review August 18, 2006. Revised manuscript received February 12, 2007. Accepted February 15, 2007. ES0619962