Open Access

Abstract. Multi-year objective analyses (OA) on a high spatiotemporal resolution for the warm season period (1 May to 31 October) for ground-level ozone and for fine particulate matter (diameter less than 2.5 microns (PM2.5)) are presented. The OA used in this study combines model outputs from the Canadian air quality forecast suite with US and Canadian observations from various air quality surface monitoring networks. The analyses are based on an optimal interpolation (OI) with capabilities for adaptive error statistics for ozone and PM2.5 and an explicit bias correction scheme for the PM2.5 analyses. The estimation of error statistics has been computed using a modified version of the Hollingsworth–Lonnberg (H–L) method. The error statistics are "tuned" using a χ2 (chi-square) diagnostic, a semi-empirical procedure that provides significantly better verification than without tuning. Successful cross-validation experiments were performed with an OA setup using 90% of data observations to build the objective analyses and with the remainder left out as an independent set of data for verification purposes. Furthermore, comparisons with other external sources of information (global models and PM2.5 satellite surface-derived or ground-based measurements) show reasonable agreement. The multi-year analyses obtained provide relatively high precision with an absolute yearly averaged systematic error of less than 0.6 ppbv (parts per billion by volume) and 0.7 μg m−3 (micrograms per cubic meter) for ozone and PM2.5, respectively, and a random error generally less than 9 ppbv for ozone and under 12 μg m−3 for PM2.5. This paper focuses on two applications: (1) presenting long-term averages of OA and analysis increments as a form of summer climatology; and (2) analyzing long-term (decadal) trends and inter-annual fluctuations using OA outputs. The results show that high percentiles of ozone and PM2.5 were both following a general decreasing trend in North America, with the eastern part of the United States showing the most widespread decrease, likely due to more effective pollution controls. Some locations, however, exhibited an increasing trend in the mean ozone and PM2.5, such as the northwestern part of North America (northwest US and Alberta). Conversely, the low percentiles are generally rising for ozone, which may be linked to the intercontinental transport of increased emissions from emerging countries. After removing the decadal trend, the inter-annual fluctuations of the high percentiles are largely explained by the temperature fluctuations for ozone and to a lesser extent by precipitation fluctuations for PM2.5. More interesting is the economic short-term change (as expressed by the variation of the US gross domestic product growth rate), which explains 37% of the total variance of inter-annual fluctuations of PM2.5 and 15% in the case of ozone.

Abstract.Multi-year objective analyses (OA) on a high spatiotemporal resolution for the warm season period (1 May to 31 October) for ground-level ozone and for fine particulate matter (diameter less than 2.5 microns (PM 2.5 )) are presented.The OA used in this study combines model outputs from the Canadian air quality forecast suite with US and Canadian observations from various air quality surface monitoring networks.The analyses are based on an optimal interpolation (OI) with capabilities for adaptive error statistics for ozone and PM 2.5 and an explicit bias correction scheme for the PM 2.5 analyses.The estimation of error statistics has been computed using a modified version of the Hollingsworth-Lönnberg (H-L) method.The error statistics are "tuned" using a χ 2 (chi-square) diagnostic, a semiempirical procedure that provides significantly better verification than without tuning.Successful cross-validation experiments were performed with an OA setup using 90 % of data observations to build the objective analyses and with the remainder left out as an independent set of data for verification purposes.Furthermore, comparisons with other external sources of information (global models and PM 2.5 satellite surface-derived or ground-based measurements) show reasonable agreement.The multi-year analyses obtained provide relatively high precision with an absolute yearly averaged systematic error of less than 0.6 ppbv (parts per billion by volume) and 0.7 µg m −3 (micrograms per cubic meter) for ozone and PM 2.5 , respectively, and a random error generally less than 9 ppbv for ozone and under 12 µg m −3 for PM 2.5 .This paper focuses on two applications: (1) presenting longterm averages of OA and analysis increments as a form of summer climatology; and (2) analyzing long-term (decadal) trends and inter-annual fluctuations using OA outputs.The results show that high percentiles of ozone and PM 2.5 were both following a general decreasing trend in North America, with the eastern part of the United States showing the most widespread decrease, likely due to more effective pollution controls.Some locations, however, exhibited an increasing trend in the mean ozone and PM 2.5 , such as the northwestern part of North America (northwest US and Alberta).Conversely, the low percentiles are generally rising for ozone, which may be linked to the intercontinental transport of increased emissions from emerging countries.After removing the decadal trend, the inter-annual fluctuations of the high percentiles are largely explained by the temperature fluctuations for ozone and to a lesser extent by precipitation fluctuations for PM 2.5 .More interesting is the economic short-term change (as expressed by the variation of the US gross domestic product growth rate), which explains 37 % of the total variance of inter-annual fluctuations of PM 2.5 and 15 % in the case of ozone.

