Skip to main content

ATBD Ocean Algorithm

4. Algorithm description: ocean

4.1 Strategy

The basic strategy of the C6 dark target aerosol retrieval over the ocean (C6-O) is unchanged from that described in ATBD-96 and ATBD-09 (also in Tanré et al., 1997, Levy et al., 2003; Remer et al., 2005). The mechanics of the ocean algorithm are illustrated below in Figure 4. The algorithm is based on a ‘look-up table’ (LUT) approach, i.e., radiative transfer calculations are used to generate pre-computed theoretical outputs for a set of aerosol and surface parameters and compared with the observed radiation field. The algorithm assumes that one fine and one coarse lognormal aerosol modes can be combined with proper weighting η, to represent the ambient aerosol properties over the target, i.e.,


Spectral reflectance from the LUT is compared with MODIS-measured spectral reflectance to find the ‘best’ (least-squares) fit. The best fit, or an ‘average’ of a set of the best fits is used as the solution to the inversion. Although the core inversion remains similar to the process described in Tanré, et al. (1997), and updates including sediment masking and the special handling of heavy dust including dust retrievals over glint included in ATBD-09, minor updates to masking of clouds and the inclusions of multiple wind speed look-up tables are new. In this section, we describe the over-ocean LUTs, the cloud masking and pixel selection logic, the retrieval strategy details, and over-ocean product description.

4.2 Aerosol models (modes) and lookup tables

The aerosol optical models for C6 are unchanged from C5. They were derived mainly from data gleaned from sun/sky photometers (e.g. AERONET; Holben et al., 1998), and from analysis of errors in the products from previous versions of the MODIS algorithm. For C6-O, there are four fine modes and five coarse modes, with spectral refractive indices and lognormal size parameters listed in Table 2A. Based on analyzing Dubovik et al. 2000, coarse modes 7, 8 and 9 were modified from those used in C4 (Remer et al., 2005). Based on Mie code, tables 2B-D describe spectral properties of the aerosol modes, for unit AOD=1.0 (at 0.55 μm). 






Like all previous versions, the DT LUT employs the MODRAD (Ahmad et al., 1982) radiative transfer (RT) code to simulate TOA reflectance for a coupled ocean/atmosphere. The spectral reflectance at the satellite level is assumed to be from a combination of radiation from the surface and the atmosphere. The ocean surface calculation includes sun glint reflection off the surface waves (Cox and Munk, 1954), reflection by whitecaps (Koepke, 1984) and Lambertian reflectance from underwater scattering (sediments, chlorophyll, etc). In previous collections (through C5), the surface wind speed was assumed at 6.0 m/s. However, work by Zhang and Reid (2010), noted that near surface wind patterns significantly influence ocean wave and glint patterns, and incorrect assumptions about these patterns can bias the aerosol retrieval. Concurrently, Kleidman et al. (2012) compared MODIS C5 DT-ocean data with hand-held Microtops sun photometer data from the Marine Aerosol Network (MAN) (Smirnov et al., 2009) and found that there were residual MODIS errors related to wind speed. Sensitivity studies suggested that the problem would be enhanced closer to glint. Following other algorithm teams (e.g., Sayer et al., 2012a; Herman et al., 2005; Kahn et al., 2007), we now introduce wind speed dependence to the MODIS DT-ocean aerosol retrieval. This takes on the form as an additional step in interpolation of the DT LUT.

Embedded within MODRAD are wind speed dependent models to account for the “roughness” of the sea surface (waves and whitecaps, Cox and Munk, 1954) and the foam fraction (Koepke, 1984). In addition to the standard 6 m/s wind speed with 0.16%, the DT LUT includes simulations for three additional wind speeds, 2 m/s, 10 m/s  and 14 m/s, having foam fractions of 0.01%, 1% and 3%, respectively (e.g. Monahan and Muircheartaigh, 1980). The LUT is interpolated with respect to the 2m wind speed from the NCEP-GDAS 1˚ × 1˚ analysis. Wind speeds less than 2 m/s are assumed to be 2 m/s, and greater than 14 m/s are assumed to be 14 m/s; otherwise the LUT is linearly interpolated between the nearest two indices. Zero water leaving radiance is assumed for all compared wavelengths, except for at 0.55 μm, where a fixed reflectance of 0.005 is used. The atmospheric contribution includes multiple scattering by gas and aerosol, as well as reflection of the atmosphere by the sea surface. Thus, spectral reflectance was computed for each of the nine aerosol models described in the Table 2. Six values of  τa, (normalized to 0.55 μm) were considered for each mode, ranging from a pure molecular (Rayleigh) atmosphere (τa = 0.0) to a highly turbid atmosphere (τa = 3.0), with intermediate values of 0.2, 0.5, 1.0 and 2.0. For each model and aerosol optical depth at 0.55 μm, the associated aerosol optical depths were stored for the other six wavelengths, including the blue (at 0.47 μm). Computations are performed for combinations of 11 solar zenith angles (6°, 12°, 24°, 36°, 48°, 54°, 60°, 66°, 72°, 78° and 84°), 16 satellite view zenith angles (0° to 72°, increments of 6°) and 16 relative sun/satellite azimuth angles (0° to 180°, increments of 12°) for a total of 2816 angular combinations.

