Accounting for the effects of surface BRDF on satellite cloud and trace-gas retrievals : a new approach based on geometry-dependent Lambertian equivalent reflectivity applied to OMI algorithms

Most satellite nadir ultraviolet and visible cloud, aerosol, and trace-gas algorithms make use of climatological surface reflectivity databases. For example, cloud and NO2 retrievals for the Ozone Monitoring Instrument (OMI) use monthly gridded surface reflectivity climatologies that do not depend upon the observation geometry. In reality, reflection of incoming direct and diffuse solar light from land or ocean surfaces is sensitive to the sun–sensor geometry. This dependence is described by the bidirectional reflectance distribution function (BRDF). To account for the BRDF, we propose to use a new concept of geometry-dependent Lambertian equivalent reflectivity (LER). Implementation within the existing OMI cloud and NO2 retrieval infrastructure requires changes only to the input surface reflectivity database. The geometry-dependent LER is calculated using a vector radiative transfer model with high spatial resolution BRDF information from the Moderate Resolution Imaging Spectroradiometer (MODIS) over land and the Cox–Munk slope distribution over ocean with a contribution from water-leaving radiance. We compare the geometry-dependent and climatological LERs for two wavelengths, 354 and 466 nm, that are used in OMI cloud algorithms to derive cloud fractions. A detailed comparison of the cloud fractions and pressures derived with climatological and geometry-dependent LERs is carried out. Geometry-dependent LER and corresponding retrieved cloud products are then used as inputs to our OMI NO2 algorithm. We find that replacing the climatological OMI-based LERs with geometry-dependent LERs can increase NO2 vertical columns by up to 50 % in highly polluted areas; the differences include both BRDF effects and biases between the MODIS and OMI-based surface reflectance data sets. Only minor changes to NO2 columns (within 5 %) are found over unpolluted and overcast areas.