Introduction
Long-term series of surface objective analyses (OA) of chemical species are valuable tools for understanding the historical evolution of pollution, providing long-term comparisons with models, building a climatology of surface pollutants, evaluating the efficiency of existing pollution control and abatement measures and regulations, computing pollution trends, and supporting epidemiological studies.Among the most important surface pollutants are ground-level ozone and PM 2.5 (particulate matter with diameter less than 2.5 microns), which are primary contributors to poor air quality in the US (EPA, 2012).These pollutants are also the main constituents of smog and, together with NO 2 , form the basis of the Canadian Air Quality Health Index (AQHI; Stieb et al., 2008).
Ozone is not directly emitted but produced by complex photochemical reactions of nitrogen oxides (NO x ) and volatile organic compounds (VOCs) from natural and anthropogenic sources.In urban areas, over one hundred chemical reactions are involved in ozone production (Jacobson, 2002;Seinfeld and Pandis, 2006).Ozone is an oxidant and is detrimental to health (Hazucha et al., 1989;Berglund et al., 1991).For example, it exacerbates asthma, especially with simultaneous presence of allergenic pollen (White et al., 1994;Newman-Taylor, 1995;Cashel et al., 2003).It also impacts agricultural productivity (Skärby and Selldén, 1984;Tingey et al., 1991), causes injuries and additional stress to forest ecosystems (Reich, 1987;Badot, 1989;Chappelka and Flagler, 1991) and to materials by cracking rubber and polymers (Cass, 1991).Ozone is also an important source of the hydroxyl radical, which breaks down many pollutants and certain greenhouse gases (GHG) and acts itself as an effective GHG (Jacobson, 2002;IPCC, 2007;Houghton, 2009).The atmospheric lifetime of tropospheric ozone is from a few weeks above the boundary layer (Tarasick et al., 2010) to a few hours near the surface at night (IPCC, 2007).
Another hazardous pollutant is fine particulate matter, or PM 2.5 .Primary emissions of PM 2.5 can be manmade (mostly by burning of fossil fuels in power plants or vehicles and various industrial processes) or produced naturally (volcanoes, dust storms, forest or grass fires and sea spray).Secondary formation of PM 2.5 is also an important source, originating basically from the atmospheric transformation of inorganic and organic species (Seinfeld and Pandis, 2006).PM 2.5 aerosols can cool or warm the atmosphere via interaction with incoming solar radiation (aerosol direct effect) or via their ability to act as cloud condensation or ice nuclei, and thus play a role in cloud formation (indirect effect) (Hobbs, 1993;Jacobson, 2002;IPCC, 2007;Houghton, 2009).Health impacts of PM 2.5 are also numerous, including the promotion of asthma (Newman-Taylor, 1995;Cashel et al., 2003) and other respiratory problems (ALA, 2012), the stimulation of high plaque deposits in arteries producing vascular inflammation, oxidative stress and atherosclerosis (a hardening of the arteries that reduces their elasticity), leading to heart attacks and related cardiovascular problems and therefore an increase in morbidity and mortality (Pope et al., 2002;Sun et al., 2005;Reeves, 2011).According to Pope et al. (2002), a 10 µg m −3 increase in PM 2.5 has been associated with a 6 % increase in death rate.Exposure to ambient fine particles was recently estimated to have contributed 3.2 million premature deaths worldwide in 2010 due largely to cardiovascular disease, and 223 000 deaths from lung cancer (IARC, 2013).Together, ozone and PM 2.5 can also trigger bronchial micro-lesions, which facilitate the penetration of macromolecules such as pollen augmenting the allergenic reaction (Gervais, 1994).Table 1 summarizes the main environmental and health issues related to both pollutants.Given the above evidence, it is therefore of paramount importance to provide the public and health specialists with the best information about these pollutants.
This study uses an optimal interpolation (OI) technique to produce a long-term series of ground-level ozone and PM 2.5 concentrations over a large region at a relatively low cost.Models are generally characterized by known deficiencies whereas measurement systems suffer from representativeness problems and lack of sufficient coverage, therefore providing often only local information.The OI technique, used in operational meteorology for decades, provides an optimal framework to extract the maximum information of both model and observations (Rutherford, 1972;Daley, 1991;Kalnay, 2003;Brasnett, 2008).Producing maps of objective analyses based on OI on a regular basis has numerous applications in air quality (AQ): 1. initializing numerical AQ models at regular time intervals (usually every 6 or 12 h) with appropriate fields having overall bias and error variance that are ideally minimum (Blond et al., 2004;Tombette et al., 2009;Wu et al., 2008), 2. providing users with a more accurate picture of the true state of a given chemical species by using an appropriate optimal blend of model fields together with observations so that it produces the best possible analysis (given available data), not only in the vicinity of observation points but also elsewhere in a given domain, even where the observation network does not have an optimal density, 3. building potentially useful maps of health indices (Air Quality Health Index), environmental indices or pollutant burden on ecosystems, 4. producing surface pollutant climatology, 5. computing temporal trends.
One of the key components of data assimilation, or objective analysis, is error statistics but the prescription of adequate error statistics for air quality can be challenging.Unlike the free troposphere or the stratosphere, boundary layer problems and complex topography make it difficult to produce error covariance statistics for ground pollutants such as ozone (Tilmes, 1999(Tilmes, , 2001)).Moreover, models are often imprecise over complex boundary layer surfaces.However, the relatively flat topography found over eastern and central North America and the importance of transport of ozone and  (Jacobson, 2002;IPCC, 2007).
PM 2.5 above and within the boundary layer make these pollutants excellent candidates for objective analysis and data assimilation since the correlation length is much larger than model resolution.Consequently, information can be spread around efficiently over more than one model grid point.Producing a long series of multi-year analyses retrospectively may also pose many other technical problems, causing discontinuities or inconsistencies in time series: for example, changes of model version, set of imprecise or out-of-date emissions inventories or changes in instrumentation through the years.However, in this study, efforts were made to eliminate biases in the analyses likely caused by the above factors.This study focuses on (1) obtaining optimal (adaptive) error statistics that can follow diurnal and seasonal cycles and adapt, to a certain degree, to changes over time; and (2) reducing, as much as possible, systematic errors so that the analyses form an unbiased, consistent and coherent data set throughout the whole period.Multi-year analyses presented here combine the information provided by a long series of air quality model outputs from the Canadian Air Quality Regional Deterministic Prediction System (AQRDPS), that is, the CHRONOS (Canadian and Hemispheric Regional Ozone and NO x System) model for the period 2002-2009 and GEM-MACH (Global Environmental Multi-scale coupled with Model of Air quality and CHemistry) model for the period 2010-2012.The observations used during the same period are taken from US and Canadian air quality mon-itoring surface networks.One of the main applications of OA, which is presented in this paper, is the summer climatology.Other methods to derive a multi-year climatology within the troposphere exist, such as the traditional spatial domain-filling techniques using observations and trajectories (Tarasick et al., 2010).However, uncertainties of up to 30-40 % are noted with the meteorological inputs of trajectory models (Harris et al., 2005).Moreover, in the boundary layer, complex dispersion and turbulence tend to render trajectories near the surface less precise than that of higher levels, as reported by Tarasick et al. (2010).The climatology presented here avoids the process of high uncertainties associated with back trajectories near the surface.The spatial interpolation used in this study is accomplished naturally through exponential functions (see Sect. 2.1) so that meteorological and chemical patterns from the model are preserved.Numerous studies on stratospheric ozone climatology, based on satellite observation combined with various mapping techniques, have appeared in the last decade or two (Fortuin and Kelder, 1998;McPeters et al., 2007;Ziemke et al., 2011).There have also been numerous studies on tropospheric climatology, based mostly on ozonesondes (Logan, 1999;Logan et al., 1999;Tarasick et al., 2010).However, comparatively little research has been done focusing on multi-year analyses or climatology specifically for groundlevel ozone and PM 2.5 for the whole domain of North America.The MACC (Monitoring Atmospheric Composition and Climate) reanalysis project of the global tropospheric composition (http://www.gmes-atmosphere.eu/about/project/) is one of the few initiatives to produce a long series of surface and tropospheric pollutant analyses but has a focus in Europe.RETRO-40 (REanalysis of the TROpospheric chemical composition over the past 40 yr, http://retro.enes.org) is also a project of global reanalyses initiated in Europe.In the US, efforts to combine model results and observations of surface ozone and PM 2.5 have been attempted (Berrocal et al., 2010(Berrocal et al., , 2012;;McMillan et al., 2010) using a hierarchical Bayesian space-time modeling.However, these studies have focused only on the eastern US.
Multi-year analyses of ground-level ozone and PM 2.5 presented in this paper are delivered on a comparatively high spatiotemporal resolution (15 km or 21 km -on an hourly basis) for both Canada and the US.The rest of our manuscript is organized as follows: Sect. 2 describes the methodology and theory of OA, models and observation used.Section 3 presents a cross-validation and comparison of OA results with external sources.Section 4 presents results of long-term averages of OA during the CHRONOS era (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) and the GEM-MACH era (2010-2012), from which summer climatology is derived for both periods.Section 5 introduces pollution trends and analyzes inter-annual fluctuations obtained using OA.Finally, Sect.6 provides a discussion about certain aspects of the methodology and special issues related to OI, and Sect.7 contains our summary and conclusions.

The analysis scheme
The analysis scheme used in this study is based on an optimal interpolation (OI) that utilizes two independent pieces of information: surface air quality short-term forecasts and surface observations.The OI method adopted is standard and has been described elsewhere (for a review, see Daley, 1992, chap. 4;Ménard, 2000;Kalnay, 2003, Sect. 5.4 andMénard andRobichaud, 2005).In this section, the analysis scheme is extended to include a semi-empirical adaptive scheme for error statistics designed to obtain optimal weights for model and observations.A residual bias correction for PM 2.5 is also proposed.The basic goal of an objective analysis scheme is to find an expression that minimizes the error variance of the combined field of model and observation.The background information x n f at a time n provided by the model can be written as: where x n t , ε n f are, respectively, the true value and the background error at time n.Similarly, the observation provides the following information: with ε n o including instrumental and interpolation error.The background error variance is then defined as and the observation error variance as It can be shown that the optimization problem yields the following form for the analysis vector x n a (e.g., Daley, 1991;Kalnay, 2003): where x n f is the background field vector (with length as the total number of model grid points, N, for both background and analysis vectors).The background field vector is obtained from a short-term forecast, and H is an operator that performs an interpolation from the model grid point space to the observation space (here we use a bilinear interpolation), y n o is the vector that contains all the observations at a given time n (length of this vector is equal to the number of total valid observations, e.g., nobs), and K is the optimal Kalman gain matrix to be defined below (with dimension is N*nobs) and contains information about optimal error variances (Eqs.3 and 4).The second term on the right-hand side of Eq. ( 5) is called analysis increment (Daley, 1991;Kalnay, 2003) and could be viewed as the correction to the model due to observations that brings the analyses closer to the true value.In OI, error statistics are stationary and prescribed from past experiments (rather than as time-evolving as in a Kalman filter), and error correlations are given as functions of space (rather than matrices defined on a specific grid), so there is no need to interpolate the error correlations onto the observation locations.The computation of the Kalman gain does, however, involve the inversion of a matrix (see Eq. 6).In meteorology, because of the large number of observations, this inversion is calculated in batches in smaller domains using either data selection or a compactly supported covariance function (Daley, 1991;Houtekamer and Mitchell, 1998).In air quality, the number of surface observations at a given time is limited (in North America approximately 1300 observations or less per species) and hence the inversion of the matrix can be computed directly.The analysis scheme presented in this paper is therefore equivalent to a 2-D-VAR (two-dimensional variational analysis).The derivation of the analysis equation (Eq.5) can be obtained as a best linear unbiased estimator (BLUE) or by assuming that the errors are Gaussian distributed.Other assumptions associated with Eq. ( 5) are (i) the observation errors are uncorrelated with the background errors, and (ii) interpolated observations are linearly related to the model state (e.g., the observation operator is linear).Furthermore, for OI the background error correlation is modeled as a function, generally assumed to be isotropic and homogeneous.It could be shown that in Eq. ( 5), the optimal gain matrix K has the following form (see Kalnay, 2003): where B is the background error covariance matrix defined on the model grid.But in OI, each term in Eq. ( 6) is computed as a function between pair of points.For example, for a pair of observation sites k 1 and k 2 , a typical form of the optimum interpolation has the following form (Daley, 1991;Ménard, 2000): (assuming here that the background error correlation is a homogeneous isotropic first-order auto-regressive model).The error covariance between a given site k 1 and a particular analysis grid point (i, j ) is also given explicitly as which represents the background error covariance between a given station k 1 and the nearest model grid point.B is the background error covariance matrix itself; R the observation error covariance matrix; x(k 1 ) and x(k 2 ) the position in space of the corresponding stations k 1 and k 2 ; and x(i, j ) the grid point position.Finally, σ f(j,i) and L c represent respectively the background error standard deviation and the correlation length at a specific grid point (i, j ) and are assumed in this study to be constant throughout the domain.However, σ f (k) and σ o (k) for all observation locations k are defined locally and obtained with an autocorrelation model fitting the observation-minus-forecast residuals, a procedure known as the Hollingsworth-Lönnberg (H-L) method (see Sect. 2.2).Equations ( 7) and (8) describe a first order autoregressive model which is positive definite (the product of all the terms on the right-hand side of the equation is always positive; see Daley 1991;Ménard, 2000), so that the right-hand side of Eq. ( 6) is also positive definite.In this study it should be noted that the background error covariance is univariate (e.g., no correlation between ozone and PM 2.5 ).
Unlike meteorology, where the objective analysis or data assimilation cycle of 6 h are most commonly used (Houtekamer and Mitchell, 1998;Gauthier et al., 1999;Brasnett, 2008), in air quality a cycle of 1 h is more appropriate (Blond et al., 2004;Tombette et al., 2009;Wu et al., 2008) since there is a significant diurnal variation of surface pollutants and care must be taken to resolve short or intermittent episodes.Therefore, over the entire study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) analyses were produced hourly.Only warm season (1 May-31 October) analyses are presented due to unresolved biases issues for the cold season for PM 2.5 and due to the fact that ozone and related photochemical species are less of an environmental threat in winter in most of North America.
To verify the consistency of error statistics, the innovation vector (d = y n o − Hx n f ) is computed to obtain the chi-square diagnostic in real time as follows (Ménard, 2000): where p is the number of ingested observations, and is called the innovation matrix, d is the innovation vector, and d T its transpose.In theory, the value of chi-square divided by the number of observations should be close to unity (Ménard, 2000;Ménard and Lang-Ping, 2000) for optimal error statistics.The matrix S given by Eq. ( 10) needs to be inverted only one time per analysis for the non-adaptive scheme, and, in the case of the adaptive scheme, several times until convergence is achieved.In both cases, the matrix inversion is performed using a Cholesky decomposition (Golub and Van Loan, 1996, Sect. 4.2).A potential problem with this method may exist when the matrix S is of moderate to large size (typically when the dimension is greater than 1500 × 1500).
In this case, the inversion may become inaccurate and alternative methods should be used.In this study, the maximum number of monitoring sites used for analyses is always inferior to that number (around 1100 for the period 2002-2009, progressively increasing to around 1300 in 2012).At a specific grid point, the analysis error variance σ 2 a is always smaller than both the background error variance σ 2 f and the observation error variance σ 2 o at a specific grid point (e.g., Kalnay, 2003 for a derivation), e.g., with σ 2 a defined as follows: where ε n a is the analysis error and σ 2 f and σ 2 o defined by Eqs.(3) and ( 4), respectively at a particular grid point.Note that since Eq. ( 11) was derived using the least square method theory, there is no assumption whether the distribution needs to be Gaussian or not (Kalnay, 2003).According to Eq. ( 11), mapping historical evolution of pollutants is therefore more appropriate, with an objective analysis being more accurate than model and observations, each of them taken separately.Following the same theory, the analysis error matrix A could also be derived and has the following form (see Kalnay, 2003 for a demonstration):