4.3  Pixel selection: cloud glint and sediment masking

The masking of clouds and sediments and the selection of pixels are described in Remer et al., (2005) and Levy et al., (2013).  Prior to ocean retrieval, the reflectance is organized into 10 km (at nadir) boxes of 20 by 20 pixels (at nominal 500 m resolution) for the 10 km product, and 3 km boxes of 6 by 6 pixels for the 3 km product, and are corrected for gas (water, carbon dioxide and ozone) absorption (see Appendix 1). This requires degrading the resolution of the 250 m channels (ρ0.65 and ρ0.86) to match the 500 m resolution of the other channels used in the retrieval. Note that if land is encountered in any of the 400 pixels, the entire box is processed using the land algorithm.

Much attention has been paid in the algorithm to the difficult task of separating usable “clear” pixels from cloudy or cloud contaminated pixels.  The standard MxD35 cloud mask uses many tests including the brightness in the visible channels to identify clouds.  This procedure can incorrectly identify heavy aerosol as clouds, and will thereby miss retrieving significant aerosol events over ocean.  On the other hand, relying on IR-tests alone permits low altitude, warm clouds to escape and be misidentified as 'clear', introducing cloud contamination to the aerosol products. Therefore, our cloud mask over ocean combines spatial variability tests (e.g. Martins et al., 2002) along with tests of brightness in visible and infrared channels. Appendix 2 describes the cloud mask over ocean. To identify underwater sediments in shallow water (near coastlines) the masking procedure of Li et al. (2003) is used in addition to the cloud mask.

The algorithm sorts the pixels that remain after cloud and sediment masking according to their ρ0.86 value, discarding the darkest and brightest 25%, and thereby leaves the middle 50% of the data.  This filter is used to eliminate residual cloud contamination, cloud shadows, or other unusual extreme conditions in the box.  Because the ocean cloud mask and the ocean surface are expected to be less problematic than their counterparts over land, the filter is less restrictive than the one used in the land retrieval. For the 10 km retrieval at least 10 of the 400 pixels in the original box must remain after masking and filtering. For the 3 km retrieval, at least 5 of the 36 pixels in the original box must remain. If the threshold for the minimum number of pixels is not met, no retrieval is attempted and fill values are given for the aerosol products in the retrieval box.  If there are at least 10 good pixels in the 0.86µm channel and at least 30 good pixels in the remaining 5 channels, the mean reflectance and standard deviation are calculated for the remaining 'good' pixels at the six pertinent wavelengths.

Ocean Glint and Internal Consistency:

The ocean algorithm was designed to retrieve only over dark ocean (i.e. away from glint). There is a special case when we retrieve over glint which is described below. The algorithm calculates the glint angle, which denotes the angle of reflection, compared with the specular reflection angle. The glint angle is defined as

  Θglint=cos–1((cosθscosθv)+(sinθssinθvcosΦ))                            (25)

where θs, θv, and Φ are the solar zenith, the satellite zenith and the relative azimuth angles (between the sun and satellite), respectively. Note that Fresnel reflection corresponds to Θglint = 0.  If Θglint > 40˚, we can avoid glint contamination and proceed with the retrieval. The algorithm performs several consistency checks of the spectral reflectances. Depending on the outcome of these consistency checks, the algorithm may either declare the reflectances to be beyond the range necessary for a successful inversion and exit the procedure, or continue onto the inversion after assigning quality flags (Quality Assurance – QA) to each wavelength.

4.4 Retrieval Algorithm

The ocean inversion procedure described in ATBD-96, Levy et al., (2003) and Remer et al., (2005) remains the same for the C6-O algorithm. Following Tanré et al. (1996), we know that the MODIS-measured spectral radiance (0.55 – 2.13 µm) contains almost three pieces of independent information about the aerosol loading and size properties. With some assumptions, the algorithm can derive three parameters: the τ at one wavelength (), the ‘reflectance weighting parameter’ (the over-ocean definition of Fine Weighting - η) at one wavelength (ηλ=0.55 ) and the ‘effective radius’ (re), which is the ratio of the 3rd and 2nd moments of the aerosol size distribution. The effective radius is represented by choosing a single ‘fine’ (f) and single ‘coarse’ (c) aerosol mode for weighting with the η parameter. The inversion is based on the look-up table (LUT) of four fine modes and five coarse modes (Table 2).  Although the LUT is defined in terms of a single wavelength of optical thickness, the parameters of each of the single mode models define a unique spectral dependence for that model, which is applied to the retrieved value of , to determine optical thickness at the other reported wavelengths.

