Article pubs.acs.org/est
Emulation and Sensitivity Analysis of the Community Multiscale Air Quality Model for a UK Ozone Pollution Episode Andrew V. Beddows,* Nutthida Kitwiroon, Martin L. Williams, and Sean D. Beevers King’s College London, Waterloo, London, SE1 8WA, United Kingdom S Supporting Information *
ABSTRACT: Gaussian process emulation techniques have been used with the Community Multiscale Air Quality model, simulating the effects of input uncertainties on ozone and NO2 output, to allow robust global sensitivity analysis (SA). A screening process ranked the effect of perturbations in 223 inputs, isolating the 30 most influential from emissions, boundary conditions (BCs), and reaction rates. Community Multiscale Air Quality (CMAQ) simulations of a July 2006 ozone pollution episode in the UK were made with input values for these variables plus ozone dry deposition velocity chosen according to a 576 point Latin hypercube design. Emulators trained on the output of these runs were used in variancebased SA of the model output to input uncertainties. Performing these analyses for every hour of a 21 day period spanning the episode and several days on either side allowed the results to be presented as a time series of sensitivity coefficients, showing how the influence of different input uncertainties changed during the episode. This is one of the most complex models to which these methods have been applied, and here, they reveal detailed spatiotemporal patterns of model sensitivities, with NO and isoprene emissions, NO2 photolysis, ozone BCs, and deposition velocity being among the most influential input uncertainties.
■
the “local” sensitivity to small variations in the baseline value of an input, equivalent to taking a partial derivative of the output with respect to that input.13 Such local sensitivity coefficients can be computed directly with a set of coupled auxiliary equations or separately as the Decoupled Direct Method (DDM).14 Extending this to compute second derivatives takes into account nonlinearity in the model sensitivity, and “cross sensitivities” may be calculated to two or more inputs perturbed together.15 Cohan et al.16 applied this to CMAQ version 4.3, extrapolating further from the base case model run by using first and second derivatives in a second order Taylor series expansion, thus extending direct methods into somewhat more than local SA. In the context of sensitivity to input uncertainty, however, it is more informative to perform a “global” SA, one in which all inputs are perturbed simultaneously but independently over their full uncertainty ranges.13 The prohibitively large number of runs required for the global SA, and related tasks like Monte Carlo uncertainty analysis, of models with many uncertain inputs means that some kind of meta-modeling technique is needed to produce a quick to run model surrogate. This idea is not new; Fang et al.17 describe a number of methods for producing meta-models, including polynomial regression, splines, kriging, and neural networks. Other methods which have been used in moderately
INTRODUCTION The adverse health effects of elevated ozone and NO2 concentrations are well documented.1 Consequently, legislation requires that concentrations of pollutants are measured at fixed site monitoring stations, and models such as the Community Multiscale Air Quality (CMAQ) model2 are used to estimate pollutant concentrations between sparsely distributed monitoring locations. The output of such models is subject to uncertainty, caused both by the model’s imperfect representation of reality and by uncertainties in inputs which are propagated through to the outputs.3 This work concerns methods to ascertain the relative contributions of different input uncertainties to the total output uncertainty, revealing where efforts to improve the accuracy of input data will yield the greatest improvement in model performance. The response of CMAQ output to input variability has been widely studied, for example, to assess the effect of changing emissions scenarios,4,5 the effect of changing emissions in combination with different chemical mechanisms,6 and the effect of varying ozone deposition to vegetation.7 The validity of such studies may be assessed by evaluating the model’s ability to reproduce known real-world behavior in a dynamic evaluation.8−10 A closely related task is to change the components of the modeling system or its inputs, for example, the meteorological model, vertical resolution, boundary conditions (BCs), and chemical mechanism, in order to find the best configuration for prediction in a particular domain.11,12 The above examples can be described as informal sensitivity analysis (SA) methods. An alternative approach is to compute © XXXX American Chemical Society
Received: Revised: Accepted: Published: A
November 21, 2016 April 24, 2017 April 26, 2017 April 26, 2017 DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology complex applications include “Stochastic Response Surface Methods”18 and “High Dimensional Model Representations”.19 All of these methods require the computer model to be run a number of times at different input settings to provide training data to build the meta-model. The Taylor series expansions of DDM sensitivity coefficients mentioned above have been used as meta-models of CMAQ20,21 but differ in that only one or two runs are required. This, however, inherently means that considerable extrapolations are involved in their construction. Kriging metamodels of deterministic computer models were first proposed by Sacks et al.,22 who observed that the Kriging predictor can be made to exactly interpolate the output of the training runs, mirroring the way the computer model produces identical output if run again at the same input settings. This “Gaussian process emulation” method was developed further in a Bayesian framework, with the aim of incorporating observational data to calibrate the model input factors by Kennedy and O’Hagan.23 The same statistical framework can be used for global sensitivity analysis24 and, in theory, allows properties of the model output uncertainty distribution to be derived analytically, but Oakley and O’Hagan25 concede that this is best done by using the emulator as a model surrogate in a Monte Carlo analysis. Most examples of the use of emulation in the environmental sciences are found in climate modeling, with a number of case studies described by Kennedy et al.26 Gaussian process emulators of several climate models have been used in a history matching process with observational data,27−29 which calibrates model inputs in a conservative manner to reduce their uncertainty ranges, and Sexton et al.30 used emulators to attempt the more ambitious task of producing calibrated probabilistic climate projections. When perturbing inputs simultaneously in a global SA, the sensitivity to each can be quantified by estimating how much the associated output variance would be reduced if that input were fixed.13 Various techniques for calculating this are described by Saltelli et al.,13 all requiring considerable numbers of model runs, so the use of emulators provides obvious benefits. One such method, the Fourier Amplitude Sensitivity Test (FAST), was used by Lee et al.31 with a global aerosol model. Modeled concentration of cloud condensation nuclei was emulated with respect to eight uncertain parameters with a separate emulator and FAST for each grid cell. Lee et al.32 extended this study to include 28 input variables, this time also using the emulators for Monte Carlo uncertainty estimation. The work presented here is similar to that carried out by Lee et al.,32 in that emulation is used in conjunction with FAST to allow global SA, but differs in that a screening method was used prior to emulation in order to select which of the hundreds of uncertain inputs in the CMAQ model would be included in the analysis.
Figure 1. 81 km (green), 27 km (red), and 9 km (blue) grid cell domains used in the CMAQ model run.
km, an intermediate domain with a grid size of 27 km, and an inner domain producing the output presented here, covering the UK with a grid size of 9 km. The model was used in off-line mode, with the Weather Research and Forecasting model34 providing meteorological inputs. BCs for the outer domain were derived from data provided as part of the Air Quality Model Evaluation International Initiative,35 emissions for the outer and middle domains from the European Monitoring and Evaluation Programme (EMEP, http://www.ceip.at/), and emissions for the inner domain from the UK National Atmospheric Emissions Inventory (NAEI).36 A full description of this base-case run, along with a comprehensive performance evaluation, is contained in the DEFRA regional model evaluation report.37 Input Variable Screening. Typically, in complex models, only a small subset of inputs have significant influence on a particular output,13 so screening the inputs to identify this subset makes subsequent analyses more tractable. Following the method described by Morris,38 the inputs were perturbed one at a time within a range of one-half to double their baseline values, with perturbations of magnitude 0.75 times that of the baseline value. This was repeated several times for each input with a different starting point on each occasion, and the remainder of the inputs were held at different constant values on each occasion, producing a number of “elementary effects” on the model output for each input, the means and standard deviations of which together indicate the overall importance and linearity/monotonicity of that input. More detail on this method is given in pages S2−S3 of the Supporting Information. For the emissions inputs, which take the form of spatiotemporally varying fields, the whole field was perturbed by the same amount in order to keep the analysis tractable. Reaction rate constants were perturbed by changing the preexponential factor, allowing normal spatiotemporal variation due to temperature fluctuations across the domain, and photolysis rates were perturbed before the model modified them according to cloud cover.
■
METHODS Base-Case Model Run. July 2006 was modeled with CMAQ 4.7.1, using the Carbon Bond 5 (CB-05) chemical mechanism. 33 Newer CMAQ releases exist, but ozone chemistry is largely unchanged and the CB-05 mechanism is still widely used. A 10 day spin-up allowed the SA to be performed for 21 days covering a significant ozone pollution episode and several days on either side. This was a typical UK summer pollution episode with hot dry weather and anticyclonic conditions bringing air masses from central Europe. Figure 1 shows the outer spatial domain with a grid size of 81 B
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology
equating to 18.6 runs per variable. Accuracy of emulation was evaluated with leave-one-out cross validation (LOOCV), where the output of each of the training runs is predicted using an emulator trained on the remaining runs. The mean absolute error (MAE) of the predictions was calculated as follows,
A useful addition to the method, suggested by Campolongo et al.,39 uses the mean of the absolute values of the elementary effects, μ*, as a single metric to rank the inputs in order of importance. In the same paper, Campolongo et al.39 enhance the method to improve its space filling properties. This space filling method is implemented in the function morris of the package sensitivity40 in the R statistical computing language41 and is used in the screening exercises described here. It was expected that the results of this screening would differ across the domain, which contains 9360 grid cells, but performing the screening for all of these cells would have been prohibitive; so, a subset was chosen comprising the 22 grid cells containing the monitoring sites used in the DEFRA regional model evaluation.37 A total of 223 inputs were screened: 40 gas-phase species from the BCs and 27 from the emissions and the 156 reactions from the CB-05 core chemical mechanism. This was performed at 4 p.m. and midnight, on both the 10th of May and 19th of July (during the same study, some SA was also performed for May 2006) with a ten day spin-up. In this total of 188 screening exercises, the highest μ* value attained by each input factor was used to discard those factors which fell below a threshold of 2 ppb influence on both ozone and NO2 output. In such a screening exercise, it is fortuitous if a natural partitioning appears between influential and noninfluential inputs. Unfortunately, this was not the case, and the 2 ppb threshold was chosen somewhat arbitrarily, although it turned out to be justified, as we show in the Results section. At the end of this process, 12 reactions were retained, along with CO and ozone BCs, and eight emissions variables. In the analyses which follow, the emissions inputs for the 9 km domain were perturbed separately from those for the 81 and 27 km domains, and ozone deposition velocity was added to make a total of 31 variables. They are listed in the section describing the FAST methods, where their uncertainty distributions are also assigned. Gaussian Process Emulation. Emulators were constructed using the R DiceKriging package, described by Roustant et al.,42 who also give advice on its use. A great deal of information is also available in the references cited in the Introduction to this paper and in pages S4−S8 of the Supporting Information, which include brief mathematical details. The package was used with default settings which implement universal Kriging with a Matérn 5/2 covariance function. The input values at which to make the CMAQ training runs were chosen using a Latin hypercube sampling43 (LHS) plan. LHS is a stratified sampling technique, having desirable properties for Monte Carlo uncertainty analysis, but here is simply used to efficiently sample the input variable space without introducing correlations between inputs. This provided factors by which to multiply input variables, sampling each from one-half to double its baseline value. This range was larger than the ranges of the input uncertainties to be used in the sensitivity analyses and so provided scope to widen them or to perform other experiments as required. A separate emulator was made for ozone or NO2 concentration at every location and time step used in the analyes which follow. Emulation gives a point prediction of the model output at untried input values and a variance which describes the uncertainty in the emulation. The aim was to construct emulators which were accurate enough to allow the point predictions alone to be used with confidence in place of CMAQ. All of the results presented use emulators with 31 input variables, made using Latin hypercubes with 576 points,
MAE =
1 n
n
∑ |yi ̂ − yi | i=1
(1)
where y is the output of the left-out model run and ŷ is the value of that output predicted by the emulator. This gives an indication of the average amount by which an emulator will be in error when predicting CMAQ output at untried input values. The MAEs for emulators of ozone and NO2 concentrations were calculated for 125 randomly chosen time steps and locations. The distributions of these MAE statistics are shown in Figure S2. The median MAE for ozone emulators is around 1.3 ppb, and three-quarters have MAEs below 1.9 ppb, while for NO2 over three-quarters of the emulators have a MAE of below 0.5 ppb. These results were considered adequate to proceed with using the emulators for the analyses which follow. FAST. A variance-based global SA entails decomposing the model output variance into a series of terms of increasing dimensionality and dividing each of those terms by the total variance.44 Such a series for a model with k inputs would have k terms giving the first-order effects, or “main effects”, for each input, the next (k2) terms giving the second-order effects, those involving combinations of two inputs, and so on. Each of these effects can be thought of as the amount by which the output variance would be reduced if the inputs involved in that term were fixed. The number of terms in this decomposition is 2k − 1 and hence becomes prohibitive to calculate. Instead, a “total effect” index may be used, which is the sum of all the terms in the series involving a particular factor, equating to the main effect plus interactions. A clear exposition of the mathematical details can be found in Saltelli et al.,45 and further information is also given in pages S8−S10 of the Supporting Information. The Fourier Amplitude Sensitivity Test46 (FAST) provides a computationally efficient procedure for estimating the main effects by relating the probability distribution representing the uncertainty in each input to a different frequency. The input variable space is then explored by a search curve which traverses each dimension at a rate proportional to the frequency assigned to the corresponding input. This curve defines a set of values of the inputs at which to run the model, which induce a periodicity in the output. Fourier analysis of this output then gives Fourier coefficients from which the main effect for each factor can be derived. The extended FAST47 estimates main and total effects indices at the same time by choosing two frequencies for each input, one of which is assigned to that input, as in the classic FAST, and the other to all the remaining inputs. This induces two periodicities in the model output, the first of which can be used to calculate the main effect, and the second to calculate the total effect. More detail is given on pages S11−S12 of the Supporting Information, and the extended FAST is implemented in the function fast99 of the R package sensitivity.40 In order for these indices to correctly represent the proportion of output variance attributable to each input, those inputs should be independent, and while an analysis which takes correlations between them into account is theoretically possible, the extra complexity involved means that a first SA is generally performed with the assumption of C
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology independent inputs.48 While there is a likelihood that some of the CMAQ inputs are not independent, for example, NO and NO2 emissions from road transport, this does not mean that there are easily quantifiable correlations between the uncertainties in those inputs which could be included in an analysis. The variables retained following the screening process are given in Table 1, along with their uncertainty distributions,
which are centered on 1.0, so that values sampled from them equate to a factor by which to multiply the input. The reactions are labeled following the CB05 definition,33 and rate uncertainties are taken from the NASA/Jet Propulsion Laboratory Atmospheric Chemical Data Evaluation.49 Emissions uncertainties are from a NAEI uncertainty estimation50 which gives two estimates created using different methods. The larger estimate of ±30% for NOx and VOCs is used here, as there are major assumptions in the smaller estimate of ±7− 10%, such as the conjecture that official statistics are “subject to very limited uncertainty” and that missing sources of emissions would not contribute significantly to the total. Each of the listed emissions species was treated as a separate input for the 9 km domain (labeled “UK” in the results) and as another input for the 81 and 27 km domains taken together (labeled “EU” in the results). In the absence of available information, the same level of uncertainty was assigned to BCs and to ozone deposition velocity as to emissions. These uncertainty distributions were used in a separate FAST, with a separate emulator of either ozone or NO2 output, for each of the 504 hourly time steps in the 21 days from the 11th to the 31st of July 2006 at a number of locations. Brief descriptions of these locations are given in pages S13−S15 of the Supporting Information.
Table 1. Distributions Used to Characterize Input Uncertaintiesa
a
reaction
distribution
R1: NO2 photolysis R3: O3 + NO R7: NO2 + O3 R9: O3 photolysis R10: O1D + M R11: O1D + H2O R28: NO2 + OH R30: HO2 + NO R66: OH + CH4 R74: HCHO photolysis R87: C2O3 + NO R112: PAR + OH NO emissions NO2 emissions ISOP emissions PAR emissions XYL emissions ETH emissions CO emissions OLE emissions O3 BCs CO BCs O3 deposition velocity
normal, sd = 0.2 normal, sd = 0.1 normal, sd = 0.15 normal, sd = 0.3 normal, sd = 0.5 normal, sd = 0.08 normal, sd = 0.4 normal, sd = 0.3 normal, sd = 0.1 normal, sd = 0.4 normal, sd = 0.5 normal, sd = 0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3 uniform, ±0.3
■
RESULTS The base case modeled and measured ozone concentrations51 at the grid squares containing the Harwell monitoring site in Oxfordshire and the North Kensington site in Greater London for the 11th to the 31st of July are shown in Figure 2. Overall, the model captures the observed behavior quite well but fails to reproduce some of the afternoon peaks in both magnitude and timing. The model performance appears slightly better at Harwell than London, possibly due to the fact that Harwell is a rural background monitoring location, sited so as to be representative of a large area, and so more comparable to the CMAQ 9 km grid squares than an urban monitoring site. Figure 3 shows total effects at the above locations as time series of stacked bars, as suggested by Saltelli et al.,52 who
All centered on 1.0.
Figure 2. Measured and modeled ozone concentrations at (a) Harwell and (b) London, the 11th to the 31st of July 2006. D
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology
Figure 3. Time series of total effects on modeled ozone concentrations at (a) Harwell and (b) London, the 11th to the 31st of July.
more important, and these are likely to be dominated by local traffic emissions. What is also clear is that deposition velocity is less important in this urban location than in the rural background at Harwell. Dry deposition to vegetation is an important ozone sink in the CMAQ model;53 thus, an urban grid square with less vegetation than a rural area has a smaller absolute amount of deposition, and uncertainty in the deposition velocity of this smaller amount is less influential on the model output. The only VOC emission appearing on these plots (and hence the only one accounting for more than 1% of the variance) is isoprene, of which biogenic emissions increase markedly in the kind of hot weather conditions54 experienced during this episode. Isoprene emissions are widespread across rural areas of the UK, and isoprene has strong ozone creation potential.55 Ethene, xylene, and carbon bond species, of which anthropogenic emissions comprise a large part, did come through the screening process but were not signinficant in the final analyisis, providing an illustration of the nature of global SA. It does not necessarily mean that the model is insensitive to these inputs in isolation, but when many inputs are perturbed at the same time to which the model is more sensitive, their influence is effectively “drowned out”, indicating that the accuracy of those other inputs is more important in determining the accuracy of the model output. A further noteworthy feature in Figure 3b is the more “spiky” nature of the top surface of the plot, indicating that interactions between input uncertainties become important for short, isolated periods of time. Closer inspection reveals that these points are shortly before midnight in every case, and given that this would also be shortly after sunset at this time of year, these interactions could be related to the transition from the daytime photochemical regime to night time chemistry. One of these periods is examined in more detail in Figure 4, where the overall height of the bars represents the total effect as the sum
normalize the sensitivity indices so the height of the bars represent a true fraction of the total variance. Here, the indices are not normalized, meaning that the total height of the bars gives an indication of the size of the input interaction effects, as the effect of an interaction between several inputs is included in the total effect of each of them. The FASTs were performed with all 31 variables identified in the screening process, but only those accounting for more than 1% of the variance at any time step are shown on the time series plots. For example, CO boundary conditions were retained after screening but accounted for less than 1% of the variance in the FASTs. These plots and those for five additional locations in the Supporting Information have between 16 and 18 variables, demonstrating that the screening threshold was not too high; it is better to allow some variables through the screening, which later turn out to be unimportant, than to set the threshold too high and risk ruling out potentially important inputs. Of particular note on the plot for Harwell in Figure 3a is the dominance of ozone BCs, shown in blue, at the start and end of the period, when wind speeds were relatively high in comparison to those in the middle of the period when deposition velocity, in green, becomes more important. Also, as the weather became warmer in this middle period chemical processes start to dominate; BCs are important when there is little photochemical activity and vice versa. Diurnal variation in sensitivity to NO2 photolysis, in yellow at the bottom of the plot, is clearly evident. It is important to note that this does not go to zero, even at night, because although the value of the total effect gives the sensitivity at just that particular time step, this sensitivity is to uncertainty in the rate for the whole period, indicating that modeled ozone concentration at night is still sensitive to photolysis occurring earlier in the day. Figure 3b shows the results of the SA for London. In comparison to Harwell, uncertainty in NO emissions is much E
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology
NO, NO2, and ozone is such that those emissions rapidly decline in influence with distance from the source. In common with the ozone plot for London in Figure 3b, the London NO2 plot shows more interactions than the plots for the other areas but this time with a less strict temporal pattern. There are some similarities, however, with both plots showing peaks in interactions shortly before midnight on the 11th and 12th which are associated with sharp peaks in the influence of reaction rates.
■
DISCUSSION This analysis has revealed a complex spatiotemporal pattern of the sensitivity of modeled ozone concentrations to input uncertainties, showing that the most effective strategy for improving prediction accuracy varies greatly by time and place. Rural background predictions may be improved for significant periods of time with more accurate BCs; but, during ozone episodes, particularly in urban areas, this will have little effect, and improved NO emissions and reaction rate data will be of greater benefit. Photolysis reactions are among the most important; thus, accurate modeling of solar radiation and cloud cover will also help to improve predictions. Studies for US domains using the DDM to calculate model sensitivities have also shown that CMAQ ozone concentrations are more sensitive to emissions of NOx than VOCs.16,21 Also using DDM-based methods in a US domain, Napelenok et al.56 find as we do that ozone sensitivity to NOx and VOC emissions and ozone BCs has high temporal variability throughout their study period. They find greater sensitivity to VOC emissions than we do, despite using the same chemical mechanism. This could be due to differences in absolute values of emissions in different inventories but also may be attributable to the fact that the global SA perturbs many more variables at the same time; the VOC sensitivity is effectively drowned out by the fact that the model is more sensitive to so many other variables. Studies of sensitivity to the modeled chemistry have tended to focus on switching chemical mechanisms, for example,11,57 rather than perturbing individual reaction rates, so are not comparable to the results presented here. The key difference in the methods used here to other SA methods used in air quality modeling is that perturbing many inputs simultaneously gives a more reliable estimate of which has the most influence on the model output at any given time. The NO2 analysis is dominated by uncertainty in NO emissions, except for remote locations where NO2 pollution is not generally a problem, so improving the accuracy of emissions inventories would appear to be the only way of significantly improving NO2 predictions, at least when only the variables which have been included in this analysis are considered. However, it is likely that errors in meteorology play a part in driving uncertainty here and are also a significant driver of ozone prediction errors. Meteorological inputs have not been included and to do so would represent significant progress. This presents specific difficulties, however, such as violating conservation of mass between the model grid cells when perturbing variables such as wind speed and direction. It has been suggested that varying wind speed uniformly across the whole domain can alleviate this problem.58 However, there are other variables which are deterministically linked to wind speed, such as planetary boundary layer height, so simply changing the wind speed without a means of modifying those variables accordingly would lead to modeling a physically unrealistic scenario.
Figure 4. Main effects and interactions for ozone sensitivity in London, 10 pm on the 11th of July.
of the main effect and interactions. Uncertainty in UK NO emissions can be seen to be dominant, and the high degree of interaction with the other factors is clearly evident. For some inputs, NO2 photolysis and ozone BCs, for example, the main effect is a small proportion of the total effect. This means that uncertainty in that input by itself is not causing a great deal of uncertainty in the model output; it is only when there is also uncertainty in other inputs with which it interacts that it is able to induce considerable uncertainty in modeled ozone concentration. Five other locations are shown in the Supporting Information, and the progression through the Figures S4 to S8 is from urban to gradually more rural locations, finishing at Strath Vaich, a remote site in Northern Scotland. Here, uncertainty in BCs is clearly important, and emissions have little influence on the model output. The quenching of O(1D), shown in red, is an influential input, as is ozone photolysis immediately below and in combination with NO2 photolysis; it can be seen that the chemistry normally associated with the tropospheric background is driving the model behavior at this location. These inputs are also influential at the next most remote location, Yarner wood in the Dartmoor national park, and moving through the other less remote rural sites and then into Manchster Picadilly and London, the transition from the tropospheric background regime to the urban regime is apparent. Nitrogen Dioxide. Plots showing FAST time series for NO2 output at London, Harwell, and Strath Vaich are shown in Figures S9−S11. In all cases, the analyses are dominated by UK NO emissions and appear less interesting than the ozone results, but there are a few points worthy of note. Again, Strath Vaich is greatly influenced by the chemistry we would expect in the tropospheric background, and emissions become more important moving toward urban areas. The influence of NO2 emissions is generally small but is largest in London plot and nonexistent at Strath Vaich. The emissions inputs only have significant amounts of NO2 in major conurbations and transport routes, and the speed of the reactions which cycle F
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology
(4) Mueller, S. F.; Mallard, J. W. Contributions of Natural Emissions to Ozone and PM(2.5) as Simulated by the Community Multiscale Air Quality (CMAQ) Model (vol 45, pg 4817, 2011). Environ. Sci. Technol. 2011, 45, 7950. (5) Yu, Y.; Sokhi, R. S.; Kitwiroon, N.; Middleton, D. R.; Fisher, B. Performance characteristics of MM5-SMOKE-CMAQ for a summer photochemical episode in southeast England, United Kingdom. Atmos. Environ. 2008, 42, 4870−4883. (6) Arnold, J. R.; Dennis, R. L. Testing CMAQ chemistry sensitivities in base case and emissions control runs at SEARCH and SOS99 surface sites in the southeastern US. Atmos. Environ. 2006, 40, 5027− 5040. (7) Emberson, L. D.; Kitwiroon, N.; Beevers, S.; Büker, P.; Cinderby, S. Scorched Earth: how will changes in the strength of the vegetation sink to ozone deposition affect human health and ecosystems? Atmos. Chem. Phys. 2013, 13, 6741−6755. (8) Pierce, T.; Hogrefe, C.; Trivikrama Rao, S.; Porter, P. S.; Ku, J.-Y. Dynamic evaluation of a regional air quality model: Assessing the emissions-induced weekly ozone cycle. Atmos. Environ. 2010, 44, 3583−3596. (9) Gilliland, A. B.; Hogrefe, C.; Pinder, R. W.; Godowitch, J. M.; Foley, K. L.; Rao, S. Dynamic evaluation of regional air quality models: Assessing changes in O3 stemming from changes in emissions and meteorology. Atmos. Environ. 2008, 42, 5110−5123. (10) Foley, K. M.; Hogrefe, C.; Pouliot, G.; Possiel, N.; Roselle, S. J.; Simon, H.; Timin, B. Dynamic evaluation of CMAQ part I: Separating the effects of changing emissions and changing meteorology on ozone levels between 2002 and 2005 in the eastern US. Atmos. Environ. 2015, 103, 247−255. (11) Appel, K. W.; Gilliland, A. B.; Sarwar, G.; Gilliam, R. C. Evaluation of the Community Multiscale Air Quality (CMAQ) model version 4.5: Sensitivities impacting model performance Part I - Ozone. Atmos. Environ. 2007, 41, 9603−9615. (12) Appel, K. W.; Roselle, S. J.; Gilliam, R. C.; Pleim, J. E. Sensitivity of the Community Multiscale Air Quality (CMAQ) model v4.7 results for the eastern United States to MM5 and WRF meteorological drivers. Geosci. Model Dev. 2010, 3, 169−188. (13) Saltelli, A.; Chan, K.; Scott, E. Sensitivity Analysis; Wiley: Chichester, UK, 2008. (14) Dunker, A. M. The decoupled direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 1984, 81, 2385−2393. (15) Hakami, A.; Odman, M. T.; Russell, A. G. High-order, direct sensitivity analysis of multidimensional air quality models. Environ. Sci. Technol. 2003, 37, 2442−52. (16) 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−48. (17) Fang, K.; Li, R.; Sudjianto, A. Design and modeling for computer experiments; Chapman & Hall/CRC: Boca Raton, FL, 2006. (18) Isukapalli, S. S.; Roy, a.; Georgopoulos, P. G. Stochastic response surface methods (SRSMs) for uncertainty propagation: application to environmental and biological systems. Risk analysis: an official publication of the Society for Risk Analysis 1998, 18, 351−63. (19) Li, G. Y.; Rosenthal, C.; Rabitz, H. High dimensional model representations. J. Phys. Chem. A 2001, 105, 7765−7777. (20) Pinder, R. W.; Gilliam, R. C.; Appel, K. W.; Napelenok, S. L.; Foley, K. M.; Gilliland, A. B. Efficient Probabilistic Estimates of Surface Ozone Concentration Using an Ensemble of Model Configurations and Direct Sensitivity Calculations. Environ. Sci. Technol. 2009, 43, 2388−2393. (21) Tian, D.; Cohan, D. S.; Napelenok, S.; Bergin, M.; Hu, Y.; Chang, M.; Russell, A. G. Uncertainty Analysis of Ozone Formation and Response to Emission Controls Using Higher-Order Sensitivities. J. Air Waste Manage. Assoc. 2010, 60, 797−804. (22) Sacks, J.; Welch, W.; Mitchell, T.; Wynn, H. Design and analysis of computer experiments. Statistical science 1989, 4, 409−423.
An alternative proposition would be to remove meteorological variability from the analysis to some degree by grouping the results at each location by back trajectory and general synoptic conditions and only including time steps when modeled meteorology was in good agreement with observations. It may also be useful to group some inputs, such as anthropogenic VOCs, to make the results easier to interpret and more policy relevant. This analysis could not have been achieved without the use of emulation as CMAQ is simply too computationally expensive to allow sufficient numbers of runs to be made; over 2.5 million emulated “runs” were required to produce each of the time series plots in the Results and Supporting Information. The model is one of the most complex of those which have been emulated in the literature, and the emulators are among the largest in terms of number of inputs. Despite this, the LOOCV demonstrates good accuracy in emulation, showing that global SA of complex air qualty models is a feasible proposition. Readers interested in using these tehniques in their own work should consult references cited earlier, in particular Campolongo et al.,39 Pujol et al.,40 Roustant et al.42 and Saltelli et al.45,47 The extra information on the methods which can be found in pages S2−S11 of the Supporting Information is intended as a readable entry point into this literature.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.6b05873. Additional details of the methods; descriptions of locations referred to in the Article; time series plots for NO2 SA; additional locations for ozone (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]; phone: +44 (0)20 7848 4009; fax: +44 (0)20 7848 4047. ORCID
Andrew V. Beddows: 0000-0002-7455-4087 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by a MRC-PHE Centre for Environment and Health PhD studentship awarded to the first author who also thanks Dr. Marta Blangiardo of Imperial College London for advice during the course of the project.
■
REFERENCES
(1) Review of evidence on health aspects of air pollution; WHO European Centre for Environment and Health: Copenhagen, Denmark, 2013; http://www.euro.who.int/en/health-topics/environmentand-health/air-quality. (2) Byun, D. W.; Schere, K. L. Review of the governing equations, computational algorithms, and other components of the models-3 Community Multiscale Air Quality (CMAQ) modeling system. Appl. Mech. Rev. 2006, 59, 51−77. (3) Denby, B.; Dudek, A.; Walker, S. E.; Costa, A. M.; Monteiro, A.; van den Elshout, S.; Fisher, B. Towards uncertainty mapping in airquality modelling and assessment. Int. J. Environ. Pollut. 2011, 44, 14− 23. G
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology
Transl.) 1993, 1 407−414; Matematicheskoe Modelirovanie 1990, 2, 112−118. (45) Saltelli, A.; Annoni, P.; Azzini, I.; Campolongo, F.; Ratto, M.; Tarantola, S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput. Phys. Commun. 2010, 181, 259−270. (46) Cukier, R. I.; Levine, H. B.; Shuler, K. E. Nonlinear sensitivity analysis of multiparameter model systems. J. Comput. Phys. 1978, 26, 1−42. (47) Saltelli, A.; Tarantola, S.; Chan, K. A quantitative modelindependent method for global sensitivity analysis of model output. Technometrics 1999, 41, 39−56. (48) Saltelli, A.; Tarantola, S.; Campolongo, F.; Ratto, M. Sensitivity Analysis in Practice; Wiley: Chichester, UK, 2004. (49) Sander, S.; Abbatt, J.; Barker, J.; Burkholder, J.; Friedl, R.; Golden, D.; Huie, R.; Kolb, C.; Kurylo, M.; Moortgat, G.; Orkin, V.; Wine, P. Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies, Evaluation No. 17; JPL Publication 10-6; NASA Jet Propulsion Laboratory: Pasadena, CA, 2011. (50) Estimation of Uncertainties in the National Atmospheric Emissions Inventory; AEA Technology: Abingdon, Oxfordshire, UK, 2003; http://naei.defra.gov.uk/reports. (51) DEFRA. UK-Air Data Archive; Available at: http://uk-air.defra. gov.uk/data/ [accessed 20/01/2015]. (52) Saltelli, A.; Tarantola, S.; Chad, K. Presenting Results from Model Based Studies to Decision-Makers: Can Sensitivity Analysis Be a Defogging Agent? Risk Anal. 1998, 18, 799−803. (53) Pleim, J.; Ran, L. Surface Flux Modeling for Air Quality Applications. Atmosphere 2011, 2, 271−302. (54) Laffineur, Q.; Aubinet, M.; Schoon, N.; Amelynck, C.; Müller, J.F.; Dewulf, J.; Van Langenhove, H.; Steppe, K.; Šimpraga, M.; Heinesch, B. Isoprene and monoterpene emissions from a mixed temperate forest. Atmos. Environ. 2011, 45, 3157−3168. (55) Derwent, R.; Hov, Ø. Application of sensitivity and uncertainty analysis techniques to a photochemical ozone model. J. Geophys. Res. 1988, 93, 5185. (56) Napelenok, S. L.; Foley, K. M.; Kang, D.; Mathur, R.; Pierce, T.; Rao, S. T. Dynamic evaluation of regional air quality model’s response to emission reductions in the presence of uncertain emission inventories. Atmos. Environ. 2011, 45, 4091−4098. (57) Luecken, D.; Phillips, S.; Sarwar, G.; Jang, C. Effects of using the CB05 vs. SAPRC99 vs. CB4 chemical mechanism on model predictions: Ozone and gas-phase photochemical precursor concentrations. Atmos. Environ. 2008, 42, 5805−5820. (58) Hanna, S. R.; Lu, Z.; Frey, H. C.; Wheeler, N.; Vukovich, 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.
(23) Kennedy, M. C.; O’Hagan, A. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2001, 63, 425−464. (24) Oakley, J. E.; O’Hagan, A. Probabilistic sensitivity analysis of complex models: a Bayesian approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2004, 66, 751−769. (25) Oakley, J.; O’Hagan, A. Bayesian inference for the uncertainty distribution of computer model outputs. Biometrika 2002, 89, 769− 784. (26) Kennedy, M. C.; Anderson, C. W.; Conti, S.; O’Hagan, A. Case studies in Gaussian process modelling of computer codes. Reliability Engineering & System Safety 2006, 91, 1301−1309. (27) Edwards, N. R.; Cameron, D.; Rougier, J. Precalibrating an intermediate complexity climate model. Climate Dynamics 2011, 37, 1469−1482. (28) McNeall, D. J.; Challenor, P. G.; Gattiker, J. R.; Stone, E. J. The potential of an observational data set for calibration of a computationally expensive computer model. Geosci. Model Dev. 2013, 6, 1715− 1728. (29) Williamson, D.; Goldstein, M.; Allison, L.; Blaker, A.; Challenor, P.; Jackson, L.; Yamazaki, K. History matching for exploring and reducing climate model parameter space using observations and a large perturbed physics ensemble. Climate Dynamics 2013, 41, 1703−1729. (30) Sexton, D. M. H.; Murphy, J. M.; Collins, M.; Webb, M. J. Multivariate probabilistic projections using imperfect climate models part I: outline of methodology. Climate Dynamics 2012, 38, 2513− 2542. (31) Lee, L. A.; Carslaw, K. S.; Pringle, K. J.; Mann, G. W. Mapping the uncertainty in global CCN using emulation. Atmos. Chem. Phys. 2012, 12, 9739−9751. (32) Lee, L. A.; Pringle, K. J.; Reddington, C. L.; Mann, G. W.; Stier, P.; Spracklen, D. V.; Pierce, J. R.; Carslaw, K. S. The magnitude and causes of uncertainty in global model simulations of cloud condensation nuclei. Atmos. Chem. Phys. 2013, 13, 8879−8914. (33) Updates to the Carbon Bond Chemical Mechanism: CB05; ENVIRON: Novato, CA, 2005; http://www.camx.com/publications/ model-development.aspx. (34) Skamarock, W. C.; Klemp, J. B. A time-split nonhydrostatic atmospheric model for weather research and forecasting applications. J. Comput. Phys. 2008, 227, 3465−3485. (35) Rao, S. T.; Galmarini, S.; Puckett, K. Air Quality Model Evaluation International Initiative (AQMEII): Advancing the State of the Science in Regional Photochemical Modeling and Its Applications. Bull. Am. Meteorol. Soc. 2011, 92, 23−30. (36) UK Emissions of Air Pollutants 1970 to 2008; AEA Technology: Didcot, Oxfordshire, UK, 2010; http://naei.defra.gov.uk/reports. (37) DEFRA Phase 2 Regional Model Evaluation; King’s College London: London, UK, 2013; http://uk-air.defra.gov.uk/library. (38) Morris, M. D. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics 1991, 33, 161−174. (39) Campolongo, F.; Cariboni, J.; Saltelli, A. An effective screening design for sensitivity analysis of large models. Environmental Modelling & Software 2007, 22, 1509−1518. (40) Pujol, G.; Iooss, B. with contributions from Laurent Gilquin, A. J.; Gratiet, L. L.; Lemaitre, P. Sensitivity: Sensitivity Analysis; 2014; R package version 1.8-1; http://CRAN.R-project.org/. (41) R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2014; http://www.R-project.org/. (42) Roustant, O.; Ginsbourger, D.; Deville, Y. DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization. Journal of Statistical Software 2012, 51, 1−55. (43) McKay, M.; Beckman, R.; Conover, W. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 1979, 21, 239−245. (44) Sobol’, I. On sensitivity estimation for nonlinear mathematical models. Mathematical Modeling and Computational Experiments (Engl. H
DOI: 10.1021/acs.est.6b05873 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