Adaptive error statistics
Great care should be given to the production of error statistics.Not doing so may destroy the effective optimality of an assimilation scheme (Daley, 1992;Tilmes, 2001;Sandu and Chai, 2011).Moreover, inappropriate error statistics could cause numerous rejection of valid observations in the quality control step (Gauthier et al., 2003; see also Robichaud et al., 2010).Innovations are the best sources of information to compute error statistics (Daley, 1991;Blond et al., 2004).The technique used here consists of pairing up different monitoring sites, calculating the covariance of OmP (observation minus model prediction) between the paired stations, plotting the result as a function of distance, and fitting an autoregressive correlation model as a function of distance but excluding the data at the origin.This method was originally adopted in meteorology by Gandin (1963) using climatology and by Rutherford (1972) utilizing a short-term forecast.This technique was further developed by Hollingsworth and Lönnberg (1986) and became a standard in the design of optimum interpolation.For this paper's focus on surface air quality, not all stations could be adjusted with a correlation function due to insufficient data or too much noise in the data, likely caused by the proximity of emission sources or surface and boundary layer effects.Consequently, a regional average of error statistics obtained for other monitoring stations was provided as a replacement for these particular stations where statistical estimation could not be obtained.Figure 1 shows an example of the application of the Hollingsworth and Lönnberg method for a typical site (here the Goddard Space Flight Center air quality monitoring station): σ 2 f is the intercept of the fitted first-order autoregressive model, and σ 2 o is the residual (or nugget) error variance at zero distance.As a result of the fitting, an estimation of the local isotropic correlation length, L ci , at site i is also obtained.Since the correlation model does not allow for nonhomogeneous background error correlations, a spatially averaged uniform correlation length is used in the optimum interpolation computer code (see Sect. 6 for a discussion of the implications of the homogeneous assumption).A sensitivity analysis of the two parameters revealed that a smaller analysis error can be obtained whenever Eq. ( 9) was satisfied (e.g., normalized chi-square χ 2 p is one).On this basis, an adaptive scheme has been developed by utilizing the chi-square diagnostic to scale the error statistics.The scheme to tune error statistics parameters (correlation length scale and background error variance σ f ) is presented below.The two tunable parameters are automatically adjusted "on-line" using a recalculation of the K matrix for each iteration.We will first describe the case when the initial value of χ 2 p is greater than one.The algorithm proceeds as follows: 1. Let n be the (first) iterate for which the Kalman gain is recalculated (at a given time step).First L c is reiterated, until there is convergence to a minimum of L c and also to a minimum Chi-square.If this minimum value is one, that is, then the convergent value of L c that is associated with the condition of Eq. ( 15) is the correlation length scale used in the optimum interpolation code.
2. If Eq. ( 15) is not fulfilled, e.g., whenever χ 2 p does not reach the value one, then an adjustment on σ 2 f is performed as well: until the chi-square condition Eq. ( 15) is reached.
Appendix A gives more mathematical details for the reason why the formulation for Eq. ( 14) works; that is, because (a) it allows the value of chi-square equals 1 as a possible fixedpoint solution and (b) Eq. ( 14) is a contracting transformation for an initial chi-square greater than one.For such transformation, Banach (1922) demonstrated that it has a unique fixed point so that the transformation converges in that case.
Finally, for initial chi-square less than one, the correlation length is inflated until Eq. ( 15) is verified.Although the latter situation occurs less frequently, it is also convergent, as shown by numerical experiment.In that case, note that the correlation length undergoes inflation (rather than contracting) according to Eq. ( 14).Note also that the rationale for Eq. ( 15) is based on the findings of Menard (2000) where it is demonstrated that this leads to optimal error statistics.
The adaptive scheme proposed in this study requires inverting Eq. ( 10) several times until convergence.This procedure avoids tedious work of constructing new error statistics set for each hour, season and year, as would be required otherwise (e.g., without the adaptive scheme).Nevertheless, a set of basic errors statistics were constructed for one month during summer for both 2004 (CHRONOS era) and 2012 (GEM-MACH era), while for other periods the adaptive scheme was relied upon to adjust for changing conditions.Typically less than 10 iterations are required to achieve the minimization procedure (satisfying the criteria of convergence within less than 1 %).This methodology is used to produce a long series of multi-year analyses.Since several million hourly analyses were required for this study, care was taken to limit the CPU time.However, the solver of the Optimal Interpolation scheme was computationally optimized so that an hourly map of objective analysis could be produced within a minute by a typical Linux station with the adaptive scheme.
In Sect. 3 it is demonstrated that the adaptive scheme shows better skill than the non-adaptive scenario (e.g., original un-tuned error statistics obtained from H-L, as in Fig. 1).Contraction of the correlation length has been done by other authors in a similar context, such as Frydendall et al. (2009) Fig. 1.Determining error statistics from the Hollingsworth and Lönnberg (H-L) method.Fitting model follows a FOAR (first order autoregressive) model for the error covariance (COV).For a particular monitoring station, the intercept of the curve represents the background error variance (indicated by an X on the vertical axis) and Lc i , the correlation length (value at 1/e).Note that averages (prior to fit) are calculated in bins of 30 km.The cutoff distance is taken to be 900 km.
(their Eq. 4), who scaled the correlation length obtained from the Hollingsworth and Lönnberg method (1986) according to the density of monitoring in a surface pollutant measurement network.The procedure used in this study, although very different from Frydendall et al. (2009), follows the same rationale that whenever the initial chi-square is greater than one, for instance, results improve significantly when the correlation length from H-L is reduced by a factor from 1 to 5 (as in Frydendall et al., 2009).More research needs to be done to obtain a general theoretical framework on how to optimally tune error statistics parameters.
The procedure presented in this paper (Eqs.14, 15 and 16) is shown to be a practical solution to get a better analysis but the authors do not claim that this method is a general solution.The adaptive procedure reduces the correlation length, thereby significantly reducing the systematic error and, to a lesser extent, the random error (as will be shown in Sect.3).The mean value of the correlation length, L c , is contracted to a range between 30-100 km, as opposed to the 75-300 km range obtained from the H-L method, or non-adaptive scheme.This reduced value for the correlation length is more in agreement with the correlation length used in the CMAQ (Community Multiscale Air Quality Modeling System) data assimilation algorithm (e.g., correlation length of 60 km; see Sandu and Chai, 2011).In the free troposphere, the correlation length is about one order of magnitude higher (500-1000 km), according to the literature (Liu et al., 2009;Tarasick et al., 2010;van der A et al., 2010 and others).These values are closer to the original H-L results found in this study when not using the adaptive scheme.To explain this difference, the authors suggest that the boundary layer and the topography or the proximity to emission sources results in a lower the correlation length required for surface assimilation.

Bias correction
Bias correction for analysis is a difficult problem and some hypotheses have to be made in order to obtain a solution (Dee and da Silva, 1998;Ménard, 2010).Two cases are discussed, both of which have a tractable solution: either the observation systematic errors are small with respect to the systematic model error, or vice versa.In the first case, a model bias correction can be developed, as seen below; in the latter case, an observation bias correction needs to be calculated.When both model and observations have significant biases, a bias correction scheme (for both model and observations) can still be developed, provided that statistics of the bias errors are known and have different characteristic length scales (see Ménard, 2010).In this study, the observation bias is assumed to be small compared to the model bias.The diagnostic model bias correction would subsequently be calculated as follows.Assuming the model bias ê is known, an unbiased analysis x is obtained by the following equation (i.e., the sequential form; see Eq. ( 41) in Ménard, 2010): The Kalman gain matrix is standard here, i.e., this is the same as the one used in Eq. ( 5).Grouping terms in Eq. ( 17) yields to: The last term on the right-hand side of Eq. ( 18) can be identified as being the analysis bias.As the observation bias is assumed to be small, the residuals O − A = y o − Hx a can be used as a source of information of the analysis bias ε a .Since < O−A > is only defined at the observation locations, we extend it to the whole model domain surface as follows: with O − A region is time and spatial average of O−A over a certain region and where 2 is defined by Equations ( 19) and ( 20) are applied to four regions having a form of ellipses (instead of a rectangular form to avoid "corner effects"), where a and b are the ellipse semi-axes, x and y the horizontal distance from the center of the ellipse.
The time average ε a is computed over a whole season and for four different regions, as depicted in Fig. 2b.The values chosen for a and b in Eq. ( 20) were obtained through trial and error so that they reduce the bias as much as possible in the cross-validation mode.Values of the seasonal average bias correction are constant through the region (within a given ellipse) but vary with time of day and changes with season.The O − A itself are evaluated a posteriori using a previous set of objective analyses.The bias correction is uniform inside a given region with its value equal to the region average.Outside of a region, the bias correction decays as a function of the square distance 2 .Outside of all the four regions, the bias correction is the sum of the individual distance-decaying region contributions (bottom part of Eq. ( 19), e.g., when > 1, the grid point is outside the elliptical region).Two criteria were used to define the ellipse position and parameters (e.g., values of aand b in Eq. 20): (1) the ellipses were centered where the density of stations are higher, and (2) the highest model biases (e.g., time-averaged analysis increments) were located inside the ellipses so that modeling of bias correction outside the ellipses were not critical and had little impact.The numerical value for a is taken as 1000-1200 km and b in the range 250-500 km, depending on the particular region.This is an empirical bias correction that was tailored for the problem under study.The authors do not claim that this is a general bias correction algorithm but a simple semi-empirical method that turns out to verify slightly better than without using bias correction (see the verification described below).Finally, it is important to point out that the adaptive scheme (Eqs.14 and 16 above) is applied before the bias correction scheme in this methodology because it was found that the adaptive scheme more efficiently corrects the bias than the bias correction scheme itself (see the verification below).Therefore, there is a need to use the bias correction scheme only for PM 2.5 , (Eqs.18 and 19) to remove the residual bias.

