OMI air-quality monitoring over the Middle East

Using Ozone Monitoring Instrument (OMI) trace gas vertical column observations of nitrogen dioxide (NO2), formaldehyde (HCHO), sulfur dioxide (SO2), and glyoxal (CHOCHO), we have conducted a robust and detailed time series analysis to assess changes in local air quality for over 1000 locations (focussing on urban, oil refinery, oil port, and power plant targets) over the Middle East for 2005–2014. Apart from NO2, which is highest over urban locations, average tropospheric column levels of these trace gases are highest over oil ports and refineries. The highest average pollution levels over urban settlements are typically in Bahrain, Kuwait, Qatar, and the United Arab Emirates. We detect 278 statistically significant and real linear NO2 trends in total. Over urban areas NO2 increased by up to 12 % yr−1, with only two locations showing a decreasing trend. Over oil refineries, oil ports, and power plants, NO2 increased by about 2–9 % yr−1. For HCHO, 70 significant and real trends were detected, with HCHO increasing by 2– 7 % yr−1 over urban settlements and power plants and by about 2–4 % yr−1 over refineries and oil ports. Very few SO2 trends were detected, which varied in direction and magnitude (23 increasing and 9 decreasing). Apart from two locations where CHOCHO is decreasing, we find that glyoxal tropospheric column levels are not changing over the Middle East. Hence, for many locations in the Middle East, OMI observes a degradation in air quality over 2005–2014. This study therefore demonstrates the capability of OMI to generate long-term air-quality monitoring at local scales over this region.


Introduction
It is well established that poor air quality can significantly impact human health, ecosystems and agriculture, the built environment, and regional climate (Monks et al., 2009).The human and monetary costs associated with increasing levels of air pollution are substantial.For example, the World Health Organization estimates that globally nearly 7 million premature deaths were attributed to household and ambient air pollution during 2012 (WHO, 2014).Similarly, the Organisation for Economic Co-operation and Development estimated that outdoor air pollution is costing its 34 member states, plus the People's Republic of China and India, an estimated USD 3.5 trillion a year in terms of the value of lives lost and ill health (OCED, 2014).
Regulatory control of urban air quality is typically most effective when important target areas are extensively monitored, so pollutant emissions can be quantified and their atmospheric chemical processing well understood.In addition to in situ measurements made by ground stations and aircraft campaigns, satellite observations also form an important component of air-quality monitoring (Martin, 2008;Duncan et al., 2014).Satellite measurements of pollutants, such as nitrogen dioxide (NO 2 ), carbon monoxide (CO), tropospheric ozone O 3 , formaldehyde (HCHO), sulfur dioxide (SO 2 ), and aerosol particulate matter (PM), have been widely used to characterise their global atmospheric distributions, to quantify surface precursor emissions, and to evaluate local air quality (see e.g.Martin, 2008;Streets et al., 2013;Duncan et al., 2014, and references therein).Furthermore, with the advent of successive and longer-duration satellite missions, particularly the sequence of ultraviolet-visible (UV-VIS) instruments of GOME (Global Ozone Monitoring Experiment) (Burrows et al., 1999), SCIAMACHY (Scanning Imaging Absorption Spectrometer for Atmospheric CHartographY) (Bovensmann et al., 1999), and GOME-2 (Callies et al., 2006), together with the OMI (Ozone Monitoring Instrument) sensor (Levelt et al., 2006), there is a growing ability to study long-term changes in air quality from space.Some notable studies that have used satellite trace gas measurements to examine air-pollution trends include Richter et al. (2005), van der A et al. (2006, 2008), Ghude et al. (2008), De Smedt et al. (2010, 2015), Russell et al. (2012), Schneider and van der A (2012), Hilboll et al. (2013), Jin and Holloway (2015), Krotkov et al. (2016), Lamsal et al. (2015), Lelieveld et al. (2015), Schneider et al. (2015), and Duncan et al. (2016).
The Middle East is a region where long-term changes in air quality have probably been less well studied, in comparison to Asia, Europe, and North America.Besides frequent dust storms (Furman and Hadar, 2003), the region's air quality is characterised by year-long high ozone levels (Lelieveld et al., 2009), with an observed summer maximum (Liu et al., 2009;Zanis et al., 2014).The high ozone levels are in part due to long-range transport, strong local emissions, and favourable conditions for ozone photochemistry (Lelieveld et al., 2009).In situ instruments frequently record high pollutant concentrations in urban areas which exceed recommended guidelines (e.g.Modarres and Dehkordi, 2005;Nasralla and Seroji, 2007;Abdul-Wahab, 2009;Munir et al., 2013;Rashki et al., 2013).
The severity and variability of the region's air pollution can be directly observed from space.For example, several studies have reported appreciable trends in NO 2 vertical column over Middle Eastern cities prior to 2011, with increases of the order of 2-9 % yr −1 found over Tehran, 4-5 % yr −1 over Jeddah, 6-7 % yr −1 over Riyadh, and 10-20 % yr −1 over Baghdad, depending on the satellite instrument and the averaging period (van der A et al., 2008;Schneider and van der A, 2012;Hilboll et al., 2013).Similarly, De Smedt et al. (2010) also found an increasing trend of 1-3 % yr −1 in GOME and SCIAMACHY tropospheric HCHO columns over Tehran and Baghdad during 1997-2009 but a decrease of about 2 % yr −1 over Riyadh.
More recently, Lelieveld et al. (2015) examined annual changes in OMI NO 2 and SO 2 columns over the Middle East over 2005-2014, reporting that, after increases in the period 2005-2010, there was a reduction in these gases, either due to new regulatory legislation or due to falls in economic output associated with regional conflicts and geopolitical controls.In particular, decreases in NO 2 tropospheric columns of the order of 40-50 % were noted over Damascus and Aleppo since 2011, coinciding with the start of Syria's civil war.A similar regional-scale study by Krotkov et al. (2016) also observed that during 2005-2008 OMI NO 2 column increased by 20 % but remained approximately constant thereafter, whereas OMI SO 2 columns dropped by 20 % after 2010, only recovering to 2005 levels in 2014.A further OMI NO 2 column trend analysis over the region's major cities by Duncan et al. (2016) likewise reported decreases over Damascus and Aleppo, but of about 3-4 % yr −1 , with increases of about 2-6 % yr −1 elsewhere.However, besides these valuable studies, to the best of our knowledge, long-term changes in local air quality for many smaller Middle Eastern cities and towns have not been reported.In this study we aim to remedy this situation by determining the changes in air pollution over local population centres and also oil/energy infrastructure from space using a decade's worth of observations from the OMI instrument.Our target areas are therefore (1) cities and towns, (2) largescale oil refineries and ports, and (3) coal, gas and oil-fuelled electricity generating power plants.
To track the air quality over our specified targets, we use OMI tropospheric vertical column observations of NO 2 , HCHO, and glyoxal (CHOCHO), together with retrieved boundary layer column concentrations of SO 2 .Although each of these reactive species has different sources and sinks, they are all established key indicators of anthropogenic emissions, active photochemistry, and air pollution.For example, the dominant sources of NO x (=NO + NO 2 ) in the troposphere are the combustion of fossil fuels, biomass burning, emissions from soil, and lightning.Boundary layer SO 2 is predominantly generated by the burning of sulfur-laden fossil fuels and the refinement of sulfur ores; volcanic SO 2 emissions are typically injected high into the atmosphere, well above the boundary layer.The chemical reactions of NO 2 and SO 2 lead to the formation of nitrate and sulfate aerosols, which contribute to PM 2.5 (particulate matter with diameters < 2.5 µm), another critical pollutant (Kim et al., 2015).HCHO and CHOCHO are reaction products from the oxidation of anthropogenic, biogenic, and pyrogenic volatile organic compounds (VOCs); they can also be directly emitted from fires.Observed HCHO and CHOCHO distributions, therefore, contain the signature of underlying VOC emissions.Quantifying VOC emissions is important, as the oxidation of VOCs in the presence of high NO x and sunlight leads to the formation of tropospheric ozone, which is a major air pollutant contributing to photochemical smog as well as a key greenhouse gas and atmospheric oxidant (Monks et al., 2015).Apart from SO 2 , which has a lifetime of about 1 week, the trace gases are also relatively short lived (of the order of hours to a day), so the spatial displacements from emission sources are often small.
In this study, we have two broad goals: (1) to establish which locations have the highest pollution levels and (2) to use time series analysis in order to determine which locations have statistically significant trends and to quantify the trend magnitudes.The structure of this paper is as follows.In Sect. 2 we introduce the OMI data products.In Sect. 3 we discuss how the OMI data is gridded and describe the time series analysis.We present our results in Sect.4, with a discussion on their validity in Sect. 5. Finally, we conclude the paper in Sect.6.

