A method to represent ozone response to large changes in precursor emissions using high-order sensitivity analysis in photochemical models

Abstract. Photochemical grid models (PGMs) are used to simulate tropospheric ozone and quantify its response to emission changes. PGMs are often applied for annual simulations to provide both maximum concentrations for assessing compliance with air quality standards and frequency distributions for assessing human exposure. Efficient methods for computing ozone at different emission levels can improve the quality of ozone air quality management efforts. This study demonstrates the feasibility of using the decoupled direct method (DDM) to calculate first- and second-order sensitivity of ozone to anthropogenic NOx and VOC emissions in annual PGM simulations at continental scale. Algebraic models are developed that use Taylor series to produce complete annual frequency distributions of hourly ozone at any location and any anthropogenic emission level between zero and 100%, adjusted independently for NOx and VOC. We recommend computing the sensitivity coefficients at the midpoint of the emissions range over which they are intended to be applied, in this case with 50% anthropogenic emissions. The algebraic model predictions can be improved by combining sensitivity coefficients computed at 10 and 50% anthropogenic emissions. Compared to brute force simulations, algebraic model predictions tend to be more accurate in summer than winter, at rural than urban locations, and with 100% than zero anthropogenic emissions. Equations developed to combine sensitivity coefficients computed with 10 and 50% anthropogenic emissions are able to reproduce brute force simulation results with zero and 100% anthropogenic emissions with a mean bias of less than 2 ppb and mean error of less than 3 ppb averaged over 22 US cities.


Introduction
Ozone is a natural trace constituent of the troposphere that is influenced by emissions from human activities (Warneck, 2000).Precursor emissions that influence tropospheric ozone formation are nitrogen oxides (NO x ), volatile organic compounds (VOC) and other compounds such as CO with atmospheric chemical reactions similar to VOCs.Government agencies regulate ozone precursor emissions to achieve air quality objectives and to reduce adverse effects on human health and welfare.Ozone air quality management strategies are complex to develop because (1) the chemical reactions that produce and destroy ozone constitute a nonlinear system, and (2) tropospheric ozone has a lifetime of days to weeks, permitting transport across jurisdictional boundaries by prevailing winds.
Computer models are used to simulate tropospheric ozone and to quantify the effects from emission reduction strategies (Rao et al., 2011).The most widely used ozone models are photochemical grid models (PGMs) that represent the atmosphere as a three-dimensional (3-D) grid of cells, describe physical processes (emission, transport, deposition) that change the mass of individual compounds within grid cells, and describe photochemical reactions that transform chemicals one to another within each grid cell.Tropospheric ozone concentrations have wide spatial and temporal variability.PGMs are applied at spatial scales from local (a few km) to global and with sub-hourly time resolution over time periods of a year or more.
Developing more efficient methods for applying PGMs to evaluate many different emission reduction strategies can improve the quality of ozone air quality management efforts.
As applied in the US, ozone models must be able to represent maximum values, such as daily maximum 8 h average (DMA8), to assess compliance with ambient standards and to provide multi-year ozone frequency distributions for use in human exposure assessments (US EPA, 2012).US regulators are interested in evaluating the range of emissions extending from current conditions (100 % anthropogenic emissions) to zero US anthropogenic emissions.PGMs allow regulators to ask, for example, what is the ozone frequency distribution after emissions are reduced by control strategies?The answer to this question varies from one city to another, as well as within each city.A traditional "brute force" approach of repeating model simulations with altered emissions requires many such simulations to provide answers.
This study applies an existing mathematical method called the decoupled direct method (DDM; Dunker, 1981) from which to conduct PGM sensitivity analyses for annual simulations at continental scale.The goal is to produce algebraic models based on the sensitivity analysis that can accurately represent complete annual frequency distributions of hourly ozone at any location and any anthropogenic emission level between zero and 100 %, adjusted independently for NO x and VOC.In this context, accuracy refers to how well the algebraic model can reproduce the brute force approach of re-running the PGM with altered emission inputs.Challenges include the computational burden of computing emission sensitivity in a continental scale, year-long model application.Strategies for applying DDM efficiently are proposed and evaluated.The methods proposed here are compared to a similar method described recently by Simon et al. (2013).