The retrieval requires a single fine mode and a single coarse mode. The difficulty is in determining which of the (4 x 5 =) twenty combinations of fine and coarse modes and their relative optical contributions that best mimics the MODIS-observed spectral reflectance. To accomplish this the reflectance from each mode is combined using η as the weighting parameter,



where  is a weighted average reflectance of an atmosphere with a pure fine mode 'f' and optical thickness  and the reflectance of an atmosphere with a pure coarse mode 'c' also with the same.  For each of the twenty combinations of one fine mode and one coarse mode, the inversion finds the pair of  and η0.55  that minimizes the ‘fitting error’ (ε) defined as




where Nλ is the sum of good pixels at wavelength λ,  ρλm  is the measured MODIS reflectance at wavelength λ, ρλray  is the reflectance contributed by Rayleigh scattering, and ρλLUT   is  calculated from the combination of modes in the look-up table and  defined by Equation (26). The 0.01 is to prevent a division by zero for the longer wavelengths under clean conditions (Tanré et al. 1997).  Note the inclusion of the Rayleigh reflectance scattering in Equation (27), as compared with the formula in C4-O (Remer et al., 2005). The inversion requiresρ0.87 LUT  to exactly fit the MODIS observations at that wavelength and then finds the best fits to the other five wavelengths via Equation (27).  The 0.87 µm channel was chosen to be the primary wavelength because it is expected to be less affected by variability in water leaving radiances than the shorter wavelengths, yet still exhibit a strong aerosol signal, even for aerosols dominated by the fine mode. By emphasizing accuracy in this channel variability in chlorophyll will have negligible effect on the optical thickness retrieval and minimal effect on ηλ=0.55.

The twenty solutions are then sorted according to values of e. The ‘best’ solution is the combination of modes with accompanying   and ηλ=0.55 that minimizes ε.  The solution may not be unique.  The ‘average’ solution is the average of all solutions with ε< 3% or if no solution has ε < 3%, then the average of the 3 best solutions. Once the solutions are found, then the chosen combination of modes is the de facto derived aerosol model and a variety of parameters can be inferred from the chosen size distribution including spectral optical depth, effective radius, spectral flux, mass concentration, etc.

      Final Checking

Before the final results are output, additional consistency checks are employed.  In general, if the retrievedτ at 0.55 µm is greater than –0.01 and less than 5, then the results are output. Negative optical depths are rounded to exactly zero. There are exceptions and further checking for heavy dust retrievals made over the glint.  The final QA-confidence flag may be adjusted during this final checking phase.

    Special case: Heavy dust over glint.

If Θglint ≤ 40˚, a case where we normally do not retrieve, we check for heavy dust in the glint.  We use a similar technique as before during the masking operations when we noted that heavy dust has a distinctive spectral signature because of light absorption at blue wavelengths.  In the situation of identifying heavy dust over glint we designate all values of ρ0.47/ ρ0.65  < 0.95 to be heavy dust.  If heavy dust is identified in the glint, the algorithm continues with the retrieval, although it sets QAC=0.  This permits the retrieval, but prohibits the values from being included in the Level 3 statistics (Remer et al., 2005). If heavy dust is not identified in the glint, then the algorithm writes fill values to the aerosol product arrays and exits the procedure.

4.5 Algorithm Flow Chart

Figure 4: Flowchart of the over-ocean aerosol retrieval algorithm (adapted from Remer et al., 2005, updated for C6).     

4.6 Retrieved ocean products

As discussed earlier, the primary retrieved products of the ocean algorithm are the total τ at 0.55 μm (), the Fine (Mode) Weighting (η or η 0.55 )  and which Fine (f) and which Coarse (c) modes were used in the retrieval. The fitting error (ε) of the simulated spectral reflectance is also retrieved. Both the ‘best’ and ‘average’ solutions are reported. From these primary products a number of other parameters can be easily derived. Examples include the effective radius (reff) of the combined size distribution, the spectral total, fine and coarse t values (, and ), and a measure of the columnar aerosol mass concentration (MC). Table 3 lists the products retrieved and derived during the ocean retrieval algorithm. Details about the ‘Quality_Assurance_Ocean’ (QA) SDS are given in the Appendix. Note that a parameter’s type does not indicate whether a parameter should be used in a quantitative way.  Each parameter must be independently validated via comparison with ground-truth measurements, where validation means that the product can be shown to be comparable with a set of ground-truth observations on a global scale, within some defined expected uncertainty. Prior to operational C6 processing, a testbed of 8 months of Aqua-MODIS granules were used to evaluate a few of the derived parameters. This provisional evaluation is described in Section 5. This ATBD document will be updated as the MODIS products are evaluated and reach validation status


Some of the ocean products are combined with products from land (discussed in the next section) as the Joint products. For τ , two joint products are reported, the ‘Optical_Depth_Land_And_Ocean’ and the ‘Image_Optical_Depth_Land_And_Ocean’.  The first product is constrained by QAC whereas the “Image” product includes all QAC values in order to provide the most complete visulization of the AOD.