The Ozone Monitoring Instrument (OMI)
The Dutch-Finnish Ozone Monitoring Instrument (Levelt et al., 2006) is a nadir-viewing UV-VIS 2-dimensional charge-coupled device (CCD) spectrometer, launched on board NASA's Aura satellite in July 2004.OMI orbits the Earth in a Sun-synchronous polar orbit, crossing the Equator at 13:30 local time (LT) in its ascending mode.The instrument has a 114 • field-of-view producing a 2600 km wide swath which contains 60 cross-track pixels that range in size from 14 × 26 km 2 at nadir to 28 × 160 km 2 at the swath edges.With these viewing geometry and orbital characteristics, OMI achieves global coverage daily (in nominal operational mode).However, from 2007 onwards, OMI's coverage has been considerably reduced due to problems with certain rows of its CCD detector (see Sect. 3.1 for further details).

SAO formaldehyde
The official NASA HCHO product is provided by the updated Smithsonian Astrophysical Observatory (SAO) retrieval, as described in González Abad et al. (2015) and evaluated in Zhu et al. (2016).HCHO slant columns are retrieved through a direct non-linear least-squares fitting of spectral radiances within the interval 328.5-356.5 nm.The retrieval algorithm includes dynamic calibration of solar and radiance wavelengths, the use of a daily radiance reference spectra, an under-sampling correction, and computation of commonmode residual spectrum.The cross sections of HCHO and other absorbers are fitted, together with a Ring effect correction, scaling and closure polynomials, and a spectral shift parameter.The retrieved slant columns are converted to vertical columns using air mass factors (AMFs) taken from lookup tables pre-computed using the VLIDORT radiative transfer model (Spurr, 2008), which uses a priori HCHO profiles from a global 2.0 • × 2.5 • GEOS-Chem chemistry transport model simulation (originally described in Bey et al., 2001).Effective cloud fraction (CFR) and cloud-top pressure (CTP) are taken from the OMI O 2 -O 2 cloud product (Acarreta et al., 2004), whilst the surface reflectivity for clear-sky scenes is extracted from the OMI mode Lambertian equivalent reflectivity dataset created by Kleipool et al. (2008).A daily post-processing normalisation correction (a function of latitude and detector row) is applied to reduce retrieval biases, minimise noise, and reduce cross-track striping.Observations with effective cloud fractions > 20 % are rejected.Over our domain of interest (20-80 • E, 10-50 • N) we find that the median uncertainty of a single measurement is about 60-70 %, with an interquartile range (IQR) of about 175 %.

DOMINO NO 2 vertical columns
Tropospheric NO 2 vertical columns are from the KNMI DOMINO v2 product (Boersma et al., 2011a, b), which is a well-established dataset (see e.g.Hains et al., 2010;Wang et al., 2016).The NO 2 slant columns are retrieved using the differential optical absorption spectroscopic (DOAS) technique by fitting the absorption cross sections of NO 2 , O 3 , and H 2 O, along with a synthetic Ring spectrum, to observed reflectance spectra in the 405-465 nm interval (Boersma et al., 2007).The retrieved total slant columns are then assimilated into the TM4 chemistry model (Dentener et al., 2003), to estimate and remove the stratospheric NO 2 component.Finally, tropospheric vertical columns are then derived by applying altitude-resolved AMFs from a pre-computed lookup generated with the KNMI DAK radiative transfer model using TM4 NO 2 profiles from a 2.0 • × 3.0 • global simulation.Surface reflectivity and cloud parameters are taken from Kleipool et al. (2008) and Acarreta et al. (2004), respectively.As with HCHO, observations with effective cloud fractions > 20 % are rejected.The median uncertainty of an individual measurement over our target region is about 60-70 %, with an IQR of about 54 %.

BIRA CHOCHO vertical columns
Glyoxal vertical columns are retrieved from OMI spectra using an updated DOAS retrieval originally devised for GOME-2 (Lerot et al., 2010).The retrieval uses a daily mean Earthshine Pacific Ocean radiance spectrum and applies a row-dependent wavelength calibration.An initial pre-fit of liquid water optical depth (in 405-490 nm) is performed, before the absorption cross sections of CHOCHO and other interfering absorbers are fitted in the spectral window 435-490 nm to retrieve the slant columns.Tropospheric vertical columns are obtained using pre-computed look-up tables of altitude-resolved AMFs, which use a priori profiles over land from a global 2.0 • × 2.5 • IMAGES simulation (Müller and Stavrakou, 2005).Over the oceans a single oceanic profile is used for the AMFs, derived from airborne MAX-DOAS (multi axis differential optical absorption spectroscopy) measurements over the Pacific Ocean (Volkamer et al., 2015).Unlike the retrievals of HCHO and NO 2 , the cloud effects are not accounted for in the AMF calculation (via the independent pixel approximation); instead, only clear-sky observations are retained by rejecting scenes with effective cloud fractions > 20 %.The median uncertainty of an individual measurement is about 50-60 %, with an IQR of about 250 %.

NASA sulfur dioxide
To assess changes in anthropogenic SO 2 , we use planetary boundary layer (PBL) columns (NASA product OMSO2 v1.2.0) determined using a principal component analysis (PCA) retrieval that is sensitive to surface concentrations (Li et al., 2013).The algorithm applies the PCA technique to OMI UV Sun-normalised radiances (310.5-340nm) over an SO 2 -free region (in the equatorial Pacific) to establish the principal components of the main physical and measurement spectral features that do not correspond to SO 2 absorption.A set of principal components and SO 2 radiance Jacobians (i.e. that describe the radiance sensitivities to changes in the SO 2 column) are then iteratively fitted to observed OMI radiances to obtain the slant column density.An estimate of the boundary layer SO 2 is then calculated based on the assumptions of the vertical SO 2 distribution (Krotkov, 2014).The algorithm applies the PCA technique to each OMI detector row individually and uses the VLIDORT radiative transfer code to compute the SO 2 Jacobians, with RT model inputs based on a fixed atmospheric profile and climatological SO 2 profile over the summertime eastern US, the latter corresponding to an effective AMF of 0.36 (Fioletov et al., 2016).Scenes with strong ozone absorption of > 1500 Dobson units (1 DU = 2.69 × 10 16 molecules cm −2 ) are excluded due to spectral interference in the retrieval.For this study we only use SO 2 data from OMI rows 4-54 (0-based) and with a cloud radiance fraction < 0.3 (Krotkov, 2014), which corresponds to an effective cloud fraction of about 15 % (Krotkov et al., 2016), broadly consistent with the other species.Although SO 2 over the Middle East is mostly unaffected by volcanic emissions over 2005-2014(Krotkov et al., 2016)), transient volcanic SO 2 enhancements were removed if they exceeded a threshold of 5 DU (Fioletov et al., 2011).The estimated uncertainty of the SO 2 PBL is about 0.5 DU in the tropics and ∼ 0.7-0.9 at higher latitudes, based on the analysis of the root mean square (RMS) and standard deviation values for instantaneous field-of-view (IFOV) observations in different latitudinal bands (see Li et al., 2013;Krotkov, 2014).