Models (trial fields)
In this study, we use CHRONOS (Canadian Hemispheric and Regional Ozone and NO x System), a chemical transport model (CTM) for air quality prediction of oxidants on both regional and hemispheric scales in Canada (Pudykiewicz et al., 1997).This model has been adopted as the trial field for the first period, 2002-2009.CTMs or any air quality models solve the mass balance equation for chemical species.They also have a photochemical module, which is the only tool available to simulate chemical transformations and capable to reproduce, in an approximate way, the chemistry of lower tropospheric pollution at every grid point of a given domain (Jacobson, 2002;Pudykiewicz et al., 1997;Seinfeld and Pandis, 2006;Pagowski et al., 2010).Archived operational model outputs have been used and the algorithm of OI applied in an off-line fashion for every hour during the period.For the period 2010-2012, an on-line model GEM-MACH replaced CHRONOS as the main model of the AQRDPS (Air Quality Regional Deterministic Prognos- tic System) suite at the Canadian Meteorological Centre (CMC).GEM-MACH is a limited area air quality model with the same gas-phase chemistry as CHRONOS (e.g., ADOM-II) but its meteorological driver is now accessible on-line.Its boundary conditions are driven by the operational version of the regional GEM model (Moran et al., 2012).This new operational model is a technical improvement over the CHRONOS model, primarily due to its better boundary conditions and that it also brings the possibility of full coupling between air quality and meteorology.The domain of both models covers most of North America and the model resolution is 21 km for CHRONOS and 15 km for GEM-MACH.The reader is referred to Côté et al. (1998) for further information on the meteorological driver GEM.Main characteristics and more details about CHRONOS and GEM-MACH are to be found in Appendix B.

Observation system
Figure 2a shows the location of surface observations for ozone used by the OA scheme (valid for summer 2010).The site density is high over the eastern US, the US west coast and the US Gulf States, lower elsewhere in the US and southern Canada, and almost vanishing in northern Canada and Alaska.For the PM 2.5 network (Fig. 2b), the number of sites is approximately two times less, although the geographical Ozone is usually measured by ultraviolet absorption with instrument requirements specified under the US National Ambient Air quality Standards (NAAQS) (www.epa.gov/air/ozonepollution/actions.html).Instrument noise error is assumed to be 1 ppbv (one part per billion in volume).However, the standard deviation of the observation error, which includes the representativeness error, is believed to be higher than 5 ppbv (Fleming et al., 2003).For PM 2.5 , TEOM (Tapered Element Oscillating Microbalance) and beta attenuation monitors have been accepted under NAAQS since 1990 (www.epa.gov/particles/actions.html).One of the most commonly used PM 2.5 monitors (TEOM-SES), however, largely underestimates concentrations in winter.The correction needed to account for that depends mostly on temperature, especially when the daily temperature is below 10 • C, due to the volatilization of the semi-volatile component of particulate mass (Allen and Sioutas, 1997).In this study, the analyses focused on the warm season (1 May-31 October), which experiences much less this instrument bias problem since temperature is normally above 10 • C in the US and southern Canada.Overall uncertainties due to PM 2.5 instrument noise are evaluated to 2 µg m −3 (Pagowski et al., 2010).The effectiveness of monitors to represent the pollution concentration in a given area depends largely on local sources and sinks, topography and meteorology, monitor location and the spatiotemporal variability (Brauer et al., 2011).Problems related to spatial representativeness and other specials issues will be addressed in future studies.

Validation of results
In this section, validation tests are presented to assess the performance of the basic objective analysis compared to different options: adaptive vs. non-adaptive, and with or without bias correction.First, cross-validation will be shown, which consists of reprocessing the objective analysis but with 90 % of the data to produce OA outputs, leaving out 10 % of the data to perform the verification itself.This group of 10 % of observations has never been seen by the analysis and is hence considered as a set of independent data.Three sets of additional similar verifying experiments are then performed and put together for the final verification.Second, comparisons with data from other sources will also be presented.
Objective analyses should be unbiased, have low random error and high reliability.Three metrics are proposed to evaluate the performance of the multi-year objective analyses (OA), and also to compare with the model results.The three metrics for performance evaluation used in the crossvalidation are: (1) average O−P (observation minus prediction) and O−A (observation minus analysis), (2) standard deviation of O−P and O−A, and (3) frequency of being correct within a factor two (FC2).Together, these metrics are nonredundant (Chang and Hanna, 2004) and were used as a set throughout this study.In fact, these metrics respectively measure the systematic bias, random error (accuracy) and reliability.The latter is a more robust measure of the performance, which is not sensitive to "outliers" or "compensating errors" (Chang and Hanna, 2004).

Cross-validation tests
Since performing cross-validation involves reprocessing of objective analyses several times to obtain a sufficient amount of data for verification purposes, only specific years were selected for this evaluation.The warm seasons of 2005 and 2007 during the CHRONOS era were chosen for the validation process of ozone and PM 2.5 respectively, and the year 2011 during the GEM-MACH era for both ozone and PM 2.5 .This 3 yr sample provides enough information so that a high degree of statistical confidence (e.g., p value < 0.05) exists for the results obtained.The verification examined four different regions for ozone and PM 2.5 .Figure 3 shows the results for ozone: the top left panel is for North America during the warm season (1 May-31 October), top right panel is for Canada, and bottom left and right for eastern and western US, respectively.The orange and blue navy curves are in all panels associated with the systematic and random errors, respectively, for the basic (non-adaptive) objective analysis scheme, and the green and cyan curves for the tuned (adaptive scheme) objective analysis.The red curves depict the model systematic errors and the black curves depict the model random errors.The latter curves are shown for purpose of reference comparison with OA.A very significant reduction of both errors (systematic and random) is obtained with the objective analyses at almost any time of the day as compared to the model forecast.For ozone, the adaptive scheme (iteratively adjusted error statistics) shows the smallest errors (random and systematic: i.e., green and cyan curves) in any region.Whenever a green dot appears on top (Fisher's test for the variance) and/or bottom (T test for average) for a specific hour, it means that the adaptive vs.
non-adaptive are significantly different at the level of confidence exceeding 95% (p value < 0.05).
The performance of the analyses for PM 2.5 during the warm season of 2007 in the cross-validation mode is shown in Fig. 4. In this particular verification, a new experiment is introduced that is adaptive with bias correction (Adap/BC) since it was established that an attempt to explicitly correct the bias for PM 2.5 was achievable and successful using Eq. ( 19).For the CHRONOS era, it is demonstrated that the adaptive scheme with bias correction (gray and pink curves for systematic and random error respectively) yields the best results overall, especially during the daytime period (e.g., more obvious for the 15Z-24Z period; see top and bottom left panels).Finally, Fig. 5 indicates that for 2011 independent verification is also excellent for OA during the GEM-MACH era (green and cyan curves, respectively) as compared to the model (red and black curves), the latter revealing relatively large biases during the daytime (12:00-24:00 UTC).In both eras, the verification shows nearly unbiased analyses as compared to the model as well as much lower random errors.
Values of FC2 for both Canada and the US for both ozone and PM 2.5 for different hours of the day (00:00 UTC, 06:00 UTC, 12:00 UTC and 18:00 UTC) were also computed (Tables 3 and 4).In principle, FC2 values can vary between 0 (absolutely unreliable) up to 1 (absolutely reliable).It is shown that, overall, the best FC2 scores are obtained with the OA adaptive scheme for ozone and the OA adaptive with bias correction (BC) in the case of PM 2.5 .At any time and anywhere, FC2 scores for OA are largely superior (e.g., more reliable) to that of the model, as should be expected from Eqs. ( 11) and ( 13).
Success of the cross-validation suggests not only that the methodology presented in Sect. 2 is sound but also implies that OA yields reasonably precise values in areas where there are no observations (but not too far away from other surrounding monitoring sites).This builds a case for using OA in several applications whenever observations are often missing (e.g., producing a climatology of pollutants with OA, or calculating trends from OA rather than directly using observations).

Comparison with other sources (global models, satellite and IMPROVE)
It is also interesting for validation purposes to compare results with external and totally independent information from various sources, such as model and satellite climatology.ozone concentration mostly over northern and central US compared to OA while GEM-AQ seems to be halfway between MOZART and CHRONOS, although GEM-AQ underestimates ozone over most of the US.The CHRONOS model (upper left panel), significantly underestimates ozone over many regions.The average of the three models provides a much better agreement with OA than each model taken individually (results not shown).Thus, OA could serve as a point of comparison and verification for global or regional models or with global surface climatologies such as provided, for (2) 1 ppbm = 1.66 ppbv.
example, by RETRO-4 or MACC or from other sources such as the CDC/EPA2 project in the US. Figure 7 compares a surface global climatology obtained from the satellite instrument MODIS (Moderate Resolution Imaging Spectroradiometer) for PM 2.5 for the period 2001-2006(van Donkelaar et al., 2010) ) with a climatology obtained from OA for the period 2004-2009 (near 18:00 UTC, which is roughly the time of satellite overpass).For the purpose of comparison, the same methodology described for the warm season was extended to the whole year (both warm and cold season) for the period 2004-2009.Although the years are different, the comparison is again instructive and could indicate flaws or weaknesses in both OA and satellite.The result is that although the respective climatology roughly agrees, important differences appear over some areas such as the Rocky Mountains (southwest US) and northern Mexico.These differences could be caused by satellite retrieval artefacts over higher elevation (R. Martin, personal communication, January 2012) or by imprecise Mexican emissions not represented correctly in the CHRONOS model, or due to meteorological conditions.Other likely explanations include a lack of dense monitoring in these areas, making OA less reliable over western higher elevation regions, and complex terrain meteorology not well simulated by the trial field (model), affecting the quality of OA.Other satellite climatologies exist for PM 2.5 such as produced by van Donkelaar et al. (2006) for the period January 2001 to October 2002 but the result was found similar to that of Fig. 7.
Comparisons were also made between the annual reconstructed fine mass concentrations (PM 2.5 ) obtained from the IMPROVE (Interagency Monitoring of Protected Visual Environments) and from the US EPA's Chemical Speciation Network (CSN) with selected results obtained in this study.For example, for the annual mean 2005-2008 (see IM-PROVE report V, 2011, Fig. 2.2.8b), the comparison with this study's results (obtained also from averaging years 2005-2008 for PM 2.5 ; figure not shown) indicates the same spatial patterns (maxima in southern California and eastern US) with mean maximum values that are also in agreement (in the range 15-21 µg m −3 ).