Background
The nonlinear relationship between ozone formed and initial precursors is illustrated using box model simulations in Fig. 1.This figure shows a response surface of daily maximum 1 h average ozone (ppb) constructed from 121 simulations with varying initial NO x (ppb) and VOC (ppbC) using the 2005 version of the Carbon Bond chemical mechanism (Yarwood et al., 2005).The ozone response surface is curved throughout indicating that ozone should generally be expected to respond nonlinearly to changes in NO x and/or VOC.When the ratio of VOC / NO x is high, ozone is reduced most effectively by reducing NO x (NO x -limited condition).In contrast, when the ratio of VOC / NO x is low, ozone is reduced most effectively by reducing VOC (VOC-limited condition).
The line through points A and B in Fig. 1 defines a section of changing NO x at constant VOC (200 ppbC).If the derivatives of ozone with respect to NO x (S (1)  = ∂O 3 /∂NO x ; S (2)  = ∂ 2 O 3 /∂ 2 NO x ; etc.) are known, then the response of ozone to a relative NO x perturbation ( N = NO x / NO x ) can be approximated using a second-order Taylor series Figure 2 shows the section through points A and B in Fig. 1, including first-and second-order approximations to the ozone response from changing NO x at both A and B using Eq.(1).The example chosen for Fig. 2 is particularly demanding in selecting a highly curved section of the ozone response surface that traverses a local maximum.This example shows that for purposes of representing ozone in the range of NO x from zero to A, (1) second-order representations are generally, but not always, more accurate than first order; (2) second-order representations do not necessarily provide good accuracy (within a factor of two) for response to large perturbations (factor of two change in NO x ); and (3) accuracy is improved by applying Eq. (1) at the center of the range of interest rather than at one extreme.These findings are not surprising but are presented here to explain choices made in this study for applying sensitivity methods to the more complex problem of ozone sensitivity to emissions in a 3-D PGM.

Decoupled direct method
Two widely used methods for computing emissions sensitivity within a PGM are the DDM (Dunker, 1981) and the adjoint method (Menut et al., 2000).Both methods can efficiently compute first-order sensitivity coefficients (S (1)  = ∂C/∂E) for a model output such as concentration (C) with respect to an input such as emission (E).The adjoint method calculates only first-order sensitivity coefficients (Menut et al., 2000).However, Hakami et al. (2003) demonstrated the application of DDM to second-and third-order sensitivity coefficients (S (2) = ∂ 2 C/∂E 2 ; S (3) = ∂ 3 C/∂E 3 ) and called the application high-order DDM (HDDM).The sensitivity coefficients needed for this study were computed using HDDM following Dunker et al. (2002) and Cohan et al. (2010).