Data gridding
We monthly average the OMI observations onto a highresolution 0.05 • × 0.05 • grid using an area-weighting tessellation algorithm (Spurr, 2003;Hewson et al., 2015).The gridding algorithm properly accounts for the areal proportions of grid cells underlying the satellite footprint and inversely weights each observation according to the measurement uncertainty and OMI ground pixel size.Spatial zoom orbits are not included in our analyses, nor are scenes with solar zenith angles > 70 • .We also fol-low the data recommendations provided with each product to reject non-optimum observations.The quality of the level 1B radiance data from certain rows of OMI's CCD detector is known to be affected by blockage effects, wavelength shifts, and stray light originating from outside the nominal field of view.This dynamic behaviour is the well-known OMI row anomaly (http://www.knmi.nl/omi/research/product/rowanomaly-background.php), which impacts atmospheric retrievals from the affected rows.The temporal variability of the row anomaly is also known to compromise the derivation of long-term trace gas trends (De Smedt et al., 2015).For this reason, we follow a similar approach to De Smedt et al. (2015), by using the OMI XTrackQualityFlags (Dutch Space, 2009) to construct a static mask based on the most affected rows at the end of 2013, to then discard row-anomaly observations over the entire 2005-2014 period.This quality filtering largely removes any statistically significant sampling trends in the generated monthly datasets, except for the SAO HCHO and DOMINO NO 2 data.Therefore, for these gases we further use the XTrackQualityFlags to identify affected observations not flagged by the static mask, whose measurement uncertainties we then increase by a factor of 1000, so that these observations are included in the analysis without affecting the monthly averaged fields (i.e. they are assigned a very low weight in the averaging).On average there are typically 20-35 samples per grid cell per month.To reduce noise in the monthly gridded data, we smooth the data with a 0.15 • × 0.15 • Gaussian filter of 1σ width.For the noisier CHOCHO fields, a 2σ Gaussian is used.The spatial smoothing enhances localised "hot spots" and is preferable to averaging the data onto a coarser-resolution grid where the atmospheric signatures of target features can be lost.The median uncertainties of the gridded data are about 4 % for NO 2 , 6 % for HCHO, 27 % for CHOCHO, and 40 % for SO 2 .
Figure 1 shows the 2005 annual distributions of each species over the broad Middle East region.Clearly visible in the NO 2 , HCHO, and CHOCHO maps are the major urban areas of Riyadh, Baghdad, Tehran, Jeddah, and pollution enhancements along the eastern and northern coasts of the Persian Gulf.Intense SO 2 hotspots are found over the Jeddah and Mecca region, Kharg Island (in the Gulf), and near Kerman (Iran), where the Sarcheshmeh Copper Complex smelter is found.