Objective analyses and analysis increments for
surface ozone and PM 2.5 over North America

A climatology of summer pollution
One of the main applications of multi-year series of analyses is to produce a climatology of surface pollutants.In this study, a climatology (ground-level ozone and PM 2.5 ) is created by averaging objective analyses produced using the methodology described in Sect. 2 during the summer months (June-July-August) for all available years in this study.Mapping a climatology with OA is more appropriate and more precise than either using models or observations alone (according to Eqs. 11 and 13) as long as the OA biases are relatively small throughout the whole study period.Moreover, in OA, missing observations are not critical for the quality of OA if there are other observations available in the neighborhood (according to independent verification, as previously discussed in Sect.3). Figure 8 shows the monitoring of OA systematic (bias) and random (standard deviation) errors for the period 2002-2012 and compares with the model in use at the time for (a) ozone and (b) PM 2.5 .Throughout the study period, the OA errors (bias and systematic) are rather stable as compared to the model.In fact, the OA systematic error (bias) is near zero for both ozone (Fig. 8a) and PM 2.5 (Fig. 8b), which is not the case with the model systematic error, which is much higher than that of OA.The latter is due to changes of model version, improvement of the parameterizations, and changes of biogenic and anthropogenic emissions through the years.As previously discussed, the very low bias of OA is due mostly to the impact of the adaptive scheme and to some extent to the bias correction in the case of PM 2.5 .The random error for OA is approximately 2 times less than that for model ozone (Fig. 8a) and approximately 1.5 times less for OA-PM 2.5 as compared to the model (Fig. 8b). Figure 8 shows a relatively high accuracy for OA: e.g., an average absolute systematic error less than 0.6 pbbv and 0.7 µg m −3 respectively for ozone and PM 2.5 and a random error generally less than 9 ppbv for ozone and under 12 ug m −3 for PM 2.5 during the warm season.The fact that OA for ozone and PM 2.5 is virtually unbiased with relatively low systematic errors provides confidence in using it for different applications.Figure 9a and b show average OA outputs for the two main eras for ozone and for PM 2.5 , respectively: prior to and including 2009 using the CHRONOS model, and after 2009 using the GEM-MACH model.The top panels of both figures are computed averages of all the objective analyses for all hours during the summer (June-July-August), for all available years.The bottom panels are also for summer but valid at midday (2 p.m. LT: local time).The left panels are for ozone and the right panels for PM 2.5 .One can observe that the highest levels of ozone and PM 2.5 , key components of smog, in both periods tend to be observed in the eastern US (south of the Great Lakes) and southern California, as expected (according to US/EPA; www.epa.gov/airquality/greenbook/; those regions also correspond to the highest frequency of NAAQS non-attainment).The diurnal variation is stronger for ozone than that for PM 2.5 since the maps at 2 p.m. LT (bottom panels) are quite different than the average computed for all hours (top panels).This particular time of day is of interest because (1) during the warm season, the planetary boundary layer is often well mixed so that the pollutant values are more representative of the whole boundary layer rather than just the surface values, and (2) it is roughly coincident with geo-synchronous satellite passage such as MODIS.The analyses of the second era (2010-2012), using the GEM-MACH model (Fig. 9b), roughly indicate the same situation as for the CHRONOS era (Fig. 9a), except for an increase of average values of background ozone, especially in northern Canada.However, the reliability of OA is uncertain in northern Canada.Together both Fig, 9a and b (multi-year averaged objective analyses) satisfactorily represent a summer climatology for the past decade for ground-level ozone and PM 2.5 over North America.To our best knowledge, this is the first peer-review manuscript of ground-level climatology (ozone and PM 2.5 ) for both Canada and the US based on a series of multi-year analyses on a high spatiotemporal resolution (15-21 km; 1 h).Similar analyses in the US covered only the eastern part of the country (Berrocal et al., 2010(Berrocal et al., , 2012; CDC/EPA project).

Long-term averages of analysis increments
In principle, a long-term average of analysis increment (correction to the model due to observations) reveals, among other things, how much the model is different from the analysis for various time of day and for various regions and chemical species.In this study the mean increment is used to analyze and monitor the change between the CHRONOS era to the GEM-MACH era. Figure 10a and b depict the aver-age analysis increment (AI) for both eras (CHRONOS and GEM-MACH, respectively).A dipolar structure in the zonal direction of the AI (indicated by + and − in the figures) is noted for ozone during the CHRONOS era (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), meaning that the model had the tendency to overestimate in the eastern part of North America and to underestimate in the western part (Fig. 10a, left panels).This tendency is also present for the GEM-MACH model (Fig. 10b, left panels) but negative increments seem to be augmented for this latter model in the east while positive increments have dimin-ished in the west, so overall the zonal gradient remains almost unchanged.PM 2.5 analysis increment climatology reveals that the CHRONOS model generally underestimated PM 2.5 (positive AI, Fig. 10a, right panels) whereas the GEM-MACH model overestimated (negative AI) in the eastern US and underestimated (positive AI) in the west, leading to the presence of a noticeable zonal dipolar pattern appearing for PM 2.5 during the GEM-MACH era (Fig. 10b, right panels).Observations of persistent AI may provide insight on the AQ forecasting suite behavior and possible model weaknesses.In principle, patterns in analysis increment or a presence of a dipole (persistent negative values in eastern US and positive in western US) may reveal structures of compensating model errors that are corrected by the OA.Therefore, a climatology of analysis increment is also necessary in the monitoring of forecasting systems because it reveals how much the correction is needed to get close to the true value.In fact, the random error augmented for PM 2.5 when the switch from CHRONOS to GEM-MACH took place in November 2009 (see Fig. 8b), which is consistent with the above findings regarding analysis increment patterns for PM 2.5 .Note that despite the models' larger systematic and random the OA adaptive scheme naturally dampens erratic model behavior, as revealed by Fig. 8.More specifically, it shows a low and steady bias near zero through time.

A study of decadal trend and mapping of warm season pollution
Another important application of a long series of analyses is to calculate temporal trends over different regions in order to evaluate the effectiveness of pollution control measures and regulation, and also to monitor any rise of the background pollution related to intercontinental transport.As previously explained, missing observations are not critical for OA.Moreover the multi-year analyses presented here offer low biases and random error (see Fig. 8 and Sect.3), they can be used to compute trends.To complement this analysis of trend, differences in the model grid space between two similar years (2005 and 2012) are also mapped.Using OA for trends or for mapping purposes is more advantageous than performing computations of only local observations at specific sites.This is because, according to Eq. ( 11), the overall objective analysis error variance (σ 2 a ) is always smaller than the observation variance error (σ 2 o ) or the model/background error variance (σ 2 b ), no matter if we have Gaussian distributions errors or not for both observations and model forecasts.Furthermore, OA has a series of quality controls such as background checks so that outliers are filtered out more efficiently.Particularly, the background check step (comparison of observations with model prediction, e.g., O−P) provides a powerful quality control test that rejects data which are approximately five (ten) times the standard deviation of O−P for ozone (PM 2.5 ).For example, in this study, the background check was found effective in eliminating the zerospan test of the ozone analyzer, which is spurious and which sometimes is not filtered out from Canadian observation raw data.The background check was also able to remove, at times, data influenced too heavily by the immediate prox-imity of local strong sources of PM 2.5 , such as acute spikes observations produced by too close proximity to local fires or fireworks, making data non-representative at the scale of the OA, i.e., 15-21 km.Moreover, if an observation was missing, there was no hole in the spatiotemporal sequence, since the analysis provided a likely value at a specific site where otherwise a break in the observation sequence would have been.
Finally, OA provides maps and permits the study of the interannual evolution across geographical regions for the whole of North America at once, not only at a single site.

Decadal trends and inter-annual fluctuations for ozone and PM 2.5
Inter-annual fluctuations over the past decade or so are driven by the five main driving mechanisms: (1) impact of better regulation and anti-pollution measures (EPA, 2010(EPA, , 2012;;Oltmans et al., 2013), (2) meteorological fluctuations (Thompson et al., 2001;Camalier et al., 2007;Zheng et al., 2007;Davis et al., 2011), (3) socio-economic changes (e.g., recession; Friedlingstein et al., 2010;Granados et al., 2012;Castellanos and Folkert Boersma, 2012), ( 4) increase of background levels due to intercontinental transport (Cooper, 2010(Cooper, , 2012)), and ( 5) inter-annual variability of stratospherictropospheric exchanges (Tarasick et al., 2010;Lin et al., 2012a, b).Following Cooper et al. (2010Cooper et al. ( , 2012) ) and Vautard et al. (2006), a general trend analysis is computed by using the percentiles to take into account the changes of the full distribution through time.High percentiles (e.g., 95th, 98th or 99th) changes are more likely to indicate local changes in emission whereas low percentiles changes rather indicate global or background changes (Cooper at al., 2010(Cooper at al., , 2012;;Lelieveld et al., 2004).However, increases of low percentiles could also be caused by lesser urban nocturnal titration of ozone by NO x (Vautard et al., 2006).According to Fig. 11a and Table 5a, ozone high percentiles are decreasing overall whereas low percentile as well as the median and the mean are all increasing.This implies that the standard deviation of the distribution is becoming smaller with time.Note that high percentiles of ozone values mostly occur in the afternoon.For PM 2.5 (Fig. 11b), high percentiles are also decreasing with time but low percentiles are neither increasing nor decreasing (see also Table 5a).Note that Table 5 is also broken down into individual regions (eastern US: Table 5b, western US: Table 5c, eastern Canada: Table 5d, and western Canada: Table 5e).Overall, similar results are found on the regional scale.As previously mentioned, the decrease of high percentile is likely associated with successful anti-pollution measures and regulations, which act to lower the peak values through time.The authors believe that the general downward decadal trend of the 95th, 98th and 99th percentile for both ozone and PM 2.5 (black dashed line in Fig. 11a that summertime extreme ozone events in many US urban areas have indeed decreased (Lefohn et al., 2008), which agrees with the results obtained here.In general, Fig. 11a, b and Table 5 also agree with results of Cooper et al. (2010Cooper et al. ( , 2012)).Background ozone is increasing, likely due to intercontinental transport from emerging countries (Stohl et al., 2002;Cooper et al., 2010Cooper et al., , 2012)).One exception is the region of western Canada where the decrease of high percentile is not statistically significant and even reverses to show an increase at the percentile 95th.Therefore, in that region, it is likely that the effect of the growing oil business overcame the decreasing trends associated with better regulations seen in other regions.Further analysis was conducted to focus on the causes of inter-annual fluctuations and not only on long-term trends.To study this, the decadal trend (dash black line in Fig. 11a, b) is removed and the inter-annual fluctuations correlated with selected predictors such as the US temperature and the mean monthly precipitation of each summer (June-July-August) for ozone (N = 11 yr, which is 2002-2012) and PM 2.5 (N = 9, which is 2004-2012).Correlations were also computed between those fluctuations with known economical indices such as the Industrial Dow Jones (http://www.djaverages.com)and various forms of the gross domestic product growth rate (http://www.tradingeconomics.com/united-states/gdp-growth) in order to analyze the influences of economic fluctuations on pollution levels.A multiple regression model using a stepwise-like procedure establishes the explained variance of the main predictors to the overall model (stepwise procedure of SAS, Statistical Analysis Software version 9.3).Stepwise regression is an automatic procedure to select the best predictors for a regression model using a sequence of F-tests in a stepwise manner (for more details, see Weisberg, 1985 or Draper andSmith, 1981).Table 6 presents the regression equation along with the percentage of the statistical model explained by each predictor.For PM 2.5 , gdpmol (previous year warm season gross domestic product growth rate) explains 37 % of the variance of the 98th percentile fluctuations (devp98-PM 25 ) whereas the mean summer precipitation (pjjaus) explains 27.5 % of the total variance of the inter-annual fluctuations.The statistical model itself explains 64.5 % (R 2 = 0.645) of the total variance.For ozone, as expected, the mean temperature of June-July-August (ttjjaus) explains most of the fluctuations for the 98th percentile (R 2 = 76 %) whereas the warm season gdpmo (gross domestic product from May through October of the current year) explains 11 % of the fluctuations.Finally, the US gross domestic product growth rate from July to December of the previous year (e.g., gdpjdl) explains the remainder (4 %) for a total of 91 % for the whole model.These links between pollution and recent economic fluctuations established here have also been observed elsewhere, notably Castellanos and Folkert Boersma (2012), who attributed the acceleration of the downward trend of tropospheric NO 2 column over Europe in 2008-2009 to distinct changes in anthropogenic activity in Europe linked to sharp downturns in gross domestic product caused by the global economic recession.Finally, downward trends shown in Fig. 11a, b compare well to national US trends for ozone and PM 2.5 (EPA, 2012, Fig. 4).

Mapping differences between 2005 and 2012 for ozone and PM 2.5 using OA
The goal of this section is to complement the trend computation performed in Sect.5.1 by mapping geographic changes.Differences are mapped in the analysis grid space for surface ozone and PM 2.5 between two similar years.The methodology to select these two similar years is given below.Beside the reduction of emissions due to regulations, the factors that most influence the formation of a photochemical pollutant such as ozone are the weather (more specifically the temperature; see for example, Camalier et al., 2007), the short-term fluctuations of economy (Friedlingstein et al., 2010;Castellanos and Folkert Boersma, 2012;Granados et al., 2012) 3 The theory about principal component analysis and loading factors can be found elsewhere in textbooks such as Morrison (1976).
(see Fig. C1 in Appendix C).Note that the weather fluctuations control only the third principal component axis (results not shown).This means that the first two components (first and second principal axes on the graph) could be interpreted in terms of the dominant variables (gross domestic product and area burned by wildfires).The proximity of the two years 2005 and 2012 in the principal component analysis plot (depicted by the arrow on Fig. C1) indicates similarities of the above-mentioned factors.These two years (2005 and 2012) were therefore selected because they both show roughly similar dominant loading factors: e.g., similar gross domestic product growth rate in both the US and Canada that are positive and above average in both countries for both years, and also similar fire regime, which is above average in the US for both years and below average in Canada for both years (see Table 7).Computing differences between these two "similar" years minimizes biases caused by the inter-annual fluctuations of the above-mentioned dominant factors during the study period, allowing the impact of anthropogenic emissions changes to be isolated.As shown in Appendix C, other years also showing similarities with 2012 (2006, 2007 and 2011) could also have been selected but 2005 was chosen because (1) 2005 OA was fully tested with an in-depth cross-validation and shown to verify very well against independent observations, and (2) the difference with 2012 gives the longest period difference (7 yr).Nevertheless, results obtained with other pairs of years (e.g., 2012 minus 2006, 2012 minus 2007) were found to show similar patterns to that of 2012-2005 (figures not shown).

Ozone
Maps depicting summer averages (June-July-August) for 2005 (Fig. 12a) and for 2012 (Fig. 12b) have been produced using the OA adaptive scheme described in Sect. 2. The difference between the two selected years (2012 minus 2005) is presented in Fig. 12c.The details for each region (in % change for different percentiles, average and standard deviation trends) are given numerically and computed separately in the observation space in Table 8a.The gray areas over oceans and over continental northern Canada are considered artefacts (unreliable zones in Fig. 12c) and therefore not included in computations in Table 8.The reasons for this are that very few observations are available in those locations, and a problem with boundary conditions was present with the CHRONOS model such that the objective analysis could not correct the model due to lack of available observations in northern Canada and over oceans.Therefore, in those regions where the analysis error is equal to that of the model, the analysis has no skill.In Sect.6 we will discuss this issue in more detail.Increases of average ground-level ozone are noticeable in several parts of the US from 2005 to 2012, but especially obvious over high plains, foothills, and also in western Canada (Alberta), where positive differences become meaningful and present over large areas (e.g., regions within the black ellipse in Fig. 12c).Part of those changes could be attributed to increasing pollution brought by the intercontinental transport originating from emerging countries (Stohl et al., 2002;Cooper et al., 2005Cooper et al., , 2010)), a change of local patterns of emissions linked with forest fire activity or a possible change of vertical transport due to deep stratospheric intrusions inter-annual variability and, finally, to growing regional socio-economic factors.For the latter, oil and gas drilling in the US and oil sands activities in western Canada have both increased during the recent decade.For example, US crude oil production increased from 5 million barrels production per day in 2008 to 6.  , 2012).However, the western parts of the US and Canada are more likely to be impacted by free tropospheric air masses containing enhanced ozone from upwind emission sources, such as Asia or from the stratosphere, than the eastern parts of North America because of higher topography (Cooper et al., 2012 and ref-erences therein;Lin et al., 2012a, b).This could partially explain the absence of a negative trend of high ozone percentiles and a strong increase of low percentiles in western Canada (according to Table 5e).Finally, differences between the two years could also be partly explained in terms of temperature and possibly other meteorological factors conducive to ozone formation, which are different for 2005 and 2012 locally and regionally (see EPA website for ozone corrected due to weather, http://www.epa.gov/airtrends/weather.html and Thompson et al., 2001 for a review of statistical correction methods for year-to-year meteorological changes impacting of ozone formation).

PM 2.5
A similar computation is presented for PM 2.5 in Fig. 13.Table 8b gives the details in percent change for all percentiles and statistical first two moments (mean and standard deviation) for each region in the observation space.As observed for ozone, PM 2.5 has significantly decreased in the eastern US (particularly in regions near and south of the Great Lakes, e.g., blue areas in Fig. 13c) where most of the industrial US activities usually takes place.Differences in Fig. 13c suggest decreases of 5 to 20 µg m −3 during the whole period (from 2005 to 2012).Since it is known that a change of 10 µg m −3 in PM 2.5 is associated with an approximate 6 % change of death rate (Pope et al., 2002), this improvement in air quality should have a significant positive impact on the health of the inhabitants of these regions.On the other hand, a significant increase of PM 2.5 up to about 10 µg m −3 in the western part of North America, as seen in Fig. 13c, is expected to significantly reduce the health of exposed people in those regions.The spatial scale of changes in eastern North America (meso-scale) contrasts with that in the west, which seems to be more localized (smaller spatial scale).In the east, the decrease suggests a generalized positive effect of anti-pollution measurements adopted through the years, which in turn has had a successful impact in reducing domestic emissions.For ozone and PM 2.5 , especially in the western part of North America, positive increase from 2005 to 2012 noted locally are more symptomatic of growing local socio-economic and industrial activities as mentioned above (e.g., increasing oil and gas drilling activities).
But it could also be linked with increase of fire occurrence at specific locations, which generates locally great amounts of PM 2.5 .In other words, even if the total area burned by wildfires in the US was roughly similar in 2005 compared to 2012, as mentioned above, at the local scale, local fires and wildfires are not likely to have occurred in the exact same locations.This could partially explain the very local scale nature of the differences in some areas.Note that in the future, the increase of PM 2.5 and ozone is likely to occur in the western US and Canada.According to the US National Council study, 200-400 % increases in burned area are expected per degree of warming in the western US (National Academy of Sciences, 2011).The impact of anti-pollution measures (negative trends) is obvious in eastern parts of the US and southern Ontario (Canada) (as indicated by the black ellipse on Fig. 13c) as well as in some local parts of southern California, while in other geographical areas other factors, as mentioned above, seem to overcome these reduction measures such that it reverts to positive differences over a broad area inside the black ellipse in the western US and Canada (e.g., Montana, Idaho, Utah, Colorado, North and South Dakota, Alberta and Saskatchewan).Finally, both Figs.12c (ozone) and 13c (PM 2.5 ) show differences of roughly the same pattern (some increase in the west and a general decrease in the east), indicating some consistency between the two pollutants.The results obtained in this section for the US are broadly consistent with those produced by other sources (EPA, 2010, 2012 for ozone and PM 2.5 ; IMPROVE, 2011 for PM 2.5 ).

Implication of homogeneous and isotropic assumptions
Homogeneous and isotropic assumptions were used in deriving the analysis scheme in Sect. 2. Whenever the urban-scale pollution gradient is high, the topography is irregular or the density of stations is weak, homogeneous and isotropic assumptions can be questionable.In this study, the fact that the density of the data over the US and southern Canada is high (Fig. 2), at least in urban centers and for some cases near coastline stations (California and US eastern seaboard), the assumption can be made that homogeneity and isotropy are not critical, at least for these locations.Moreover, the verification with independent data in different regions and with different types of sites (rural, semi-urban and urban; results not shown) does not indicate major differences or any serious problem with the above hypotheses.Simulations using nonhomogeneous background error statistics were performed in  -2005).Note that the areas indicated unreliable are due to model artefacts and unavailable observations so that differences are meaningless and therefore not shown.These areas are also where the analysis error is too high (see Fig. 14a, b).Within the ellipse is a region of significant changes over a widespread region (see text).Note: units are in ppbv.
the context of this study but were found to be inconclusive.More studies are needed to examine the impact of using the homogeneous and isotropic assumptions, which is beyond the scope of the present study.Reference is made to the work of Frydendall et al. (2009) for the treatment of anisotropy and to Blond et al. (2004) for non-homogeneous treatment of error statistics for air quality assimilation.

Biases, accuracy and limitations of multi-year analyses
Multi-year analyses of high spatiotemporal resolution (15-21 km; every hour) were built in this study instead of reanalyses4 .Standard chemical reanalyses, which use the same model setup to generate all the historical analyses, tend to require an enormous amount of human resources and very tedious work (e.g., RETRO-40 or MACC projects in Europe or CDC project in the US).Conversely, multi-year objective analyses for surface, as produced in this study, are simply off-line objective analyses reprocessed from the archived operational model outputs and are less demanding on resources.For example, the computing cost of integration is significantly lower and could be done on a Linux machine since model outputs are precalculated.However, the main advantages of reanalysis are retained only if care is taken to eliminate the systematic biases at any time in the analyses, as was done in this study.Otherwise, incorrect trends could be produced due to various changes of model versions, out-of-date same resolution throughout the whole study period, whereas multiyear analyses can allow for some change during the period.
emissions and/or increases in resolution if the biases are not eliminated in the objective analysis scheme.
The long-term analyses presented here are unbiased, or have very small biases (see Figs. 3,4,5 and 8), and are available on a 21 km grid prior to 2009 and 15 km after 2009.They are also available in terms of hourly, daily, monthly, seasonally, yearly and multi-year averages.The model domain covers most of North America but since surface observations are only dense over the continental US and southern Canada, it is only in these regions that observations can constrain the model and that the confidence in the results is high.A typical map of analysis error using Eq. ( 13) is presented in Fig. 14a for ozone and Fig. 14b for PM 2.5 .In areas where the density of stations is high (see Fig. 2), the analysis error could be 2-4 times lower than in those locations where the density is low.Note that values above a certain threshold are not plotted on the maps of Fig. 14 because there are insufficient observational data in these regions.Therefore, there is no added value of analyses compared to the model in these regions and the analysis is thus considered unreliable.A useful simple application of the analysis error map could be in the assessment of an optimal network density.In regions where the analysis error is high (low), it is necessary (unnecessary) to increase the network density.In locations where there are no monitoring stations but other available observations, the OA still yields good results.This is considered an advantage for some applications (climatology, trend computation) using OA instead of observations or modeling alone.For example, missing observations or model biases are not an issue for OA.

Comparison of trend analyses with other studies
The summertime trends in ozone found in this study are consistent with observations made in other studies.A review of reported ozone trends by Chan and Vet (2010) found mean positive trends ranging between 0.3-1.0ppbv yr −1 for Canada and the US.A study by Vautard et al. (2006) for Europe also established a positive trend in ozone levels of 0.65 ppbv yr −1 , which is within the same range of values.The above results are also consistent with Table 5a of this study, that is, an increase of 0.47 ppbv yr −1 for the median (percentile 50th) and a mean upward trend of 0.3 ppbv yr −1 for southern Canada and the US.As global temperatures continue to rise, more favorable conditions for ozone formation are likely to occur due to factors such as increases in biogenic and soil emissions (IPCC, 2007;Zheng et al., 2008) and increases in wildfires (IPCC, 2007;Jaffe et al., 2003, Houghton, 2009;National Academy of Science, 2011).These factors, which in the recent past exhibited a downward trend during the warm season, may result in higher future percentiles of ozone and PM 2.5 (see US EPA, 2010, Fig. 35 for an outlook for the eastern US).Moreover, expanding gas and oil drilling could also contribute more to ozone formation in the short-term future.

Spatial scale representativeness
It is also important to discuss the spatial scale representativeness of OA.Spatial variability is important and pollutant concentration gradients could be high at times, raising the issue of the local representativeness in the objective analysis, especially in urban environments.For ozone, the urban-scale gradient could be strongly dependent on the distribution of sinks (mostly NO x titration associated with urban traffic) and certainly not picked up correctly by the OA.In urban areas, the spatial variation of NO x is important in determining the surface ozone urban-scale gradient.It is possible for future work to utilize land use regression models, which are highly correlated with NO x , to capture the urban-scale gradient as a means of producing better data fusion of information at higher resolution.The spatial scale represented by the OA is the scale of the analysis grid (15-21 km).More research is needed to capture adequately local scale features.
For PM 2.5 , Brauer et al. (2011) indicate that in some Canadian cities significant spatial variation exists, while in others PM 2.5 mass is spatially homogeneous.Therefore, within-city spatial variation of ozone and fine particulate matter (PM 2.5 ) is case specific.NO x and ultrafine particles (diameter less than 1 micron) are known to be highly correlated but this association does not necessarily hold true for larger particles.Difficulties also arise because primary and secondary PM 2.5 observations are considered together and treated as one single family.However, information on the chemical composition of PM is required in order to elucidate the processes governing production, transport and deposition at different scales, as well as chemical composition, which differs between and within primary and secondary PM (Hobbs, 1993;Jacobson, 2002;Seinfeld and Pandis, 2006).Furthermore, primary PM 2.5 exhibits spatial variability over small scales while secondary particles tend to be more uniformly distributed (Blanchard et al., 1999;Pinto et al., 2004).These issues will be addressed by future work within the context of OA further development.

Summary and conclusions
The purpose of this study is twofold: (1) to present multi-year analyses of ground-level ozone and PM 2.5 over North America during the warm season (1 May-31 October), from which a climatology could be built; and (2) to apply the results for computation of a summer (JJA) decadal trend and map differences of summer average of two similar years in order to isolate the impact of reduction of anthropogenic emissions over the past decade or so.The multi-year analyses themselves form a coherent and continuous data set for the period 2002-2012 for ozone and 2004-2012 for PM 2.5 .The analyses are freely available upon request.To the best of the authors' knowledge, no such multi-year analyses have ever been presented for both ground-level ozone and PM 2.5 for both the US and Canada covering a large portion of North America over a decadal period and at such a high spatiotemporal resolution (15-21 km; every hour) over a decade or so.The analyses are based on a methodology that utilizes a modified optimal interpolation scheme adapted for air quality.Analyses have been obtained through an optimal combination of the CMC's Air Quality Regional Deterministic Prediction System (ARDQPS: composed of CHRONOS 2002-2009and GEM-MACH 2010-2012) for which model outputs are available for every hour and AIRNow surface observations database (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), supplemented with extra Canadian stations (not part of the AIRNow program added during the GEM-MACH era).A semi-empirical tuning procedure (referred to here as the adaptive scheme) using the Chi-square statistic was developed and applied on-line to adjust some sensitive parameters of the error statistics set.The methodology was tested successfully and verification with independent data has shown excellent results.The impact of the adaptive scheme was shown to significantly reduce both the bias and random error.An explicit bias correction scheme was also used for PM 2.5 to further reduce the residual biases.By running this optimal interpolation over a decade or so, a high-quality integrated estimate of two of the main components of human health-threatening smog has been produced.The multi-year analyses presented here are at high spatiotemporal resolution and show a relatively high accuracy with an average absolute systematic error less than 0.6 pbbv and 0.7 µg m −3 respectively for ozone and PM 2.5 and a random error generally less than 9 ppbv for ozone and under 12 µg m −3 for PM 2.5 during the warm season.
Long-term averages are presented as summer climatology maps (June-July-August) for ground-level ozone and PM 2.5 .The objective analyses obtained are also used to compute trends for ozone and PM 2.5 .Low percentiles of ozone exhibit an upward trend (southern Canada and US together) while high percentiles of ozone show a downward trend overall for North America during the warm season.Some local exceptions to the overall downward trend in high percentile ozone are found in the northwestern part of North America (northwest US and Alberta).The results presented in this study are comparable with other studies that have examined long-term trends in ozone (Vautard et al., 2006;Chan and Vet, 2010;Cooper et al., 2010Cooper et al., , 2012;;EPA, 2010EPA, , 2012)).The decreasing trends in the high percentiles of ozone and PM 2.5 strongly suggest that domestic emissions reduction has been effective; this is especially obvious for the eastern parts of North America.The reduction in high percentile concentrations of these pollutants implies that human and environmental health risks associated with air pollutant exposure have decreased over the last decade, at least in the eastern US.However, global (background) transport of ozone is increasing and, combined with climate warming, could produce a further increase in ozone (high and low percentiles) in the future.Moreover, oil and gas industries are still developing in North America at a high rate, which could also contribute to increases in ozone and PM 2.5 in the future.
Finally, a multiple linear regression (stepwise-like procedure) confirms that a significant part of the variance of interannual fluctuations of high percentiles for both ozone and PM 2.5 is linked with economical fluctuations.For PM 2.5 , the gross domestic product growth rate of the preceding warm season (gdpmol) explains 37 % of the total variance of the inter-annual variation of the 98th percentile (after removing the linear trend tendency), e.g., devp98.For ozone, although the summer temperature explains most of the variance (76 %), various forms of the gross domestic product growth rate explain up to 15 % of the variance.Economic recession can trigger noticeable short-term changes in anthropogenic emissions that can reduce pollution.Sharp downturns in the economy are linked to decreases in industrial, construction, transportation and in other human activities in North America during the recession of 2008-2010.Presumably, this has lead to an analogous decrease in the high percentiles for ozone and PM 2.5 during that period.
The multi-year analyses presented here were intended mainly for model evaluation, computation of regional pollution trends and for epidemiological studies.Unresolved issues include the treatment of random high pollution such as forest fires.Monitoring stations in the vicinity or downwind generally record high levels of PM 2.5 , since forest fire emissions are not captured by the operational ARDQPS suite.Consequently, the OA quality control could reject part of the data associated with such an extreme event.Another unresolved issue is the inability of the long-term average or climatology to correctly capture fine-scale pollution gradients.These topics will be addressed in future work.

Fig. 2 .
Fig. 2. Map of available observation sites (circa 2010) used for multi-year analyses for (A) ozone and (B) PM 2.5 .The four ellipses superimposed on the bottom part of panel (B) are indicated for the four regions where a bias correction is performed for PM 2.5 .Equations for the bias correction (inside and outside the ellipses) are given in the text (see Eqs. 19 and 20).

Fig. 3 .
Fig. 3. Warm season cross-validation for the year 2005 (CHRONOS era) for ozone vs. time of day (UTC).The metrics used are mean (systematic error) and standard deviation (random error) of OmP (observation minus precipitation) and OmA (observation minus analysis).The diurnal variation of systematic and random errors, respectively, for model (red and black curves), non-adaptive objective analysis (orange and blue navy) and adaptive OA (green and cyan) are presented for four different regions.Top left panel: North America, top right: Canada, bottom left: eastern US, and bottom right: western US. Green dots at the top (bottom) of each panel indicate a successful F test variance (T test bias) for statistical significance of the difference between two selected experiments (i.e., adaptive vs. non-adaptive scheme).Units are in ppbv.
Figure 6 presents objective analyses averaged during the whole year of 2005 (bottom left panel), CHRONOS model output (upper left panel) with compatible yearly outputs from the MOZART model (Horowitz et al., 2003, Model for OZone And Related Tracer -version 2: horizontal resolution of 2.8 • : upper right panel) and GEM-AQ (Global Environmental Multiscale coupled with Air Quality) model (Kaminski et al., 2008, resolution of 1.5 • : bottom right panel).Although the global MOZART (Model for Ozone and Related Tracers)and GEM-AQ models have lower horizontal resolution (few hundreds kilometers) than that of the Canadian AQ model suite (resolution of 21 km for OA-CHRONOS), the comparison is nonetheless instructive.In fact, the general pattern of the two global models is roughly in agreement with the objective analyses (OA) especially near coastlines and the Gulf of Mexico where the uncertainty of OA is higher due to less observations available there.MOZART overestimates

Fig. 4 .Fig. 5 .
Fig. 4. As Fig. 3 but with cross-validation for the year 2007 for PM 2.5 .The diurnal variation of systematic and random errors, respectively, for model (red and black), non-adaptive OA (orange and blue navy) and adaptive OA (green and cyan) and adaptive OA with an explicit bias correction (gray and pink) are presented for four different regions.Top left panel: North America, top right: Canada, bottom left: eastern US, and bottom right: western US.Significance tests for difference are as in Fig. 3.The differences tested are between the adaptive with bias correction vs. the adaptive scheme with no bias correction.NA: adaptive scheme OFF, Adap: Adaptive scheme ON, NBC: no bias correction, BC: with bias correction.Units are in µg m −3 .

Fig. 6 .
Fig. 6.Comparison of surface ozone OA annual average for 2005 (all hours, all seasons combined) with external sources.Top left panel: CHRONOS model 2005 (no observation ingested).Bottom left: OA -average for 2005.Top right: MOZART model annual average (version 2).Bottom right: GEM-AQ annual average for 2005.High ozone values are in red, low values are in blue.Note that (1) for the sake of the comparison only for this figure, units are given in mass mixing ratio (ppbm) instead of volume mixing ratio (ppbv);(2) 1 ppbm = 1.66 ppbv.

Fig. 7 .
Fig. 7. Comparison of surface PM 2.5 climatology obtained from (A) satellite surface-derived PM 2.5 2001-2006 (MODIS), van Donkelaar et al. (2010); (B) OA PM 2.5 average for all hours, all seasons 2004-2009 near 18:00 UTC (at approximately the time of satellite overpass).Note that both figures (A) and (B) have the same color bar.High values are in red, low values in blue (units are in µg m −3 ).

Fig. 8 .
Fig. 8. Evolution of the systematic and random error for model and OA suite over a decade for (A) ozone (2002-2012) in ppbv, (B) PM 2.5 (2004-2012) in µg m −3 .The model in use is indicated at the bottom of the figure (e.g., either CHRONOS or GEM-MACH).

Fig. 9 .
Fig. 9. (A) Long-term average OA (CHRONOS era; i.e., before 2010) for summer months (JJA) for surface ozone and PM 2.5 .Top left panel: all hours ozone analysis.Bottom left: ozone analysis at 2 p.m. LT.Top right: all hours PM 2.5 analysis.Bottom right: PM 2.5 analysis at 2 p.m. LT.High ozone values are in red and low values are in blue.(B) as (A) but for the GEM-MACH era (2010-2012).Units are ppbv for ozone and µg m −3 for PM 2.5 .

Fig. 10 .
Fig. 10.(A)Long-term average analysis increment (CHRONOS era: i.e, before 2010) for summer months June-July-August (JJA) for surface ozone (in ppbv) and PM 2.5 (in µg m −3 ).Top left panel: all hours ozone analysis increment.Bottom left: ozone analysis increment at 18:00 UTC.Top right: all hours PM 2.5 analysis increment.Bottom right: PM 2.5 analysis increment at 18:00 UTC.Red values are strong positive corrections to the model, blue values are strong negative corrections.(B) as (A) but for theGEM-MACH era (i.e., 2010-2012).

Fig. 11 .
Fig. 11.Long-term trend of percentiles (5th, 25th, 50th, 75th, 95th, 98th and 99th) and standard deviation (STD DEV) for (A) ozone distribution (units are in ppbv) and (B) PM 2.5 distribution (units are in µg m −3 ).Computations were done in the observation space for summer months using all hours available.Procedure REG of SAS(SAS, 1988  for more details) was used to draw the dash line trend (decade tendency).

Fig. 12 .
Fig. 12.Comparison of surface average ozone summer months (JJA) in 2012 vs. 2005: (A) OA ozone 2005 (CHRONOS era), (B) OA ozone 2012 (GEM-MACH era), (C) difference (OA-2012 minus OA-2005).Note that the areas indicated unreliable are due to model artefacts and unavailable observations so that differences are meaningless and therefore not shown.These areas are also where the analysis error is too high (see Fig.14a, b).Within the ellipse is a region of significant changes over a widespread region (see text).Note: units are in ppbv.

Fig. 14 .
Fig. 14.Analysis errors based on Eq. (13) for (A) ozone and (B) PM 2.5 .Deep blue corresponds to very small analysis errors whereas red is associated to higher errors.The locations where there are no values plotted are where the analysis has no skill (i.e., OA is unreliable due to model erratic behavior and/or no observations available to correct model values).

Table 1 .
Summary of environmental and health impacts of ozone and PM 2.5 .
1 Ozone in presence of sunlight decomposes into O 2 and O.The oxygen molecule combines with water vapor to give OH

Table 2 .
Number of stations available from US/EPA AIRNow database and Canadian stations for 2005 (CHRONOS era) and 2012 (GEM-MACH era).
trieval Now) program.Raw data are provided by numerous US local air quality agencies (between 150 and 200 agencies in US) as well as Canadian agencies 1 .The AIRNow US EPA ozone and PM 2.5 real-time data base has been available to us for surface ozone observations since 2002 for a large part of North America and since 2003 for PM 2.5 .However, the multi-year analysis retrospective examines data since 2004 for PM 2.5, at which time the observation network became more stable.By 2012, approximately 1200 AIRNow sites provided data on an hourly basis, plus an extra 100 stations originating from Canadian provinces and territories that are not part of the AIRNow program.

