Accounting for surface reﬂectance anisotropy in satellite retrievals of tropospheric NO 2

. Surface reﬂectance is a key parameter in satellite trace gas retrievals in the UV/visible range and in particular for the retrieval of nitrogen dioxide (NO 2 ) vertical tropospheric columns (VTCs). Current operational retrievals rely on coarse-resolution reﬂectance data and do not account for the generally anisotropic properties of surface reﬂectance. Here we present a NO 2 VTC retrieval that uses MODIS bidirectional reﬂectance distribution function (BRDF) data at high temporal (8 days) and spatial (1 km × 1 km) resolution in combination with the LIDORT radiative transfer model to account for the dependence of surface reﬂectance on viewing and illumination geometry. The method was applied to two years of NO 2 observations from the Ozone Monitoring Instrument (OMI) over Europe. Due to its wide swath, OMI is particularly sensitive to BRDF effects. Using representative BRDF parameters for various land surfaces, we found that in July (low solar zenith angles) and November (high solar


Introduction
Since the first satellite observations of tropospheric NO 2 from the Global Ozone Monitoring Experiment (GOME)  launched in 1995, the spatial resolution of space-borne UV/VIS instruments has been gradually improved. The pixel size of the OMI sensor  on the Aura satellite launched in 2004 is up to 13 × 24 km 2 at nadir, which is much smaller than the pixel size of earlier instruments such as GOME (40 × 320 km 2 ) and SCIAMACHY (30 × 60 km 2 ) (Bovensmann et al., 1999). The improvement in spatial resolution increasingly allows the sensors to detect NO 2 pollution features on a regional scale, and retrieval algorithms should take full advantage of this capability. For satellite NO 2 retrievals, Published by Copernicus Publications on behalf of the European Geosciences Union. measurement precision and uncertainty depend on a number of factors. A detailed general error analysis was presented by Boersma et al. (2004). It showed that the retrieval errors are dominated by the uncertainty in the tropospheric air mass factor (AMF trop ), estimated to be of the order of 20-50% for polluted-scene pixels.
One of the key input parameters for the calculation of the AMF trop is the surface reflectance. It affects retrievals directly through the clear-sky AMF trop and indirectly through the cloud retrievals. Reflectance of light from the terrestrial surface is generally an anisotropic phenomenon and the angular pattern is controlled by spectral and structural features of the surface cover (Kimes, 1983;Li and Strahler, 1986). Depending on the given viewing and illumination geometry, surfaces may appear brighter or darker. This effect is described mathematically by the bidirectional reflectance distribution function (BRDF) (Nicodemus et al., 1977). The BRDF represents an intrinsic property of the surface and describes the scattering of a parallel beam of incident light to a reflected direction in the hemisphere. Since it is defined as a ratio of infinitesimals, it cannot be measured directly. For BRDF estimation from satellite remote sensing, observations over sufficiently large angular ranges are first atmospherically corrected and then fitted to a semi-empirical BRDF model (Engelsen et al., 1998;Lucht et al., 2000). Multiangular instruments such as the Multiangle Imaging Spec-troRadiometer (MISR) (Diner et al., 1998(Diner et al., , 2005 and the POLarization and Directionality of the Earth's Reflectances (POLDER) (Leroy et al., 1997;Lallart et al., 2008) measure multiple-angle views over a short time span. In contrast, sensors with a single field of view such as the MODerate resolution Imaging Spectroradiometer (MODIS) (Justice et al., 1998) must accumulate sequential observations of the same scene under different viewing geometries over a specified time period. Given the BRDF, several associated reflectance-related quantities can be derived, as described in Schaepman-Strub et al. (2006). In this paper, we use the bidirectional reflectance factor (BRF), defined as the ratio of the radiance reflected by this surface to the radiance reflected by a lossless Lambertian (isotropic) surface under the same irradiance, and the directional-hemispheric reflectance (or black-sky albedo), defined as the integral of the BRDF over all viewing geometries.
In current NO 2 remote sensing retrievals, the assumption of an isotropic (Lambertian) reflecting surface is used. The Lambertian equivalent reflectance (LER) is defined as the reflectance of an isotropic surface, for which the modeled and measured reflectivity at the top of the atmosphere (TOA) are equal, assuming a pure Rayleigh scattering atmosphere without clouds or aerosols in the radiative transfer model (Koelemeijer et al., 2003). LER data sets used in previous operational NO 2 retrievals (Herman and Celarier, 1997;Koelemeijer et al., 2003) were constructed from older satellite instruments with coarser spatial resolutions and mapped onto a grid that is much coarser than the pixel sizes of the more recent instruments. Recently, a new LER data set of Kleipool et al. (2008) with an improved resolution of 0.5 • × 0.5 • and generated from high resolution OMI observations has become available and has been introduced in the operational Dutch OMI NO 2 (DOMINO) product for all data after 17 February 2009 (Boersma et al., 2009a). Apart from their coarse resolution, these LER climatologies also do not account for interannual and short-term variability. In the DOMINO product actual snow and ice are taken into account based on the NISE ice and snow cover data set  built on passive microwave observations. Moreover, assuming a constant reflectance irrespective of viewing geometry is expected to affect the accuracy of the retrieval, especially for instruments with a wide off-nadir viewing range such as OMI (2600 km swath) and GOME-2 (1920 km swath).
In this study we focus on data from the OMI instrument. Figure 1 shows that the Earth's curvature increases the endswath off-nadir viewing zenith angle (VZA) of OMI from 57.5 • at the satellite to 70 • at the surface, which is the relevant angle for calculating BRDF effects. The figure also shows how the at-surface VZA varies across the 60 pixels of an OMI swath.
Our prime motivation for this study is to develop a retrieval with a more accurate treatment of surface reflectance, in order to obtain more reliable NO 2 column estimates. We first present our NO 2 retrieval which is based on the DOMINO product but with considerable improvements (accurate terrain height, look-up table error corrected) as described in Zhou et al. (2009). Here we add the option to consider the angular dependence of surface reflectance. We take advantage of BRDF estimations from MODIS and we use the Linearized Discrete Ordinate Radiative Transfer model (LI-DORT) for accurate BRDF modeling. We then perform sensitivity studies to investigate the effects of surface reflectance anisotropy on the satellite NO 2 retrieval for different viewing geometries and surface BRDF characteristics. OMI vertical tropospheric columns (VTCs) of NO 2 retrieved with our method are compared with results based on MODIS blacksky albedos and TOMS/GOME LER data (Boersma et al., 2004) assumed for isotropic surface reflectance. The comparison is made for the months July and November for a domain covering most of Europe.

MODIS BRDF/albedo algorithm and products
In this study we used the operational MODIS BRDF/albedo algorithm (Lucht et al., 2000) and standard data products (MOD43B, collection 5), and developed a methodology to map this information onto the OMI pixels for the NO 2 trace gas retrieval. The MODIS BRDF/albedo data sets have high spatial resolution (500 m for observations at nadir), a high temporal resolution (retrieved every 8 days based on all clearsky observations over a 16-day interval), and an atmospheric correction that accounts for trace gas absorption, molecular and aerosol scattering, and coupling between atmospheric and surface BRDF (Vermote et al., 2002). The operational MODIS BRDF model characterizes the surface anisotropy with a linear combination of pre-set BRDF kernels (see Eq. 1 below), which are derived from detailed modeling of surface reflectance. All the MODIS BRDF/albedo products are provided with quality flags, and these products have been thoroughly validated against a variety of surface measurements taken at different locations world-wide. The validation studies suggest that the overall accuracy of the MODIS albedo (broadband, integrated from 0.3 to 5 µm) is of the order of 10% with an increasing uncertainty in winter months especially as solar zenith angle increases beyond 70 • -75 • Knobelspiesse et al., 2008;Liang et al., 2002;Liu et al., 2009;Salomon et al., 2006). Very little information is available, however, on the uncertainties of the spectral albedo in different wavelength bands, but it can be larger, especially at the short wavelengths (470 nm) relevant for the present study (Vermote and Kotchenova, 2008). The operational MODIS BRDF/albedo algorithm uses a weighted linear sum of an isotropic parameter plus two BRDF kernels, to characterize the complete surface BRDF: where K vol and K geo are the volumetric and geometric scattering kernels (Roujean et al., 1992), respectively, f iso , f vol and f geo are the isotropic, volumetric and geometric kernel coefficients, and θ , υ, and φ are the solar zenith, viewing zenith and relative azimuth angles, which are defined at the ground (see Fig. 1). By definition, the BRF of a surface is expressed as its BRDF times π (Schaepman-Strub et al., 2006). Note that in the MODIS BRDF/albedo algorithm, the scale factor π is neglected (Lucht et al., 2000), and the BRF can be derived directly from Eq. (1). In practice, the kernel coefficients are determined by an optimization procedure that identifies the best fit of the modeled reflectance from Eq. (1) to a set of atmospheric-corrected reflectance measurements (Lucht et al., 2000). Volumetric scattering is applicable to a horizontally homogeneous leaf canopy. Roujean et al. (1992) derived an expression for kernel K vol (called the RossThick kernel) for a dense leaf canopy. This kernel has a minimum near the backscatter direction and is brighter along the limbs. Geometric scattering, in contrast, expresses effects caused by the larger (inter-crown) gaps in a canopy, as from scenes containing 3-D objects that cast shadows and are mutually obscured from view at off-nadir angles. K geo used in the MODIS data processing is a reciprocal form called LiSparse-R (Lucht et al., 2000) based on the work of Wanner et al. (1995) and Li and Strahler (1992). It is derived from surface scattering and geometric shadow casting theory with an assumption of a sparse ensemble of surface objects casting shadows on the background. It has been shown that this "RossThick-LiSparse-R" model is well suited to describe BRDFs for a wide variety of land covers (Wanner et al., 1995;Lucht et al., 2000;Bicheron and Leroy, 2000).
The MODIS BRDF/albedo standard products MOD43B (collection 5, available from NASA's Land Processes Distributed Active Archive Center, https://lpdaac.usgs.gov/ lpdaac/get data) are produced by combining cloud-free, atmospherically corrected surface reflectance observations (MOD09) from both the Terra and Aqua satellites, and are provided in an Integerized Sinusoidal Grid (ISG) projection with standard tiles representing 1200 × 1200 one kilometer pixels. In this study, we make use of the first three of the four standard products for all land and coastal areas and shallow water regions (within 5 km of land and less than 50 m deep). For each pixel, the first product (MOD43B1) provides the best fit RossThick-LiSparse-R model parameters f k (λ) for the first seven spectral bands (0.47-2.1 µm) of MODIS and three additional broadbands when there are seven or more high-quality observations well distributed over the viewing hemisphere (full inversion). A backup inversion algorithm (Strugnell et al., 2001) is used for cases with insufficient or poor sampling, and for cases where the standard model fitting is of poor quality. In the backup algorithm, a global land cover classification derived from the Olson classification (Olson, 1994) and a seasonal model is employed, and archetypal BRDFs compiled from various field measurements are assigned to each land cover. For each pixel, the corresponding archetypal BRDF is assumed as an a priori guess and its shape is then constrained by the available observations. Jin et al. (2003) and Salomon et al. (2006) demonstrated that the albedo changed only slightly when the MODIS BRDF/albedo algorithm switched from the backup algorithm to the full inversion, indicating that data quality is only little reduced when the backup algorithm has to be used. A fill value is stored if the number of good observations is less than three. The second product (MOD43B2) provides per-pixel quality flags indicating first if the algorithm has produced a result for that pixel, and if so, a quality value for that pixel. The third product (MOD43B3) provides black-sky and white-sky albedos based on coefficients f k (λ) from MOD43B1. As noted already, the black-sky albedo a bs is the ratio of the hemispherically integrated total radiance to a plane parallel incident beam flux, and it is a function only of solar zenith angle (SZA). For a given solar zenith angle θ it can be determined by integrating Eq. (1) over all angles (υ, φ) of the hemisphere. The black-sky albedo of MOD43B3 is computed for the local noon solar zenith angle for each location based on the following polynomial fit, which was found to capture very well the SZA-dependence from the computationally expensive integral when θ is smaller than 80 • (Lucht et al., 2000): −0.007574 − 0.070987 θ 2 + 0.307588 θ 3 + f geo (λ) −1.284909 − 0.166314 θ 2 + 0.041840 θ 3

Surface reflectance and BRDF parameter datasets for OMI NO 2 retrieval
Channel 3 (459 nm-479 nm) MODIS BRDF/albedo products are used here because this channel is closest to the window centered at 440 nm used in our NO 2 retrieval (Zhou et al., 2009). Systematic errors induced by the wavelength inconsistency are expected to be minor. As reported by Kleipool et al. (2008) the albedo differences between 470 nm (center of MODIS channel 3) and 440 nm (NO 2 VTCs) are very small for most of the land types. Obvious differences exist only over water (the average LER decreases from about 0.058 at 440 nm to 0.05 at 470 nm) and bare land (the average LER increases from about 0.135 at 440 nm to 0.15 at 470 nm). The MOD43B products are produced every 8 days based on observations over a 16-day period; this is an appropriate tradeoff between the availability of sufficient angular samples and the temporal stability of surface properties (Wanner et al., 1997). The 8-day MODIS datasets were then applied to all OMI observations from days 5 to 12 within the corresponding 16-day MODIS observation interval. First, a four-step pre-processing of the coefficients f k from MOD43B1 is performed in order to recover missing pixels in the datasets due to poor or insufficient input observations. For each missing pixel, an interpolation between the previous and following 8-day datasets, then an interpolation between the neighboring 5 × 5 pixels, and then an interpolation between corresponding datasets from the previous and subsequent years are attempted in sequence; the pre-processing stops once a value is filled in. For the first and third interpolations (temporal), only snow-free pixels with valid BRDF inversions (MOD43B2 quality flag < 4) are taken. For each such pixel, the coefficient value is taken to be the average of the two temporal neighbours, and the quality flag is assigned the lower quality of the two. A low-pass filtering with a 3 by 3 point kernel is applied to the data set before the second interpolation to avoid the spread of information from potentially noisy adjacent pixels. For the second interpolation (spatial), the pixel is marked as snow-covered if more than half of the available 25 good-quality neighboring pixels are so marked. For a snow-covered (snow-free) pixel, the coefficient is again taken as the average value of its snowcovered (snow-free) neighbours, and the quality flag is set to the worst (highest) value of them. The type of interpolation applied to filled-in values is recorded as a processing flag.
The pre-processed MOD43B products are then mapped onto the OMI pixels. For each OMI pixel, all MODIS pixels (1 km resolution) with centers located inside the OMI pixel are identified. With the geometry parameters (θ , υ, φ) known for each OMI pixel, the BRF and the black-sky albedo a bs are computed following Eqs. (1) and (2), respectively, for each of the identified MODIS pixels. Then the values of BRF and a bs as well as the coefficients f iso , f vol and f geo are averaged over the OMI pixel and stored in HDF5-EOS format together with the original OMI data of the DOMINO product. Note that for the retrieval with full BRDF treatment, we need only the coefficients for the cloud-free part and a bs for the cloudy part of a pixel, and BRF is stored only for the sensitivity studies discussed later. Further parameters mapped onto the OMI pixels are black-sky albedo from MOD43B3 calculated for the local noon SZA, snow-cover, quality and processing flags, and the percentage of valid MODIS pixels within an OMI pixel. Figure 2 is an example of the processing of the coefficient f vol for one OMI orbit on 1 December over central Europe. Panel (a) shows the original MODIS f vol data, panel b the same values after gap-filling, and panel c the coefficients mapped onto the OMI pixels. Data constructed at OMI resolution captures the fine structure in the original data quite well. The original data from MODIS typically has higher noise and more missing values in winter months due to snow and cloud contamination. On the other hand, OMI is less likely to deliver a clear-sky observation at those locations within the corresponding period.

Tropospheric NO 2 retrieval
The Dutch-Finnish OMI instrument is part of the payload of the Earth Observing System (EOS) Aura satellite launched in July 2004. The Aura satellite (Schoeberl et al., 2006) passes over the equator in a sun-synchronous ascending polar orbit at 13:45 LT (local time). In this study, we base our tropospheric NO 2 retrieval on the approach described in Zhou et al. (2009), which uses tropospheric slant columns (SCD trop ) from the Dutch OMI NO 2 (DOMINO) product data (Boersma et al., 2009a, version 1.0.2) available from ESA's TEMIS project (Tropospheric Emission Monitoring Internet Service, www.temis.nl), and calculates tropospheric air mass factors (AMF trop ) with a high-resolution topography data set and a priori vertical profiles from the TM4 chemistry-transport model provided by the DOMINO data product.
In this paper, the AMF trop calculation is supplemented with accurate modeling of surface reflectance anisotropy. Instead of calculating AMF trop with the TOMS/GOME or the Kleipool et al. (2008) albedo data set used in DOMINO (and in the original OMI cloud algorithm), MODIS BRDF parameters as described in Sect. 2.2 are used to characterize the surface BRDF, and this has some important considerations for the retrieval algorithm. Firstly, we must use a radiative transfer model that can deal accurately with bidirectionally reflecting surfaces. Secondly, since cloud parameters (cloud fraction, cloud pressure) from the ancillary cloud preprocessing algorithm depend on the choice of albedo data set, it becomes necessary to retrieve these parameters again.
The AMF trop is defined as the ratio of the SCD trop of the absorber along an average backscattered path of the photons observed by a satellite instrument to the tropospheric vertical column density (VCD trop ). The AMF trop depends on the a priori trace gas profile x a and a set of forward model parametersb which includes cloud parameters, surface albedo and surface pressure. For small optical thickness, the altitude dependence of the measurement sensitivity to the atmospheric species of interest (calculated with a radiative transfer model) can be decoupled from the shape of the vertical trace gas profile (calculated e.g. with an atmospheric chemistry transport model). The AMF trop can then be written as follows (Palmer et al., 2001;Boersma et al., 2004): where "L" is an index denoting the atmospheric layer, m L are the altitude-dependent box air mass factors, and x a,L the layer subcolumns (molecules cm −2 ) of the a priori NO 2 profile. The coefficients c L are layer-specific correction terms that describe the temperature dependence of the NO 2 absorption cross-sections. The AMF for a partly cloudy scene is based on the IPA (independent pixel approximation), and is determined as a linear combination of the AMFs calculated separately for the clear-sky and cloudy fractions of a pixel : where AMF cloud is the AMF for a completely cloudy pixel, and AMF clear the AMF for a completely cloud-free pixel, p c is the cloud pressure, p eff the surface pressure, f cl the OMI effective cloud fraction, and I cl and I cr are the radiances for cloudy and clear scenes, respectively.
We determine cloud fraction f cl and cloud pressure p c from results of the OMI cloud retrieval algorithm based on the O 2 -O 2 absorption band at 477 nm (Acarreta et al., 2004). The slant column density of O 2 -O 2 (N s ) and the measurement-derived continuum reflectance (R c ) are obtained from the OMCLDO2 Level 2 data product available from NASA's mirador earth science data search tool (http://mirador.gsfc.nasa.gov/cgi-bin/ mirador/collectionlist.pl?keyword=omcldo2). The look-up table from the operational OMCLDO2 algorithm is then used to convert the quantities N s and R c into the cloud pressure and the effective cloud fraction (Sneep et al., 2008). For each pixel considered, the pixel-averaged MODIS black-sky albedo a bs and the ground altitude derived from the global digital elevation model GTOPO30 (http://eros.usgs.gov/ #/Find Data/Products and Data Available/gtopo30 info) are used as additional inputs for the look-up table. The backscattered radiance from the clear-sky (I cr ) and cloudy fractions (I cl ) are obtained from the LIDORT model discussed below. The AMF cloud is obtained from Eq. (3), with m L = 0 for all layers below cloud. For consistency with assumptions used in the cloud retrieval algorithm, the cloud is assumed to be a Lambertian surface with albedo 0.8.

BRDF treatment in NO 2 retrieval
The box air mass factors m L are calculated using the Linearized Discrete Ordinate Radiative Transfer model (LI-DORT, version 3.3; Spurr, 2008). This is a multiplescattering model with the capability to generate simultaneous fields of radiances and weighting functions in a multilayer atmosphere. LIDORT can deal with both Lambertian and bidirectionally reflecting surfaces (Spurr, 2004), which makes it especially attractive for our study. Instead of interpolating from a pre-computed lookup-table as in Zhou et al. (2009), we improved the retrieval by running LIDORT on-the-fly for each pixel. Vertical profiles of temperature and NO 2 are taken from the TM4 model (Dentener et al., 2003) but scaled to the high resolution topography as described in Zhou et al. (2009). The box air mass factors are derived from LIDORT radiances and profile weighting functions according to: where the weighting function is defined as the analytic derivative of the intensity field with respect to the optical depth τ of layer L (L = 1. . . 34), and I is the intensity of the backscattered radiance. LIDORT includes a pseudospherical correction for the multiple-scattering contribution which treats the solar beam attenuation in a curved atmosphere, and an exact treatment of the single scattering contribution based on curved-atmosphere attenuation for both the solar and line-of-sight paths. This is important for nadir-geometry satellite instruments with wide-angle offnadir viewings such as OMI and GOME-2. LIDORT 3.3 includes nine possible BRDF kernel functions, and the surface reflectance is specified as a linear combination of (up to) three semi-empirical kernel functions. In our calculations, we selected the kernel functions RossThick and LiSparse that are used in the MODIS BRDF model as described in Sect. 2.1. The nonreciprocal LiSparse kernel in the original LIDORT package (Spurr, 2004) was modified to be consistent with the LiSparse-R kernel used in the MODIS BRDF/albedo algorithm by adding the factor 1/cos θ assumed for the sunlit component (Lucht et al., 2000). For the BRDF surface treatment, the pixel-averaged coefficients (f iso , f vol , f geo ) calculated in Sect. 2.2 are provided as basic inputs for the box AMF calculations.
To solve the radiative transfer equation in an anisotropically scattering medium using the discrete ordinate method, the dependence on azimuth angle is separated using a Fourier series expansion of the radiation field in terms of the cosine of the relative azimuth angle. For each BRDF kernel k, the m-th Fourier component is calculated as: The integration over the azimuth angle from 0 to 2 π is performed by double numerical quadrature over the ranges [0, π] and [−π, 0]. The number of terms of the BRDF azimuth quadrature is set to 50 to assure a numerical accuracy better than 10 −4 (Spurr, 2004). For bidirectionally reflecting surfaces, the reflected radiation field is the sum of the diffuse and direct components for each Fourier term. In LIDORT, the diffuse-field surface contributions are based on components in Eq. (6), while the direct-beam contributions are based on a precise specification of the solar beam BRDF rather than with their truncated forms based on a finite Fourier series expansion. For comparison, box AMFs are also calculated with the Lambertian surface assumption. In this case, we input either the pixel-averaged BRF or the black-sky albedo a bs (calculated in Sect. 2.2) as the Lambertian albedo to be used in LI-DORT. Note that, in contrast with the black-sky albedo, the BRF will account for viewing geometry dependence in the surface reflectance. Despite this, the underlying Lambertian assumption in the radiative transfer model will not account fully for BRDF. Figure 3 shows an example of box AMF profiles calculated with these three different surface treatments. It can be seen that the three profiles are very similar in shape, with decreasing m L towards the surface, illustrating the diminished sensitivity of the satellite instrument at lower levels due to increased scattering of light. The effect of surface treatments is most strongly felt near the surface, where the box AMFs differ by up to 10% in this example. In this case, the lowestlayer box AMF calculated with the full BRDF treatment lies almost half-way between values calculated with BRF and a bs . The differences between the box AMF curves depend on both the BRDF parameters and the geometry parameters, and we study these dependencies in the following section.

Spatial and temporal distributions of BRDF parameters
Before analyzing the impact of BRDF on the NO 2 retrieval, the general BRDF characteristics of the surface over Europe are described in this section and contrasted with the Lambertian assumption in the next section. The BRDF describes the intrinsic reflectance characteristics of the surface which is determined by the specific land surface type and its optical properties. Coefficients of a semi-empirical BRDF model can not be interpreted directly in terms of measurable biophysical variables such as leaf area index, but the latter can be derived from an empirical formula based on the BRDF parameters to distinguish different land cover types or to detect structural changes . To study temporal and spatial distributions of BRDF parameters from MOD43B1 within our domain of interest, we combined six MODIS tiles (horizontal 17-19, vertical 3-4) which together cover a major part of western and central Europe. We calculated monthly mean maps and frequency distributions of the coefficients f iso , f vol and f geo for July and November 2006, as shown in Fig. 4. Snow-covered pixels were excluded since no reliable retrieval of NO 2 can be achieved in these cases due to an incorrect estimation of effective cloud fraction (Boersma et al., 2009a). Spatial variations over land are large and many interesting features can be seen, such as the high values of f iso over Spain due to dry soils, and similarly high values over the Apulia region in southern Italy. Seasonal differences can also be identified. For example, f vol values become smaller and more homogeneous across Europe in November, a phenomenon which corresponds to the autumnal decrease of dense vegetation cover, and f geo values become smaller over northwestern France in November, this time corresponding to the enhanced shadowing of sparse vegetation types in this region. Normalized frequency distributions show the range of values across Europe in the two months. The peak of the f iso frequency distribution shifts from around 0.03 in July to 0.04 in November, with the majority of the data within a range 0.01-0.12 in both cases. For f vol , the majority of values lie in the range 0-0.05 in November but in a wider range of 0-0.12 in July. Aside from a population of zero values in both months that occur mostly over the ocean, the median value of f vol is 0.02 in July and 0.015 in November. The higher July value is due to increased multiple scattering by green-leaf facets. For f geo over land, the peak of the frequency distribution shifts from 0.003 in July to 0.006 in November, with similar ranges 0.001-0.02. The small difference in f geo corresponds to the small seasonal variation in geometric scattering over desert, evergreen needleleaf forest, urban or built-up areas (Bicheron and Leroy, 2000;. Better correlation of BRDF coefficients with the land type and vegetation structure can be found in the nearinfrared band and red band . The anomalous features over the North Sea may be due to the fact that the specular (glitter) BRDF model (Cox and Munk, 1954) characterized by non-linear parameters such as wind speed and refractive index of water is often needed for more accurate modeling of the BRDF over water surfaces. However, this is outside the scope of our study focusing on NO 2 over land.

Geometry dependence of bidirectional reflectance factor
Typical values of BRDF kernel coefficients and SZAs for July and November are summarized in Table 1 (cases A1 and A2). For these cases, BRFs were calculated as a function of viewing zenith angle υ (0 • -70 • ) and relative azimuth angle φ (0 • -360 • ) and presented as polar plots in Fig. 5, where the radius corresponds to υ and the polar angle to φ. In our convention, φ = 0 • corresponds to backward scattering conditions when the observer is on the same side of the local vertical as the sun. Curves of BRFs are also plotted as a function of υ, where the left part with negative viewing angles corresponds to φ = 240 • and the right part corresponds to φ = 60 • in the polar plot; this is representative for an OMI swath. A number of solar zenith angles are considered corresponding to typical values at different latitudes within the domain of interest.
In November, the values of BRF are generally larger and vary more strongly with υ than in July. From the BRF curves for the selected OMI swath in Fig. 5, we can see that the solar zenith angle has an important impact on the sensitivity of BRF to the viewing zenith angle change, especially in November. For example, the BRF difference between υ = 0 • and υ = 70 • increases from about 0.02 at θ = 62 • to 0.04 at θ = 74 • . In Fig. 5b we can see the impact of the "hot-spot" characteristic of the geometric kernel around υ = 30 • in the    Table 1, respectively. Right panels: variation of BRF along an OMI swath for the same cases in July (b) and November (d). The OMI swath is marked by dark blue lines in the left panels.
backward scattering in July, while this disappears at the large solar zenith angles in November (Fig. 5d). In Fig. 6 relative differences between BRF and black-sky albedo are plotted as a function of υ for typical solar zenith angles and three sets of BRDF coefficients for each of the two months. The first set of BRDF coefficients for each month is the same as that used in Fig. 5, and the other two sets represent typical values over northern Poland and northern Italy, referred to as cases B and C in Table 1. These two contrasting areas were chosen since the coefficients f geo and f vol differ significantly, which implies different BRDF characteristics. The SZA differences between these two areas are also considered. The relative difference between BRF and black-sky albedo is a measure of the difference induced by ignoring viewing zenith angle dependence in the reflectance. Trends and values from the same month with different BRDF coefficients are comparable. In November, the reflectance difference increases very fast with υ and can become as large as 50% for the outermost pixels. Note also that the difference is asymmetric with respect to the two different sides of the swath.

Sensitivity of NO 2 VTCs to the surface reflectance treatment
To evaluate the sensitivity of NO 2 VTCs to the surface reflectance treatment, we first calculated a set of box AMFs with the full BRDF treatment, and then generated two more AMF sets based on the Lambertian surface assumption, taking BRF and a bs as the input Lambertian albedo, respectively. The BRDF/BRF comparison reflects the difference induced by radiative transfer modeling without BRDF treatment, while the BRDF/a bs comparison characterizes the difference induced by ignoring the viewing angle dependence. For  Table 1. given a priori profiles, we calculated the clear-sky AMF trop for the above scenarios (AMF BRDF , AMF BRF and AMF bs ) as a function of viewing zenith angle and relative azimuth angle. The a priori profiles used in the study are shown in Fig. 7; these are taken from TM4 model output at OMI overpass times over Germany in summer (11 July) and late autumn (11 November) 2006. These are respectively, typical scenarios for a well-mixed boundary layer in summer and a much more pronounced NO 2 maximum located close to the surface in winter. A recent study by Huijnen et al. (2010) showed that compared to other (regional) air quality models, TM4 a priori partial columns tend to be too large in the boundary layer and to peak at lower levels; this is due to an implementation error for the NO 2 tracer field. Therefore, we modified the lower levels of the winter TM4 a priori profile according to the shape of the EURAD-IM profile in Huijnen et al. (2010), and compared the results with that calculated with the original TM4 a priori profile. In Fig. 8, polar plots of AMFs are shown for three scenarios, referred to as cases A1, A2b and A2a in Table 1. Patterns of AMFs in the polar plots reflect the dependence of air mass factors on viewing geometries given a typical a priori profile for that month. (a) and (c) represent typical scenarios for July and November, respectively, calculated with typical TM4 NO 2 profiles, BRDF coefficients and solar zenith angles of that month. (b) was calculated with the same BRDF coefficients and solar zenith angles as (c), but with the same July NO 2 profile as (a). Comparing (a) and (b) shows the effect of different solar illumination and BRDF parameters in winter as opposed to summer, while comparing (b) and (c) shows the additional effect of different a priori profiles. Most of the differences between summer and winter are explained by the different SZA and BRDF settings. The choice of a priori profile only plays a secondary role by modulating the pattern to some extent and changing its magnitude. For the July case a very similar pattern is obtained when the a priori is replaced by the winter profile (not shown), indicating that the influence of the a priori profile is comparatively small in summer. For July scenario (a), with well-mixed a priori profiles and a constant surface reflectance a bs , it is easy to see that the minimal AMF bs occurs near υ = θ and φ = 0 • , since the shortest average photon path leads to less absorption and scattering, and hence to an AMF minimum. In November (panels b and c) the situation is more complicated and the locations of minimas and maximas vary between the different surface treatments and are additionally modified by the selected a priori profile. This is probably due to the enhanced importance of scattering processes due to the longer photon pathways in November rendering the results more sensitive to BRDF effects and the a priori profile shape. Comparing the three sets of AMFs, we can see that in July AMF BRDF is much closer to AMF BRF and AMF bs than in November, which implies a much smaller relative NO 2 VTCs difference under the assumption of a Lambertian surface (note that since VTC trop = SCD trop /AMF trop the relative differences in NO 2 VTCs are identical, and opposite in sign, to relative differences in AMFs).
To study further the sensitivity of NO 2 VTCs differences to the input parameters, we show the relative difference of NO BRF 2 and NO bs 2 compared to NO BRDF 2 in Fig. 9 for the OMI swath marked in Fig. 5. BRDF coefficients and solar zenith angles for the six scenarios in July and November (Table 1) are considered. In addition, for November all the three a priori profiles in Fig. 7 are used, but we only plot the results with the two TM4 profiles since the results for the EURAD IM and TM4 profiles in Fig. 7b are nearly identical. NO 2 VTCs differences can be as high as 20% in November, but are mostly below 5% in July. Consistent with the result in Fig. 8, the NO 2 VTCs difference is sensitive to both the specific set of BRDF coefficients and the choice of a priori profile. Comparing scenarios B and C with the respective scenarios A shows that different BRDF coefficients and solar zenith angles over different regions in Europe lead to significant variations between the NO BRF 2 difference curves. All curves show a certain degree of asymmetry with respect to the relative azimuth angle. This is caused by the surface anisotropy which makes the surface appear brighter or darker depending on whether the observer is on the same or opposite side of the local vertical as the sun. The results may therefore differ significantly for pixels with similar viewing zenith angles but located on opposite sides of the swath. For an OMI orbit in November, maximum differences of both the NO BRF 2 and NO bs 2 tend to occur for the outermost pixel, on the opposite side of the swath, and it can be seen that the difference of NO BRF 2 is smaller than NO bs 2 for most of the pixels.

Comparison of OMI NO 2 from different surface treatments
We applied our NO 2 retrieval to all OMI observations from 2006 and 2007 and using different surface reflectance treatments as well as different albedo data sets. The following sections will show comparisons of corresponding monthly mean NO 2 VTCs averaged over the two years and mapped onto a 0.05 • × 0.05 • grid. Each grid cell was assigned a weighted mean of all OMI pixels covering the cell. The weighting was done according to OMI pixel size with smaller pixels in the centre of the swath given more weight than the larger pixels at the sides. Only OMI pixels with a cloud radiance fraction smaller than 50% and not contaminated by snow (based on the NISE data set) were considered. For the 2007 data, pixels affected by the row anomalies beginning in June 2007 (see http://www.knmi.nl/omi/research/validation/ cama/badrows.txt) were screened out.

Effect of different surface treatments
In Fig. 10a and b, monthly mean maps for July and November of NO 2 VTCs retrieved with the full BRDF treatment are contrasted with the values obtained with the two different Lambertian surface assumptions (BRF and black-sky albedo). The maps reveal a lot of detail, such as high tropospheric NO 2 columns over densely polluted regions including the Benelux region, the Po Valley, and industrial areas in Germany and Poland, and correspondingly low values over the Alps and other rural areas. Elevated values are also seen along ship tracks over the English Channel and west of Spain. The significant differences of NO 2 VTCs between July and November are mainly due to the increased NO 2 lifetime in winter (Schaub et al., 2007). Figure 10c to f show the relative differences between full-BRDF NO 2 VTCs and those retrieved with the Lambertian surface assumptions (NO BRF 2 and NO bs 2 ). The maps show a smooth spatial variation of the relative monthly mean differences in both seasons. The noisy values at high latitudes in November are due to the very limited number of cloud-and snow-free pixels over these areas. Relative differences are smaller than 12% for most of the domain. Since the differences are a function of geometry parameters as seen in Fig. 9, averaging over all pixels over the same location results in a smaller difference than obtained for individual pixels. Difference maps of NO BRF 2 show larger spatial variation than those of NO bs 2 , which corresponds well with the trend in Fig. 9 showing a larger sensitivity of NO BRF 2 to the differences in BRDF characteristics and solar zenith angles between northern and southern areas. The retrieval with a bs results in an underestimation of NO 2 VTCs over the whole domain in November, which can be explained by the fact that most pixels of the OMI swaths have negative relative differences as seen in Fig. 9b.

Effect of different surface reflectance data sets
Different surface reflectance data sets available today for satellite trace gas retrievals show substantial differences which may result in corresponding differences in retrieved NO 2 VTCs. Figure 11 shows a comparison of monthly mean maps for July and November (averaged over the two years) of MODIS BRF, MODIS black-sky albedo, TOMS/GOME LER and OMI LER (Kleipool et al., 2008). The figures were calculated by first mapping the reflectance data sets onto the OMI pixels and then applying the same data screening and gridding as for NO 2 VTCs in the previous section. The TOMS/GOME LER data set uses the spectral dependence of the GOME database of Koelemeijer et al. (2003), but scales the albedo itself to match the TOMS 340/380 nm database (Herman and Celarier, 1997). We used TOMS/GOME LER to generate another OMI NO 2 VTC data set for comparison with the MODIS-based data (see below) while OMI LER (440 nm) is shown in Fig. 11 only for reference. It can be seen that the differences between the TOMS/GOME LER and the MODIS reflectance data sets are more significant in winter when the snow and cloud contamination is expected to affect the coarse resolution observations  Table 1). much more strongly. The improved resolution of MODIS and OMI results in smoother distributions as compared to TOMS/GOME LER. The OMI LER values are generally closer to the MODIS data sets over the arid areas of Spain and over the plains north and south of the Alps in winter. However, the significantly higher values of OMI LER over the northeastern part of Europe point to some problems with residual snow and/or cloud contamination in this data set. NO 2 VTCs retrieved with the TOMS/GOME LER (NO LER 2 ) are shown in Fig. 12 and compared with the BRDF surface treatment. The maps of NO LER 2 show a general pattern very similar to that of NO BRDF 2 in Fig. 10. However, the relative differences ( Fig. 12c and d) are significant and generally larger in November than in July. In November, the contrast between polluted and rural areas is smaller than  Table 1 (A1, B1 and C1 are summer cases; A2a, A2b, B2 and C2 are winter cases). that seen in Fig. 10b, and some hot spots (such as the Swiss plateau) are absent. The mean relative differences of NO LER 2 are much larger than those of NO BRF 2 and NO bs 2 , especially in November, suggesting that the differences between the TOMS/GOME LER and the generally lower values of MODIS black-sky albedo have a more profound impact on the retrieved NO 2 VTCs than the effects of surface anisotropy. The patchy structure of the relative NO 2 VTCs differences ( Fig. 12c and d) is due to sharp transitions between adjacent grid cells in the TOMS/GOME LER data set provided on a grid of 1 • × 1 • resolution.
In July, the relative NO 2 VTCs differences are mostly below 15% over land. However in November, NO LER 2 is lower than NO BRDF 2 by 20%-60%. This is most likely due to snow and cloud contamination which is expected to affect the TOMS/GOME LER data more strongly than the MODIS data due to the coarse spatial resolution of the GOME sensor. Snow or cloud contamination leads to a high bias in surface reflectance (compare Fig. 10f with Fig. 10b) and therefore an overestimation of air mass factors. Furthermore, NO 2 VTCs over polluted areas (where the NO 2 a priori profile has a stronger peak at low altitudes) are more sensitive to variations in surface albedo. The comparatively high NO LER 2 values over the arid areas of Spain are consistent with too low TOMS/GOME LER data values compared with MODISderived reflectances.
Comparisons of the cloud fractions and cloud pressures derived from MODIS black-sky albedo and the TOMS/GOME LER show that the change of albedo affects the cloud fraction more strongly (correlations in f cl between the two data sets are only 0.59 and 0.73 in July and November, respectively) than the cloud pressure (correlations in p c are 0.96 and 0.92 for July and November, respectively). The generally lower values of a bs lead to higher mean values of f cl and p c (f cl : 4.87% and 8.22%; p c : 859.8 hPa and 881.4 hPa for July and November, respectively) compared to those from the TOMS/GOME LER (f cl : 3.73% and 6.39%; p c : 855 hPa and 864.3 hPa for July and November, respectively). Note that these numbers are based on pixels that have been pre-selected for low cloud radiance fractions. Given the significant differences between albedo data sets we think it is mandatory to re-run the OMI cloud algorithm with any new albedo database as is done in our study.

Comparison of monthly mean NO 2 VTCs from different parts of the swath
For OMI, viewing geometry varies considerably across the swath, but remains relatively constant for the same pixels in subsequent swaths. Hence, we expect that different parts of the swath are affected differently by BRDF effects (cf. Fig. 9). If BRDF effects are ignored in the retrieval, however, NO 2 VTCs obtained from pixels near the left-hand limit of the swath may differ systematically from values obtained at the centre or the right-hand limit. To test this hypothesis, we binned NO 2 VTCs according to their location within an OMI swath and computed monthly mean fields for three different bins (left, centre, right) separately. From the sixty pixels of each OMI swath, we selected eight pixels for each of the three bins (pixels 3-11 for the left, 26-34 for the centre, and 52-58 for the right bin), discarding the two outermost pixels on each side. Figure 13 shows NO BRDF 2 and NO bs 2 results for the three bins in November, averaged over 2006 and 2007 (cloud radiance fraction < 50%), as well as the relative differences for NO bs 2 . The patterns of NO 2 VTCs have more similarity in the maps of the left and rightmost pixels while the values of the center pixels tend to be significantly lower. As suggested by the relative difference maps (bottom panel), the use of BRDF surfaces clearly has a positive effect, bringing the July November ) in mean values when the surface is treated as Lambertian, assuming the BRF as albedo for July (c) and November (d), and the black-sky albedo (as Lambertian input) for July (e) and November (f). three bins into closer agreement with each other. However, it does not fully correct for the significant differences between the centre and edge bins. The unambiguous identification of BRDF effects is complicated by several factors. First, the same location at the surface is seen at different times of the day (approx. 2 h difference in local time between left and right limiting bins). The diurnal cycle of NO 2 emissions and photochemistry may therefore contribute to the differences of NO 2 VTCs across the swath, though in November, the diurnal cycle in NO 2 VTCs is not very strong as demonstrated in Boersma et al. (2009b). Secondly, bins cover different days in November 2006 and 2007, for which meteorological conditions may not be equal. Another, probably dominant, factor is the tendency for cloud fractions to be smaller for the center pixels (mean cloud radiance fraction 20%); larger pixels at swath edges are less likely to be cloud free (mean cloud radiance fraction 30%) (Krijger et al., 2007). Since it is not possible (yet) to model complex cloud-related effects in the retrieval algorithm, the tropospheric NO 2 retrieval for the cloud-contaminated pixels has a higher uncertainty than that for clear-sky scenes. High NO 2 VTCs for the side pixels over some areas in the northern part of Europe should therefore be treated with caution. In general, the pattern in the relative differences shown in the bottom row of Fig. 13, with an overestimation of NO 2 VTCs for the left-hand pixels and a marked underestimation in the center, correlates well with the corresponding viewing angle ranges (45 • , 63 • for the left and −8 • , +8 • for the center pixels) in Fig. 9b. The center and leftmost pixels exhibit the largest differences, which can be higher than 15% in absolute value over some areas such as Poland and northern Germany, further confirming the results of Fig. 9b. Although not fully conclusive, the analysis here demonstrates the potential benefits of accounting for surface BRDF effects in the retrieval. As also demonstrated by this analysis, a quantitative proof is difficult and will require an extensive statistical analysis applied to multiple years of observations; we aim to address this issue in a follow-up paper.

Conclusions and outlook
A new satellite tropospheric NO 2 retrieval accounting for the dependence of surface reflectance on the illumination and viewing geometry was presented and applied to two years of OMI observations over the major part of western and central Europe. We developed a methodology which, for each OMI satellite pixel, calculates pixel-averaged BRDF parameters based on high temporal and spatial resolution BRDF data from the MODIS instrument. These parameters were then used as input for the air mass factor calculations with the radiative transfer code LIDORT. In this way we fully account for surface BRDF effects and the surface-atmosphere coupling due to multiple scattering and reflection. Cloud parameters (cloud fraction, cloud pressure) were recalculated for each pixel based on the OMI cloud retrieval algorithm using the MODIS black-sky albedo for the surface reflectance in order to be consistent with the NO 2 retrieval.
We studied the spatial and temporal variation of the isotropic (f iso ), volumetric (f vol ) and geometric (f geo ) BRDF coefficients, and for the corresponding BRDF contributions with representative solar zenith angles for July and November 2006, we studied the BRDF dependence on geometry parameters. An accurate surface treatment of BRDF effect LEFT CENTER RIGHT Fig. 13. Binning of mean NO 2 VTCs in November averaged over 2006 and 2007 retrieved with full BRDF surface treatment (first row) and black-sky albedo (second row), and mean relative differences of NO bs 2 (third row) in regard to the location of pixels in OMI swaths. Results in the left panels are from pixels located at the 3rd to 10th left-most positions in the swath, those in the middle panels from pixels at the eight center positions, and those in the right panels from the pixels located at the 3rd to 10th right-most positions.
was found to be more important in winter, when variations in BRDF with land type and latitude-dependent SZA across Europe can strongly affect the BRDF characteristics. To evaluate the effect of a full BRDF treatment versus the traditional Lambertian surface approximation on the NO 2 retrieval, we compared the NO 2 VTCs of the new approach with two sets of results using the BRF and black-sky albedos as Lambertian inputs. With high SZA and enhanced NO 2 profile loading in winter, the polar plots of NO 2 VTCs exhibit a more complicated pattern. NO 2 VTCs are more sensitive to surface reflectance treatment in November than in July; retrieval differences between NO BRF 2 or NO bs 2 and NO BRDF 2 for an OMI swath can be up to about 20% (15%) for the outermost (inner half of) pixels, and are sensitive to specific choices of BRDF coefficients, SZA values and a priori profile.
To analyze the influence of the new treatment of surface anisotropic reflectance on the OMI NO 2 retrieval, we studied not only the mean NO 2 VTCs in July and November averaged over all clear-sky pixels, but also the binned NO 2 VTCs according to the location of pixels within OMI swaths. Patterns in these NO 2 VTCs correspond closely with trends seen in the sensitivity study above, and this demonstrates that the accurate treatment of surface anisotropic reflectance is especially important when individual pixels are analyzed, since the NO 2 VTCs difference with a Lambertian surface assumption depends strongly on geometry parameters and BRDF characteristics.
Furthermore, we demonstrated the potential improvement of our MODIS BRDF-based retrieval over other available retrievals based on the TOMS/GOME LER data set and blacksky albedos. Benefiting from the higher spatial and temporal resolution, the contrast between the polluted and clean areas is enhanced with the BRDF-based results. Moreover, with a more accurate BRDF-based calculation of AMFs, retrieved NO 2 VTCs for the same location and time period tend to agree better between the different subsections of the swath.
In our future work, the quality of the tropospheric NO 2 columns will be assessed by comparison with ground-based NO 2 measurements and the method will be applied to several years of OMI observations over Europe to study the temporal and spatial variations of the NO 2 columns. This study also suggests the special need in further studies on accounting for surface anisotropic reflectance effect in tropospheric NO 2 retrieval, in particular for satellites with wide swaths (e.g. GOME-2), and future geostationary instruments, for which changing solar zenith angles during the measurements will contribute to surface anisotropic reflectance effects.