Introduction
Satellite ultraviolet and visible (UV-vis) nadir backscattered sunlight trace-gas, aerosol, and cloud retrieval algorithms must accurately estimate the reflection by the Earth's surface in order to produce high-quality data sets.Surface reflectivity climatologies used in most current algorithms are typically gridded monthly Lambertian equivalent reflectivities (LERs) that have been derived from satellite observations (e.g., Herman and Celarier, 1997;Kleipool et al., 2008;Russell et al., 2011;Popp et al., 2011).These climatologies generally have no dependence on the observation geometry.However, it is well known that both ocean and land reflectivity depend upon viewing and illumination geometry.
Here, we give some basic definitions that have been used in the literature to provide context to our problem and for clar-Published by Copernicus Publications on behalf of the European Geosciences Union.ity because sometimes different definitions have been used for similar or the same quantities.This dependence is described by the bidirectional reflectance distribution function (BRDF), mathematically expressed as where dI (ω r ) is the portion of total radiance reflected in the direction defined by the vector ω r due to the illuminating irradiance, dF , from the direction defined by the vector ω i : dF (ω i ) = I (ω i )cos(θ i )dω i .I (ω i ) is the radiance incident on the surface from the direction of ω i , θ i is the angle between the normal to the surface and the direction of illuminating light, and dω i is the element of the solid angle (Nicodemus, 1965;Schaepman-Strub et al., 2006;Martonchik et al., 2000).The reflected radiance I (ω r ) is calculated by integrating the product of BRDF and dF over all directions of the incident radiation.I (ω r ) provides a boundary condition at the surface for computations of the top-of-atmosphere radiance.When the surface is illuminated by a parallel beam of light, the integral over the solid angle of reflected light is where θ r is the zenith angle of reflected light (subscript "r" is omitted in the remainder of the paper for simplicity) provides the so-called black sky albedo (BSA) of the surface.Integration in Eq. ( 2) is carried out over the solid angle of 2π for the upper hemisphere.This equation is general, but it follows from Eq. (2) that the BRDF for a special case of a perfect Lambertian surface is equal to 1/π .The frequently used dimensionless bidirectional reflectance factor (BRF) is defined as "the ratio of the radiant flux reflected by a sample surface to the radiant flux reflected into the identical beam geometry by an ideal (lossless) and diffuse (Lambertian) standard surface, irradiated under the same conditions as the sample surface" (Schaepman-Strub et al., 2006).In general, the relationship between BRF and BRDF for an arbitrary surface can be obtained from Eq. (1) by using BRDF = 1/π for an ideal Lambertian surface, i.e., BRF(ω i , ω r ) = dI (ω i , ω r )/dI Lam (ω i ) = πBRDF(ω i , ω r ).(3) BRF and BRDF are both inherent properties of the surface that do not depend on the illumination conditions (Schaepman-Strub et al., 2006).While BRDF is a function describing a surface for all possible illuminating and reflected directions, the BRF refers to a specific illumination and observational geometry for a given measurement.BRF from satellite observations can therefore differ significantly for the same area over different days due to variations in sun-satellite geometries.In other words, for a given surface BRDF is always the same (neglecting seasonal changes), but BRF changes from day to day depending on observational conditions.
Many satellite UV-vis algorithms are based on the socalled mixed Lambert equivalent reflectivity (MLER) model, first introduced by Seftor et al. (1994).For example, the MLER concept is currently used in most trace-gas (Boersma et al., 2011;Bucsela et al., 2013) and cloud (Acarreta et al., 2004;Joiner and Vasilkov, 2006) retrieval algorithms for the Ozone Monitoring Instrument (OMI), a Dutch-Finnish UVvis sensor (Levelt et al., 2006) onboard the NASA Aura satellite.The MLER model treats cloud and ground as horizontally homogeneous Lambertian surfaces and mixes them using the independent pixel approximation (IPA).According to the IPA, the measured top-of-atmosphere (TOA) radiance is a sum of the clear sky and overcast subpixel radiances that are weighted with an effective cloud fraction (ECF).The ECF is calculated by inverting (4) at a wavelength not substantially affected by rotational Raman scattering (RRS) or atmospheric absorption, where I m is the measured TOA radiance, I g and I c are the precomputed clear sky (ground) and overcast (cloudy) subpixel TOA radiances, and R g and R c are the corresponding ground and cloud LERs, respectively.The MLER model typically assumes R c = 0.8.This value of R c was used by McPeters et al. (1996) for a UV total column O 3 algorithm and independently derived by Koelemeijer et al. (2001) for use in near-infrared O 2 A-band cloud pressure retrievals.The assumption of R c = 0.8 effectively accounts for Rayleigh scattering in partially cloudy scenes (Ahmad et al., 2004;Stammes et al., 2008).This approach also accounts for scattering/absorption that occurs below a thin cloud.In this paper we also assume R c = 0.8 for the OMI cloud and NO 2 algorithms.
The MLER model compensates for photon transport within a cloud by placing the Lambertian surface somewhere in the middle of the cloud instead of at the top (Vasilkov et al., 2008).As clouds are vertically inhomogeneous, the pressure of this surface corresponds not necessarily to the geometrical center of the cloud but rather to the so-called optical centroid pressure (OCP) (Vasilkov et al., 2008;Sneep et al., 2008;Joiner et al., 2012).The cloud OCP can be thought of and modeled as a reflectance-averaged pressure level reached by back-scattered photons (Joiner et al., 2012).Cloud OCPs are the appropriate quantity for use in trace-gas retrievals from satellite instruments (Vasilkov et al., 2004;Joiner et al., 2006Joiner et al., , 2009)).
In one of the early studies to explore the effects of surface BRDF on satellite trace-gas retrievals, Zhou et al. (2010) showed how various treatments of surface reflectance, including BRDF, affect OMI tropospheric NO 2 retrievals over Europe.Their study, which covered the months of July and November, suggested that accounting for surface BRDF effects can change NO 2 retrievals by up to 20 % with the largest effects at high view angles.Ignoring the surface BRDF can also introduce NO 2 retrieval errors that vary with land type (Noguchi et al., 2014;Kuhlmann et al., 2015).Russell et al. (2011) studied the effect of using different surface albedo products on the NO 2 columns and found that the impact of the surface albedo can be up to ±40 % for land.In an effort to improve NO 2 retrievals over China, Lin et al. (2014Lin et al. ( , 2015) ) revised the calculation of tropospheric air mass factor (AMF) in the Dutch OMI NO 2 (DOMINO) product using improved information for cloud, aerosols, and BRDF from the Moderate Resolution Imaging Spectroradiometer (MODIS); they reported better agreement with independent NO 2 observations.Similarly, McLinden et al. (2014) improved OMI NO 2 standard product (Bucsela et al., 2013) for the Canadian oil sands region using high-resolution MODIS retrievals.Our motivations for this work follow from these studies that offered valuable insights into the effects of the surface BRDF on NO 2 retrievals.We continue in this line of investigation by (1) examining in detail the BRDF effect on retrieved cloud parameters that are important inputs for trace-gas retrievals including NO 2 ; (2) additionally investigating BRDF impact on cloud and NO 2 retrievals over ocean; and (3) providing a computationally efficient method of accounting for BRDF effects over both land and water that can be incorporated into existing retrieval algorithms with minimal changes.
To account for surface BRDF in the existing MLER cloud and trace-gas algorithms, we introduce the concept of a geometry-dependent surface LER.The geometry-dependent LER is derived from TOA radiance computed with Rayleigh scattering and BRDF for the exact geometry of a satellitebased pixel.This approach does not require any major changes to existing MLER trace-gas and cloud algorithms.The main revision to the algorithms requires replacement of the existing static LER climatologies with LERs calculated for specific field-of-view (FOV) sun-satellite geometries.The geometry-dependent surface LER approach can be applied to any current and future satellite algorithms that use the MLER concept.
The main goal of this paper is to document a new global surface reflectivity product that will be publicly available and could be easily used within several existing operational satellite trace-gas and cloud algorithms.We implement the geometry-dependent LERs based on a MODIS BRDF product and use these LERs within OMI cloud and NO 2 algorithms.Henceforth, when we refer to geometry-dependent LERs, this refers to a MODIS-based data set.We compare the cloud and NO 2 retrievals based on the geometrydependent LER with the retrievals based on the climatological LER derived from TOMS and OMI measurements.Henceforth, climatological LERs refer to products derived from OMI and TOMS.The differences between those retrievals include both BRDF effects and possible biases between the MODIS and other instrument (OMI and TOMS) reflectance data sets.The existing operational algorithms make use of climatological LER products.By comparing the products retrieved with the geometry-dependent LER with those retrieved with the climatological LER, we address a practical question of how large the differences in various satellite products would be if the climatological LERs were replaced with the geometry-dependent LERs.
It should be noted that the MODIS BRDF product is derived from the atmospherically corrected TOA reflectances (i.e., aerosol and Rayleigh scattering effects are removed at the high spatial resolution of MODIS).In contrast, the climatological LERs currently used in OMI algorithms, from either the Total Ozone Mapping Spectrometer (TOMS) or OMI, are derived by correcting only for Rayleigh scattering and thus include aerosol effects (see details in Herman and Celarier, 1997;Kleipool et al., 2008).Therefore, the use of the geometry-dependent LER product in trace-gas algorithms over heavily polluted regions may also require an explicit account of aerosols (Lin et al., 2015).In this study we do not consider aerosol effects.For all radiative transfer calculations, we use the VLIDORT code (Spurr, 2006).VLIDORT computes the Stokes vector in a plane-parallel atmosphere with a non-Lambertian underlying surface.It has the ability to deal with attenuation of solar and line-of-sight paths in a spherical atmosphere, which is important for large solar zenith angles (SZA) and viewing zenith angles (VZA).We account for polarization at the ocean surface using a full Fresnel reflection matrix as suggested by Mishchenko and Travis (1997).Unlike Lin et al. (2014Lin et al. ( , 2015)), we use a vector code because neglecting polarization can lead to considerable errors for modeling backscatter spectra in UV-vis.This is particularly the case for modeling backscatter spectra over the ocean where reflection of unpolarized light from the flat ocean surface at the Brewster angle leads to perfect linear polarization (Vasilkov et al., 1990a, b).
Both the RRS and O 2 -O 2 algorithms utilize the MLER concept.We use 354 and 466 nm in the RRS and O 2 -O 2 algorithms, respectively, to compute ECF.It should be noted that the ECF implicitly accounts for non-absorbing aerosols, treating them as clouds and this increases cloud fraction.However, the increase of cloud fraction due to the presence of aerosols cannot correctly reproduce an increase of diffuse solar light at the surface caused by aerosol scattering.This may introduce some error in the calculation of the clear-sky subpixel radiance because the BRDF effect depends on a ratio of diffuse to direct solar light.
The OMI RRS cloud algorithm is detailed in Joiner et al. (2004), Joiner andVasilkov (2006), andVasilkov et al. (2008).OCP is derived from the high-frequency structure in the TOA reflectance caused by RRS in the atmosphere.The OCP is retrieved by a minimum-variance technique that spectrally fits the observed TOA reflectance within the spectral window of 345.5-354.5 nm.The RRS algorithm does not report the cloud OCP for ECF < 0.05 due to large retrieval errors at small values of ECF (Vasilkov et al., 2008).
Our O 2 -O 2 cloud algorithm retrieves OCP from OMIderived oxygen dimer slant column densities (SCDs) at 477 nm.Our algorithm spectral fitting differs from KNMI's in that it utilizes temperature-dependent O 2 -O 2 cross sections (Thalman and Volkamer, 2013) and incorporates a new fitting technique similar to that developed by Marchenko et al. (2015) for NO 2 SCD retrieval.The fitting procedure derives the O 2 -O 2 SCD using retrieved O 3 and NO 2 slant column estimates from independent OMI algorithms.This is an implementation choice that is designed to minimize potential errors due to cross talk between O 3 , NO 2 , and O 2 -O 2 cross sections during the fitting procedure.
The OCP is estimated using the MLER method to compute the appropriate AMFs (Yang et al., 2015).To solve for OCP, we invert where VCD is the vertical column density of O 2 -O 2 (VCD = SCD/AMF), AMF g and AMF c are the precomputed (at 477 nm) clear sky (ground) and overcast (cloudy) subpixel AMFs, P s is the surface pressure, and f r is the cloud radiance fraction (CRF) given by f r = ECF • I c /I m .Lookup tables of the TOA radiances and AMFs were generated using VLIDORT.Temperature profiles needed for computation of VCD and AMF are taken from the Global Modeling Initiative (GMI) chemistry transport model (Strahan et al., 2007) driven by the NASA GEOS-5 global data assimilation system (Rienecker et al., 2011).Comparisons of the retrieved OCPs with those from the operational KNMI OMI O 2 -O 2 algorithm, OMCLDO2, have shown good agreement with a correlation coefficient of ∼ 0.99 for ECF > 0.2 when identical surface climatological LERs are used.