Methods
Modeling was performed with version 5.40 of the Comprehensive Air quality Model with extensions (CAMx; ENVI-RON, 2012) with model configuration and inputs developed for a previous study of North American background ozone (Emery et al., 2012).Calendar year 2006 was modeled for an outer domain with 36 km grid cells covering the conterminous US and portions of Canada and Mexico (Fig. 3).Two nested inner domains with 12 km grid cells covered the eastern and western US.The nested 12 km domains communicated with the 36 km grid and with each other by twoway nesting.Meteorological and emissions data were developed by the US EPA for the Air Quality Model Evaluation International Initiative (AQMEII) program (Rao et al., 2011).Meteorology was modeled using the Weather Research and Forecast (WRF) model (Skamarock et al., 2008) at 12 km resolution (Vautard et al., 2012).Emissions data were also at 12 km resolution (Pouliot et al., 2012) and included the US EPA 2005 national emissions inventory grown to 2006 (EPA, 2010), the Environment Canada 2006inventory (Environment Canada, 2011), a 1999 inventory for Mexico grown to 2006, biogenic emissions from BEIS version 3.14 (Vukovich and Pierce, 2002), and fire emissions (Coe-Sullivan, 2008).Chemical boundary conditions for the 36 km grid were down-scaled from 2006 GEOS-Chem version 8-03-01 global model output (Zhang et al., 2011).
Two annual CAMx simulations were conducted with HDDM to compute ozone sensitivity to US anthropogenic emissions of NO x (N) and VOC (V ) at first-and secondorder, yielding sets of 5 ozone sensitivity coefficients from each simulation: N V .Each simulation ran for roughly one month using 48 cores on eight Intel X5675 CPUs.The first HDDM simulation computed sensitivity coefficients with anthropogenic NO x and VOC emissions reduced to 50 %, corresponding to the example discussed above for point B in Fig. 1.The resulting Taylor series expansion to represent ozone at any NO x and VOC level as a function of relative changes in NO x ( N = NO x /NO x ) and VOC ( V = VOC/VOC) from the 50 % emissions case is where the subscript "50" indicates a parameter calculated with 50 % anthropogenic emissions.The second HDDM simulation computed sensitivity coefficients with anthropogenic NO x and VOC emissions reduced to 10 % in order to better represent ozone near zero anthropogenic emissions.Sensitivity coefficients from the 10 % case were combined with sensitivity coefficients from the 50 % case to yield a three-equation algebraic model to estimate ozone for any emissions level between 0 and 100 % NO x and VOC: (2) N( 10) Below 15 % NO x (at any VOC) Eq. ( 3a) is applied exclusively, while above 25 % NO x (at any VOC) Eq. ( 3b) is applied exclusively.Within the 15-25 % NO x range (at any VOC), sensitivity coefficients from both simulations are linearly interpolated across the NO x dimension (N) using Eq.(3c).The transition points between Eqs. (3a), ( 3b) and (3c) (i.e., N = 15 % and N = 25 %) were selected for this application based on results of performance tests with the 3-D PGM, described below.Similarly, defining transition points using N alone performed better than using both N and V because ozone is predominantly limited by NO x emissions in the 3-D PGM simulations.The transition points for Eqs.(3a)-(3c) could be adjusted for different applications.
Month-long brute force CAMx simulations were performed to evaluate the accuracy of Eqs. ( 2) and (3).Simulations were performed for January and July with 100 % NO x and VOC, 25 % NO x (100 % VOC), and zero NO x and VOC US anthropogenic emissions.The 100 % and zero anthropogenic emission cases bound our range of interest for applying Eqs. ( 2) and (3).The 25 % anthropogenic NO x emissions case was useful in developing Eq.(3).
CAMx computes sensitivity coefficients with HDDM for the entire model domain enabling application of Eqs. ( 2) and (3) at any location within the domain.Results were evaluated statistically for all monitoring sites in 22 cities, marked in Fig. 3, by computing the mean bias (MB = (HDDM -BF)/n) and mean error (ME = |HDDM−BF|/n) for the n hourly ozone values predicted by the HDDM and brute force (BF) methods.HDDM and BF model results for grid cells containing monitors were paired to calculate the MB and ME statistics.
In addition, results are shown graphically at the locations of an urban and a rural monitor near Dallas, Texas.Dallas was selected because (1) the area is out of compliance with EPA's ambient ozone standard set at 75 ppb; (2) ozone in Dallas is influenced both by local production and regional transport (Kemball-Cook et al., 2009); and (3) the selected locations exhibit HDDM performance attributes compared to

