Environ. Sci. Technol. 2006, 40, 2663-2670
Adjoint Sensitivity Analysis for a Three-Dimensional Photochemical Model: Implementation and Method Comparison PHILIP T. MARTIEN* AND ROBERT A. HARLEY Department of Civil and Environmental Engineering, University of California, Berkeley, CA, USA 94720-1710 DAN G. CACUCI Department of Nuclear Engineering, University of California, Berkeley, California 94720-1730 and Institute for Nuclear Technology and Reactor Safety, University of Karlsruhe, 76131 Karlsruhe, Germany
Photochemical air pollution forms when emissions of nitrogen oxides (NOx) and volatile organic compounds (VOC) react in the atmosphere in the presence of sunlight. The goal of applying three-dimensional photochemical air quality models is usually to conduct sensitivity analysis: for example, to predict changes in an ozone response due to changes in NOx and VOC emissions or other model data. Forward sensitivity analysis methods are best suited to investigating sensitivities of many model responses to changes in a few inputs or parameters. Here we develop a continuous adjoint model and demonstrate an adjoint sensitivity analysis procedure that is well-suited to the complementary case of determining sensitivity of a small number of model responses to many parameters. Sensitivities generated using the adjoint method agree with those generated using other methods. Compared to the forward method, the adjoint method had large disk storage requirements but was more efficient in terms of computer processor time for receptor-based investigations focused on a single response at a specified site and time. The adjoint method also generates sensitivity apportionment fields, which reveal when and where model data are important to the target response.
Introduction Near the earth’s surface, air pollutants such as ozone and particulate matter (PM) pose concerns for human health, natural ecosystems, agricultural crop yields, and the preservation of building materials (1). Ozone and secondary PM form by complex nonlinear chemical processes that take place in the atmosphere. As a consequence, these air pollution problems have proved to be difficult to understand and control. Three-dimensional air quality models (AQMs) have been developed to represent the relevant coupled processes of emissions, advection, diffusion, photochemical reactions, and surface deposition (2). Such models require thousands of inputs and model parameters, and it is impossible to know a priori exactly which among the thousands of inputs are most important. Sensitivity analysis addresses the question * Corresponding author e-mail:
[email protected]. 10.1021/es0510257 CCC: $33.50 Published on Web 03/17/2006
2006 American Chemical Society
of how the various model responses, quantities of interest derived from the model output, change with perturbations to the model data, the inputs and parameters that define the model. Sensitivity analysis produces first-order derivatives, or sensitivity coefficients, that can be used to rank the importance of model data to a defined response. A number of studies have applied 3-D photochemical models to rank the reactivity of volatile organic compounds (VOC) (3-5). Other studies have investigated the sensitivity of key responses to perturbations in model parameters, such as reaction rate coefficients, to improve our understanding of the chemical mechanisms in AQMs (6). If sensitivities to all data items are available, then from estimates of the data variances and covariances one can estimate the uncertainty of a given response (7). A complete uncertainty analysis based on the use of sensitivity coefficients has not been made for a threedimensional AQM. Previous work on estimating uncertainty in air quality models has applied simplified models and used statistical methods to generate uncertainties (8). Other studies have conducted an uncertainty analysis examining a limited set of inputs (9, 10) or using surrogates to approximate the propagation of errors across the modeling domain (11). A computationally efficient method for the case of many parameters and relativity few responses is the adjoint sensitivity analysis procedure (ASAP) (7, 12). Like the Green function (13-15), the adjoint function, also called the importance function (16), is independent of perturbations in the model data and perturbations in the resultant concentrations. Unlike the Green function, the adjoint function depends on the specific model response selected, and thus must be recalculated for each target response. The compensation for this dependence on the response, however, is that the adjoint function is an N-vector instead of an N × N matrix like the Green function kernel. This greatly reduces the necessary computer memory and storage requirements of the method. Adjoint models have found wide application in the atmospheric sciences (17), particularly within meteorological models for variational data assimilation (18-21). Cho et al. (22) used the adjoint method for calculating sensitivities in an atmospheric chemistry model. Variational data assimilation has also been used to adjust emissions to minimize differences between predicted concentrations and measurements (23, 24). An adjoint model was applied to a simplified multi-box AQM for generating sensitivities (25). Sirkes and Tziperman (26) discussed the distinction between discrete and continuous formulations of adjoint models and showed that discrete models, if not properly formulated, can display numerical instabilities. Schmidt and Martin (27) used an adjoint method in a three-dimensional model to assess the importance of long-range transport of ozone and its precursors to the Paris region. Hess and Vukic´evic´ (28) used a coarsegrid chemical transport model of the free troposphere and the adjoint to the model to examine the transport of a pollutant plume between Asia and Hawaii. Sandu et al. (29, 30) produced and demonstrated automatic code-generation software to facilitate direct and adjoint sensitivity analysis studies of chemical kinetic systems. Recently, Sandu et al. (31) applied these computational tools in a three-dimensional air quality model to conduct sensitivity analysis and data assimilation for an air pollution study in East Asia. The objective of this research is to implement a continuous adjoint sensitivity method in a 3-D photochemical AQM and compare the resource requirements and products obtained from the adjoint procedure to that of other sensitivity VOL. 40, NO. 8, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2663
methods. We distinguish this work from that of previous studies in that (1) we make direct comparison to other sensitivity methods in a single model, emphasizing the strengths of each; (2) we explicitly demonstrate how sensitivities to many parameters, including emissions, initial and boundary conditions, deposition velocities, and chemical reaction coefficients, may be efficiently calculated using the adjoint function; and (3) we discuss the unique properties of the adjoint function as a sensitivity apportionment tool. The Methods section of this paper briefly summarizes the derivation of the adjoint sensitivity analysis procedure (ASAP) and then outlines the numerical methods developed. An evaluation of the new method, its efficiency and usefulness, is provided in the Results and Discussion section.
Methods Theory. The mixing ratio ci of a trace gas species i can be predicted by solving an equation of mass conservation for a system of N coupled trace gases in the atmosphere:
( )
∂ci ∂(ujci) ∂ci ∂ ) K + + ∂t ∂xj ∂xj jj∂xj ri(c1, ..., cN, p1, ..., pM, x, t) + qi (1) where uj represents components of the mean wind; Kjj represents components of the diagonal turbulent diffusivity tensor; ri represents sources and sinks of species i due to chemical reactions; and qi represents elevated point-source emissions. The reaction term ri introduces nonlinear coupling among the species and includes a number, M, of additional model parameters pl, such as chemical reaction rate coefficients. In this paper, eq 1 is applied for times t ∈(0, T] and for a limited domain such that x1 ∈(0, L1), x2 ∈ (0, L2), and x3 ∈ (0, L3), where x1 and x2 are curvilinear (contravariant) horizontal components of a terrain-following space coordinate x and x3 is its vertical component. In eq 1 and subsequent equations, we use an index summation notation whereby a repeated index represents summation over that index. Equation 1 is subject to both initial and boundary conditions:
ci(x, t ) 0) ) ci0(x); (initial condition)
(
∂ci uj bi ) uj ci - Kjj ∂xj
( ) ∂ci )0 ∂xj
)
T
0
L1
0
L2
0
L3
0
c(x, t, r)‚f(x, t) dx3dx2dx1 dt (3)
If I is the index of a particular pollutant of interest, the response function may be selected such that f ) (0,...,0, fI,0...,0), where fI is a scalar function peaked around a particular location and time of interest and approaches or equals zero otherwise. For example, if the target response is a 1-hour average concentration in a single model grid cell (with dimensions ∆x1, ∆x2, ∆x3) and there are nt time steps in the averaging hour (nt∆t ) 1 hour), then the response function is the product of normalized boxcar functions in time and space:
fI ) ftfx1fx2fx3
(4)
where ft ) (nt∆t)-1 for times within the target hour, fxj ) (∆xj)-1 within the target grid cell, and ft ) fxj ) 0 for all other times and grid cells. It is possible to have more complex responses that depend on multiple pollutants and grid cells. For example, if the target response was the 1-hour average NOx (NO + NO2) concentration, then the vector response function, f, would require two nonzero components: one corresponding to NO and the other to NO2. If the response was total population exposure to a pollutant cI, then the nonzero component I of the response function would be the population distribution. For nonlinear equations, such as eq 1, the derivation of the adjoint equation begins by forming the forward sensitivity equations, also called the linear tangent equations of the model. For details, see refs 7, 31, and 33. The forward sensitivity equations, solved for the change in concentration δc given a change in data δr with respect to a base case r0, take the general form
F(δc, c0, r0) ) δQ(δr, c0, r0)
(5)
where the zero superscript denotes the unperturbed value. For the specific case of eq 1, the forward equations are
( )
∂c0i ∂ δKjj + Dilδpl + δqi (6) ∂xj ∂xj
; (lateral outflow)
j)1, 2
∂ci -K33 ) ei - vdici; (bottom) ∂x3 The term ci0 represents the initial concentrations, bi the lateral boundary concentrations, ei the surface emissions, and vdi the dry deposition velocity. The function ci depends on the independent variables t and x. It also depends on the model inputs and parameters, also called the model data, represented by a vector
(2)
A model response, a scalar quantity of interest derived from the predicted concentrations, may be expressed as the inner 9
∫∫ ∫ ∫
( )
; (lateral inflow)
∂ci ) 0; (top) ∂x3
2664
R(r) ) 〈c(x, t, r), f(x, t)〉 ≡
∂δci ∂(uj0δci) ∂(δujc0i ) ∂ 0 ∂δci + Kjj - Jikδck ) + ∂t ∂xj ∂xj ∂xj ∂xj
j)1, 2
r ) (ci0, bi, ei, qi, u1, u2, u3, K11, K22, K33, vdi, pl)
product of the pollutant vector with the N-vector response function, f.
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 8, 2006
The term Jij represents the N × N Jacobian matrix that contains the rate of change of each element in the chemical reaction vector with respect to each pollutant (Jij ) ∂ri/∂cj); the term Dil represents the N × M matrix of “direct derivatives” that contains the rate of change of each element in the reaction vector with respect to each of the M parameters in r (Dil ) ∂ri/∂pl). Thus, F and δQ (in eq 5) are, respectively, defined by the left-hand side and right-hand side of eq 6. Equation 6 relates perturbations in the model data to the variation in pollutant concentrations. The required initial and boundary conditions are supplied by forming the differential of the model initial and boundary conditions (31, 33). The differential of the desired response (as defined by eq 3) is
δR ) 〈δc, f〉
(7)
To apply the forward method prescribed by eq 7 to calculate the change in response due to a perturbation in data, one must first apply eq 6 to determine δc.
The ASAP proceeds by forming the inner product of eq 5 with an, as yet unspecified, adjoint function ψ:
+
∫∫ ∫ T
L1
0
0
L2
0
[(
∂c0i ψi δei + δK33 - δvdic0i ∂x3
)]
x3)0
(8)
dx2 dx1dt (12)
Because F(δc) is linear in δc, an adjoint operator A(ψ) exists (7) such that
The utility of eq 12 arises from the fact that once one evaluates ψ, which is independent of any perturbations in the model data (δr), it specifies the change in response due to any perturbation in terms of integrals of known quantities. For a given response, R, one can calculate the unnormalized and semi-normalized sensitivities (respectively SI and sI) of R to a change in any particular data element δRI using
〈ψ, F(δc)〉 ) 〈ψ, δQ〉
〈ψ, F(δc)〉 ) 〈δc, A(ψ)〉 + P(ψ, c0, δr, r0)
(9)
where P is a term evaluated at the boundaries of the integration defined by the inner product. The final conditions and boundary conditions on ψ are set such that the boundary term P is independent of δc (33). If we define the adjoint function by requiring that A(ψ) ) f, then the first term on the right-hand side of eq 9 equals δR (see eq 7). We can combine eqs 8 and 9 to relate the adjoint function directly to the response, that is,
δR ) 〈ψ, δQ〉 - P(ψ, c0, δr, r0)
(10)
The specific adjoint operator for the model equation of this study is given by
( )
∂ψi ∂(u0j ψi) ∂ 0 ∂ψi K - JTik ψk ) fi ∂t ∂xj ∂xj jj ∂xj
(11)
Details of the derivation are provided elsewhere (7, 31, 33). Equation 11 coupled with the final and boundary conditions listed below completely specifies ψ.
δR and sI ) R0I SI δRI
SI )
where I is a fixed index (no summation on I) and R0I is an unperturbed reference data value. Often, a given data element is not a real number but is instead a real-valued function RI ) RI(x, t). In this case, one can still calculate sensitivities by considering either multiplicative or additive perturbations to the data. Multiplicative perturbations take the form δRI ) δR0I (x, t), where is a real number. Sensitivities are then formed with respect to changes in the averaged data element. As a specific example, consider the case of the sensitivity of the response to emissions for a uniform scaling of the emissions.
SI )
∫∫ ∫
δ
T
L1
0
0
L2
0
[ψIe0I ]x3)0dx2 dx1dt δe0I
∫∫ ∫ T
ψi(x, t ) T) ) 0; (final condition)
(
u0j ψi + K0jj
∂ψi )0 ∂xj
(
∂ψi )0 ∂xj
)
)
L1
L2
0
0
-
∫∫ ∫ T
0
-
L2
0
∫∫ ∫ T
0
L1
0
[(
L3
0
L3
0
T
0
L1
0
L2
0
[ψIe0I ]x3)0dx2 dx1dt
(15)
where e0I ) (∫T0 ∫L01 ∫L02e0I dx2 dx1dt)(TL1L2)-1. If one applies an additive perturbation that is constant in time and uniform
∂ψi ) 0; (top) ∂x3
spatially, such as δRI ) δR0I , one can determine a potential sensitivity, SI. For example, the potential sensitivity of the response generated by a constant additive perturbation to emissions is given by
∂ψi ) 0; (bottom). ∂x3
SI )
i
L2
0
3
2
L3
0
i
∂c0i δK11 ∂x1
)]
i 0(x)
dx3dx2dx1
dx3 dx2dt (inflow only)
ψi δu2(b0i - c0i ) + u02δbi + ∂c0i δK22 ∂x2
)]
L2 0
∫∫ ∫ T
0
L1
0
L2
0
[ψI]x3)0dx2 dx1dt
δe0I
∫∫ ∫ 0
L1 0
δe0I
T
1
ψi δu1(b0i - c0i ) + u01δbi +
[(
(14)
j)1, 2
i
L1
∫∫ ∫
sI ) e0I SI )
L3
0
[ψIe0I ]x3)0dx2 dx1dt
and
∫ ∫ ∫ ∫ ψ δQ dx dx dx dt + ∫ ∫ ∫ ψ (x, t ) 0)δc 0
L2
0
; (lateral outflow)
The specific change in response is given by T
0
)
e0I
; (lateral inflow)
v0diψi - K033
0
L1
0
j)1,2
u03ψi + K033
δR )
(13)
dx3 dx1dt (inflow only)
L1
0
)
L2
0
[ψI]x3)0dx2 dx1dt (16)
The potential sensitivity reveals, not the sensitivity to a scaling of the actual emissions, but the sensitivity to potential emissions anywhere in the domain. Implementation. The numerical solution of eq 11, solved backward in time from time t ) T to t ) 0, is inherently unstable. Unlike the cases of the model equation and the forward sensitivity equation, the adjoint equation effectively includes negative diffusion, which enhances numerical perturbations instead of smoothing them, leading to an unstable solution. To avoid this numerical instability, we make a change in the independent time variable: τ ≡ T t. Another useful transformation is a change in sign of the wind input. If we define wj ≡ -uj, and apply the change in the time variable, then eq 11 becomes VOL. 40, NO. 8, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2665
( )
∂ψi ∂(w0j ψi) ∂ 0 ∂ψi ) K + JTikψk + fi + ∂τ ∂xj ∂xj jj ∂xj
(17)
The initial and boundary conditions for eq 17 become
ψi(x, τ ) 0) ) 0; (initial condition)
(
w0j ψi - K0jj
(
∂ψi )0 ∂xj
)
∂ψi )0 ∂xj
)
; (lateral inflow)
j)1,2
; (lateral outflow)
j)1,2
w03ψi - K033 v0diψi - K033
∂ψi ) 0; (top) ∂x3
∂ψi ) 0; (bottom) ∂x3
It is important to note that all data that depend on t in the original adjoint equation, still depend on t ) T - τ in this modified adjoint equation. Because of the change in the time variable, what was the final condition has become the initial condition. Because of the change in sign of the velocity components, what was the inflow lateral boundary condition has become the outflow boundary condition and vice versa. With these changes, the adjoint equations become similar to those of the forward method, and thus nearly identical numerical methods can be applied. This study applied the California Institute of Technology (CIT) model (34, 35) for the numerical solution of the model equations and the 3-D version of the decoupled direct method (DDM-3D) (36), based on CIT, for the numerical solution of the sensitivity equations. CIT was selected because one of the objectives of the study was to make a comparison of sensitivities generated with both DDM and ASAP, and CIT had implemented DDM for many different types of model parameters. The chemical mechanism used in this study is an extended version of the SAPRC99 mechanism (37). As for CIT and DDM-3D, the numerical solution of the adjoint equations (eq 17) employed a first-order accurate (in time) operator-splitting method. For the horizontal advection and horizontal diffusion operations, the numerical algorithms used were identical to what is used for DDM-3D (36) except that the boundary value for ψi was set to zero. The advection solver is a linear finite-element method with a nonlinear filter (34, 38). The numerical solution method applied for ASAP, like that of DDM-3D, does not account for the nonlinearity of the advection solver in the differentiation of the advection scheme. (See ref 39 for a discussion of differentiation of nonlinear advection solvers.) For the vertical advection step, no bottom boundary condition is required other than w03 ) 0. At the model top, for inflow, the equations ψtop ) 0 and ψN3 ) ψtop satisfy the boundary conditions; no boundary condition is required for outflow. The quantity ψN3 is defined at the center of the top cell, where N3 is the number of vertical layers. The numerical solution of the model equations couples the vertical diffusion and chemistry operators (34). In DDM3D, these operators were decoupled, increasing efficiency, but creating an inconsistency between the numerical solution of the model equations and the sensitivity equations. For ASAP, the coupling of the vertical diffusion and chemistry steps was maintained by solving the block tridiagonal system of simultaneous equations:
Bψn+1 ) Cψn + f 2666
9
(18)
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 8, 2006
where n and n + 1 denote time steps and where matrixes B and C have dimension (N × N3) × (N × N3) and depend on both the transposed Jacobian matrix and a centereddifference vertical diffusion operator. The adjoint vector ψ and the response function vector f have dimension (N × N3) × 1 corresponding to each species for each vertical layer. The top boundary condition for the diffusion operator is satisfied by setting the diffusion flux equal to the advection flux at the model top, noting that for inflow the advective flux is zero, since ψtop ) 0. The components of B, C, and f and the times for which they were evaluated differ for active versus steady-state species, as defined in the SAPRC99 chemical mechanism. This is because for steady-state species eq 1 reduces to a set of algebraic equations and, subsequently, both the forward sensitivity equations and the adjoint equations differ for active versus steady-state species. Details of the numerical solution method can be found elsewhere (33). In the implementation of ASAP, automatic code generation was applied for the Jacobian and the direct derivative matrixes. However, the model formulation was continuous, as opposed to discrete (26), and no automated code generation was applied for the formulation of derivatives. Instead, existing code used for DDM-3D was applied directly, with little modification. Updating the adjoint function with respect to vertical diffusion, chemistry, and the response function required factoring the sparse, but stiff, leading matrix of eq 18. To solve this system of equations, we used the UMFPACK routines (40) and took advantage of the sparsity due to the block tridiagonal structure. Once the adjoint function was calculated, the integration of known quantities supplied the sensitivities (eqs 12, 13). To perform these integrations, the trapezoidal rule was used to integrate over time, where instantaneous values are known at the time interval boundaries. Simple summation was used for spatial integration since the cell values approximate volume-average quantities. ASAP was implemented in three steps. A preliminary step involved modification of the CIT model to output all concentration fields and data needed for subsequent steps. Next, program ASAP1 was written to implement the numerical methods described above to compute the adjoint function. Then, program ASAP2 was written to compute sensitivities from the perturbed data and the adjoint function. One of the of main challenges in implementing ASAP was that the adjoint function could not be computed as the air quality model equations were being solved because the adjoint equation starts at the model ending time. Rather than trying to store all the data required by the adjoint equation in computer memory, these data were simply output to disk. The modified CIT model generated three additional output files: (1) all the hourly data read in by the model itself, (2) the average of the species concentration fields before and after the chemistry time step, and (3) instantaneous concentration fields output at the initial time and at the end of each time step. The time interval for output was selected to be less than the minimum advection time step and to be short enough that the concentration fields did not vary significantly from one output interval to the next. This required that the maximum CIT model time step be set not to exceed the output interval. Experimentation suggested that a 3-minute output interval was sufficient. ASAP1 reads the hourly data and the averaged concentration fields; it computes the adjoint function for a given response; and it outputs the instantaneous adjoint function fields, also at 3-minute intervals. ASAP2 reads the hourly data, the instantaneous concentrations, instantaneous adjoint fields, and, optionally, an hourly perturbed data input. ASAP2 generates sensitivities, semi-normalized sensitivities, and optionally a change in the response. It is only necessary to
FIGURE 1. Comparison of ASAP and DDM semi-normalized sensitivities to those generated using the brute-force method by (a) sensitivity type and (b) location. rerun the air quality model if the base data changes. For a given base case, it is only necessary to rerun ASAP1 if a different model response is being studied. ASAP2 can generate sensitivities of a given response to any model data; ASAP2 runs very quickly compared to the air quality model and ASAP1.
Results and Discussion Sensitivity Comparisons. Sensitivities generated with ASAP were compared (Figure 1) to those generated using two other methods: (1) the “brute-force” method (BFM), whereby an approximation to the local derivative was formed as a centered difference of two model runs, and (2) the decoupled direct method, DDM-3D (36). The test case selected for making sensitivity comparisons was a three-day ozone episode from the Southern California Air Quality Study (41). Model set up and evaluation have been described elsewhere (5, 42). The sensitivity of 1-hour and 8-hour peak ozone on day three of the simulation at each of four sites were derived and compared using the three methods for the following parameters (and brute-force parameter perturbations): surface emissions of lumped aromatics ((20%), ethene ((15%),
toluene ((25%), and NO ((5%); boundary conditions of ozone ((10%), and NO ((10%); reaction rate coefficients for NO2 + hν ((1%), HO + NO2 ((1%), HCHO + hν f 2 HO2 + CO ((2%), and O3 + NO ((2%); and deposition velocity for NO ((10%), NO2 ((10%), and ozone ((10%). Sensitivities of 1-hour ozone to the initial conditions for NO ((5%) and ozone ((5%) were also compared for times 1 h, 2 h, and 3 h after the simulation start time. The four sites selected for comparisons included coastal Hawthorne, upwind of the core urban area and most of the emissions in the modeling domain (see Figure 2), central Los Angeles, and inland sites Azusa and Rubidoux, which experienced some of the highest ozone levels during the episode. Figure 1a compares semi-normalized sensitivities by type from both DDM and ASAP versus BFM. Figure 1b shows the same comparison by site. Semi-normalized sensitivities from both ASAP and DDM showed good agreement with those from BFM: R2 ) 0.977 for ASAP and R2 ) 0.985 for DDM. The largest negative sensitivities were the sensitivities of (day three) peak ozone at Azusa and Central LA to NO surface emissions and the sensitivities of peak ozone at Azusa to the rate coefficient for the reaction HO + NO2, which also showed the largest deviations from the 1:1 line for both DDM and ASAP. The largest positive sensitivities were those of ozone at Rubidoux (1 h after the start of the simulation) to ozone initial conditions and the sensitivities of peak ozone at Azusa to the rate coefficient for the NO2 + hν reaction. Some ASAP sensitivities of peak ozone at Hawthorne to rate coefficients deviated significantly from the 1:1 line while, relative to BFM, DDM underpredicted sensitivities of peak ozone at Central LA to boundary conditions. Nonlinearities in the ozone responses, which are not produced by either DDM or ASAP, are likely the cause of some of the differences in this comparison. Earlier investigations (43, 44) have explored the importance of nonlinearities in such comparisons. Decreasing the data perturbations for the BFM to reduce the effects of nonlinearities often improved the comparisons until the perturbations resulted in ozone differences that were within the range of numerical roundoff errors. Another possible source of discrepancy may stem from the fact that the nonlinear components of the advection solver were not differentiated in either ASAP or DDM-3D. (See also ref 43 for discussion of this source of discrepancy.) The comparisons of Figure 1 did not reveal any benefit to having preserved the coupling of vertical diffusion and chemistry solvers in ASAP, relative to DDM in which the coupling of these solvers was not preserved. A more thorough investigation of this point is warranted using numerical methods with a higher-order accuracy in time. Resource Requirements. After evaluating the sensitivities generated by ASAP, we compared its resource requirements with those of the other two methods. The CPU use, memory requirements, and disk storage space requirements were all compared as the number of sensitivities generated for a single response was increased. Performance tests were run on a Dell PC with a single Intel Pentium-4 2.8 GHz processor and 2 GB of system memory. The CPU requirements of BFM scale linearly with the number of sensitivities required. For each centered-difference calculation of a sensitivity, the air quality model was run twice. Two simulations require approximately 0.5 CPU seconds per model grid cell per simulated day. The CPU requirements of DDM were about the same as those of CIT for a single sensitivity: one sensitivity required about 0.5 CPU seconds per cell per simulated day. For additional sensitivities, however, DDM-3D was more efficient; each additional sensitivity cost only about 10% of the cost of the first. This does not continue indefinitely however. After about 100 data sensitivities, for any realistic problem, computer memory limits the number of sensitivities that can be VOL. 40, NO. 8, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2667
FIGURE 2. Estimated spatial distribution of the semi-normalized sensitivity (ppb) of ozone to NO emissions, using DDM-3D. The distribution shown is for the lowest model layer at noon on 25 June. simultaneously run with DDM. Beyond 100 sensitivities, depending on the problem size and system memory, it is necessary to start another simulation. The CPU requirements of ASAP are almost entirely (more than 90%) due to the inversion of the coupled vertical diffusion-Jacobian matrix in ASAP1. After the adjoint function is computed in ASAP1, computing the sensitivities for a given response is extremely efficient. The cost of computing one sensitivity with ASAP is about 3 CPU seconds per cell per simulated day. The cost of computing 1000 is about the same. In fact, computing hundreds of sensitivities becomes routine with ASAP. The ASAP becomes more efficient than brute force after about five sensitivities and more efficient than DDM after about 40 sensitivities for a single response. Recall, however, that to generate sensitivities to a different response, for example at a different site or time, requires an additional ASAP1 simulation. The maximum memory requirements are about the same, roughly 2 KB per grid cell, for all three methods for the sensitivity of one response to a single parameter. The DDM, while somewhat less demanding of maximum memory than BFM, will exceed the memory resources of most systems after about 100 sensitivities. For the grid configuration of this study, there were 12 000 cells. For this configuration, 100 KB memory per cell is equal to about 1 GB of system memory. For ASAP, additional sensitivities do not affect the maximum memory use, since maximum memory is used in computing the adjoint function in ASAP1. Once the adjoint function is computed, the memory requirements to compute sensitivities are minor. The disk storage requirements for BFM and DDM are roughly the same. For both methods, the base case concentration fields represent about 100 KB of storage per cell per simulated day. Each additional sensitivity need not generate a large amount of additional output. For BFM, usually only one additional field for the relevant response for each of the two simulations required using the centereddifference method are needed. For DDM, usually only one additional field for the relevant response is needed. The ASAP is a more disk intensive process than the other methods. Generating the adjoint function takes nearly 10 times the storage capacity of the base model output. This is because, in order to generate the adjoint function, ASAP requires all instantaneous concentration fields to be output at high time resolution. In this case, two sets of instantaneous concentrations were output every 3 min, 20 times the frequency of the hourly output usually generated. The adjoint function itself is also output at the same time interval to calculate sensitivities. For the grid configuration of this study and a 60 h simulation, ASAP required about 18 GB of disk space. Once the adjoint function is generated, however, generating sensitivities to many parameters for a single response does not require any appreciable additional storage capacity. Many applications of air quality models involve more days and 2668
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 8, 2006
larger domains, than those used for this study. It would not be unusual, for example, to have a study domain that contained 20 times the number of grid cells and 2-3 times the number of days simulated here. For such an application, the ASAP as currently implemented would require about 1 TB of disk space, which, while not prohibitive, imposes special considerations in the computer system design for the analysis of multiple responses. The ASAP appears to be most efficient when the number of responses to investigate is small compared to the number of desired sensitivities. Because of the efficiency of the method, it is not necessary to predefine the data with respect to which one wishes to calculate sensitivities. The method can be used to calculate all sensitivities; once calculated, the sensitivities may then be sorted to determine which are most important. Sensitivity Analysis Products. This section presents analysis products derived from ASAP and compares these to products available from DDM. With ASAP, for a given response, one does not need to guess which sensitivities to examine a priori because the extra effort to generate all sensitivities is low. With DDM, for a given parameter, the analysis of multiple responses requires little additional effort. In this sense the ASAP and DDM tools are complementary. It is often the case, however, that the number of interesting responses is small compared to the number of parameters. While for a single response ASAP can generate many sensitivities more efficiently than either BFM or DDM, one might argue that ASAP, while efficient, does not provide anything really new. There are, however, products that the ASAP provides that would be quite difficult to generate using forward sensitivity techniques: the adjoint function can be used to determine when and where perturbations to model data contribute to the change in model response. In other words, ASAP can be used as a linear source apportionment tool, or more accurately, as a sensitivity apportionment tool. The adjoint function, appropriately integrated, yields the firstorder potential sensitivity of the response to each model input and parameter. Similarly, the adjoint function yields the spatial and temporal contribution of model inputs and parameters to the target response or sensitivity. As a source apportionment method, ASAP has the advantage that unlike other source apportionment tools (39), with ASAP one need not define the source regions a priori: ASAP maps them as part of the analysis. However, ASAP does not produce a complete source attribution since it provides only a firstorder analysis. DDM generates fields of sensitivities of various possible responses to a single parameter. Figure 2 presents the case where the responses are ozone at many locations, all model grid cells, and where the single parameter is the surface emissions of NO. The estimated spatial distribution of seminormalized sensitivity is shown in Figure 2, near the time of peak ozone at Azusa and Rubidoux. NO emissions suppress
FIGURE 3. Estimated spatial distribution of the potential sensitivity (ppb/(kg/km2 hr)) of peak ozone at Azusa on 25 June to (uniform) perturbations in NO emissions, using ASAP. Potential sensitivity was calculated as the sum over each hour of ψNO∆t∆x1∆x2, where ψNO is the component of the adjoint function corresponding to NO, ∆t is the time step, and ∆x1 and ∆x2 are the horizontal dimensions of a grid cell. Fields are shown on 24 June at (a) 1200 PST, (b) 1500 PST, and (c) 1800 PST.
FIGURE 4. Estimated spatial apportionment of the semi-normalized sensitivity (ppb) of the 1-hour peak ozone at Azusa on 25 June to (actual) NO emissions on 24 June at (a) 1200 PST, (b) 1500 PST, and (c) 1800 PST. The sensitivity contribution was calculated using ASAP as the sum over each hour of ψNOeNO∆t∆x1∆x2 where ψNO is the component of the adjoint function corresponding to NO, eNO represents surface emissions of NO, ∆t is the time step, and ∆x1 and ∆x2 are the horizontal dimensions of a grid cell.
maximum ozone at Azusa, where semi-normalized sensitivities are less than -140 ppb, suggesting VOC-sensitive chemistry at this site. This is useful information if one wants to understand how changes to NO emissions influence many areas and if one wants to examine how that influence varies with time.
4a). Thus, while offshore NO emissions have the potential to affect the next day’s peak ozone as shown in Figure 3a, there simply were not many sources of NO emissions over the ocean. As the air mass that affects Azusa on 25 June moves onshore, the NO emissions entrained into the air mass increase, as does the magnitude of their contribution to Azusa’s peak ozone (Figures 4b,c). DDM provides a source-based sensitivity analysis procedure, whereby the sensitivity of a large number of receptors to a single perturbation of sources can be analyzed. ASAP provides a receptor-based sensitivity analysis procedure, whereby the sensitivity of a single receptor to a large number of source perturbations can be analyzed. While, these two methods are complementary, to date, most sensitivity analyses of photochemical models have been source-based, using forward methods such as DDM. Emerging tools, such as ASAP, based on the adjoint function, facilitate receptorbased analyses.
The ASAP generates fields of a completely different nature (Figure 3). The NO component of the adjoint function when integrated over each surface-level model grid cell and over all time steps provides the first-order potential sensitivity of peak ozone at Azusa to each point in a uniformly perturbed NO emissions field. Figure 3a shows that at noon on 24 June, increased NO emissions from a region offshore would contribute positively to peak ozone at Azusa on 25 June. A few hours later (Figure 3b) the region of positive contribution is reduced and, as the flow pushes onshore, a region of negative contribution develops near the coast. By the evening of 24 June, NO emissions reduce peak ozone at Azusa (Figure 3c). For comparison to Figure 3, Figure 4 shows the actual apportionment of the semi-normalized sensitivity of the 1-hour peak ozone at Azusa on 25 June to NO emissions on the previous day. Since the amount of NO emissions offshore (from shipping) was small relative to that from other sources, its contribution to Azusa peak ozone was also small (Figure
Acknowledgments This research is an extension of work originally funded by the California Air Resources Board under contract no. 98309. The statements herein are those of the authors and not necessarily those of the California Air Resources Board. We are grateful to the Bay Area Air Quality Management District VOL. 40, NO. 8, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2669
for leaves of absence granted to P. T. Martien to complete this research as part of his doctoral dissertation. We also acknowledge the helpful comments provided by anonymous reviewers.
Literature Cited (1) National Research Council. Rethinking the Ozone Problem in Urban and Regional Air Pollution; National Academy Press: Washington, DC, 1991. (2) Russell, A.; Dennis, R. NARSTO critical review of photochemical models and modeling. Atmos. Environ. 2000, 34, 22832324. (3) Bergin, M. S.; Russell, A. G.; Milford, J. B. Effects of chemical mechanism uncertainties on the reactivity quantification of volatile organic compounds using a three-dimensional air quality model. Environ. Sci. Technol. 1998, 32, 694-703. (4) Khan, M.; Yang, Y.-J.; Russell, A. G. Photochemical reactivities of common solvents: comparison between urban and regional domains. Atmos. Environ. 1999, 33, 1085-1092. (5) Martien, P. T.; Harley, R. A.; Milford, J. B.; Russell, A. G. Evaluation of incremental reactivity and its uncertainty in Southern California. Environ. Sci. Technol. 2003, 37, 1598-1608. (6) Milford, J. B.; Gao, D.; Russell, A. G.; McRae, G. J. Use of sensitivity analysis to compare chemical mechanisms for air quality modeling. Environ. Sci. Technol. 1992, 26, 1179-1189. (7) Cacuci, D. G. Sensitivity and Uncertainty Analysis. Volume I. Theory; Chapman & Hall/CRC: Boca Raton, FL, 2003. (8) Bergin, M. S.; Noblet, G. S.; Petrini, K.; Dhieux, J. R.; Milford, J. B.; Harley, R. A. Formal uncertainty analysis of a Lagrangian photochemical air pollution model. Environ. Sci. Technol. 1999, 33, 1116-1126. (9) Reynolds, S.; Micheals, H.; Roth, P.; Tesche, T.; McNally, D.; Gardner, L.; Yarwood, G. Alternative base cases in photochemical modeling: their construction, role, and value. Atmos. Environ. 2001, 30, 1977-1988. (10) Hanna, S. R.; Lu, Z.; Frey, H. C.; Vukovich, N. W. J.; Arunachalam, S.; Fernau, M.; Hansen, D. A. Uncertainties in predicted ozone concentrations due to input uncertainties for the UAM-V photochemical grid model applied to the July 1995 OTAG domain. Atmos. Environ. 2001, 35, 891-903. (11) Moore, G. E.; Londergan, R. J. Sampled Monte Carlo uncertainty analysis for photochemical grid models. Atmos. Environ. 2001, 35, 4863-4876. (12) Cacuci, D. G. The Forward and the Adjoint Methods of Sensitivity Analysis, Uncertainty Analysis; Ronen, Y., Ed.; CRC Press: Boca Raton, FL, pp. 71-144, 1988. (13) Hwang, J.; Dougherty, E. P.; Rabitz, S.; Rabitz, H. The Green’s function method of sensitivity analysis in chemical kinetics. J. Chem. Phys. 1978, 69, 5180-5191. (14) Cho, S.-Y.; Carmichael, G.; Rabitz, H. Sensitivity analysis of the atmospheric reaction-diffusion equation. Atmos. Environ. 1987, 21, 2589-2598. (15) Vuilleumier, L.; Harley, R. A.; Brown, N. J. First- and secondorder sensitivity analysis of a photochemically reactive system (a Green’s function approach). Environ. Sci. Technol. 1997, 31, 1206-1217. (16) Marchuk, G. I. Adjoint Equations and Analysis of Complex Systems; Kluwer Academic Publishers: Boston, MA, 1995. (17) Hall, M.; Cacuci, D. Physical interpretation of the adjoint functions for sensitivity analysis of atmospheric models. J. Atmos. Sci. 1983, 40, 2537-2546. (18) Lewis, J. M.; Derber, J. C. The use of adjoint equations to solve a variational adjustment problem with advective constraints. Tellus 1985, 37A, 309-322. (19) Daley, R. Atmospheric Data Analysis; Cambridge University Press: Cambridge, UK, 1991; pp 242-262. (20) Navon, I. M.; Zou, X.; Derber, J.; Sela, J. Variational data assimilation with an adiabatic version of the NMC spectral model. Mon. Weather Rev. 1992, 120, 1433-1446. (21) Errico, R. M. What is an adjoint model? Bull. Am. Met. Soc. 1997, 78, 2577-2591. (22) Cho, S.-Y.; Carmichael, G.; Rabitz, H. Relationships between primary emissions and acid deposition in Eulerian models determined by sensitivity analysis. Water Air Soil Pollut. 1988, 40, 9-31.
2670
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 8, 2006
(23) Elbern, H.; Schmidt, H.; Ebel, A. Variational data assimilation for tropospheric chemistry modeling. J. Geophys. Res. 1997, 102(D13), 15967-15986, doi: 10.1029/97JD01213. (24) Elbern, H.; Schmidt, H.; Talagrand, O.; Ebel, A. 4D-variational data assimilation with an adjoint air quality model for emission analysis. Environ. Model. Software 2000, 15, 539-548. (25) Menut, L.; Vautard, R.; Beekmann, M.; Honore´, C. Sensitivity of photochemical pollution using the adjoint of a simplified chemistry-transport model. J. Geophys. Res. 2000, 105(D12), 15379-15402, doi: 10.1029/1999JD900953. (26) Sirkes, Z.; Tziperman, E. Finite difference of adjoint or adjoint of finite difference? Mon. Weather Rev. 1997, 125(14), 33733378. (27) Schmidt, H.; Martin, D. Adjoint sensitivity of episodic ozone in the Paris area, to emissions on the continental scale. J. Geophys. Res. 2003, 108(D17), doi: 10.1029/2001JD001583. (28) Hess, P. G.; Vukic´evic´, T. Intercontinental transport, chemical transformations, and baroclinic systems. J. Geophys. Res. 2003, 108(D12), doi: 10.1029/2001JD001583. (29) Sandu, A.; Daescu, D. N.; Carmichael, G. R. Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part I - theory and software tools. Atmos. Environ. 2003, 37, 50835096. (30) Daescu, D. N.; Sandu, A.; Carmichael, G. R. Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part II - numerical validation and applications. Atmos. Environ. 2003, 37, 5097-5114. (31) Sandu, A.; Daescu, D. N.; Carmichael, G. R.; Chai, T. Adjoint sensitivity analysis of regional air quality models. J. Comput. Phys. 2005, 204, 222-252. (32) Dunker, A. M. The direct decoupled method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 1984, 81, 2385-2393. (33) Martien, P. T. Forward and Adjoint Sensitivity Analysis in Eulerian Photochemical Air Quality Models. PhD dissertation, University of California, Berkeley, CA, 2004. (34) McRae, G. J.; Goodin, W. R.; Seinfeld, J. H. Numerical solution of the atmospheric diffusion equation for chemically reacting flows. J. Comput. Phys. 1982, 45, 1-42. (35) Harley, R. A.; Russell, A. G.; McRae, G. J.; Cass, G. R.; Seinfeld, J. H. Photochemical modeling of the Southern California air quality study. Environ. Sci. Technol. 1993, 27, 378-388. (36) Yang, Y.-J.; Wilkinson, J. G.; Russell, A. G. Fast, direct sensitivity analysis of multidimensional photochemical models. Environ. Sci. Technol. 1997, 31, 2859-2868. (37) Carter, W. P. L. Documentation of the SAPRC99 Chemical Mechanism for VOC Reactivity Assessment, Final Report to the California Air Resources Board under contract no. 92-329 and 95-308; Center for Environmental Research and Technology: University of California, Riverside, CA, 2000. (38) Forester, C. K. Higher order monotonic convective difference schemes. J. Comput. Phys. 1977, 23, 1-22. (39) 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. (40) Davis, T. A.; Duff, I. S. An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM J. Matrix Anal. Appl. 1997, 19, 140-158. (41) Lawson, D. R. The Southern California air quality study. J. Air Waste Manage. Assoc. 1990, 40, 156-165. (42) McNair, L. A.; Harley, R. A.; Russell, A. G. Spatial inhomogeneity in pollutant concentrations, and their implications for air quality model evaluation. Atmos. Environ. 1996, 30, 4291-4301. (43) Hakami, A.; Odman, M. T.; Russell, A. G. Nonlinearity in atmospheric response: A direct sensitivity analysis approach. J. Geophys. Res. 2004, 109, D15303, doi: 10.1029/2003JD004502. (44) 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.
Received for review May 31, 2005. Revised manuscript received December 19, 2005. Accepted February 9, 2006. ES0510257