Time series construction
To determine the temporal variability of the OMI vertical column data over urban areas we use the Global Rural-Urban Mapping Project Version 1 (GRUMP v1) settlement points database (Balk et al., 2006;SEDAC, 2015) to identify the geolocation coordinates of 818 cities and towns that have populations ranging from over 50 to nearly 7 million inhabitants (as determined for the year 2000).We then construct a 10-year time series (of 120 months) by averaging those monthly gridded data that lie within ±2 grid cells of the city or town location, which corresponds to a radial distance of approximately 10 km from the urban centre.Visual inspection of these spatial masks overlaid onto Google Earth imagery shows this filtering criteria is well suited to capture the extent of most urban areas (see e.g.Fig. S1 in the Supplement).For Baghdad, Riyadh, and Tehran, which have larger urban spread, the average of ±4 grid cells is used instead, consistent with an approximate 20 km radial distance.In addition to the trace gas vertical columns, we also construct coincident time series for the associated cloud fraction, cloud-top pressure (or height), AMFs (where appropriate), and number of grid-cell samples.This helps clarify whether any observed trends in the trace gases are real, by applying the same time series analysis to all the retrieved parameters.
We use this same approach to examine air-quality variability over oil refineries, oil ports, and power plants.The 2010 Oil Refining Survey (Kootungal, 2010) was used to identify 41 major crude oil refineries, whereas the Global Energy Observatory (GEO) free online resource (http:// globalenergyobservatory.org) was used to locate 18 oil ports and 155 power plants.Where different types of power plants (e.g. oil versus gas fuelled) were closely co-located, a single geolocation coordinate was used to mark that target.The geographical distribution of the selected targets is shown in Fig. 2 (left panel), and, as an example, Fig. 2 (right panel) shows the urban spatial filtering masks, applied to the observed 2005 annual OMI NO 2 distributions over northern Iran.The latter figure shows that this approach treats the larger cities, such as Tehran, as individual urban regions, with intention of resolving trends of separate districts, and that in some cases the spatial masks can partially overlap for targets which are close to one another.
Considering all locations, we find that by applying a weighted average to the data within ±2 grid cells of the target, the median uncertainties of the trace gas vertical column time series data points reduce to 0.37 ± 0.26 % for NO 2 , 0.73 ± 0.64 % for HCHO, 1.82 ± 1.90 % for CHOCHO, and 1.63 ± 2.58 % for SO 2 .

Time series analysis
Each time series consists of monthly mean OMI observations, y(t), over a given target site, where time t is in fractional years.We analyse each time series individually for each location following a consistent procedure.First, we filter the time series data for outliers, rejecting observations that lie beyond 2.5 median absolute standard deviations (Leys  et al., 2013).This helps to reduce difficulties in the analysis and interpretation of noisy vertical column data when encountered; the filtering is not required for AMFs, cloud parameters, or the number of samples.Second, we linearly interpolate across missing data points, applying a quality filter, so that if > 20 % of (randomly scattered) values are missing from the time series, then the dataset is rejected.Typically 0-25 % of data points may be missing from a time series, but, on average, we find that less than 6 % are missing.Third, we then fit to the data a model function, F (t), consisting of a linear component plus a four-term harmonic Fourier series, defined as follows: where µ is the mean value of the time series data at time t = 0, ω is the linear trend of the variable (per year), and the parameters A n and B n are the Fourier series coefficients that essentially model seasonal and inter-seasonal variability.Many other trend studies have fitted very similar functions (see e.g.van der A et al., 2006;Gardiner et al., 2008;De Smedt et al., 2010;Hilboll et al., 2013;Jin and Holloway, 2015).We fit this model to the data using a non-linear least-squares Levenberg-Marquardt algorithm, which generates an estimate of the fit parameters µ, ω, A n , and B n plus their uncertainties and covariance.Fourth, we check that the trend ω is real, i.e. significant at the 95 % confidence level.The generally accepted rule for determining whether a trend is statistically significant is |ω/σ ω | > 2, provided the lag-one autocorrelation of the fit residual is small (Weatherhead et al., 1998;van der A et al., 2006).To determine the precision (σ ω ) of the trend, we follow the approach of Gardiner et al. (2008) by bootstrap resampling (with replacement) the initial fit residuals to reconstruct the fitted function with representative noise.The model function F (t) is then refitted to this data, and the fit parameters are recalculated.This process is repeated 2000 times to build a sampling distribution for each of the fit coefficients, with the difference between the 2.5th and 97.5th percentiles representing the coefficient's associated 2σ uncertainty (De Smedt et al., 2010).The advantage of the bootstrap resampling method is that it enables non-normally distributed data to be treated robustly (Gardiner et al., 2008).
Fifth, we additionally filter the initial fit residuals with short-and long-term filters, to derive several other quantities of interest.For this we follow closely the curve fitting routines adopted by the NOAA's Earth System Research Laboratory (ERSL) Global Monitoring Division (GMD), which are based on the original study of Thoning et al. (1989) and are fully described online (at http://www.esrl.noaa.gov/gmd/ccgg/mbl/crvfit/crvfit.html).The general approach is to transform the residuals into the frequency domain using a Fast Fourier Transform (FFT), apply a low-pass filter function to the frequency data, and then transform the filtered data to the time domain using an inverse FFT.We use a low-pass frequency filter H (f ) of the form where f is the frequency (cycle per year) and f c is the frequency response of the filter.We filter the initial fit residual twice, once with a short-term cut-off value for smooth-ing the data and once with a long-term value to remove any remaining seasonal oscillation and to track inter-annual variability.Once the residual has been filtered, several curves can be constructed.
1.A smoothed function fit F S (t), which is the function fit F (t) plus the residual filtered using the short-term cutoff value.Its uncertainty is given by σ 2 , where σ S is the uncertainty of the short-term filtered residual and σ F is the uncertainty of the model function estimated from the covariance matrix of the fit parameters using error propagation.
2. A long-term trend fit F T (t), which is the linear component of F (t) plus the residual filtered using the longterm filter.It represents the long-term trend with the seasonal cycle removed.Its uncertainty is σ 2 , where σ L is the uncertainty of the long-term filtered residual.
3. A de-trended seasonal cycle F C (t), which is computed by subtracting the long-term trend fit from the smoothed function fit, that is, It represents the annual seasonal oscillation with any long-term trend removed.Its uncertainty is given by σ 2 The average seasonal amplitude is defined as the mean peak to trough difference of F C (t).
4. A growth rate curve F G (t), which is the rate of change of the long-term trend F T (t), determined using a threepoint (quadratic) Lagrangian interpolation to compute the first derivative.Its uncertainty is given by σ 2 G = 2 × σ 2 F T since the derivative is approximately equal to the taking of the difference of two data points 1 year apart and plotting this difference midway between the two points.The average growth rate, G, is then simply the median value of F G (t).
The statistical uncertainty of the residual filters (σ S and σ L ) are calculated following Thoning et al. (1989), via where σ rsd is the residual standard deviation, r (k−j ) are the lags in a first-order auto-regressive process (defined as r(k) = r(1) k for k = 1, 2, . ..), and n c is the number of filter weights.The filter weights, c i , are the values of the impulse response functions of the filter transfer functions, computed by applying the filters to a delta function impulse in the time domain (Thoning et al., 1989).We compute the filter uncertainties for the short-and long-term cut-off values individually and apply them to calculate the uncertainties in the derived curves.In this study, we use short-term and long-term filters for 200 and 667 days, respectively.Sensitivity tests using filters of 100, 150, 500, and 720 days indicate that these two filter values offer the best compromise for slightly higher correlations between the data and fitted curves versus slightly smaller curve uncertainties.Furthermore, we also evaluated the long-term filter values (500, 667, and 720 days) on the DOMINO NO 2 and SAO HCHO products over urban areas (818 targets) to assess their impact on the growth rate (note the filters do not affect the fitted linear trend µ).We found differences in the growth rates are small when using the 500 and 720 day filters, compared with the default 667-day filter, e.g. for NO 2 the differences are less than 1 %, and for HCHO they are 1-2 %. Figure 3 shows an example of a time series fit to observed NO 2 data over Dahuk in Iraq (43.00 • E, 36.87 • N, population: 65 683), where a statistically significant large upward linear trend of 2.77 ± 0.35 × 10 14 molecules cm −2 yr −1 is found.This corresponds to a linear growth of 12.23 ± 1.54 % relative to the observed 2005-2014 median VCD.In this example, |ω/σ ω | = 7.9, and the uncertainties of the trend (F T ) and smoothed curves (F S ) are about 5 % and 7.5 %, respectively.The median growth rate G is 12.44 ± 8.18 % yr −1 , whilst the mean seasonal amplitude is 1.78 ± 0.32 × 10 15 molecules cm −2 (about 79 ± 14 % relative to the median column).A similar analysis of the coincident time series of the NO 2 AMF, cloud fraction, cloud-top pressure, and number of samples reveals no other significant trend.This indicates that the upward growth in NO 2 is not caused by a trend in any other retrieval parameter and is real at the 95 % confidence level.This observed increase in NO 2 over the city is likely linked to its recent population growth (estimated at 280 400 people in 2012, see e.g.MOP-KRG, 2012) and rapid urban expansion (Mustafa et al., 2012).Although no underlying trends have been reported in the OMI data (Boersma et al., 2011b;Li et al., 2013;González Abad et al., 2015), as a precaution we performed the time series analysis on gridded OMI data over the remote Pacific Ocean (60 • N-60 • S, 90-170 • W).No statistically significant trends were found for any species.

Results
In this section we present our main results, with Tables 1-4 providing a concise summary of the analysis.Tables S1 to S2 (both excel files) in the Supplement provide a more complete classification of the analysis, where the highest ranked median levels and absolute linear trends for each target category are tabulated.Only the top 50 ranked locations are given for urban and power plants categories.

Average pollution levels
We use the observed median vertical columns determined for all locations within each category (i.e.not the fitted µ in Eq. 1) to assess average pollution levels (see Table S1).The overall 2005-2014 average median columns are 26.98 ± 1.85 × 10 14 molecules cm −2 for NO 2 , 3.97 ± 0.30 × 10 15 molecules cm −2 for HCHO, 18.71 ± 2.09 × 10 13 molecules cm −2 for CHOCHO, and 0.21 ± 0.09 DU for SO 2 .Relative to the overall median column, the highest mean NO 2 values are found over urban areas and oil ports, which are about 5 and 15 % higher than refineries and power plants, respectively.Whereas for HCHO, we find that the highest average columns are found over oil ports, with values about 5 % higher than those for refineries, 15 % higher than urban areas, and 16 % higher than power plants.The highest average CHOCHO levels are also found over oil ports, with refineries, power plants, and urban areas being about 5, 16, and 25 % less, respectively.For SO 2 the highest average values are found over oil ports and refineries, which are 57 % higher than power plants and 81 % higher than urban regions.Hence, the general conclusion we can make is that average trace gas levels tend to be highest over oil ports and refineries (with the exception of NO 2 ).
The interquartile range (IQR) of the observed median columns gives some indication of the spread of the trace gas levels within a given target category.For NO 2 the IQR varies between 68 % for urban areas, 59 % for refineries, 45 % for oil ports, and 53 % for power plants (as a percentage of the overall average median values).The higher IQR over urban areas may reflect the larger variability of underlying NO x emissions between different towns and cities.For HCHO, the IQR only varies between 13 % for urban areas and ports, 20 % for refineries, and 26 % for power plants.
The similar IQR values likely indicate a lower variation in HCHO sources and sinks over different locations.For CHO-CHO, the IQR varies between 25 % for urban areas to 51 % over refineries, with the IQR over ports and plants about 46-47 %.For SO 2 , the IQR is lowest over urban areas (42 %) but increases substantially over ports (100 %), power plants (114 %), and to a maximum over refineries (148 %).
Figure S4 shows the statistical box-and-whisker plots of the observed median vertical columns over urban targets for each country.The highest average pollution levels are typically found in Bahrain, Kuwait, Qatar, and the United Arab Emirates (UAE).For example, NO 2 columns over Bahrain, Kuwait, and the UAE are 52, 66, and 53 % above the median urban level of 28.13 × 10 14 molecules cm −2 , respec- Considering average NO 2 levels over individual urban locations, we find Kuwait has 25 settlements in the top 50 highest ranked places, whereas Saudi Arabia has 9 locations, Bahrain 5, Iran 5, and the UAE 5.The highest average NO 2 columns were found over Funtas (in Kuwait), with a median value of 63 × 10 14 molecules cm −2 .Funtas is one of several coastal settlements adjacent to the Mina Al-Ahmadi and Mina Abdullah refineries, which all have high average NO 2 levels.These two refineries had median vertical column levels of 55.07 × 10 14 molecules cm −2 , which were the highest of all refineries.Similarly, the adjacent Mina Al-Ahmadi port had the highest columns over all ports (of 64.11 × 10 14 molecules cm −2 ).Thus, this general area is a particularly intense NO 2 hotspot, as indi- cated in Fig. 1, due to the close proximity of several emission sources.For comparison, the highest NO 2 columns over power plants were found at the Tarasht Shahid Firouzi power station near Tehran, which had a median level of 74.18 × 10 14 molecules cm −2 .For HCHO, Iran with 16 locations, Yemen 13, and Saudi Arabia 6 dominated the top-50 ranked urban settlements.The five locations with the highest median values were located in Yemen, along its western coast, with the highest columns found over al-Marawiah, which had an average value of 5.28 × 10 15 molecules cm −2 .Over refineries, the highest median HCHO columns of 4.86 × 10 15 molecules cm −2 were detected at the Umm Said refinery in Qatar.In addition, 6 of the highest 10 ranked power plants were also located in Qatar, particularly those near Doha and Mesaieed.However, higher median column values of 4.99 and 5.03 × 10 15 molecules cm −2 were found over the Bandar-e Khomeini port and the Petroshimi and Bandar Imam power plants, respectively.Both of these sites are located close to the Bandar Imam Petrochemical facility in south Iran, which may indicate a potentially strong VOC source.
The top-50 ranked urban locations for SO 2 were dominated by Kuwait (26 locations) and Iran (13 locations).The highest median level of 0.73 DU was found at Rafsanjan (Iran), which is located close to the Sarcheshmeh copper mine and smelter facility.The coastal settlements adjacent near the Shuaiba, Mina Al-Ahmadi, and Mina Abdullah refineries (in Kuwait) again had high pollutant levels, typi- The highest average CHOCHO values over urban areas were found mostly in Iran (17 locations), Saudi Arabia (13), and Kuwait (17).The cities of Mecca and Tehran have the largest average levels of 44.32 and 39.44 × 10 13 molecules cm −2 , respectively.The highest levels over refineries are found at Tehran (34.87 × 10 13 molecules cm −2 ) and the three main coastal Kuwait refineries of Shuaiba, Mina Al-Ahmadi, and Mina Abdullah (30.34-30.80× 10 13 molecules cm −2 ).
The oil port of Mina Al-Ahmadi also has the highest median level of 30.48 × 10 13 molecules cm −2 .The Mecca OCGT power plant recorded the highest level of 43.81 × 10 13 molecules cm −2 ; the other notable CHO-CHO levels over power facilities were at Tarasht Shahid Firouzi and Besat (in Tehran), with columns of 42.72 and 39.85 × 10 13 molecules cm −2 , respectively.

Seasonal variability
To assess the the seasonal variability over each target, we determined the average peak-to-peak difference of its cor-responding time series and used these values to compute the median seasonal amplitude over all locations, within each target category.For NO 2 , the median seasonal peak-topeak amplitudes are very similar irrespective of target categories.Typical amplitudes of 10-12 × 10 14 molecules cm −2 are observed (about 39-46 % relative to the overall median).The highest seasonal peak-to-peak amplitude of 66 × 10 14 molecules cm −2 was found over the Besat Thermal Power Plant in Iran, whereas the lowest seasonal peakto-peak amplitude of 1.49 × 10 14 molecules cm −2 was determined over Ataq in Yemen.Unlike NO 2 , HCHO exhibits greater seasonal variability (relative to the overall median) with seasonal amplitudes that vary from 3.02-3.43× 10 15 molecules cm −2 (or about 60-86 %), with a high of 5.65 × 10 15 molecules cm −2 over Dahuk in Iran and a low of 1.20 × 10 15 molecules cm −2 over the Salalah OCGT power plant in Dhofar (Oman).Similarly, CHOCHO also exhibits large seasonal variability as median amplitudes vary from 14.38-16.48× 10 13 molecules cm −2 , or 77-88 %, with the highest amplitude of 44.53 × 10 13 molecules cm −2 over Qaemshahr (Iran) and lowest amplitude over Bilin (Palestine) of 8.89 × 10 13 molecules cm −2 .The seasonal variations in HCHO and CHOCHO are predominately driven by seasonal variations in biogenic VOC emissions (Müller et al., 2015).Even larger seasonal variability is demonstrated for SO 2 , as amplitudes range from 0.21-0.44DU (100-210 %), with the lowest amplitude of 0.09 DU over Al Ghaydah (Yemen) and highest over 1.15 DU over Shahreza (Iran).We find that the SO 2 and CHOCHO retrievals are generally noisier than their NO 2 and HCHO counterparts, and thus their corresponding time series, even after outlier filtering and residual smoothing, are noisier too (see e.g.Figs.S3-S4).Thus, the higher seasonal amplitudes found in the CHOCHO and SO 2 data may simply reflect larger point-to-point variability in their time series data.

Linear trends µ
For NO 2 , statistically significant real linear trends were determined for 198 of 818 urban locations (a detection rate of 24 %).However, in Palestine only 1 out of a potential 274 urban targets had a real trend.Neglecting the Palestine results, in this instance, increases the overall detection rate to 36 %.The corresponding "detection percentages" for oil refineries, oil ports, and power plants are 42, 33, and 37 % respectively (Table 1).Urban trends ranged from −3.05 ± 0.93 to 12.23 ± 1.54 % per year (relative to the observed median column), with the highest trend in Dahuk (Sect.3.3, Fig. 3).Generally, Iraq and Iran have the highest linear NO 2 trends.For example, 5 of the top 7 highest linear trends were found in Iraq, whereas 21 Iranian cities appeared in the top 50 (see Table S2).Overall, the median linear trend is about 3 % yr −1 , although, for the top 50 highest ranked locations, the trends were of the order of 2-12 % yr Trends over refineries ranged from about 2-6 % yr −1 , with a maximum found over the Daura refinery, near Baghdad; its capacity is smaller than some of the larger refineries, which may indicate the influence of other nearby NO 2 sources.For Iran, Iraq, and the UAE, about 60 % of the refineries studied showed an NO 2 increase.Trends over oil ports ranged from about 2-9 % yr −1 , with the highest trend found over Umm Qasr in Iraq.Trends over power plants ranged from about 2-8 % yr −1 , with 5 of the top highest 10 trends found in Iran.The highest trends for power stations were detected over the two Sabiya plants in Kuwait which had trends of 8.13 ± 0.14 % yr −1 .Figure 4, which shows the geographical distribution of NO 2 trends (panel a), indicates that there are several local and regional locations where consistent increases of NO 2 are observed.For example, these are several targets with increasing trends situated close to Riyadh, Tehran, Baghdad, and Muscat (in Oman), as well as Isfahan and Yazid (both in Iran).Regional trend hot spots include the areas south-west of Tehran and a thin corridor stretching from northern Jordan to Lebanon, passing through southwest Syria.
For HCHO, statistically significant real linear trends were determined for only 4 % of urban locations.The corresponding trend detection percentages are 15 % for oil refineries, 22 % for oil ports, and 17 % for power plants (Table 2).Urban trends ranged from 2-7 % yr −1 .The highest absolute linear urban trend was found at al-Wakrah in Qatar (Fig. S2), although this was only about 5 % yr −1 in relative terms; the highest percentage trend was found over Attaif in Saudi Arabia (6.95 ± 2.42 % yr −1 ; see Table S2).Generally, Saudi Arabia has the highest number of settlements with increasing HCHO trends (13 of the corresponding 34 trends found).Elsewhere there were 6 trends found in Oman, 5 in Qatar, 3 in Iran, 3 in the UAE, 1 Iraq, and 1 Israel.Only six trends were detected over oil refineries, which ranged from 2-3.5 % yr −1 , with a maximum found over the Ruwais refinery in the UAE.Only four trends over oil ports were found, which ranged from 2-4 % yr −1 , with the highest trend found over the Ras Laffan port at Al Khawr (Qatar).Trends over power plants ranged from 2-7 % yr −1 but were only found over Saudi Arabia (8 stations), Iran (5), Qatar (4), and the UAE (5). Figure 4 (panel b) shows that the target locations with trends are mostly found along the western Gulf coast, particularly along the Saudi Arabian coast near Ad-Dammam, and also near Doha in Qatar.
For SO 2 , very few trends were detected, only 2 % over urban targets, 7 % over refineries, 11 % over oil ports, and 6 % over power plants (Table 4).Over urban areas, trends ranged from −61.64 ± 28.73 % yr −1 in Hamismusayt (also known as Khamis Mushait, Saudi Arabia) to 118.49 ± 42.78 % yr −1 in Azrashahr (Iran); the latter value is inflated by a median SO 2 column of approximately zero over its location.Notably, 11 of the 18 trends were detected in Iran.Over refineries, three trends were detected of about 9-15 % yr −1 .Over oil ports, the Iranian oil ports of Bandar-e Khomeini and Kharg Island showed decreases of about 6 % yr −1 .Lastly, over power plants 5 locations had decreasing trends of −7 to −22 % yr −1 , and 4 stations had increasing trends of about 6-305 % yr −1 (again the latter value has a median SO 2 column of approximately zero).Figure 4 (panel c) shows the geographical distribution of the SO 2 trends; it is evident that several closely located targets exhibit similar trends, particularly (1) near Abba (in south-west Saudi Arabia), (2) around Neka, Sari, and Behshahr in northern Iran, which are close to the Shamid Salimi power plant at Nowzarabad, and (3) Tarbiz, Azrashahr, and Maragheh in Iran.

Median growth rates G
Figure 5 show the growth rates that correspond to sites with statistically significant linear trends (i.e.where we can be sure that the trend is not attributed to variations in other re-Atmos.Chem.Phys., 17, 4687-4709, 2017 www.atmos-chem-phys.net/17/4687/2017/trieval parameters and where there is enough signal in the data themselves to give a meaningful trend), and Fig. 6 shows their difference, here defined as the linear trend (in % yr −1 ) minus the growth rate (in % yr −1 ).Differences between the linear trend and growth rate occur due to the latter tracking inter-annual variations in the data not accounted for in the linear part of Eq. (1) (see Thoning et al., 1989).While the overall median differences are quite small (typically less than 2 % yr −1 ), for individual locations they can be much larger.Typically for NO 2 and HCHO, the differences can range between ±5 % yr −1 (Table S3).Such occurrences are nonnegligible.For example, the NO 2 linear trend and growth rates over Ibri (in Oman) were about 3 and 8 % yr −1 respectively (i.e.their difference is 5 % yr −1 ).For SO 2 the differ-ences can be even more substantial (i.e.tens of percent per year), but this is likely an effect from filtering much noisier fit residuals (e.g. as shown in Fig. S4).Furthermore, Fig. 6 shows that there are no clear spatial patterns or coherence to the geographical distribution of the linear and growth differences.This raises the interesting question for AQ trend studies: should one consider linear changes in a trace gas or use a growth rate?The former does not capture shorter-term variations in growth, but the latter is more susceptible to the choice of fit residual filter and has larger uncertainties.5 Trends: fact or fiction?
It is, perhaps, tempting to try explain the trends in terms of variations in underlying emissions and/or atmospheric chemistry.However, it seems more prudent to carefully assess the validity of our presented results instead.We do this by examining our analytical uncertainties, by performing a series of sensitivity tests, and by discussing our results in the context of other recent satellite studies.

Uncertainties
Calculated uncertainties in the linear trend and other derived curves mostly reflect the noise of the trace gases time series.For NO 2 , and considering only those locations with real statistically significant trends, we find that the median uncertainty in the linear trend is about 1 %.The median uncertainties are about 4 % for the trend curves, 6-7 % for the smoothed curves, 5-6 % in the growth curves, and about 6-8 % for the seasonal cycles, respectively.For HCHO, the corresponding values are very similar: 5 % in the trend curves, 7 % in the smoothed, 7 % in the growth, and 9 % in the seasonal curves, with the median uncertainty in the linear trend at about 1 %.For CHOCHO, uncertainties are slightly higher, being 10 % in the trend curves, 15 % in the smoothed, 14 % in the growth, and 18 % in the seasonal curves.The median uncertainty in the linear trend is about 2 %.The uncertainties are noticeably higher for SO 2 , with errors of 10-23 % in the trend curves, 15-34 % in the smoothed, 14-32 % in the growth, and 18-40 % in the seasonal curves.The median uncertainty in the SO 2 linear trend is about 2-5 %.Thus for the most part, the uncertainties in the linear trends are comparable or smaller than than the trends themselves.However, uncertainties in the growth rates are comparable or higher than the derived average growth trends.

Sensitivity tests
Subtle differences in the data analysis approach may affect the choice of locations where trends are detected or may change existing trend directions and magnitudes.Given this situation several additional tests were carried out to determine the sensitivity of our results to various parameters.These are outlined below.Tables S4-S7 record how many trends were now detected in each test and whether or not they were previously detected in the default approach (as outlined in Sect.3.3).
-Test 1: construct each 10-year time series using a mask of ±4 grid cells (∼ 20 km radius around each target) instead of the default ±2 grid cells (∼ 10 km radius around each target).
For NO 2 , there were 22 and 12 extra trends detected over urban and power plant targets, respectively, compared to our default analysis (shown in Table 1).The number of NO 2 trends detected over oil ports and refineries were unchanged.However, some trends previously detected were now missing, compensated by new trends over other locations.For example, in the case of urban targets, 170 locations in this test were also previously detected in the default scenario, but 28 of the previous trends were missing, and 50 new trends were detected (see Table S4).Nevertheless, median linear trends were of similar order (about −2 to 10 % yr −1 ).The highest urban NO 2 trends were still at Dahuk and Irbil (in Iraq) but were marginally changed from 12.23 ± 1.54 % yr −1 to 10.39 ± 1.42 % yr −1 and from 8.76 ± 1.17 % yr −1 to 8.18 ± 1.16 % yr −1 , respectively.The Umm Qasr and Zubayr oil terminals were still the highest ranked oil ports, at 8-9 % yr −1 .Similarly, the Sabiya CGGT and OCGT power plants still had the highest trend of about 8 % yr −1 .
For HCHO this resulted in an extra 25 trends being detected (spread over all categories); nevertheless, the median linear trends were still about 2-3 % yr −1 .The majority of trends previously found in the default case were still present in this test.However, notably, the number of detections over oil ports doubled from 4 to 8, with the highest trend now found over Sitrah in Bahrain at 4.03 ± 1.14 % (the trends over the Ras Laffan and Al-Ruwais oil ports were now displaced to second and third highest).That said, both al-Wakarah and Doha still had the highest (absolute) trends for cities (of about 3-4 % yr −1 ), and again the Ras Abu Fontas OCGT power plant in Doha (Qatar) also topped the highest ranked trends over power stations.
For SO 2 the highest urban trends were still over Ilam and Tabriz (Iran) and also Baghdad and Irbil (Iraq).Although 18 urban trends were still detected, only 9 were previously found in the default scenario (including those over Sari, Neka, and Maragheh), thus 9 new trend locations were found.Over oil ports the highest trend was now at Umm Qasr (23.09 ± 5.87 % yr −1 ), as a trend was not detected at Bandar-e Khomeini (BIK); Kharg Island still reported a decrease of about −6 % yr −1 .SO 2 trends were found over the same three refineries and were still about 9-12 % yr −1 .Median trends over urban areas changed from 9.8 to 12.0 % yr −1 , but refineries and power plants remained at about 9 and 6 % yr −1 , respectively.The Sabiya CGGT and OCGT power plants in Kuwait were still the highest-ranked trends.
For CHOCHO, in addition to the trend of al-Hawr, only an extra urban trend at Shushtar (Iran) of −5.33 ± 2.58 % yr −1 was detected.An additional trend at the Ras Laffan (Qatar) refinery of about −4.27 ± 1.91 % yr −1 was also found; the trend over the Ras Laffan power stations was also present as in the default scenario.
In summary, although increasing the averaging radius to ∼ 20 km around each target has a small impact on overall median trends, it can result in differences at individual locations.Generally however, most of the trends found in the default approach were still present.
-Test 2: use of different cloud fraction filters.
For this test, we used only OMI NO 2 and HCHO observations with a cloud fraction ≤ 40 % in the analysis; CHOCHO and SO 2 and were untested, as the use of a strict cloud threshold is advised (Krotkov, 2014).
For NO 2 , there were 45 fewer trends detected in total (mostly over urban locations, 165 compared to 198 previously).The majority of trends in the default analysis were still present -hence the relative median trends were still about 3-4 % yr −1 .Furthermore, the top ten urban ranked trends were also mostly unchanged, as were those over oil ports, refineries, and power plants, compared to the default scenario.
For HCHO there were two less detections in total, with the highest urban trend found over Ardebil (Iran) at 7.06 ± 2.53 % yr −1 .Again the majority of locations with trends were consistent with those detected in the default scenario; hence the median trends are still about 3 % yr −1 .
Thus, in summary, increasing the cloud fraction for NO 2 and HCHO slightly decreases the number of trends found, and, on average, has a small impact on overall median trends.Differences at individual locations can still occur, compared to the default approach.
In another test, we also applied a stricter effective cloud fraction filter of 10 % to all species -for SO 2 this correspond to a cloud radiance fraction of 20 % -although this may introduce a clear-sky bias relative to all-sky conditions (Krotkov et al., 2016).This had little impact on SO 2 , although fewer trends in total were detected for NO 2 (232) and HCHO (59); no trends were found for CHOCHO.Generally, most of the highest ranked trends for each species were very similar to the default scenario.This indicates that a 20 % cloud fraction filter is likely an optimum choice for detecting air-quality trends over this region without affecting trend magnitudes.
Despite the use of the static masks and OMI level-1B XTrackQualityFlags in the data gridding (Sect.3.1), statistically significant sampling trends are present for a number of locations where VCD trends occur (Tables 1-4).Therefore, we repeated our analysis, but using only data from OMI's unaffected rows.
However, despite using data from the unaffected detector rows, we still find sampling trends present in some time series.For example, for NO 2 , 227 statistically significant VCD trends were now detected in total, but 57 of those had a sampling trend.As a consequence, for NO 2 only 136 "real" trends in total were detected (instead of 278; see Table 1); in particular, over 100 fewer trends were found over settlements and nearly half the amount over power plants.Nevertheless, trends were found over the same locations (especially in the top 20 ranked urban) with median trends now at 3.0-6.5 % yr −1 (was 3-4 % yr −1 ).The Daura refinery (near Baghdad) and the Umm Qasr port still had the highest trends of these categories (but now of 5.78 ± 0.95 and 7.76 ± 0.91 % yr −1 , respectively).
For SO 2 there were 7 fewer trends detected in total, but notably 9 of the default 18 urban trends were still detected.Median trends were about −7 to 13 % yr −1 (formerly −7 to 10 % yr −1 ).Urban trends now ranged from about −12 to 150 % yr −1 and power plants −7 to 75 % yr −1 .A decrease of about 7 % yr −1 was found over Kharg Island port, as was an increase of 9.69 ± 2.40 % yr −1 over the Daura refinery.
For CHOCHO, only two trends over the towns of Ardebil (of about 25 % yr −1 ) and Shushtar (of about −18 % yr −1 ), both in Iran, were detected.Thus, in summary, we find that using less observations results in fewer trends being detected (since there are less observations used in the averaging, which increases the noise).Furthermore, we find that use of the unaffected rows results in slightly higher median trends overall.Selection of different detector rows also changes the pixel sizes included in the analysis, which may impact the results.
To further reduce noise in the gridded data we increase the degree of smoothing by using a larger Gaussian filter of 0.35 • × 0.35 • with a 2σ width.This had a small effect on our results, as the total number of trends detected for all species and categories was more or less the same in the default scenario, and nearly all the trends locations in the default scenario are still detected.Median trends were also unchanged.
-Test 5: no filtering for outliers in VCD time series analysis.
Whilst filtering for VCD outliers in the time series fitting is advantageous in improving the modelling fitting, it may also remove genuine data points that may affect trend detection.
For NO 2 the number of detections and their locations are approximately the same, albeit with 27 new urban trends replacing 26 default trends.Tehran is now the top-ranked city (in absolute terms), although its relative trend of 8.82 ± 2.59 % yr −1 is still less than that found at Dahuk (15.86 ± 2.33 % yr −1 ).Considering all target categories, we find that median trends range from about −2 to 16 % yr −1 (as opposed to about −3 to 12 % yr −1 ).Similarly for HCHO, the number of detections and their locations are approximately the same, and the range of trends is unchanged.
For the most part SO 2 trends are also in the same places, with a few exceptions over urban and power plant targets.Trends over urban and refinery targets cover a higher range −23 to 148 % yr −1 and 10-20 % yr −1 , respectively, than previously.This results in higher median values of 14 and 15 % yr −1 for these categories (compared to the default case).Median values over oil ports and power plants are unchanged.
Lastly, for CHOCHO, only one urban trend was now detected over Shushtar (Iran) at −8.65 ± 3.83 % yr −1 , and two extra trends were found over the closely located power plants, both at the Az Zour plant complex in Mina Said (Kuwait), which had a trend of −3.64 ± 1.75 % yr −1 .In summary, not filtering for VCD outliers only has a small effect on our results.Typically, most trends in the default scenario are detected, albeit with some changes in trend magnitudes (mostly for NO 2 and SO 2 ).
-Test 6: focus only on cities with > 500 000 people using a spatial mask of ±16 grid cells (∼ 80 km radius around each target).
In this test we focused on the integrated signals from large population centres, which correspond to 32 targets only.For NO 2 there are 22 trends detected (a detection rate of 69 %), for HCHO there were 9 trends (32 % detection rate), and for SO 2 3 trends were found (9 % detection rate).No trends were found for CHOCHO.
We find that NO 2 trends are similar to those found in the default analysis.For example, trends over Tehran, Baghdad, Riyadh, Kirkuk, Orumiyeh (Iran), Isfahan, and Ad-Dammam are approximately the same, although at Irbil the linear trend is now 6 % yr −1 ; previously it was about 9 % yr −1 .
For HCHO there were three locations which could be compared to the default case.These were: Attaif (now 3 % yr −1 was previously 7 % yr −1 ), Abu Dhabi (now 2 % yr −1 was 3 % yr −1 ), and Ad-dammam (still about 3 % yr −1 ), but there were 6 new trend locations found with trends from 2 to 3 % yr −1 , consistent with HCHO trends over other urban settlements.
For SO 2 only Tarbiz has a trend that was previously detected (61 % yr −1 ) but now with a value of 37 % yr −1 .We now find two other trends of about 15 % yr −1 at Kermanshah (Iran) and 22 % yr −1 at al-Mawsil (Mosul, Iraq).
In summary, for NO 2 and HCHO, we find that the trend detection percentages (relative to the 32 targets) increase, but generally the trends over co-existing locations found in the default analysis do not really change.
Therefore, whilst different approaches in the time series analysis may, in some cases, affect trends found over individual locations, overall there are a large number of targets (per gas species) where trends of approximately the same magnitude are consistently detected.This gives some level of confidence in the robustness of the trend analysis.

Comparative studies
Another way to corroborate our results is through comparison to other similar independent studies.For example, Duncan et al. (2016) performed NO 2 trend analyses for the world's major cities using the official NASA product, as described in Bucsela et al. (2013), available from the NASA Goddard Earth Sciences Data Active Archive Center (http://disc.sci.gsfc.nasa.gov).In that study, monthly mean values were based on the average of OMI data falling within 0.3 • × 0.3 • boxes centred over the cities.We find good agreement in trend magnitudes (within cited errors) for Tehran, Baghdad, Kirkuk, Kuwait City, Beirut, Aleppo, Mecca (not classified in our analysis as a real trend), and Jerusalem (not classified here as significant).There is disagreement at Mosul (4.43 ± 1.18 % yr −1 , here 2.14 ± 1.02 % yr −1 ), Riyadh (0.00 ± 1.19 % yr −1 , here 3.27 ± 0.85 % yr −1 ), Homs (−0.81 ± 0.99 % yr −1 , here 1.26 ± 0.74 % yr −1 and classified as real but not significant), and Damascus (−3.72 ± 1.10 % yr −1 , here −0.65 ± 0.92 % yr −1 and classified as not significant or real).These small discrepancies likely occur because of the different choices of OMI data sets and different analysis methodology.
Furthermore, the analysis of SCIAMACHY NO 2 over 2002-2012 by Schneider et al. (2015) revealed relative trends of 3.2 ± 1.1 % yr −1 over Tehran, 2.1 ± 1.1 % yr −1 over Riyadh, and 8.5 ± 1.2 % yr −1 over Baghdad, which are compatible with our NO 2 trends.Similarly, a top-down multispecies inversion involving OMI DOMINO NO 2 columns by Miyazaki et al. (2016) indicated NO x emissions trends of 3.7 % yr −1 for Tehran and 4.7 % yr −1 for Kuwait City, which are broadly consistent with our results.However, Miyazaki et al. (2016) also found an emissions trend of −6.0 % yr −1 for Dubai, whereas we find an NO 2 trend of −0.56 ± 1.05 % yr −1 (classified as real but not significant).
We have also compared our results with the global catalogue of SO 2 emission sources produced by Fioletov et al. (2016).Over the Middle East we were able to identify thirty corresponding targets, for which we compared the linear trends of the derived SO 2 emissions to our derived SO 2 trace gas trends.We find five locations where the SO 2 trends in both studies agree (within errors) and where we class the trends as real and significant.For example, there is excellent agreement at Kharg Island (here −6.89 ± 2.20 % yr −1 versus −7.84 ± 2.50 % yr −1 ) and the Besat power plant in Tehran (in this study 3.54 ± 3.35 % yr −1 , versus 3.86 ± 3.61 % yr −1 ).However, we also find good agreement at an additional 16 sites, where we class the trends as real but not significant, thus indicating the good correspondence between the two independent studies.
Lastly, both Lelieveld et al. (2015) and Duncan et al. (2016) have used independent economic and social information to attempt an interpretation of observed changes in OMI NO 2 and SO 2 over the Middle East.We have further compared our derived trends to the linear growth in the Organization of the Petroleum Exporting Countries' (OPEC) oil production and demand data, and additionally, population, GDP per capita, and energy consumption per capita (EC) data, from the World Data Bank (http://data.worldbank.org).We find such comparisons are, at the very least, difficult to www.atmos-chem-phys.net/17/4687/2017/Atmos.Chem.Phys., 17, 4687-4709, 2017 interpret.For example, over 2005-2014 Iran's GDP grew by 2.8 % yr −1 , inline with the NO 2 and HCHO median linear trends of 2-4 % yr −1 , but its oil production fell by −2.5 % yr −1 , conflicting with the increasing OMI trace gas trends found over refineries.

Summary and Outlook
We have performed a robust and detailed time series analysis on 10 years of OMI trace gas observations of NO 2 , HCHO, SO 2 , and CHOCHO to assess changes in local air quality for over 1000 urban, oil, and energy target locations over the Middle East during the period 2005-2014.We find that the highest average pollution levels of HCHO, SO 2 , and CHOCHO are over the major oil ports and refineries, compared to urban areas and power plants.For HCHO and CHOCHO, the average trace columns are about 15-25 % higher over the oil ports and refineries, whereas for SO 2 the columns are about 60-80 % higher.In contrast, NO 2 is found to be slightly higher over the urban areas and oil ports by about 5-15 %.The highest average pollution levels over urban settlements are typically in Bahrain, Kuwait, Qatar, and the UAE.Other notable pollutant hotspots include: (1) Kuwait City, Tehran, and Mecca; (2) the west coast of Yemen where HCHO levels are high; and (3) elevated SO 2 near Rasanjan and over Kharg Island.We find that the observed vertical columns can exceed average levels by about 40-320 %, depending on the trace gas species.
Our analysis shows that NO 2 linear trends over urban locations range from about −3 to 12 % yr −1 , although only two locations showed a decrease in NO 2 .Linear trends over oil refineries and oil ports are about 2-6 % yr −1 and 2-9 % yr −1 , respectively.Trends over power plants range from 2-8 % yr −1 , with 5 of the topmost 10 trends found in Iran.For HCHO, we find that urban trends are 2-7 % yr −1 .Only six trends of 2-3.5 % yr −1 were detected over oil refineries, and only four trends of 2-4 % yr −1 were detected over oil ports.Trends over power plants ranged from 2-7 % yr −1 .The increasing HCHO trends are mostly found along the western Gulf coast, particularly along the Saudi Arabian coast near Ad-Dammam, and also near Doha in Qatar.Very few SO 2 trends were detected.Over urban areas, trends ranged from about −60 to 120 % yr −1 , with 11 of the 18 trends being detected in Iran.However, some of the high relative SO 2 trends can be attributed to median levels of approximately zero over the location.Over refineries, three trends were detected of about 9-15 % yr −1 , whilst the Iranian oil ports of Bandar-e Khomeini and Kharg Island showed decreases of about 6 % yr −1 .Apart from two locations, we find that CHO-CHO levels are not changing over the Middle East.Derived growth rates are, on average, very similar to the fitted linear trends, although differences can occur for individual locations, particularly for SO 2 .We find that our derived linear trends are generally consistent with other independent OMI trend studies over this region.
Therefore, based on this analysis, we can conclude that for a large number of locations, air quality has deteriorated over 2005-2014.Whilst effective regulatory measures have been established in some countries, it is clear that effective pollutant emission controls are required to limit the health impact on the region's population.Whether this goal is achievable is an open question, especially in countries experiencing civil unrest (e.g.Iraq, Syria, Yemen, and Palestine) or increasing economic growth (e.g.Iran).
In the near future, the TROPOspheric Monitoring Instrument (TROPOMI), which is expected to be launched in 2017 and has a smaller spatial observational footprint compared to OMI, should further help evaluate changes in air quality at local scales.However, it is important that such future measurements are carefully validated and integrated with in situ measurements and chemical transport models, so that they can be utilised properly to influence air-quality policy decision making in Middle Eastern countries.

Figure 1 .
Figure 1.The OMI 2005 annual mean distributions of NO 2 , HCHO, SO 2 , and CHOCHO.The OMI data have been averaged onto a 0.05 • × 0.05 • grid using observations with cloud fractions < 20 % and solar zenith angles ≤ 70 • .Observations affected by the row anomaly are excluded, as described in Sect.3.1, and the gridded data have been smoothed with a 0.15 • × 0.15 • Gaussian filter of 1σ width (2σ for CHOCHO).

Figure 2 .
Figure 2. (Left) The geographical distributions of the target locations showing cities and towns (in red) as specified by the GRUMP v1 settlement point data base (Balk et al., 2006; SEDAC, 2015), oil refineries (in blue) (based on Kootungal, 2010), and oil ports and power plants based on the Global Energy Observatory (GEO) online resource (http://globalenergyobservatory.org).(Right) OMI 2005 NO 2 annual mean of over northern Iran, with a spatial filtering mask applied to extract observations over urban centres (see Sect. 3.2).The OMI data have been averaged onto a 0.05 • × 0.05 • grid using observations with cloud fractions < 20 % and solar zenith angles ≤ 70 • .Observations affected by the row anomaly are excluded, as described in Sect.3.1, and have been smoothed with a 0.15 • × 0.15 • Gaussian filter of 1σ width.
Figures S2 to S4 show similar fits for HCHO, CHOCHO, and SO 2 , respectively.

Figure 3 .
Figure 3.An example of a time series fit to observed NO 2 data overDahuk, Iraq, as outlined in Sect.3.3.(Top)  The monthly OMI NO 2 vertical columns are indicated by dark grey filled circles, whilst the light grey filled circles represent the fitting residual, which has been smoothed with a short-term 200-day filter (dashed blue line) and long-term 667-day filter (red dashed line).The solid red line is the longterm trend F T (t), given by the linear component of the fitted function F (t) (Eq. 1) plus the residual filtered using the long-term filter.The solid blue line is the smoothed fitted curve F S (t) given by F (t) plus the residual filtered using the short-term filter.(Middle) The NO 2 vertical column growth rate in × 10 15 molecules cm −2 yr −1 , which is the derivative of the long-term trend F T (t) shown in the toppanel.(Bottom) The de-trended seasonal cycle F C (t), which is the difference between the long-term trend and the smoothed function fit (i.e.F S (t) − F L (t)).This represents the annual seasonal oscillation with any long-term trend removed.The dark grey filled circles are the fitted harmonic components of F (t).

Figure 4 .
Figure 4. Geographical distribution of locations with statistically significant linear trends in NO 2 (a), HCHO (b), SO 2 (c), and CHOCHO (d) expressed in percent per year, relative to each location's 2005-2014 median vertical column.

Figure 5 .
Figure 5. Geographical distribution of locations with statistically significant linear trends in NO 2 (a), HCHO (b), SO 2 (c), and CHOCHO (c), but here showing the growth rate expressed in percent per year, relative to each location's 2005-2014 median vertical column.

Figure 6 .
Figure 6.Geographical distribution of locations with statistically significant linear trends in NO 2 (a), HCHO (b), SO 2 (c), and CHOCHO (d), but here showing the linear trend (in percent per year) minus the growth rate (in percent per year).

Table 1 .
Statistical summary of the DOMINO NO 2 vertical column density (VCD) analysis conducted at an approximate 10 km radius over the target locations.Statistical values are shown for those locations which have a real VCD trend, i.e. those locations without a corresponding trend in either the air mass factor (AMF), cloud fraction (CFR), cloud-top pressure (CTP), or number of samples (SAM).Values in parentheses correspond to statistics for all locations of a given category.VCD values are given in 10 14 molecules cm −2 , and IQR is the interquartile range.
tively.By comparison, Iraq and Iran were ranked only 109th and 118th.Thus, over this region we find that the more economically developed countries tend to have the highest pollution levels.

Table 2 .
Statistical summary of the SAO HCHO vertical column density (VCD) analysis conducted at an approximate 10 km radius over the target locations.Statistical values are shown for those locations which have a real VCD trend, i.e. those locations without a corresponding trend in either the air mass factor (AMF), cloud fraction (CFR), cloud-top pressure (CTP), or number of samples (SAM).Values in parentheses correspond to statistics for all locations of a given category.VCD values are given in 10 15 molecules cm −2 , and IQR is the interquartile range.

Table 3 .
Statistical summary of the NASA SO 2 vertical column density (VCD) analysis conducted at an approximate 10 km radius over the target locations.Statistical values are shown for those locations which have a real VCD trend, i.e. those locations without a corresponding trend in either the cloud fraction (CFR), cloud-top pressure (CTP), or number of samples (SAM).Values in parentheses correspond to statistics for all locations of a given category.VCD values are given in DU, and IQR is the interquartile range.

Table 4 .
Statistical summary of the BIRA CHOCHO vertical column density (VCD) analysis conducted at an approximate 10 km radius over the target locations.Statistical values are shown for those locations which have a real VCD trend, i.e. those locations without a corresponding trend in either the air mass factor (AMF), cloud fraction (CFR), cloud-top pressure (CTP), or number of samples (SAM).Values in parentheses correspond to statistics for all locations of a given category.VCD values are given in 10 13 molecules cm −2 , and IQR is the interquartile range.Values which cannot be computed are represented by NaN (not a number).
−1 , relative to each location's observed median VCD.Only two locations showed a de-crease in NO 2