Results
Figure 4 shows hourly ozone at the Dallas location for the CAMx simulations of July 2006 with 100 % and zero anthropogenic NO x and VOC emissions.Relative to background, anthropogenic emissions alter the entire diurnal distribution of ozone with photochemical ozone production causing higher daytime maxima and ozone titration by NO causing lower nighttime minima.Equation (2) was first used to compute January and July hourly ozone at Dallas and Weatherford for both the 100 % and zero anthropogenic emission cases.Equation (2) relies solely upon the ozone concentration and emission sensitivity coefficients computed with 50 % anthropogenic emissions.
Figure 5 compares the Dallas predictions from Eq. ( 2) to brute force results with each data point showing the comparison at each hour.All four cases shown in Fig. 5 have many hours with good agreement (near the 1 : 1 lines) and so discussion focuses on reasons for poorer agreement.Overall, agreement is better for the 100 % than the zero anthropogenic emissions case.The zero anthropogenic emissions case for January (Fig. 5a) has many hours with non-zero brute force ozone (15-35 ppb) but near-zero ozone predicted by Eq. (2) (i.e., points falling below the middle section of the 1 : 1 line).This type of disagreement occurs for nighttime hours when the CAMx simulation with 50 % anthropogenic emissions had near-zero ozone and consequently had near-zero ozone sensitivity to emissions.For these hours Eq. ( 2) predicts nearzero ozone at any emission level and so fails to predict that reducing NO x emissions can lead to non-zero ozone at the point where emitted NO no longer titrates all available ozone.Equation ( 3) was developed to improve predictions in these hours by incorporating sensitivity coefficients from a CAMx simulation with 10 % anthropogenic emissions that less frequently titrate ozone to near zero.The performance of Eq. ( 3) is discussed below.The July zero anthropogenic emissions case for Dallas (Fig. 5c) has fewer instances of this problem because titration of ozone to near zero is less common in July than in January.The January and July comparisons for the 100 % emissions case (Fig. 5b, d) do not exhibit this problem because when 50 % anthropogenic emissions titrate ozone to near zero, so do 100 % anthropogenic emissions and Eq. ( 2) predicts accurately.However, the January and July comparisons for the 100 % emissions case show poorer agreement for hours with low ozone (< 10 ppb) when Eq. ( 2) over-or underpredicts the brute force results.These are hours when NO x emissions strongly suppress ozone with 50 % anthropogenic emissions causing a strongly nonlinear ozone response to increasing emissions from 50 to 100 %.Predictions in these hours could be improved by incorporating sensitivity coefficients from a CAMx simulation with near 100 % anthropogenic emissions, although this was not tested.The improved performance of Eq. (3) over Eq. ( 2) for the zero anthropogenic emissions case may be seen by comparing Figs. 6 and 5 (Fig. 6 repeats Fig. 5 using Eq. 3 in place of Eq. 2).Equation (3) much improves performance for the zero anthropogenic emission cases by incorporating sensitivity coefficients from a CAMx simulation with 10 % anthropogenic emissions.There is no change in performance for the 100 % emissions case because Eq. ( 3b) is equivalent to Eq. ( 2) for anthropogenic emissions of 25 % and greater.
The performance of Eqs.
(2) and (3) at the rural Weatherford location is shown in Figs.7 and 8, respectively.Both sets of equations perform better at the rural location than at the urban location (Figs. 5, 6) because ozone is never titrated to near zero at Weatherford.The predictions of Eq. ( 2) are slightly above the 1 : 1 line at some hours for the zero anthropogenic emissions case in January.These small errors are corrected by using Eq. ( 3), which provides very close agreement with brute force results in all four cases (Fig. 8).
Review of results for monitoring sites in other cities (not shown) confirmed that Eqs. ( 2) and (3) generally perform better at rural than urban sites because ozone is less frequently titrated to near zero at rural sites.Equation (3) has transition points at anthropogenic NO x emission levels of 25 and 15 %.The upper transition point was selected by comparing the accuracy of Eqs.(3a) and (3b) for predicting ozone with 25 % anthropogenic emissions.Equation (3b) was more accurate than Eq.(3a) (Table 1) leading us to rely exclusively on Eq. (3b) for anthropogenic emission levels of 25 % and higher.The lower transition point of 15 % was selected as being between 25 %, where Eq. (3b) proved more accurate, and 10 %, where Eq. ( 3a) is exact by definition.The transition points in Eq. ( 3) are a choice that can be refined using brute force results in the region of the transition.Plots showing HDDM vs. BF for the 25 % NO x , 100 % VOC case at Dallas and Weatherford sites are shown in Figs. 9 and 10, respectively (compare to Figs. 6 and 8).
The accuracy of Eq. ( 3) over 22 cities is summarized in Table 2 for anthropogenic emission levels of 100 % NO x and VOC, 25 % NO x (100 % VOC), and zero NO x and VOC in January and July.Results for individual cities are provided in Tables S1 and S2.Note that Eq. ( 3) is exact by definition at anthropogenic emission levels of 10 and 50 %.The mean bias and mean error are smaller in January than July except for the zero emission case, where they are comparable.In  January, the greatest mean bias and mean error are 0.71 and 1.93 ppb, respectively, whereas in July the greatest mean bias and mean error are −1.77and 2.49 ppb.
A method similar to Eq. ( 3) was described recently by Simon et al. (2013).Method similarities include use of HDDM to compute second-order sensitivity coefficients and combination of sensitivity coefficients at several emission levels to improve ozone estimation accuracy.A difference is that Eq. (3) uses concentrations derived at two emission levels (10 and 50 % anthropogenic emissions) whereas Simon et al. (2013) use only concentrations derived with 100 % anthropogenic emissions.Comparing how both methods perform would be valuable for understanding their respective strengths and limitations and then developing improved methods.