Basic approach
In this section we describe our approach of generating the geometry-dependent LER.First, we average all input data over a nominal OMI pixel.The input data include MODISderived land BRDF kernel coefficients, land-water flags, terrain heights from a digital elevation model, and chlorophyll values and wind speed over water surfaces.Second, we compute the TOA radiance accounting for surface BRDF.Third, we calculate the geometry-dependent LER from the TOA radiance.Then we use this geometry-dependent LER in cloud and NO 2 algorithms to replace the climatological LERs.
The BRDF over land is calculated using the Ross-Thick Li-Sparse (RTLS) kernel model (Lucht et al., 2000): where the coefficients a iso , a vol , and a geo come from MODIS data, the isotropic kernel, a iso , describes the Lambertian part of light reflection from the surface, the volumetric kernel, k vol , describes light reflection from a dense leaf canopy, and the geometric kernel, k geo , describes light reflection from a sparse ensemble of surface objects casting shadows on the background assumed to be Lambertian.The kernels are the only angle-dependent functions, the expressions of which are given in Lucht et al. (2000).The BRDF coefficients are spatially averaged over an actual satellite FOV and used to calculate TOA radiance for its observation geometry.
The BRDF coefficients depend on wavelength.For the present study we selected two wavelengths in the UV and vis: 354 and 466 nm.These wavelengths are relatively free of atmospheric RRS and trace-gas absorption.The BRDF coefficients at 466 nm are directly taken from the MCD43GF product at 470 nm (MODIS Band 3 has a wavelength range from 459 to 479 nm and a center wavelength of 470 nm) that is provided at a spatial resolution of 30 arcsec (Schaaf et al., 2002(Schaaf et al., , 2011)).Because the MODIS product is not available at 354 nm, we plan to adjust the 470 nm LERs to account for potential spectral dependences.The adjustment applies the spectral ratio of climatological OMI-derived LERs, R g (354)/R g (470), similar to the approach of McLinden et al. (2014).In the paper we assume that the BRDF coefficients are spectrally independent to focus on the surface BRDF effects only.Using climatological data of Kleipool et al. (2008) we find that this assumption can be valid for some areas; for example, the climatological ratio R g (354)/R g (470) is close to unity (within ±5 %) over the eastern part of North America.However, this is not the case for arid and semiarid areas.We plan to release our geometry-dependent LER product computed for wavelengths other than 470 nm using a spectral correction of the BRDF coefficients.This spectral correction will be based on the ratio R g (354)/R g (470) derived from a critical analysis of different existing data sets of climatological satellite-derived LERs.
To calculate TOA radiance over water surfaces, we account for both light specularly reflected from a rough water surface and diffuse light backscattered by water bulk and transmitted through the water surface.We neglect contributions from oceanic foam that can be significant for high wind speeds.Reflection from the water surface is described by the Cox-Munk slope distribution function (Cox and Munk, 1954).We use an isotropic form of the Cox-Munk distribution in which the facet-slope variance is independent of wind direction.All computations use a wind speed of 5 m s −1 which is close to the climatological mean.
Diffuse light from the ocean is described by a Case 1 water model that has chlorophyll concentration as a single input parameter (Morel, 1988).Our Case 1 water model accounts for the anisotropic nature of light backscattered by the ocean (Morel and Gentili, 1996).A spatial distribution of chlorophyll concentration is taken from the monthly Sea-WiFS climatology.The common Case 1 water model developed for the vis (Morel, 1988) was extended to the UV using data from Vasilkov et al. (2002Vasilkov et al. ( , 2005)).To calculate waterleaving radiance, we need to know the downwelling atmospheric transmittance at the surface.The transmittance is obtained by calculating the total atmospheric direct and diffuse downwelling flux at the surface.The diffuse contribution in the transmittance will itself depend on the water-leaving radiance.To calculate the atmospheric transmittance, we introduce in VLIDORT a module for the iterative calculation of the transmittance, in which the first computation is made for a black surface, and this is then used again as input to the water-leaving contribution.This process is repeated until convergence of the transmittance is achieved (two or three iterations are sufficient).
To estimate LER over mixed surface types, we compute an area-weighted radiance for uniform land and water contributions within an OMI FOV.The LER for heterogeneous surface pixels is then calculated from this linear combination of radiances.The high spatial resolution MCD43GF product (Schaaf et al., 2011) also supplies an eight-category land water classification map at the same resolution as the BRDF parameters.We convert this map into a binary land-water mask by merging all shorelines and ephemeral water into the land category and classifying all other water subcategories simply as water.We then compute the areal fraction of land and water for each OMI FOV.For specification of the OMI pixel, we used the OMPIXCOR product that provides coordinates of OMI pixel corners (http://disc.sci.gsfc.nasa.gov/uui/datasets/OMPIXCOR_V003/summary).We used an option of overlapping pixels in the along-track direction corresponding to 75 % energy in the along-track FOV.In this option the edges of the FOV are aligned in the cross track direction but overlap in the along-track direction.
Given the computed TOA radiance, I comp , the LER is calculated by inverting where λ is wavelength, θ is the VZA, θ 0 is the SZA, φ is the relative azimuth angle, R is the LER, I 0 is the TOA radiance calculated for a black surface, T is the total (direct + diffuse) solar irradiance reaching the surface converted to the ideal Lambertian reflected radiance (by dividing by π ) and then multiplied by the transmittance of the reflected radiation between the surface and TOA in the direction of a satellite instrument, and S b is the diffuse flux reflectivity of the atmosphere for the case of its isotropic illumination from below (Eq.200 in Chandrasekhar, 1960;Dave, 1978).To speed up computations, we created lookup tables of the quantities I 0 , T , and S b for selected wavelengths.Averaging the BRDF coefficients over an OMI pixel may not be equivalent to averaging the high-resolution surface LER over the OMI pixel.We carried out a numerical experiment of calculations of TOA radiances using the highresolution BRDF coefficients and OMI geometries for the Washington-Baltimore area of the United States (Fig. 1).The TOA radiances were converted into LERs using Eq. ( 7) and then the LERs were averaged over OMI pixels.The resulting LERs were compared with that calculated from the standard procedure of averaging the BRDF coefficients first.We found that the mean LER difference was equal to 0.75 × 10 −5 with a standard deviation of 4.2 × 10 −4 , which is quite acceptable for our purposes.
It should be noted that aerosols are not included in the computation of the geometry-dependent LER.Scattering by aerosols in the atmosphere reduces the BRDF effects (Noguchi et al., 2014).Therefore, the use of the geometry- dependent LER may result in overestimation of the BRDF effects.While non-absorbing aerosols are implicitly accounted for in the cloud algorithms (see Sect. 2.3.1), the aerosols directly affect the AMF, and thus trace-gas retrievals.
Figure 2 shows a data flow diagram that summarizes the generation of the geometry-dependent LER for satellite instruments.The diagram shows input data, all the steps of processing the data, and outputs.The input data include MODIS-derived land BRDF kernel coefficients and landwater flags, chlorophyll values, terrain heights from a digital elevation model, and wind speed at a height of 10 m.All the input data are averaged over a nominal OMI pixel using the OMI pixel corner information.The averaged input data, OMI pixel geometry, and atmospheric profiles are used to compute the TOA radiance with VLIDORT.The geometry-dependent LER is calculated from the TOA radiance using Eq. ( 7) and a pre-computed lookup table of radiative transfer parameters.Values of the geometry-dependent LER for each OMI pixel along with ancillary data are written in an HDF5-EOS output file for every OMI orbit.

Geometry-dependent LER
Because reflection of incoming direct and diffuse solar light from non-Lambertian surfaces depends on satellite observational geometry, the same area observed at different geometries can have different LERs. Figure 1 shows the MODISbased high spatial resolution LER over the Baltimore-Washington area of the United States for 2 consecutive days (17 and 18 January 2005) computed using the OMI observational geometry.The SZA and VZA values are in the similar ranges for both days.However, there is a large difference in the relative azimuth angle, which varies from around 63 • for 17 January to about 118 • for 18 January.Since the land tends to have strong backward scattering, this explains the higher LER for 18 January than that for 17 January.The differences, if not accounted for, may produce errors in the trace-gas retrievals.A comparison of the computed geometry-dependent and climatological LERs at 466 nm is shown in Fig. 3 for OMI orbit 12414 of 14 November 2006.The climatological LERs (monthly) are derived from OMI observations (Kleipool et al., 2008).In general, the eastern portion of the orbital swath (that has a later Equator crossing time) has higher values of the LERs than the western part.This is an effect of the OMI observational geometry and BRDF increases in the backscattered direction over land.Over the US, the western portion of the orbital swath has higher values of the LER than the eastern portion.This is explained by the dominance of the difference in climatological surface reflectances between the east and the west (see Fig. 3b) as compared with the BRDF increase due to the OMI observational geometry.
Figure 3 shows significant differences between the geometry-dependent and climatological LERs for both land and ocean.Over land, the climatological LERs are mostly higher than the geometry-dependent LERs.This is presumably because the geometry-dependent LERs are derived from atmospherically corrected MODIS radiances while the cli- matological LERs are affected by residual aerosols.Moreover, climatological LERs are inherently contaminated by clouds due to substantially larger sizes of OMI pixels as compared with those of MODIS.This is particularly true for the Amazonian region where clouds are persistent.
Over ocean, the geometry-dependent LERs are systematically higher than the climatological LERs in areas affected by sun glint and at large VZAs.This is because the climatological LERs are based on the mode of LERs from a long time series of observations over a given area; this minimizes the impact of observations affected by sun glint and high values that occur at large VZAs.The total ocean reflectance is comprised of three components: direct and diffuse solar light reflected from the ocean surface and water-leaving light.The fraction of each component strongly depends on geometry.Reflection of direct solar light dominates in the sun-glint area.At the edges of the swath the relative contribution of reflected diffuse light increases because the sky radiance increases to the horizons and the reflection angle increases thus the Fresnel reflection increases.The higher values of LER nearer to the eastern part of the swath than at the western part are mostly due to sky light reflected from the ocean surface.The angular distribution of the sky radiance is not symmetric in the plane of satellite observations because the sun is in the western part of the swath.The sky radiance is higher in the eastern part of the swath, and it is reflected at higher angles than the light from the western part.Additionally, the higher reflection angle results in higher Fresnel reflection in the eastern part of the swath.This is confirmed by our calculations of the view angle dependence of the reflected light only, i.e., no waterleaving radiance included.
Figure 4 shows the geometry-dependent LERs computed at 466 and 354 nm and their differences for the same OMI orbit 12414 of 14 November 2006.Here, we assume that the BRDF coefficients over land are spectrally independent.The LER differences over land are thus solely due to the smoothing effect of enhanced Rayleigh scattering in UV that increases the diffuse to direct incident irradiance ratio as compared with 466 nm.Over land, LER(354) < LER(466), but the differences are relatively small (< 0.015).
Over the ocean, the LER differences additionally result from the spectral dependence of water-leaving radiance.Over the sunglint areas, the solar light reflected from the ocean surface is significantly brighter at 466 nm than at 354 nm, thus leading to higher LERs.Over areas less affected by sunglint, LER(354) > LER(466) in general due to higher amounts of water-leaving radiance.
It is interesting to note that the patterns of rivers and their tributaries are evident in the LER maps of Fig. 4 for both 354 and 466 nm.This effect is most pronounced when rivers are viewed from the OMI measurement geometry that registers the reflectance signal of Fresnel reflection from smooth river surfaces.It may be somewhat surprising that this appears at OMI spatial resolution; we can explain the effect by considering that while the LER from FOVs comprised of river areas and surrounding land is weighted linearly by the areal fraction of each, reflectance from the river surface is disproportionally high due to the Fresnel reflection in sun-  glint geometry.Outside of the regions where OMI observes glint, the LER in the Amazon basin may still be higher than expected due to the turbidity of some rivers in the Amazon floodplain that varies seasonally.

RRS algorithm
Figure 5 shows ECFs computed with geometry-dependent LERs and the differences with respect to the climatological LERs ( ECF) for OMI orbit 12414 of 14 November 2006.The largest ECFs (up to 0.05) take place over the less cloudy Amazonian areas.ECF is obviously lower for cloudy areas due to the diminished effect of surface properties on TOA radiance.The heavily cloudy areas are easily identified on the ECF map.
Figure 6a is a scatter plot of ECF retrieved with the geometry-dependent LERs vs. ECF retrieved with climatological LERs for the entire range of ECF.It shows that the scatter of data around the 1 : 1 line diminishes with increasing ECF; i.e., the difference in ECFs decreases with increasing ECF as expected.We next examine the most interesting range of ECF for trace-gas retrievals, ECF < 0.25, which corresponds to f r < 0.4-0.5.For this range, Fig. 6b shows a scatter plot of the ECFs retrieved with the geometrydependent vs. climatological LERs and Fig. 6c shows how ECF varies with ECF.Only data from 50 • S to 50 • N are used in Fig. 6 and all subsequent similar figures.This latitude range excludes areas with snow for which MODIS BRDF data are not available.On average, ECF is small and positive for the ocean (∼ 0.02).Over land ECF is even lower and ranges from ∼ −0.01 to ∼ 0.015 for ECF < 0.25.The standard deviation of ECF does not depend much on ECF.It is ∼ 0.01 over ocean and ∼ 0.015 over land.Even though ECF is small on average, it can be as large as ±0.05 for individual FOVs, which is quite substantial for the low ECF range.Figure 6d shows normalized histograms of ECFs for 0.05 < ECF < 0.25.The normalized histograms of ECF retrieved with climatological LER and ECF retrieved with BRDF are close to each other.This reflects small differences between the ECFs on average.
Figure 7 similarly shows OCPs retrieved with the geometry-dependent LER and the differences with respect to those retrieved using the climatological LERs ( OCP) for OMI orbit 12414.There are no obvious geographical pat-terns in the OCP map.OCP can be as large as ±100 hPa.The OCP differences are particularly pronounced along the edges of cloud systems.Spatial correlation between OCP (Fig. 7) and ECF (Fig. 5) is not apparent.As may be expected, OCP decreases with increasing ECF. Figure 8 is similar to Fig. 6 but for OCP.On average, OCP is small (∼ 10.0 Pa) with standard deviation of up to ∼40 Pa.As it can be seen from Fig. 8c, changes of the OCP histograms due to replacing the climatological LERs with the geometrydependent LERs are relatively small.This is a consequence of small OCP differences on average.

O 2 -O 2 algorithm
Spatial distributions of the effective cloud fraction and cloud pressure retrieved from O 2 -O 2 are quite similar to those retrieved from RRS (shown in Figs. 5 and 7).That is why we do not show maps of ECF and OCP retrieved from O 2 -O 2 .
Here we show comparisons of the cloud products retrieved from O 2 -O 2 with the geometry-dependent and climatological LERs for ECF < 0.25.Figure 9 is similar to Fig. 6 but for ECF from the O 2 -O 2 algorithm.ECF < ∼ 0.03 over land and < ∼ 0.01 over ocean.The histograms of ECF retrieved with climatological LER and geometry-dependent LER (Fig. 9d) are similar to that from the RRS cloud algorithm (Fig. 6d).OCP is noticeably lower over ocean.The standard deviation, up to 100 hPa, is also higher than that from the RRS cloud algorithm.The histograms of OCP retrieved from the O 2 -O 2 cloud algorithm (Fig. 10c) noticeably differ from those retrieved from the RRS cloud algorithm (Fig. 8c).According to Fig. 10c, lower altitude clouds (with OCP > 800 hPa) are observed more frequently over the ocean than over land.For high-altitude clouds (OCP < 450 hPa) the situation is reversed: they are observed more frequently over land than over the ocean.Both patterns in the vertical distri-bution of clouds are much less pronounced in the histograms of OCP retrieved from the RRS algorithm.
The effect of replacing the climatological surface LERs with the geometry-dependent LERs is much more pronounced for the O 2 -O 2 OCP retrievals than for the RRS retrievals.This can be explained by two physical factors.Firstly, the Rayleigh optical depth of the atmosphere in the UV (the spectral window of the RRS cloud algorithm is 345-354 nm) is much higher than in the visible (the wavelength of the O 2 -O 2 OCP retrieval is 477 nm).Higher scattering in the UV leads to a larger fraction of diffuse light illuminating the surface, which decreases BRDF effects.In the visible, the smoothing effect of Rayleigh scattering is less than in the UV, thus resulting in larger BRDF effects.Secondly, sensitivities of the OCP derived from RRS and O 2 -O 2 to surface reflectivity are different for the RRS and O 2 -O 2 algorithms.The light path of direct sunlight reflected by the surface does not contribute to the RRS signal because there is no Raman scattering involved.However, this direct light path does contribute to the O 2 -O 2 absorption.That is why the RRS algorithm is generally less sensitive to the surface and to its BRDF.For high surface reflectivity, the reflected direct solar light significantly contributes to TOA radiance and therefore causes the OCP differences related to the absence of RRS in direct solar light and the presence of O 2 -O 2 absorption in direct solar light.However, for low surface reflectivity, this mechanism becomes less significant because the fraction of the reflected direct solar light in the TOA radiance is smaller.To make the numbers characterizing the ECF and OCP differences more representative, we processed OMI data for 2 days of 14 November and 14 July 2006.Figure 11 shows the ECF and OCP differences as a function of ECF for those 2 days.The ECF differences calculated for the entire day of 14 November 2006 (Fig. 11c) are close to those calculated for orbit 12414 of that day (Fig. 9b).The OCP differences over land calculated for the entire day (Fig. 11d) are slightly lower than those calculated for orbit 12414 of that day (Fig. 10b), while the OCP differences over ocean for the entire day are quite close to those calculated for one orbit.The ECF and OCP differences are similar for different seasons.A small increase of the OCP differences in November may not be statistically significant.The data in Fig. 11 indicate that the ECF and OCP differences obtained for OMI orbit 12414 are globally representative.

BRDF effects on the OMI NO 2 retrievals
We consider the BRDF effect on the NO 2 AMFs only, because the retrieved NO 2 amount is inversely proportional to the AMF.The geometry-dependent LER approach provides an exact match of TOA radiances with the full BRDF approach but not of the photon path lengths.This simplification can lead to some biases in the calculation of AMFs and thus to biases in the retrieved NO 2 vertical columns.Zhou et al. (2010) have estimated the biases.They compared the box AMFs and NO 2 vertical columns calculated with the full BRDF with that calculated with black sky albedo and white sky albedo.According to their data, maximum differences in the box AMFs are up to 10 % at the surface and differences in the NO 2 vertical columns are smaller than 12 %.We carried out calculations of NO 2 scattering weights and AMFs with full BRDF treatment and compared them with that calculated with the corresponding geometry-dependent LER. Figure 12a shows an example of the altitude dependence of scattering weights calculated with the full BRDF treatment and the geometry-dependent LER.It can be seen that the difference between the scattering weights is small.An AMF difference for this case is 5.6 %. Figure 12b shows a scatter plot of the full BRDF AMFs vs. the geometry-dependent LER AMFs calculated for OMI measurements over the eastern US for orbit 12414 of 14 November 2006.Differences in AMFs due to different treatment of the surface are within ±6 % (at 95 % confidence interval) and always less than 10 %.
The tropospheric NO 2 AMF, AMF trop , is calculated using the MLER model with input cloud parameters from the O 2 -O 2 algorithm assuming a priori NO 2 vertical profile shapes (see Fig. 13): The effect of a surface reflectivity change, R g , of 0.01 on AMF g is shown as a function of R g in Fig. 13.The Jacobian, J = AMF g / R g , is always positive because larger surface reflectances increase satellite sensitivity to NO 2 absorption in the lowest atmosphere.J decreases with increasing R g and for unpolluted NO 2 mixing ratio profiles (Fig. 13).
An effect of replacing the climatological LERs with geometry-dependent LERs on AMF g for OMI observational geometries and ground resolution can be estimated from Figs. 3 and 13 using R g = LER(BRDF) − LER.The effect is largest over polluted regions in the eastern US, where R g is negative with values from −0.03 to −0.02 (Fig. 3), LER ∼ 0.05, and AMF g ∼ −20 to −30 %.The BRDF effect reverses over water for glint geometries and large viewing angles, but R g is large here and the effect on AMF g is reduced (i.e., small Jacobian).
To estimate the BRDF effect on AMF trop we need to account for the f r change as well.By differentiating Eq. ( 8) and assuming that AMF g and f r are independent variables and both depend on R g , we get The cloud AMF strongly depends on the OCP, since high clouds (low OCP) have a shielding effect and low clouds (high OCP), aerosols, and fog can enhance AMF c .Assuming a negligible NO 2 mixing ratio above the cloud OCP, we can neglect AMF c and Eq. ( 9) simplifies to Over land, replacing the climatological LERs with the geometry-dependent LERs reduces the surface LER on average (i.e., R g < 0), leading to smaller values of AMF g (Fig. 13).At the same time, the mean ECF increases by 0.02 (Fig. 9) and this produces even larger increases in f r ( f r ∼ 0.04).Therefore both terms in the above equation are negative meaning that switching to a geometry-dependent LER reduces AMF trop even more over land.The effect is mixed over water, since both R g or f r can change signs for certain geometries.It should be noted that we derived Eqs. ( 9) and (10) to qualitatively illustrate the effect of changing surface reflectance on AMF in cloudy conditions.The equations are not used to produce data in the figures of Sect.3.3.The data in the figures are obtained numerically using Eq. ( 8).
Figure 14 shows that the calculated impact on AMF trop arising from replacing the climatological OMI-based LERs with geometry-dependent LERs exhibits a strong spatial variation with smaller effects over ocean, unpolluted, or cloudy areas.Over land, where the geometry-dependent LER is generally lower than the climatological LER, use of the BRDF data results in lower AMFs and higher tropospheric NO 2 VCDs.The effect is enhanced over polluted areas such as eastern US, where the changes in AMF can reach up to 50 %.The effect is reduced for unpolluted and overcast conditions and mixed over oceans, because R g increases for sunglint and large VZA directions but decreases for other directions.
Figure 15 compares the clear-sky AMF trop calculated using climatological and geometry-dependent LERs for OMI orbit 12414.The use of geometry-dependent LERs generally leads to lower AMF trop by up to 29 % over land and 15 % over ocean.Differences in O 2 -O 2 cloud parameters resulting from the use of geometry-dependent LERs add additional scatter, changing AMF trop by −42 to 5 % over land and −22 to 13 % over ocean.AMF trop differences are large for low AMFs, driven by enhanced differences in LER, OCP, or f r .To make the numbers characterizing the AMF differences be more representative, we calculated the tropospheric NO 2 AMFs using the geometry-dependent LER and compared them with those calculated with the climatological LER for 2 days: 14 November and 14 July 2006.The AMF differences arising from both replacing the climatological LERs with the geometry-dependent LERs and changing the cloud parameters exhibit strong spatial variations with smaller effects over the ocean, unpolluted, or cloudy areas similar to Fig 14 .A global analysis of the AMF differences shows that the AMF differences for OMI orbit 12414 are consistent with those for both days.Figure 16 illustrates how the use of geometry-dependent LER changes NO 2 retrievals over clean and polluted areas.Consistent with previous studies by Lin et al. (2014Lin et al. ( , 2015)), AMFs are considerably lower with geometry-dependent LERs.This suggests that the current operational NO 2 products based on climatological LERs could be underestimated by up to 48 % over China.The eastern US exhibits similar but somewhat smaller differences.Minor changes are expected over unpolluted and overcast areas.

Conclusions
We developed a new concept of geometry-dependent surface LER and provided a means for computing it.Spatially averaged high-resolution MODIS BRDFs are used for computation of the geometry-dependent LER over land for OMI pixels.The Cox-Munk slope distribution function and the Case 1 water-leaving radiance model are utilized for computation of the geometry-dependent LER over ocean.This method accounts for the geometrical dependence of LER within the existing framework of MLER trace-gas and cloud algorithms with only minimal changes.It is important to note that the geometry-dependent surface LER approach can be applied to any current or future satellite algorithms that utilize MLER trace-gas and cloud algorithms.
We examined the effects of the geometry-dependent LER on OMI cloud and NO 2 algorithms.The effects on retrieved cloud parameters were relatively small on average and diminish with increasing cloud fraction.Even though the impact is small on average, it can be as large as ±0.05 for the effective cloud fraction and 100 hPa for the cloud optical centroid pressure.It should be noted that the background aerosols are included in the climatological LER; therefore, they are virtually accounted for in the ECF derived using the LER climatology.The geometry-dependent LER is calculated for aerosol-free conditions; thus the corresponding ECF should have a bias.The BRDF effects were noticeably higher for the O 2 -O 2 algorithm that uses visible wavelengths as compared with the RRS algorithm that utilizes a UV spectral range.This can be explained by the stronger smoothing effect of Rayleigh scattering in the UV as compared with the vis.
We also find that replacing the climatological OMI-based LERs with geometry-dependent LERs can increase the OMI NO 2 vertical columns by up to 50 % over highly polluted areas.Only minor changes to NO 2 columns (within 5 %) are found over unpolluted and overcast areas.It should be noted that the differences include both BRDF effects and biases between the MODIS and OMI-based surface reflectance data sets.
In the future, we plan to implement the use of geometrydependent LERs in our cloud and NO 2 OMI algorithms.Along with the use of the geometry-dependent LER product, we plan to explicitly include aerosols in the NO 2 algorithm.Further evaluation of the results with OMI data is ongoing.The proposed method of generating the geometry-dependent LERs is computationally expensive.We plan to reduce computational cost by using a neural network approach to replace VLIDORT calculations.We also plan to investigate the use of a new surface BRDF product from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm (Lyapustin et al., 2012).

Figure 1 .
Figure 1.High spatial resolution MODIS-based LERs for the Baltimore-Washington area of the United States for 17 (left) and 18 (right) January 2005 computed with the original spatial resolution of 30 arcsec but for OMI observational geometries.

Figure 2 .
Figure 2. Data flow diagram of generating the geometry-dependent LER.

Figure 3 .
Figure 3. (a) LERs computed at 466 nm for OMI orbit 12414 of 14 November 2006 using MODIS-based BRDF with OMI geometry, (b) OMI-based monthly climatology, and (c) their difference: MODIS-based minus climatological LERs.Missing MODIS BRDF data are shown in grey here and elsewhere.

Figure 5 .
Figure 5. (a) RRS-retrieved ECF computed with geometry-dependent LERs and (b) the difference between the ECFs computed with geometry-dependent and climatological LERs for OMI orbit 12414 of 14 November 2006.ECF data correspond to LER shown in Fig. 4b.

Figure 6 .
Figure 6.(a) Scatter plot of RRS-retrieved effective cloud fractions (ECFs) computed with geometry-dependent LERs vs. climatological LERs, where the 1 : 1 line is in black; (b) similar but for 0.05 < ECF < 0.25 with linear fits; (c) the mean ECF difference (diamonds) and standard deviation (error bars) as a function of ECF; (d) normalized histograms of ECF for 0.05 < ECF < 0.25.Data are for OMI orbit 12414 of 14 November 2006.

Figure 7 .
Figure 7. (a) RRS-retrieved cloud optical centroid pressure (OCP) computed with geometry-dependent LERs and (b) the difference between the OCPs computed with geometry-dependent and climatological LERs for OMI orbit 12414 of 14 November 2006.

Figure 8 .
Figure 8.Comparison of RRS-retrieved OCPs computed with geometry-dependent and climatological LERs for OMI orbit 12414 of 14 November 2006; data are for 0.05 < ECF < 0.25.(a) Scatter plot with regression lines; (b) the mean OCP difference and the standard deviation; (c) normalized histograms of OCP.
Figure 10 is similar to Fig. 8 but for OCP from the O 2 -O 2 algorithm.OCP has values up to 200 hPa.The mean OCPs are significantly larger for the O 2 -O 2 algorithm as compared with RRS.On average, OCP varies from ∼80 hPa at ECF = 0.05 to 5 hPa at ECF = 0.25 over land.

Figure 9 .
Figure 9. (a) Scatter plot of O 2 -O 2 -retrieved ECFs computed with geometry-dependent LERs vs. climatological LERs, where the 1 : 1 line is in black; (b) similar but for 0.05 < ECF < 0.25 with linear fits; (c) the mean ECF difference and the standard deviation; (d) normalized histograms of ECF.

Figure 10 .
Figure 10.(a) Scatter plot of O 2 -O 2 -retrieved OCPs computed with geometry-dependent LERs vs. climatological LERs with linear fits; data are for 0.05 < ECF < 0.25.(b) The mean ECF difference (diamonds) and standard deviation (error bars) as a function of ECF.(c) Normalized histograms of OCP.

Figure 11 .
Figure 11.Comparison of O 2 -O 2 -derived ECFs and OCPs computed with geometry-dependent and climatological LERs; data are for 0.05 < ECF < 0.25.(a) The mean ECF difference and (b) OCP difference (diamonds) and standard deviation (error bars) as a function of ECF for 14 July 2006; (c) the mean ECF difference and (d) OCP difference (diamonds) and standard deviation (error bars) as a function of ECF for 14 November 2006.

Figure 12 .
Figure 12.(a) Scattering weights calculated for an OMI pixel with the following observation angles: SZA = 53.1 • , VZA = 54.1 • , and relative azimuth angle RAZ = 60.3 • .(b) AMF calculated with full BRDF vs. BRDF-derived LER for OMI measurements over the eastern US on 14 November 2006.A regression line is shown in blue.The data are for clear skies.

Figure 13
Figure 13.(a) November mean NO 2 profiles at three locations in the eastern US from the NASA GMI model; (b) air mass factor (AMF) change due to 0.01 change in surface reflectivity as a function of surface reflectivity.Red: highly polluted profile; green: moderately polluted; blue: unpolluted profile.

Figure 14
Figure 14.(a) OMI tropospheric NO 2 air mass factor (AMF) calculated using geometry-dependent MODIS-based LER; (b) percent differences with respect to climatological LERs for OMI orbit 12414 on 14 November 2006.

Figure 15 .
Figure 15.Scatter diagrams of AMFs calculated using geometric-dependent MODIS-based LER vs. OMI-based climatological LER for the orbit 12414 for clear to moderately cloudy sky (f r < 0.5) including (a) effects of BRDF only with clouds unchanged and (b) the effects of both BRDF and O 2 -O 2 cloud parameters for land (blue) and ocean (orange).Numbers in parentheses represent percent difference at the 2nd and 98th percentile range.Percent difference in AMF with changes in surface BRDF and O 2 -O 2 cloud parameters, sorting the data by the difference with respect to (c) LER, (d) OCP, and (e) f r .