Improved satellite retrievals of NO2 and SO2 over the Canadian oil sands and comparisons with surface measurements

Satellite remote sensing is increasingly being used to monitor air quality over localized sources such as the Canadian oil sands. Following an initial study, significantly low biases have been identified in current NO2 and SO2 retrieval products from the Ozone Monitoring Instrument (OMI) satellite sensor over this location resulting from a combination of its rapid development and small spatial scale. Air mass factors (AMFs) used to convert line-of-sight “slant” columns to vertical columns were re-calculated for this region based on updated and higher resolution input information including absorber profiles from a regional-scale (15 km × 15 km resolution) air quality model, higher spatial and temporal resolution surface reflectivity, and an improved treatment of snow. The overall impact of these new Environment Canada (EC) AMFs led to substantial increases in the peak NO2 and SO2 average vertical column density (VCD), occurring over an area of intensive surface mining, by factors of 2 and 1.4, respectively, relative to estimates made with previous AMFs. Comparisons are made with long-term averages of NO2 and SO2 (2005–2011) from in situ surface monitors by using the air quality model to map the OMI VCDs to surface concentrations. This new OMI-EC product is able to capture the spatial distribution of the in situ instruments (slopes of 0.65 to 1.0, correlation coefficients of > 0.9). The concentration absolute values from surface network observations were in reasonable agreement, with OMI-EC NO2 and SO2 biased low by roughly 30 %. Several complications were addressed including correction for the interference effect in the surface NO2 instruments and smoothing and clear-sky biases in the OMI measurements. Overall these results highlight the importance of using input information that accounts for the spatial and temporal variability of the location of interest when performing retrievals.