Table 3 .
Performance of model and OA for the warm (May-October) season 2005 (cross-validation mode) evaluated using FC2 (frequency of correct value within a factor two when compared to observations) in cross-validation mode for ozone (a) for Canada (N ∼ 1440 observations) and (b) for the US (N ∼ 13200 observations).Note that Z stands for UTC (universal coordinated time); BC for bias correction.

Table 4 .
As Table3but for PM 2.5 and for the warm season 2007 (N ∼ 1200 for Canada and N ∼ 8000 for US).

Table 5 .
Trends for selected percentile (PCT), mean and standard deviation for OA-ozone (ppbv yr −1 ) and OA-PM 2.5 (µg m −3 yr −1 ) for different regions.(a) Southern Canada and US, (b) Eastern US.(c) Western US, (d) Eastern Canada, (e) Western Canada.The p value is given for statistical significance.NS indicates no statistical significance (p value > 0.15).Positive (negative) values indicate an increase (decrease) from 2005 to 2012.Computations are for all hours during summer (JJA).

Table 6 .
Multiple regression models to explain high percentile fluctuations for ozone and PM 2.5 .The % of the variance explained by each predictor is indicated below each term of the equation.tjjaus: mean US temperature (in • C) for June July and August of the current year, pjjaus: mean US precipitation (in mm) for June July and August of the current year, gdpmo: gross domestic product growth rate (in %) from May-October of the current year, gdpjdl:

Table 7 .
Numerical data for factors potentially influencing inter-annual variability of ozone and PM 2.5 for the period 2002-2012 (tacan and taus: summer temperature anomaly ( • C) for Canada and US, respectively; gdprc and gdprus: current summer gross domestic growth rate (%) for Canada and US, respectively; wfabc and wfabus: wildfire area burned in Canada and US, respectively, in millions of hectares; noajja and ensojja: summer North Atlantic oscillation index and El Niño-Southern Oscillation MEI index, respectively).

Table 8 .
Percentage changes of OA (2012 minus 2005)for ozone (a) and (b) PM 2.5 in North America.Positive (negative) values indicate an increase (decrease) from 2005 to 2012.Computations are for all hours during summer (JJA).Trends were computed using the procedure REG of SAS (SAS, SAS/STAT, 1988).