Conclusions
This study demonstrates the feasibility of applying HDDM to calculate first-and second-order sensitivity of ozone to anthropogenic NO x and VOC emissions in annual PGM simulations at continental scale.The resulting sensitivity coefficients are used to construct algebraic models using Taylor series that can accurately represent complete annual frequency distributions of hourly ozone at any location and any anthropogenic emission level between zero and 100 %, adjusted independently for NO x and VOC.We recommend computing the sensitivity coefficients at the midpoint of the emissions range over which they are intended to be applied, in this case with 50 % anthropogenic emissions.The ozone estimates at varying NO x and VOC emissions levels are only valid for the modeled time period (2006 in this case) and extrapolating results for 2006 to other years would require additional assumptions that are not considered here.The PGM itself has errors associated with inputs, algorithms, discretization, etc., although we note that these types of errors attend all model applications and are not attributable to HDDM.
When ozone predicted by the algebraic model is compared to brute force simulations with zero and 100 % anthropogenic emissions, the predictions tend to be more accurate in summer than winter, at rural than urban locations, and with 100 % than zero anthropogenic emissions.The photochemical reason for these trends is a strongly nonlinear response of ozone to NO x emissions around the point where NO is sufficient to titrate ozone to near zero.The algebraic model predictions are improved by incorporating sensitivity coefficients computed with 10 % in addition to 50 % anthropogenic emissions.Equations developed to combine sensitivity coefficients computed with 10 and 50 % anthropogenic emissions are able to reproduce brute force simulation results with zero and 100 % anthropogenic emissions with a mean bias of less than 2 ppb and mean error of less than 3 ppb averaged over 22 cities.This approach could be extended by incorporating additional sensitivity coefficients (e.g., computed at or near 100 % anthropogenic emissions) at the cost of greater computational expense to obtain the requisite sensitivity coefficients.

Fig. 1 .
Fig. 1.Representative response of maximum 1 h average ozone (ppb) to morning concentrations of emitted NO x and VOC.Section through points A and B is the ozone response to NO x at constant VOC of 200 ppbC.

Fig. 2 .
Fig. 2. Modeled response of ozone to changing NO x at VOC of 200 ppbC.Also shown are first-and second-order Taylor series approximations to the response at points A and B. Points A and B are also marked in Fig. 1.

Fig. 3 .
Fig. 3. CAMx North American modeling domain comprised of an outer 36 km resolution grid and two inner, two-way nested, 12 km resolution grids.The 22 cities used for performance evaluation are marked.

Fig. 4 .
Fig. 4. Modeled ozone for July 2006 at the urban Dallas location with 100 % (solid line) and zero (dashed line) US anthropogenic x and VOC emissions.

Fig. 5 .
Fig. 5. Hourly ozone at the urban Dallas location predicted by HDDM with Eq. (2) vs. brute force CAMx simulations for (a) zero anthropogenic NO x and VOC emissions in January, (b) 100 % anthropogenic NO x and VOC emissions in January, (c) zero anthropogenic NO x and VOC emissions in July, and (d) 100 % anthropogenic NO x and VOC emissions in July.Dashed lines show 1 : 1 agreement and solid red lines show least-squares regressions.

Fig. 6 .
Fig. 6.Hourly ozone at the urban Dallas location predicted by HDDM with Eq. (3) vs. brute force CAMx simulations for (a) zero anthropogenic NO x and VOC emissions in January, (b) 100 % anthropogenic NO x and VOC emissions in January, (c) zero anthropogenic NO x and VOC emissions in July, and (d) 100 % anthropogenic NO x and VOC emissions in July.Dashed lines show 1 : 1 agreement and solid red lines show least-squares regressions.

Table 2 .
Statistical evaluation over 22 cities of using Eq.(3) to predict ozone with 100 % NO x and VOC, 25 % NO x (100 % VOC), and zero NO x and VOC anthropogenic emissions.