Introduction
Space-based measurements of the near-surface atmospheric composition, or air quality, from near-UV-to-near-IR spectra have blossomed over the past two decades from relatively crude, research-grade products to refined, operational products suitable for monitoring and assimilation (e.g. Miyazaki et al., 2012). These products include tropospheric vertical column densities of NO 2 , SO 2 , CO, HCHO, and aerosol optical depth from nadir-viewing instruments that measure backscattered sunlight (e.g. Martin, 2008). Furthermore, through the fusion of satellite data and models, fundamental quantities crucial to air quality such as surface concentration (Lamsal et al., 2008) and emission rates (Streets This category of sensor began with the GOME (Global Ozone Monitoring Experiment) instrument  and continued with SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CartograpHY, 2002CartograpHY, -2012 (Bovensmann, 1999), OMI (Ozone Monitoring Instrument, 2004-present) (Levelt et al., 2006), and the operational GOME-2 (2006GOME-2 ( -present, 2012 instruments. Collectively these instruments have provided invaluable information on distributions and trends in NO 2 and SO 2 despite the complexities associated with the inversion of these spectra. Applications of these data to air quality issues are numerous and span wide spatial and temporal scales. More recently these data have been applied to the analysis of localized sources whose spatial extent is comparable to that of a individual resolution element, or pixel (Beirle et al., 2011;Fioletov et al., 2011). One high-profile example is the Canadian oil sands (McLinden et al., 2012). This area in the northeast corner of the province of Alberta contains a vast deposit of hydrocarbons, including an equivalent of 170 billion barrels (roughly 2.7 × 10 7 m 3 ) of oil in the form of bitumen, a viscous form of petroleum. Production has increased rapidly from about 0.6 million barrels per day (mBPD) in 1998 to 2 mBPD in 2012 and with a further doubling expected over the next decade (ERCB, 2012). Another measure of the size of the industrial complex is capital expenditures, which have fluctuated between CAD 10 and 20 billion since the mid-2000s, and are expected to remain in this range for at least the next decade. An initial study using these satellite instruments found distinct enhancements of NO 2 and SO 2 over the area of intensive surface mining (an area of about 50 km × 30 km) and with the enhancement in NO 2 increasing at a rate of about 10 % year −1 (McLinden et al., 2012). A map of the area showing the oil sands boundary and the location of the surface mining is given in Fig. 1. The rapid expansion is evident from the Landsat images. Overlaid on the Landsat images are approximate boundaries of the active mining regions (for 2011) and are used for reference in later figures. A reference location (57.1 • N, 111.6 • W) is also identified.
Even among near point sources the oil sands represent a particular challenge to satellite remote sensing due to a combination of the rapidly changing landscape and emissions, its higher latitude (57 • N), frequent snow cover, and higher boundary layer winds that disperse the emitted compounds more rapidly. These factors conspire to make signals here more difficult to observe and, perhaps more importantly, the rapid industrial evolution of the region makes some of the assumptions in the current generation of retrieval algorithms suspect.

OMI data products
Beginning with calibrated and geo-located nadir spectra, and with the exception of so-called direct inversions (e.g. Nowlan et al., 2011), UV-visible tropospheric NO 2 and SO 2 retrieval algorithms are composed of three distinct steps. The first step is a determination of the total absorption by that species (where the target species is hereafter referred to simply as the absorber) and this is quantified in terms of a slant column density (SCD), S. The SCD represents the absorber number density integrated along the path of the sunlight through the atmosphere. Since OMI measures backscattered light the path is not a direct one and includes one or more scattering events and/or surface reflections. SCDs are derived through an analysis of the measured spectra by exploiting the difference in absorption at nearby wavelengths. The second step is the removal of the stratospheric component of the total SCD, relevant only for NO 2 . This is generally done by either a filtering or statistical approach or by explicitly modelling or assimilating the stratospheric SCD and subtracting it. The final step is a conversion of the tropospheric SCD, S t , into a tropospheric vertical column density (VCD), V t , via an air mass factor (AMF), M, that accounts for the sensitivity of the instrument to the absorber. Owing to the complex path of the scattered and reflected sunlight, AMFs must be computed using a radiative transfer model. These tropospheric VCDs are the primary data product and generally serve as the starting point for scientific and monitoring applications.
Multiple OMI NO 2 VCD data products currently exist. The two global, primary products are the NASA standard product (SP) (Wenig et al., 2008;Bucsela et al., 2013; http://disc.sci.gsfc.nasa.gov/Aura/ data-holdings/OMI/omno2_v003.shtml) and the Dutch OMI NO 2 (DOMINO) (Boersma et al., 2007 http://www.temis.nl/airpollution/no2.html) processed in near real time. While there were some significant differences between their respective version 1 data products, their version 2 (v2) products are in much better agreement (Bucsela et al., 2013). The DOMINO v2 tropospheric VCD uncertainty is 25 % + 10 15 cm −2 . DOMINO and the SP use common SCDs derived using the technique of differential optical absorption spectroscopy (Platt, 1994) over the 405-465 nm wavelength range, but differ in how they remove the stratospheric SCD. The SP employs a complex high-pass-filtering approach, whereas DOMINO simulates the stratospheric NO 2 by assimilating OMI SCDs in a chemical data assimilation system . In addition there are other, offline products available such as the Empa OMI product for Europe (Zhou et al., 2009;http://temis.empa.ch/index.php) and the regional Berkeley High-Resolution (BEHR) product (Russell et al., 2011), which makes use of higher resolution input data but over a domain limited to the continental US. All of these products calculate an AMF specific to a given measurement based on a variety of factors.
In contrast, there is only one OMI SO 2 VCD planetary boundary layer (PBL) data product available. It is based on the NASA band residual method (BRM, http://disc.sci.gsfc. nasa.gov/Aura/data-holdings/OMI/omso2_v003.shtml), which uses four wavelengths between 310 and 315 nm to quantify SO 2 absorption (Krotkov et al., 2006). Unlike the case of NO 2 , the NASA SO 2 product uses a spatial and temporally invariant AMF, calculated for summertime conditions in the eastern US. Monthly average local AMFs based on the GEOS-CHEM model (Lee et al., 2009) along with other corrections (e.g. subtracting Pacific background) are only applied in the operational daily gridded (level 3) OMI SO2 product (http://disc.sci.gsfc.nasa.gov/Aura/ data-holdings/OMI/omso2e_v003.shtml).
In this work the DOMINO and SP NO 2 and the NASA SO 2 products are examined. They are collectively referred to as the "current" products.

Atmospheric chemistry models
This study also makes extensive use of two atmospheric chemistry models. The first is the Global Environmental Multi-scale -Modelling Air quality and CHemistry (GEM-MACH). GEM-MACH is the Canadian regional air quality forecast model used operationally to predict the concentrations of O 3 , NO 2 , and PM 2.5 over North America (Moran et al., 2009, Anselmo et al., 2010. The model makes use of detailed tropospheric processes for gas and particle chemistry and microphysics originating in the offline AURAMS model (A Unified Regional Air-quality Modelling System; Gong et al., 2006), but incorporated on-line into the Canadian weather forecast model (Global Environmental Multiscale model, Côté et al., 1998). A detailed description of the chemical processes found in AURAMS and GEM-MACH may also be found in Kelly et al. (2012). Both AURAMS and GEM-MACH share a sectional, speciated particle distribution -for the operational GEM-MACH forecasts used here, two bins are used: one for representing particle fine mode and the other for coarse mode. The results used here are from archived forecasts from 2010 to 2011 for a domain covering North America at 15 km × 15 km resolution. The emissions inventories for the model are from US EPA and Environment Canada data for the year 2006. Note that GEM-MACH at present does not include NO x sources for biomass burning and lightning. The cloud package used in GEM-MACH is described in Bélair et al. (2005).
The second is the GEOS-Chem chemical transport model  version v8-03-01 (http://www.geos-chem. org) driven by assimilated meteorological observations from the Goddard Earth Observing System (GEOS-5). These simulations were run at 1/2 • (latitude) by 2/3 • (longitude) resolution in a nested mode over a North American grid van Donkelaar et al., 2012). The model includes a detailed simulation of tropospheric ozone-NO xhydrocarbon chemistry as well as of aerosols and their precursors (Park et al., 2004). Canadian anthropogenic emissions are from the CAC inventory (http://www.ec.gc.ca/ inrp-npri) for 2005. Emissions from open fires for individual years are from the GFED2 inventory with monthly resolution (van der Werf et al., 2009). Lightning NO x emissions are computed as a function of cloud-top height, and are scaled globally as described by Sauvage et al. (2007) to match Optical Transient Detector/Lightning Imaging Sensor (OTD/LIS) climatological observations of lightning flashes. The global source is imposed to be 6 Tg(N) yr −1 . Higher NO x yields per flashes are used at mid-latitudes than in the tropics (Hudman et al., 2007). Different versions of the GEOS-CHEM model wave been used extensively in the retrieval of tropospheric VCDs, and have been shown capable of simulating the vertical distributions of NO 2 (e.g. Martin et al., 2002;Lamsal et al., 2008) and SO 2 (e.g. Lee et al., 2009).
Daily output from both models are sampled at the local time of the OMI overpass, roughly 13:30 LST, and from this monthly means were calculated.

OMI measurements over the oil sands
The rapid development of the oil sands raises concern over the validity of some assumptions underlying current NO 2 and SO 2 data products. Space-based nadir UV sensors are generally less sensitive to an absorber located near the surface as opposed to one aloft, and thus a key parameter in the retrieval algorithm is the assumed vertical profile of the absorber. Both NO 2 algorithms, the SP and DOMINO, rely on model-calculated NO 2 profiles: the SP uses monthly mean profiles from the Global Modeling Initiative (GMI) (Buc-sela et al., 2013) and DOMINO uses daily output from the TM4 (Tracer Model 4) chemical transport model (Boersma et al., 2007). These models, however, make use of emission inventories appropriate for the late 1990s, when oil sands' NO x emissions were significantly lower than current values (Boersma et al., 2007;B. Swartz, personal communication, 2012). The impact of this can be seen in the NO 2 profiles from these different sources of data in the vicinity of the oil sands. Figure 2a shows the annual-mean NO 2 profiles over the oil sands, sampled at the local time of the OMI overpass, from the GMI (horizontal resolution of 2 • lat × 2.5 • long) and TM4 (2 • × 3 • ) models. Also shown for comparison are profiles from the higher resolution models discussed in Sect. 2.3, GEOS-CHEM and GEM-MACH, smoothed to match the resolution of GMI. The profiles are presented as shape factors, which are volume mixing ratio (vmr) profiles normalized by their column-averaged vmr. The GMI and TM4 profiles show values that are significantly smaller in the PBL, by a factor of 2-3, than GEOS-CHEM and GEM-MACH and generally resemble shape factors for near-background profiles. In fact, their absolute value in the boundary layer is only about 0.1-0.3 ppb, and so is representative of background levels. As is discussed below, the end result of this underestimate of NO 2 in the PBL is an underestimate of tropospheric VCDs.
Beyond the outdated emissions is the related issue of model resolution. With grid boxes in excess of 200 km (latitude) × 150 km (longitude), the GMI and TM4 models cannot resolve the surface mining region whose entirety spans 50 km × 30 km. This means the same model NO 2 profile may be used at both the centre of emissions and 100 km away, where the NO 2 would be at or near background levels. This issue of profile "representativeness" has been identified as a potential error source for near point sources (Heckel et al., 2011). Figure 2b illustrates the steep horizontal gradients in the annual-mean NO 2 profiles from the GEM-MACH model by comparing vmr profiles at the maximum emissions with those at 10 and 50 km away.
Related concerns exist for SO 2 VCDs: the use of an invariant AMF, calculated for summertime conditions in the eastern US, will likely also lead to spatial and seasonal biases, although in this case the sign of the bias is less apparent. On this basis it was concluded that there is a potential for significant systematic errors in these current (DOMINO and SP NO 2 and NASA SO 2 ) products over the oil sands region. Thus, any quantitative assessment of these gases requires AMFs to be re-evaluated in this region, accounting for both the rapidly changing emissions and their small spatial scale. Furthermore, calculating NO 2 and SO 2 AMFs in a consistent manner should increase their compatibility, possibly allowing for their ratio or differences in their spatial distributions to be used to infer additional information about the nature of sources.

General concept
The general AMF framework used herein follows that of Palmer et al. (2001), Martin et al. (2002), and others. The AMF (M), defined as the ratio of the SCD (S) to the VCD (V ), or M = S/V , describes the enhancement in absorption when light traverses a slant path through a layer and represents a cornerstone of trace gas retrievals in the UV-visible portion of the spectrum. From this definition, the steps in the retrieval algorithm described in Sect. 2.2 can be expressed as where S = S t + S s , and subscripts "s" and "t" refer to stratosphere and troposphere, respectively, and M t is the tropospheric AMF. Subscripts are used to specifically reference the troposphere or stratosphere, and not used for general concepts or total (troposphere + stratosphere). In Eq. (1), S, determined from the spectral analysis, is the only quantity measured directly by OMI; S s and M t must be determined from a combination of statistical considerations, modelling, or assimilation. For SO 2 it is assumed there is no appreciable stratospheric component (S s = 0). When direct sunlight is measured, as may be the case for a ground-based spectrometer, it is straightforward to determine the path of the sunlight through the atmosphere, and hence the AMF may be deduced from purely geometrical considerations. However, when scattered sunlight is the source of information, the path is much more complex, involving, in general, multiple-scattering and surface reflection events. Moreover, the measured radiance will be a combination of many different paths. Resolving this requires radiative transfer models capable of simulating the multiple-scattering and absorption processes.
The distribution of scattered light, including nadir radiances (and thus AMFs), is controlled by many factors, including solar and instrument viewing geometry, surface reflectivity, aerosols and clouds, and the vertical distribution of the absorber (e.g. Martin et al., 2003). It is most convenient to account for the dependence on vertical distribution by using a vertically resolved AMF, m(z), such that M represents the absorber number density weighted average of m(z), where n(z) is the absorber number density vertical profile and α[T (z)) is a correction factor that accounts for the change in the absorption cross section with temperature. For both species a correction of α =1-0.003 [T (z) − T 0 ] was used, where T 0 = 220 K for NO 2 (Boersma et al., 2004) and T 0 = 273 K for SO 2 (Bucsela et al., 2013) correspond to the temperatures of the cross sections used in the spectral analysis. Values of m(z) are also referred to as box AMFs and are analogous to the scattering weights used in other studies (e.g. Martin et al., 2002) but without a normalization by a geometric AMF. Note that Eq.
(2) is quite general, but for the purpose of this work, the integration is from the surface to the tropopause. Another of the key factors in determining AMFs are clouds. Clouds behave as bright surfaces reflecting more light than that of the underlying surface. They also act to shield the sensor from any absorbers located below the cloud. The treatment of clouds used here follows the general approach used in other studies (Boersma et al., 2004;Martin et al., 2006) in which they are modelled as Lambertian reflectors located at the cloud altitude (corresponding to the optical centroid), and with an albedo of 0.8. In the case of a partially cloudy pixel, the AMF is taken as a linear combination of the cloudy and clear-sky AMF as follows (Ahmad et al., 2004): where w is the cloud radiance fraction and the subscripts "c" and "a" represent cloudy-and clear sky, respectively. The cloud radiance fraction represents the fraction of the nadir radiance that is due to the cloudy portion of the pixel. It differs from the effective cloud fraction, f , as clouds tend to be more reflective and thus contribute disproportionately to the overall radiance. They are related by where I is the nadir radiance.

Air mass factor input information
AMFs are calculated using radiative transfer models that simulate nadir radiances while accounting for all relevant physical processes, including multiple scattering by molecules and aerosols, absorption by trace gases, and reflection from the surface. Their accuracy is governed to a large extent by the accuracy of the input information with the most important being the shape of the absorber vertical profile, the surface reflectivity, cloud and aerosol information, surface pressure, and, in the case of SO 2 , the stratospheric ozone column. These are discussed below and summarized in Table 1.

Absorber profiles
Knowledge of the NO 2 and SO 2 vertical distribution is required to correctly determine the instrument sensitivity. Absorption by these species is sufficiently small that they do not impact the bulk properties of the radiation field, and thus it is only the shape of the vertical profile that is relevant. In order to better capture the large spatial gradients that these species display in the immediate vicinity of the surface mines, the GEM-MACH regional, air quality forecast model is used as a source of profile information. The 2010-2011 monthly means are used for all years as GEM-MACH output is not available prior to mid-2009. The 15 km grid size is well suited for OMI retrievals as it is comparable in size to that of the median OMI pixel size (consider only the smaller pixels) and so no additional smoothing was applied. The lack of biomass-burning and lightning NO x sources in GEM-MACH can be seen in Fig. 2a, as NO 2 profiles in the free troposphere are smaller than those of the other three models. To remedy this, the GEM-MACH profiles above the boundary layer are transitioned to monthly mean profiles from the GEOS-CHEM model. When GEM-MACH profiles are smoothed to the GEOS-CHEM resolution, there is good general agreement between the two in the boundary layer, and thus this transitioning does not generally introduce significant discontinuities. A linear weighting is used between −0.5 and +1. km, relative to the GEM-MACH boundary layer height, to ensure a smooth hybrid profile. For PBL heights below 0.5 km, the transition begins at the surface. The coarser resolution of GEOS-CHEM (by roughly a factor of 10 compared with GEM-MACH) is not an issue in the free troposphere, where horizontal gradients tend to be smaller. For consistency, hybrid SO 2 profiles are constructed in a similar manner. The use of pure GEM-MACH and GEOS-CHEM profiles is discussed as part of the sensitivity study, below. A comparison of the annual mean GEOS-CHEM and the hybrid vmr profiles (smoothed to GEOS-CHEM resolution) over the surface mining region is shown in Fig. 2c.

Surface reflectivity
For simplicity, most AMF algorithms incorporate surface reflection using a Lambertian equivalent reflectivity (LER). The often used monthly mean LER climatology derived from GOME (Koelemeijer et al., 2003) or the minimum OMI reflectivity, minLER (Kleipool et al., 2008) are, at best, available on a 0.5 • × 0.5 • grid. While in some respects the use of an OMI-derived minLER database is advantageous for OMI trace gas retrievals (same viewing conditions and wavelengths), compared with the spatial scale of interest in this work and the resolution of other input information, 0.5 • remains somewhat coarse. Furthermore, its lack of interannual variability is problematic in that the rapidly changing land cover of the oil sands may introduce a bias in any trend from an unaccounted for change in reflectivity.
An alternative source of albedo is available from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite instrument , used previously for OMI NO 2 retrievals (Zhou et al., 2009(Zhou et al., , 2010Russel et al., 2011). MODIS provides white-sky albedo (WSA) and blacksky albedo (BSA), based on 16-day averages available every 8 days, at a resolution of 0.05 • × 0.05 • (MOD43C3, collection 5). The improved spatial resolution of this product is a significant advantage when targeting near point sources, and the interannual variability accounts for the evolution of land cover changes. A limitation of this data product, however, is the spectral range -the shortest wavelength is 477 nm, similar to that required for NO 2 (∼ 440 nm) but long for the SO 2 spectral region (∼ 315 nm). The errors introduced in OMI NO 2 AMFs using this MODIS product over the more rigorous bidirectional reflection distribution function (BRDF) have been examined by Zhou et al. (2010). In the summer errors are < 5 %, while for lower sun angles errors of 10-20 % are possible depending on the shape of the absorbing profile.
The spatial patterns in the OMI albedos (Kleipool et al., 2008) at 342 nm and 477 nm were found to be generally consistent, thereby suggesting it is reasonable to simply scale the MODIS albedo by an OMI-derived ratio to arrive at a better representation of reflectivity at near-SO 2 wavelengths. As MODIS provides both BSA and WSA, where BSA is more appropriate for direct sunlight incident on the surface and WSA for diffuse light, a weighted average is used with the weighting determined by the model-calculated fraction of downwelling irradiance that is diffuse. In this work, albedo was calculated using where α bs and α ws are the MODIS monthly mean BSA and WSA at a wavelength of λ M = 477 nm and f dif is the model-calculated (see Sect. 3.3) fraction of the total downwelling irradiance at the surface that is diffuse. The averages are calculated using 100 % snow-free scenes and were smoothed to approximate the resolution of a representative OMI pixel (15 km × 30 km). The final factor in Eq. (4) adjusts the MODIS albedo to a wavelength more representative of the spectral region used in the analyses, λ = 440 nm for NO 2 or λ = 342 nm for SO 2 , using the updated Kleipool et al. (2008) OMI monthly climatology. Examples of the MODIS-derived surface albedo (via Eq. 4) for NO 2 and SO 2 are shown in Fig. 3, along with that from the OMI climatology, averaged over the summer (JJA). The MODIS maps show the for NO 2 and SO 2 albedo for 2005 and 2011. The OMI and MODIS albedos smoothed to OMI resolution ( Fig. 3a and b) appear very similar overall, but the higher resolution of the MODIS clearly shows the footprint of the surface mines ( Fig. 3c-f). They also reveal that there is an evolution with time. Over the mining area the albedo is seen to increase with time, which is consistent with a conversion of the darker forest into the brighter industrial land use. Also of note are the slightly larger values of albedo for SO 2 relative to NO 2 resulting from the wavelength adjustment.

Identification and treatment of snow
Owing to the large fraction of OMI observations that are made over snow in the oil sands region, roughly 40 %, it is important to address this aspect as well. Problems with UV-visible measurements over snow have been discussed in previous work (O'Byrne et al., 2010). In principle, snow represents an ideal surface in that its higher reflectivity means the instrument is more sensitive to the near surface and thus less sensitive to the shape of the absorber profile. In practice, however, complications arise due to issues with correctly identifying the presence of snow (either false positives or false negatives) and the subsequent choice of an appropriate surface albedo (snow or snow-free) for the determination of the cloud fraction and AMFs. The end result is often large errors in the cloud fraction and the use of an inappropriate surface reflectivity. The intent here is not to remedy the entirety of the problem as this requires an improved cloud fraction data product, which is beyond the scope of this study. Nonetheless, two elements can be readily addressed: identification of snow and the choice of snow albedo.
The presence of snow in an OMI scene is currently determined using the Near-real-time Ice and Snow Extent (NISE; http://nsidc.org/data/nise1.html), an operational, daily, global product derived from the Special Sensor Microwave Imager/Sounder (SSMIS) passive microwave sensor (Nolin et al., 1998) on a 24 km × 24 km grid. When using an independent, ground-based determination of snow cover, O'Byrne et al. (2010) found this product often missed thin snow cover, a common feature among snow products derived from passive microwave instruments. A similar evaluation was conducted here. Figure 4 shows comparisons (2000-2011) between daily snow cover from NISE compared with measurements made at two meteorological stations: Fort McMurray (30 km south of the surface mining) and Mildred Lake (in the middle of surface mining). When taking the surface-observed conditions as the truth, NISE misses snow cover roughly 30 % of the time in an annual average. Breaking this down by month (Fig. 4b) reveals that the largest discrepancies generally appear in the autumn months, when NISE missed snow more than half the time. Other snow products were examined to determine their potential suitability including the Canadian Meteorological Centre (CMC) Daily Snow Depth Analysis Data (Brown et al., 2003) and the Interactive Multisensor Snow and Ice Mapping System (IMS) (Helfrich et al., 2007; http://www. natice.noaa.gov/ims/). Both the CMC and the IMS are northern hemispheric products on a 24 km × 24 km grid and have been evaluated for their accuracy in northern Canada (Frei and Lee, 2010). IMS is an operational daily, northern hemispheric product also on a 4 km grid. Comparisons of these products with the same surface data are also shown in Fig. 4. Annually, CMC and IMS misidentified the conditions at the meteorological station 6 % of the time, compared with NISE, which was incorrect 17 % of the time. As seen in Fig. 4b CMC tended to over-predict snow in the spring, whereas IMS under-predicted snow in the fall. Delving deeper into these instances reveals that when IMS does miss snow, the snow depth is usually very thin, 1-2 cm.
When present, snow tends to increase the reflectivity of the scene. The DOMINO retrieval assumed a constant snow albedo of 0.6, and while this value may be reasonable for a global average, snow albedo depends significantly on the underying surface type with values varying between 0.2 for forests and 0.8 for open plains (Moody et al., 2007). Here, MODIS reflectivity was again used to determine snow albedo. Average BSA and WSA maps were derived from MODIS (2005MODIS ( -2012 for December, January, and February, limited to 100 % snow-covered scenes, and smoothed to 15 km × 30 km. The resultant albedo map is shown in Fig. 5 along with a snow-albedo product from OMI on a 0.5 • × 0.5 • grid (O'Byrne et al., 2010). Both are generally consistent with snow albedo values in the 0.2 to 0.4 range. However, due to its coarser resolution the OMI product misses the oil sands entirely, which, according to MODIS, is enhanced by as much as 0.2 relative to its surroundings.
The approach adopted in this work was to use the IMS product to identify snow-covered scenes. The IMS product was selected over the CMC as it is operational, whereas there is a delay before the CMC reanalysis product is available. As discussed below, snow-covered scenes are not considered in the OMI data analysis and thus the mention of snow albedo, above, may seem somewhat superfluous. However, future studies will address the issue of cloud fraction determination over snow, where an improved estimate of snow albedo becomes important.

Topography, surface pressure, and temperature profiles
Elevation, surface pressure, and temperature profiles are all taken from the GEM-MACH model with surface pressure and temperature profiles sampled at the local time of the OMI overpass and averaged monthly.

Ozone
At SO 2 wavelengths, overhead ozone impacts AMFs as it alters the depth to which the UV light can penetrate. Each OMI observation was assigned an OMI-derived total column ozone value based on the Total Ozone Mapping Spectrometer, TOMS, algorithm version 8 (Bhartia and Wellemeyer, 2002) appropriate for that day. If a value was not available for that day, the mean for that month was used instead.

Aerosols
A robust implementation of aerosol in UV-visible retrievals remains a challenge (Leitão et al., 2010). Some algorithms do not explicitly include aerosols and argue that, as the OMI cloud fraction algorithm is somewhat sensitive to them, there is some cancellation of errors . The degree to which, or even if, this applies will depend on where in the vertical the aerosols are located relative to the absorber. Other algorithms use an aerosol correction factor that can be applied to an aerosol-free AMF (Lee et al., 2009). MODIS observations of average aerosol optical depth over the oil sands suggests values are typically small, 0.1 to 0.2, and are consistent with sun photometer measurements from Fort Mc-Murray (Holben et al., 1998; http://aeronet.gsfc.nasa.gov/), roughly 30 km south of surface mining area. Previous studies, Martin et al. (2003) for NO 2 and Lee et al. (2009) for SO 2 , show that aerosol correction factors suggest aerosols would impact SO 2 AMFs in this region by a few percent at most, and thus, on this basis, aerosols are neglected. An assessment of the uncertainty introduced from this assumption is provided in Sect. 4.1, below.

Calculation of AMFs
AMFs were calculated using the VECTOR (VECTor Ordersof-scattering Radiative transfer model) (McLinden et al., 2002(McLinden et al., , 2006Wagner et al., 2007). As an initial step to verify that VECTOR can calculate AMFs that are consistent with those in the DOMINO v2.0 product, they were recalculated by VECTOR and a comparison was made. AMFs were recalculated using, where available, the same input information as used by DOMINO. This included the NO 2 profile, albedo, surface pressure, and geometry. Input information that was not available included the temperature and air density profiles, and thus a climatology was employed. For cloud-free pixels within 100 km of the reference location VECTOR AMFs were found to differ from the DOMINO AMFs on average by less than 3 %. Considering some small differences between models is expected (see, for example, Wagner et al., 2007) this agreement is acceptable (especially given that not all input parameters were identical). Using the input information described in Sect. 3.2, vertically resolved AMFs, m(z), were computed on a 0.5 km altitude grid between 0 and 16 km by successively perturbing the absorber over the 0.5 km layers and computing where I is the change in radiance for a change in the optical depth of τ due to the perturbation. It is impractical to calculate m(z) for each OMI observation. Instead, a pair of multidimensional look-up tables for m(z), I , and f dir , were generated, one for cloud-free conditions and one for cloudy conditions. Recall that f dir is used in the calculation of albedo in Eq. (4) and so is not relevant for cloud AMFs since the albedo is set to 0.8. The dependencies of the cloud-free table are surface albedo, solar zenith angle, viewing zenith angle, change in azimuthal angle, surface pressure, and column ozone. The dependencies of the cloud table are cloud-top pressure, solar zenith angle, viewing zenith angle, change in azimuthal angle, and column ozone. These are summarized in Table 2.

Environment Canada AMFs over the oil sands
AMFs were computed as described above for each OMI NO 2 and SO 2 measurement recorded within 500 km of the reference location (57.1 • N, 111.6 • W -see Fig. 1). These AMFs are referred to hereafter as Environment Canada (EC) AMFs.
A scatterplot comparing these with the DOMINO v2 AMFs is shown in Fig. 6a. The highest density of points occurs in two locations on or just below the 1 : 1 line: at an AMF of 1.5 for snow-free conditions, and around 2.5-3 for snowcovered surfaces. The points corresponding to the surface mining area, which represents a small fraction overall, lie well below the 1 : 1 as indicated. This is a result primarily of a larger fraction of the NO 2 being present in the PBL. The cluster of points representing snow-covered surfaces is generally below the 1 : 1 line as a result of the smaller values snow albedo used for the EC-AMFs, typically 0.2-0.35 as compared to 0.6 in the DOMINO product. The purpose of this scatterplot is to explore how the input information outlined in Sect. 3.2 directly impacts AMFs. As such, both TEMIS and EC AMFs use the same cloud fraction. Assessing how a cloud fraction, based on improved snow identification and snow albedo, might also impact AMFs is beyond the scope of this study. Figure 6b shows a histogram of SO 2 AMFs. Recall that the NASA algorithm uses a constant AMF of 0.36. The EC snow-free average SO 2 AMF is 0.41 with considerable variability, including minimum values of 0.25 that occur over the surface mining area. For snow scenes the average EC SO 2 AMF is 0.79. The spatial distribution of these AMFs is shown in Fig. 7, where average summer (June-July-August) NO 2 and SO 2 AMFs are shown. These maps are calculated using the pixelaveraging method of Fioletov et al. (2011) on a 1 × 1 km 2 grid and with an averaging radius of 8 km. The DOMINO v2 NO 2 AMFs (panel a) display no difference between the surface mining area and the surroundings, consistent with the background-like TM4 profiles from Fig. 2. In contrast the EC NO 2 AMFs show a distinct minimum over the surface mines but otherwise appear similar to the DOMINO AMFs if not slightly smaller. Note that all AMFs display a banding due to sampling geometry. The DOMINO : EC AMF ratio, a quantity that reflects the expected impact of VCDs, is also shown. It has a maximum of about 1.9 over the mines, but decreases rapidly with distance to about 1. The SO 2 EC-AMFs are also shown. The minimum over the mines is about 0.2-0.25 and they increase to about 0.4 in the surrounding. The minimum SO 2 AMF covers a larger area than seen for NO 2 as a result of a larger area of enhancement predicted by GEM-MACH due to its longer lifetime and/or more rapid dispersion (see Sect. 4.2, below).
To gain a better quantitative understanding of what input information is driving these changes, NO 2 AMFs were recalculated using selected TEMIS inputs. For this, snowfree OMI observations for the year 2005 and within 50 km of the reference location were considered, where the largest differences were observed (from Fig. 4c). Results from this analysis indicate that the primary driver was in fact the profile shape and that it accounted for over 90 % of the average decrease in AMF. All other changes, on average, led to small systematic changes, although increased variability was Atmos. Chem. Phys., 14, 3637-3656 observed due to their higher spatial resolution. The change in surface reflectivity, using MODIS instead of OMI, led to a small systematic decrease in AMF due to the slightly smaller average values.

Vertical column densities over the oil sands
The OMI SO 2 SCDs are known to suffer from a variable offset or bias (Lee et al., 2009;Fioletov et al., 2011Fioletov et al., , 2013 resulting in part from an imperfect removal of the stratospheric ozone absorption signal. This bias is treated as a correction that must be made to the SCD in much the same way as the stratospheric SCD must be removed. Thus, the general expression for the tropospheric NO 2 or SO 2 VCD, Eq. (1), becomes where S b is the local bias (where S b = 0 for NO 2 ). It is calculated by averaging S over all OMI pixels located between 100 and 150 km from the reference location where SO 2 is assumed to be zero (Fioletov et al., 2011(Fioletov et al., , 2013. There is a certain amount of time dependence to this correction and so it is calculated monthly. While Eq. (6) represents the methodology conceptually, it is most convenient to simply apply the calculated quantities directly to the current VCDs using where M t is the AMF for the current product (e.g. M = 0.36 for SO 2 ) and V t,EC and M t,EC are the EC tropospheric VCDs and AMFs, respectively.

Error budget and sensitivity study
Beginning with Eq. (6), the total random uncertainty in VCD, ε, can be expressed as where ε i are uncertainties in the individual terms: ε S is the uncertainty in the SCD, ε SS in the stratospheric SCD, ε S b the SCD bias, and ε M t the AMF. Values for ε S and ε SS for NO 2 are taken directly from Boersma et al. (2007) and ε S for SO 2 from Krokov et al. (2008). The uncertainty in the SO 2 bias correction is based on the standard error of the mean of the bias, averaged over all months. These values are given in Table 3. The uncertainty in AMF arises from several different sources, including uncertainties in cloud fraction, cloud pressure, surface albedo, profile shape, terrain height/surface pressure, aerosol, and stratospheric ozone. The impact of these on the AMF is assessed by varying each by an amount indicative of its own uncertainty and recalculating the AMF. For example, the surface albedo is varied by ±0.02 and the extent to which this changes the AMF is taken as its uncertainty to this parameter. This is assessed considering all OMI observations within 200 km of the reference location for the year 2005. Perturbing an input parameter in this way yields a distribution in the relative change in AMF. The standard deviation of the distribution is taken as the uncertainty to a given parameter and it is computed separately for polluted (within 50 km of the reference location) and background (between 50 and 200 km) areas. Only observations with cloud radiance fractions of 0.2 or smaller were considered. The results of this are also given in Table 3. The uncertainty due to profile shape was evaluated by recalculating AMFs but using profiles from GEOS-CHEM. Uncertainties in these parameters were assumed to be independent and so the total uncertainty in AMF was calculated by adding the individual terms in quadrature.
From Table 3 it can be seen that uncertainties in AMFs are about 20-25 % for NO 2 and 25-35 % for SO 2 , with the larger values for the polluted areas. In the case of NO 2 , the contributions from cloud fraction and profile shape are about 10 %. For SO 2 the profile shape is the largest contributor. Considering cloud fractions larger than 0.2 leads to much larger uncertainties in AMF, consistent with Lee et al. (2009). The total NO 2 VCD random error is about 0.5 × 10 15 cm −2 for background and 0.8 × 10 15 cm −2 for polluted, with the largest contributor to this being the uncertainty in SCD. The total SO 2 VCD random error is about 0.6 DU (1 DU = 2.69 × 10 16 molecules cm −2 ) for background and 0.8 DU for polluted areas, again with the largest contributor being the SCD. For locations with larger VCDs the uncertainty in AMF can become the largest contributor.
Note that these uncertainties represent an estimate of the precision of an individual measurement and do not account for systematic errors. It was the primary goal of this work to address systematic errors in the current data products, and it is believed that the approach outlined above goes a long way towards achieving this. However, it is likely some systematic errors remain, and three potential AMF effects are discussed below. One is the use of a LER as opposed to a BRDF in the treatment of surface reflection, as discussed above. Based on Zhou et al. (2010) this is estimated to be on average less than +10 %, and +5 % or smaller in the summer. It is appropriate here to discuss the impact of using either the GEOS-CHEM or the GEM-MACH profile shapes, as opposed to the hybrid. GEOS-CHEM NO 2 profiles led to change in AMFs, −7 ± 14 % (mean ± 1-σ ). While the average change is modest, the coarser resolution of GEOS-CHEM nevertheless means changes are more significant directly over the mines. GEM-MACH profiles, which are weighted more towards the PBL, had a largest impact, −31 ± 25 %.
EC-AMFs take into account the evolution of the surface albedo but have thus far neglected the change in absorber profiles as would result from increasing emissions.  Further away where there are no sources, there should be little to no change. The relatively small impact suggests that a type of saturation effect is occurring: AMFs were already sufficiently weighted towards the PBL that adding additional NO 2 there had only a modest impact. In the limiting case that all of the NO 2 resides in the PBL, further increases have no impact at all since only the profile shape is relevant. In the evaluation of the trend, an unaccounted 6 % decrease in AMF over the OMI mission, for example 2005 to 2011 or 6 years, would be equivalent to about a 1 % year −1 decrease in AMF. This would amount to an underestimate of the trend in VCD by about 1 % year −1 . As discussed in Sect. 3.2.6, aerosols were not explicitly included in the AMF calculations. While difficult to uncouple from the aerosol-biased cloud fraction effect, this is a potential source of systematic error that needs to be explored. Similar to clouds, aerosols can either enhance the AMF due to increased scattering or decrease it by shielding an absorbing layer below. To assess this, AMFs were recalculated using a single aerosol profile shape that decreases with altitude between 0 and 3 km, purely scattering, and scaled to give an optical depth of 0.1. The inclusion of aerosols acted to either increase or decrease AMF and is linked to the relative profile shapes and height of the PBL. On average, aerosols decreased NO 2 AMF within the polluted area by 6 %, with almost no impact over the background locations. Its average effect on SO 2 was about zero. A more comprehensive treatment would consider profile shapes that vary in space and time, and would also take into account a correction for aerosols in the cloud retrieval.

VCD climatologies
Long-term (2005-2011) annual average NO 2 and SO 2 VCDs were calculated to examine their spatial distribution. Only small pixels (track positions 11-50) unaffected by the row anomaly, snow-free, and observations with a cloud radiance fraction of 0.2 or less were considered here. In addition, solar zenith angles (SZAs) were limited to a maximum of 75 • for NO 2 and a maximum of 60 • for SO 2 . Finally, SO 2 VCDs were restricted to values between −5 and +15 DU with the upper limit imposed to avoid spikes from transient volcanic SO2 clouds. Averages were calculated using the oversampling, pixel-averaging method of Fioletov et al. (2011) on a 1 × 1 km 2 grid and with an averaging radius of 8 km for NO 2 and 24 km for SO 2 . The larger averaging radius for SO 2 is required due to its higher noise. Averages are shown in Fig. 8 for the DOMINO and SP NO 2 and NASA SO 2 products. In each case, the Atmos. Chem. Phys., 14, 3637-3656 corresponding maps based on EC-VCDs are also shown. The differences between the NO 2 EC-VCDs derived from the DOMINO and SP products lie mainly in how the stratosphere is removed. All maps show a clear enhancement over the region of surface mining consistent with the location of the largest emissions and in agreement with results from McLinden et al. (2012). The enhancement in NO 2 is largest over the southern mines but displays a secondary maximum over the northern cluster of mines, whereas SO 2 is enhanced primarily only over the southern mines. This is also consistent with current information on source locations: SO 2 is emitted principally from upgrading facilities (which convert bitumen to synthetic crude) in the south, whereas NO 2 also has area sources from transportation (e.g. the heavy-hauler trucks) in both the north and south. Despite the more diffuse source of NO 2 , the SO 2 enhancement is seen to cover a larger area. This is suggestive of a longer lifetime. Also worth noting is the difference in the height of the release: much of the NO x is emitted at the surface (transportation), whereas SO 2 is emitted primarily from stacks some 200-300 m above the surface, where winds are generally faster.
The DOMINO (Fig. 8a) and NASA SP (Fig. 8c) NO 2 are generally consistent with SP being 30 % larger over the mines. The differences between DOMINO and EC (Fig. 8b) NO 2 are consistent with what might have been expected from the AMFs: comparable spatial distributions with similar background values, but over the surface mines the EC VCDs are larger by a factor of up to 1.9. The same is true when comparing the NASA (Fig. 8e) and EC (Fig. 8f) SO 2 , with the EC larger by up to a factor of 1.3. Note that the EC-VCDs derived from the DOMINO (Fig. 8b) and SP (Fig. 8d) products are in very good agreement: differences are less than 10 % over the mines and about 2 × 10 14 cm −2 in the background.
The mass of these enhancements, where the enhancement is defined as the NO 2 or SO 2 above the background (presumably from oil sands operations), was determined by performing a non-linear fit of a two-dimensional Gaussian function with a constant offset to the VCDs (Fioletov et al., 2011;McLinden et al., 2012). The mass of the enhancement is relevant as it is proportional to the local emissions. In the case of SO 2 , the background is zero (as a result of the bias correction) and no offset term is used. The mass was calculated by integrating the Gaussian over all space. Note that the fit is performed using the individual VCDs and not the maps from Fig. 8, although the two approaches yield nearly identical results. The mass of the NO 2 enhancements are 4.5 (DOMINO) and 9.8 t[NO 2 ] (EC) and the mass of the SO 2 enhancements are 60 (NASA) and 80 t[SO 2 ] (EC). Thus, the EC products gives masses that are 120 % (NO 2 ) and 34 % (SO 2 ) larger than the current products.

Comparison with surface concentrations
It is difficult to validate these results directly. Ideally, coincident measurements of VCD from ground-based spectrometers or vertically integrated aircraft profiles would be used, but over the oil sands neither of these is available. The general consistency of GEM-MACH with GEOS-CHEM, a well-established model often used for AMF calculations (e.g. Lee et al., 2009) provides some confidence in the ability of GEM-MACH to simulate NO 2 and SO 2 in the PBL. Likewise, the ability of the VECTOR RT model to generally reproduce the DOMINO AMFs is also important, but these do not directly address the issue of VCD validation.
The only NO 2 and SO 2 measurements currently available in the oil sands region are from in situ ground-based (GB) instruments as part of observing network operated by the Wood Buffalo Environment Association (WBEA; http: //www.wbea.org/) (Percy et al., 2012). Their locations are shown in Fig. 9, and additional information on the sites and their instrumentation are provided in the Supplement. Note that not all of the stations measure NO 2 . A single value of NO 2 and SO 2 for each day was determined by computing their average values between 12:00 and 14:00 local time.
The correlation between these GB measurements of concentration and OMI VCDs can be assessed as has been done by several others studies (e.g. Geddes et al., 2012), yet this approach only provides information on the ability of the satellite to capture variability. A more direct method of comparison is to transform the surface vmr into a VCD, or vice versa, with either approach requiring knowledge of the vertical profile. Here the method described in Lamsal et al. (2008) is used in which tropospheric VCD is mapped to vmr using a modelled profile, where χ is the surface vmr. This assumes the model can adequately capture the spatial and temporal behaviour of this ratio. Equation (9) was applied to the EC-VCDs using the same monthly mean profiles used in the calculation of their AMFs and the resultant vmr's are referred to as EC-vmr for simplicity. The 2005-2011 average surface NO 2 and SO 2 EC-vmr maps are shown in Fig. 10, calculated using the same pixel-averaging parameters as the VCDs from Fig. 8.
In this example the EC-vmr NO 2 is based on the DOMINO v2 product. The spatial distributions generally mimic that of their VCDs, with a maximum NO 2 vmr of 2.3 ppbv and a maximum SO 2 vmr of 4.0 ppbv. Also shown here are the corresponding surface vmr maps obtained using AMFs and VCD-to-vmr mapping based on GEOS-CHEM model profiles. The GEOS-CHEM maps show similar spatial patterns to the EC maps, but GEOS-CHEM NO 2 values are typically 30 % smaller and the SO 2 are about 10 % smaller through the enhancement, with contributions to this difference from both the AMF and the column-to-surface ratio. The locations of the WBEA stations are also shown in Fig. 10, and it is clear in the case of NO 2 that the existing stations largely miss the maximum in the OMI distribution. Comparisons between average GB and EC-vmr's are shown in Fig. 11 for NO 2 and Fig. 12 for SO 2 as a function of the latitude of the WBEA station. Inset in each is the corresponding scatter plot. The EC-vmr values are calculated by averaging over measurements within 6 km for NO 2 and 12 km for SO 2 of the WBEA station. Again, the larger SO 2 radius is necessary due to its higher noise level. The GB average for a given station was computed by first determining its average as a function of wind direction (using 20 • wide bins), and then taking this average, thereby giving an equal weighting to all wind directions. This approach is favoured as satellites measure over a large area (∼ 500 km 2 for OMI) and weight the upwind and downwind directions equally. Averaging in this way, as opposed to a "simple" average that preferentially weights the prevailing wind direction, impacts the mean NO 2 values by between −15 and +15 %. For SO 2 the effect is larger, −30 % to +30 %, which reflects its more localized sources.
One additional effect was considered in advance of the GB-satellite comparisons. The GB NO 2 measurements are made by commercial, chemiluminescence analysers that rely on molybdenum converters. These instruments alternate between measuring NO x and NO, with NO 2 inferred as the difference. However, it is well known that these instruments are also sensitive to nitric acid, Fig. 11. Average (2005Average ( -2011 ground-based in situ and EC-OMIderived NO 2 surface mixing ratio measurements as a function of WBEA station latitude. Inset is the scatter plot comparing groundbased measurements with the EC-vmr's. The red represents original OMI-vmr's, while the green is after accounting for smoothing and clear-sky bias. Errors bars denote twice the standard error of the mean. The ground-based measurements have been corrected for the interference effect and include this uncertainty in their error bars (see Table S1). The lines in the inset represent the linear fits.
peroxyacetyl nitrate (PAN), and other oxidized nitrogencontaining species (Winer et al., 1974;Lamsal et al., 2008) which are mistakenly interpreted as NO 2 . A simple factor, CF, for this "interference" effect was used herein, following Lamsal et al. (2008Lamsal et al. ( , 2010, and based on model-calculated concentrations ( AN is the sum of all alkyl nitrates and PAN is peroxyacetyl nitrate). The multipliers preceding PAN and HNO 3 account for the reduced conversion efficiency of these species by the instrument. Equation (10) was evaluated separately for monthly mean GEM-MACH and GEOS-CHEM concentrations and their average CF was used to correct the GB measurements. At the relatively low NO 2 levels over the oil sands (average values of < 5 ppb), this correction factor can be considerably smaller than one due to the larger contribution of non-NO x oxygenated nitrogen species (Lamsal et al., 2008). The station-specific correction factors are given in Table S1. Values range from 0.35 to 0.7, with the smaller factors corresponding to smaller mean NO 2 levels. It is difficult to determine the accuracy of these correction factors and errors due to both the modelled species' concentration ratios and the conversion efficiencies (which can vary from instrument to instrument) may be appreciable. For simplicity, in this work the difference between the two model evaluations of CF was used as a measure of its uncertainty. From Figs. 11 and 12, the EC-vmr's are able to capture the spatial variation and gradients of the NO 2 and SO 2 displayed by the GB stations through the mining region. For NO 2 (Fig. 11), the EC-vmr's are typically smaller than the GB values (with the CF applied) by a factor of 2 (roughly 1-2 ppb). The slope of the scatterplot is 0.47 and the correlation coefficient is 0.88. The SO 2 comparison (Fig. 12) shows even better agreement with almost no bias (slope of 0.88, correlation coefficient of 0.91). A summary of the linear correlation coefficients and slopes are given in Table 4.
Given its spatial resolution, OMI is only able to provide a smoothed version of the true surface vmr distribution. Indeed, this may be the origin of some of the GB-satellite differences seen in Figs. 11 and 12. To better understand this, idealized estimates of the true NO 2 and SO 2 surface vmr distributions were constructed assuming the distribution resulting from a source "region" is reasonably represented by a 2-D Gaussian. Parameters for the Gaussians were selected so that (i) their vmr's were comparable to the average measured values at the GB stations, and (ii) after smoothing the distributions generally resembled those from Fig. 8 (although not necessarily the absolute values). For NO 2 , the sum of three Gaussians was used: one each for the north and south grouping of mines, and a smaller one for the Fort McMurray area. For SO 2 only one Gaussian was used, reflecting the lack of a significant source of SO 2 in the north or Fort McMurray. GB measurements from Fort Chipewyan (station 12) were used to define background values. These idealized distributions were then smoothed using a 15 km × 30 km 2-D boxcar to simulate the OMI pixel size, and both original and smoothed distributions (shown in the Supplement) were sampled at the location of the WBEA stations using values with 6 km for NO 2 or 12 km for SO 2 . As expected, stations near the peaks became smaller and those in the wings became larger as a result of the smoothing: NO 2 was impacted by −10 to +60 %, SO 2 by −30 to +20 %. Of particular interest is NO 2 near stations 8 and 9, which suggests that there is a local minimum between the two mining regions (see also Fig. 11). What has been ignored up to this point is the clear-sky bias in the OMI measurements resulting from the requirement that only (near) cloud-free measurements be used, with no such restriction on the GB measurements. Averaging over clearsky data only may introduce a bias from the direct effect of generally faster photochemistry due to higher levels of sunlight, including the photolysis rate of NO 2 , which is important in the NO-NO 2 partitioning. There may also be an indirect effect if cloudiness is correlated with wind direction, and hence air mass origin. The most obvious way of avoiding this bias is to sample the GB data in the same way as OMI. However, this is complicated by the method used here to compute GB station averages: first calculating averages as a function of wind direction, and then averaging over these. Another method for removing clear-sky bias is to use a measure of cloud cover at each of the surface stations to screen GB measurements, analogous to the satellite measurements, with the source of cloud information being either a surface monitor or the OMI or MODIS-Aqua cloud fraction products (e.g. Geddes et al., 2012). A third alternative, and the approach used here, is to compare NO 2 and SO 2 from GEM-MACH with and without cloud screening. For this the total column model cloud fraction, considering all cloud types, was used. The GEM-MACH cloud fraction distribution is comparable to the OMI (a peak at cloud fraction of zero and a general decrease thereafter), with the exception that the GEM-MACH has a secondary peak for cloud fraction of 1, which OMI does not. In terms of the fraction of days with clear skies, OMI has about 45 %, and GEM-MACH about 35 %. Increasing the GEM-MACH cloud threshold so that 45 % of measurements are considered cloud-free, for consistency with OMI, leads to no appreciable differences in the results below.
Limiting NO 2 to clear skies leads to a low bias, relative to all-sky conditions, in an amount varying from 5 to 50 %. This is similar to the clear-sky bias observed over the greater Toronto area (Geddes et al., 2012) and consistent with a shift in the NO x partitioning to favour NO as a result of increased photolysis. For SO 2 the opposite effect was observed: cloud screening led to a high bias, between 5 and 25 %. Scavenging of SO 2 by clouds is an efficient loss mechanism, and thus higher SO 2 concentrations are expected when considering only cloud-free conditions. Average clear-sky biases are given in the Supplement.
The impact of smoothing and clear-sky bias on the GBsatellite comparison was assessed by applying the scaling factors for each to the EC-vmr's. These are also shown in Figs. 11 and 12, and the linear corrections and slopes are given in Table 4. The EC-vmr NO 2 increased at all stations and substantially improved the agreement: a scatter plot slope of 0.64 and a correlation coefficient of 0.97. Also of note is the ability of the EC-vmr's to capture the minimum in NO 2 between the northern and southern mining regions. In contrast, there was a smaller impact on SO 2 as the two effects somewhat cancelled (slope 1.01, correlation coefficient 0.91). The combined effect of the smoothing and clear-sky bias corrections explains a large portion of the GB-satellite difference in NO 2 and highlights their importance.
To conclude this section, it is noted that the stringent SZA threshold of 60 • (combined with the snow filter) means > 99 % of OMI SO 2 measurements considered are from April to September. For NO 2 , with a SZA threshold of 75 • , the breakdown is about 75 % (April-September) to 25 % (October-March), which allows for some seasonal analysis. The comparison above was repeated but for limited to these two 6-month "seasons". For April-September, the agreement between GB and EC-vmr's improved (slope 1.00, correlation coefficient 0.94), and for October-March, the agreement decreased (slope 0.47, correlation coefficient 0.68). The relatively poor agreement in winter could simply be due to larger SZAs or a result of the shallower BL heights in the winter. Part of the underestimate by the EC-vmr's could be related to the simplified treatment of surface reflection. Zhou et al. (2010) predict a 10-20 % underestimate in the winter when using the surface modelled as a Lambertian reflector instead of the more rigorous BRDF. A more detailed investigation into these differences is beyond the scope of this study. Significant low biases have been identified in current NO 2 and SO 2 retrieval products from the Ozone Monitoring Instrument (OMI) over the Canadian oil sands arising from a combination of its rapid development and small spatial scale. Air mass factors (AMFs) were re-calculated for this region based on updated and higher resolution input information. This information includes gas absorber profiles from the high-resolution (15 km × 15 km) GEM-MACH air quality forecast model, higher spatial and temporal resolution surface reflectivity from the MODIS satellite instruments, and an improved treatment of snow via a more precise determination of snow cover and more appropriate surface albedo when snow is present. The overall impact of these new Environment Canada AMFs led to increases in the peak NO 2 and SO 2 average vertical column density (VCD), occurring over the area of intensive surface mining, by factors of roughly 2 and 1.5, respectively. Due to a lack of validation profile or VCD data in this region, comparisons were made with long-term averages of NO 2 and SO 2 from in situ surface monitors at several WBEA stations. This was achieved by using the air quality model to map the EC-vmr VCDs to surface concentrations. The OMI-EC surface concentrations were able to capture the NO 2 and SO 2 spatial distribution of the in situ instruments. The absolute values were in reasonable agreement with OMI-EC surface NO 2 and SO 2 biased low by roughly 10-30 %. This level of agreement was improved by addressing complications in these comparisons. The first was the NO 2 high bias in the chemiluminescence instruments. The other two were the effects of smoothing by OMI and the clear-sky bias.

Atmos
This work is the first attempt to homogenize the OMI NO 2 and SO 2 data products through a consistent treatment of AMFs. It also examines the impact of interannual variability in a number of AMF-dependent parameters, including profile shape, surface albedo, and ozone column. The use of output from the high-resolution air quality model GEM-MACH was valuable in several regards: input into the AMF calculations, for performing the column to surface mapping, and also in assessing the sources of bias. Overall these results highlight the importance of using input information that accounts for the spatial and temporal variability at the location of interest when performing retrievals.