MODIS Aerosol Retrieval Algorithm for Land

The surface albedo of land surfaces are usually unknown and vary with wavelength. Dark, dense vegetation has a very low reflectance in the blue and red regions. This has allowed the development of algorithms to retrieve aerosol optical depth over dark vegetation (called dark targets). Dry or unvegetated surfaces, such as deserts, and snow or ice covered surfaces reflect too much light to be used for aerosol retrieval with current methodologies. These areas are masked out in the MODIS data products showing aerosol loading over land. In the images below, global gridded results are shown for sample months of the MODIS Optical Depth over Land and Ocean product (from MOD08_M3). The values over land masses are the aerosol model corrected optical depths (explained below). Notice how the masked areas (black w. red stippling) around the globe vary in size and location with the changing of the seasons (images courtesy of Goddard Space Flight Center, NASA).

June 2001              

Sept. 2001

Dec. 2001

Mar 2002

May 2002            

The MODIS spectral channels used in the aerosol retrieval over land are:

0.47 µm (500 m resolution) 0.66 µm (250 m resolution) 2.1 µm (500 m resolution) 3.8 µm (1 km resolution)
The latter two channels are used for the determination of surface reflectance. The 11 µm window channel is used to correct for IR emission from the Earth surface. Other gaseous absorption channels, such as the near-IR (~1 µm) water vapor channel (for the correction of total precipitable water vapor) and the 9.6 µm channel (for total ozone correction) are also used. In Part I we provide a simplified description of the steps followed in the retrieval of aerosols over land. Part II consists of a more detailed explanation of the theory and methodology behind the land algorithm.

Part I. Steps of the MODIS Over Land Algorithm

1. Pixels that are dark in the blue (0.47 µm) and red (0.66 µm) channels are identified using their remotely sensed reflectance in the mid-IR channels (2.1 and 3.8 µm). This is possible owing to a correlation between land reflectance in the mid-IR and in the blue and red channels. Because aerosols do not interfere much with the transmission of these long thermal infrared wavelengths, the reflectance of mid-IR observed by the satellite can be assumed to be essentially the reflectance of the land surface, and not of particles in the atmosphere above it. 2. For pixels that are identified to have low reflectance in the mid-IR, surface reflectance (rs) in the blue and red is derived. Predictable relationships between IR reflectance (X axis) and reflectance in red (eg. 0.66 µm) or blue (eg. 0.49 µm, Y axes), have been discovered for land targets as straight lines in scatterplots. Multiplication factors are obtained as the slopes of lines of the scatterplots. 3. Preliminary optical depths (OD, or t) at 0.47 and 0.66 µm are obtained using satellite observed reflectances (r*0.47 and r*0.66) and surface reflectances (rs0.47, and rs0.66) in a look-up table generated by a continental aerosol model. Also, the single scattering path radiances (Lp) for 0.47 and 0.66µm are obtained. 4. For each grid box of pixels, an educated "guess" is made about what kind of aerosol is likely to be predominant at that location and time of the year. A decision tree is used that involves (1) OD0.47, (2) the ratio of the single scattering path radiances (Lp0.47/Lp0.66), and (3) information on global aerosol distribution. 5. Having "decided" which kinds of aerosols predominate in a given grid box, an aerosol model is selected that describes the aerosol size distribution, assigns values for refractive index and single scattering albedo, and describes the effect of nonsphericity on the phase function. The models are derived from studies at different places around the world in which in situ measurements were taken and ground instruments used to measure the ambient column aerosol size distribution. Four models are utilized (1) a continental model (same as used in step 3), (2) a biomass burning model, (3) an industrial/urban aerosol model, and (4) a dust model. 6. The chosen model is used to generate a look-up table. Satellite measured reflectance (r*) is inverted into aerosol optical depth (t), volume (or mass) concentration and spectral radiative forcing. 7. t0.47 and t0.66 are corrected for the difference between the model chosen in step 5 and the continental model. An interpolation is performed to obtain optical depth at 0.55 µm (as opposed to 0.47 or 0.66 µm). 8. Daily optical depth and mass concentration are stored on a resolution of 10x10 pixels (of 1 km nadir resolution).

Part II. A More Detailed Description of the Algorithm

Some radiative transfer theory to get us started:

The remote sensing of aerosol over the land stems from the relationship between the measured reflectance at the top of the atmosphere r* (given in apparent reflectance units - r* = pL/Foµo, where L is the radiance at the top of the atmosphere, Fo is the extraterrestrial solar flux and mo is the cosine of the solar zenith angle), and r(q,qo,f) is the surface bi-directional reflectance:

r*(q,qo,f) = ra(q,qo,f) + Fd(qo)T(q )r(q,qo,f)/(1-sp')

where q is the view direction, qo is the solar zenith angle, and f is the azimuth. ra(q,qo,f) is the path radiance, Fd(qo) is the normalized downward total flux for zero surface reflectance, equivalent to the total downward transmission. Its value is less than 1.0 due to aerosol and molecular absorption and backscattering of sunlight to space. T(q) is the upward total transmission into the direction of the satellite field of view, s the atmospheric backscattering ratio and r' the surface reflectance averaged on the view and illumination angles. In the single scattering approximation, the path radiance, ra(q,qo,f), is proportional to the aerosol optical depth, ta , the aerosol scattering phase function, Pa(q,qo,f), and single scattering albedo, wo:

ra(q,qo,f) = rm(q,qo,f) + wotaPa(q,qo,f)/(4µµo)

where rm(q,qo,f) is the path radiance due to molecular scattering, and µ and µo are cosines of the view and illumination directions, respectively. The quantity Lp = wotaPa is referred to as the single scattering path radiance.

What the satellite measures versus what we want to know:

The satellite measurement gives us (1) the observation geometry (q, qo, f, µ, and µo), and (2) the measured reflectance at the top of the atmosphere (r*). We want optical depth, ta. In order to derive the aerosol optical depth from the measured reflectance, an aerosol model that provides pre-determined values of wo and Pa for specific conditions is required. Later, we will see how an appropriate model is selected for each pixel in order to obtain the pre-determined values of wo and Pa.

Now, lets consider what's on the ground.

The radiative effect of aerosols includes backscattering and absorption of (1) the direct sunlight and (2) sunlight reflected from the surface. The more that the surface reflects light, the less the path radiance (ra) contributes to the satellite measured radiance (r*), and the greater will be the errors associated with deriving optical depth from r*. In order to reduce error, aerosol optical depth over land is only calculated over land surfaces that have low surface reflectance in the MODIS visible bandwidths. These land surfaces are called "dark targets". In fact, the remote sensing of aerosol using dark targets can be best done for surface reflectance r £ 0.06.

Which land surfaces are dark targets?

Many land covers (such as vegetation and some soils, especially wet soils) reflect little light in the red (0.60-0.68 µm) and blue (0.40-0.48 µm) wavelengths. It is said that they are "dark" in the red and the blue. Therefore, it is appropriate to use the darkest pixels in the image to estimate the aerosol optical depth (or loading) and its effect on remote sensing at these wavelengths. Before we can use the satellite measured reflectance (r*) to derive path radiance and optical depth over dark land pixels, we have to estimate the surface reflectance at the same wavelengths. In fact, the surface reflectance of these dark pixels has to be estimated within a small uncertainty of Dr = ±0.005 to ±0.01.

Selection of dark pixels, and determination of their surface reflectance (rs).

The satellite measured reflectances at 2.1 and 3.8 µm (r*2.1 and r*3.8) are used to determine which pixels are sufficiently dark to be used in the algorithm. Why? First, we should understand that reflectance of light by land surfaces across the solar spectrum is correlated to some extent. Soils usually reflect more light as wavelength increases, and the correlation between the reflectances slowly decreases as the wavelength span increases. However, parallel processes affect the reflectance in the visible (eg. 0.47 and 0.66 µm) channels and in the mid-IR channels (2.1 and 3.8 µm). Vegetation decreases reflectivity in the visible channels, due to chlorophyll absorption, but also in the mid-IR channels due to absorption by liquid water associated with the plant. Wet soil has a lower reflectance in the visible channels due to light trapping, but also in the 2.1 and 3.8 µm channels due to the liquid water absorption. Moreover, surface roughness, shadows and inclinations decrease the reflectance across the whole solar spectrum. This means that if we know the surface reflection of light in the mid-IR, we can estimate the surface reflection of light in the visible wavelengths. Second, we should understand that the longer the wavelength, the smaller the effect that aerosol scattering has on the reflectance of light measured in space. Except for dust, the aerosol effect on the radiance measured from space decreases with wavelength as l-1 to l-2. Longer wavelengths (2.1 or 3.8 µm) are less sensitive to aerosol scattering than visible wavelengths (since these wavelengths are much longer than the size of most aerosol particles) but are still sensitive to land surface characteristics. The reflectance of mid-IR observed by the satellite can be assumed to be essentially the reflectance of the land surface, and not of particles in the atmosphere above it.

How mid-IR reflectance is used to derive surface reflectance (rs) at 0.47 and 0.66 µm.

1. Pixels are grouped into grid boxes. When the spatial resolution of the channel is 1 km, there are 10x10 pixels per grid box. For channels with higher resolution (500m and 250m), more pixels are included in the grid boxes. 2. A series of criteria involving the magnitude of satellite observed mid-IR reflectance (r*2.1 and r*3.8) are applied to each grid box, and the number of pixels (N) that satisfy each criterion is found. The first criterion in the series that yields an N > 5% of the pixels in the grid box is chosen as the winning criterion. The criteria are used only over land surfaces excluding open water, clouds, ice and snow. 3. rs0.47 and rs0.66 are obtained. The criteria:

first priority : If N > 5% when 0.01 £ r*2.1 £ 0.05, then rs0.47 = r*2.1/4, and rs0.66 = r*2.1/2

second priority : If N > 5% when r*3.8 £ 0.025, then rs0.47 = 0.01, and rs0.66 = 0.02

third priority : If N > 5% when 0.01 £ r*2.1 £ 0.10, then rs0.47 = r* 2.1/4, and rs0.66 = r*2.1/2

fourth priority : If N > 5% when 0.01 £ r*2.1 £ 0.15, then rs0.47 = r*2.1/4, and rs0.66 = r*2.1/2 pixels with r*2.1 > 0.15 are not used.

Obtaining preliminary optical depths (OD) at 0.47 and 0.66 µm

Preliminary optical depths (OD) at 0.47 and 0.66 µm are obtained by assuming that the aerosol population conforms to a realistic, generalized model for aerosols over continents. This continental model assumes that the aerosols are a mixture of small and large water soluble particles, large dust like particles and small soot particles. The model assigns values for the size and volume distribution, single scattering albedo (shown in the table below), and refractive index, asymmetry factor and phase function (not shown in the table). The values obtained in the previous step (surface reflectances rs0.47 and rs0.66) plus the satellite observed reflectances (r*0.47 and r*0.66) are used with the model parameters to create a look-up table. From the look-up table, OD0.47cont, OD0.66cont and the single scattering path radiances (Lp) for 0.47 and 0.66µm are obtained.

Aerosol Properties of the Continental Model

Aerosol sub-category

mean particle radius of number distribution

rg(µm)

mean particle radius of volume distribution

rv(µm)

standard deviation of ln(r)

s

column volume of particles per column cross section

Vo (106 cm3/cm2)

single scattering albedo

wo (0.67 µm)

Water soluble * 0.005 0.176 1.090 3.050 0.96
Dust-like 0.500 17.60 1.090 7.364 0.69
Soot 0.0118 0.050 0.693 0.105 0.16
*both nucleation and accumulation modes

Deciding which category of aerosols are in the grid cells: choosing the aerosol model

For each grid box of pixels, an educated "guess" is made about what kind of aerosol is likely to be predominant at that location and time of the year. A decision tree is used that involves (1) OD0.47, (2) the ratio of the single scattering path radiances (Lp0.66/Lp0.47), and (3) information on global aerosol distribution. If the OD0.47cont < 0.15, then the Continental model should be retained. Otherwise, the Aerosol Model is chosen by following the decision tree described in the table below.
First Criterion (Q = scattering angle) Second Criterion (ratio between single scattering path radiances, Lp) Third Criterion (location, season) then use this MODEL
40° £ Q £ 150° Lp-red/Lp-blue > 0.90 globally Dust
       
40° £ Q £ 150° Lp-red/Lp-blue < 0.72 North America and Europe (year round) (100W-50E, 30N-70N) Industrial/Urban Aerosol
    South East Asia (year round) (105E-150E, 15N-45N) Industrial/Urban Aerosol
    Central America and Africa (May - Nov) (110W-50E, 0-30N) Industrial/Urban Aerosol
    South America and Africa (Dec-Apr) (110W-50E, 65S-0) Industrial/Urban Aerosol
    Central America and Africa (Dec-Apr) (110W-50E, 0-30N) Biomass Burning
    South America and Africa (May-Nov) (100W-50E, 65S-0) Biomass Burning
    Rest of the world Biomass Burning
150° £ Q £ 168° Lp-red/Lp-blue > (0.90) - (0.01)(Q - 150°) globally Dust
       
150° £ Q £ 168° Lp-red/Lp-blue < 0.72 North America and Europe (year round) (100W-50E, 30N-70N) Industrial/Urban Aerosol
    South East Asia (year round) (105E-150E, 15N-45N) Industrial/Urban Aerosol
    Central America and Africa (May - Nov) (110W-50E, 0-30N) Industrial/Urban Aerosol
    South America and Africa (Dec-Apr) (110W-50E, 65S-0) Industrial/Urban Aerosol
    Central America and Africa (Dec-Apr) (110W-50E, 0-30N) Biomass Burning
    South America and Africa (May-Nov) (100W-50E, 65S-0) Biomass Burning
    Rest of the world Biomass Burning

 

Once each grid cell has an aerosol model assigned to it, the surface reflectances (rs0.47 and rs0.66) plus the satellite observed reflectances (r*0.47 and r*0.66) are used again with the new model parameters to create a look-up table. From the look-up table, new values of optical depth (t0.47new, t0.66new) are obtained. Just like the Continental Model, the other models are based on real world measurements and assign values for the size and volume distributions, single scattering albedos, refractive indices, asymmetry factors and phase functions.

Aerosol Properties of the other Models

Aerosol category

mean particle radius of number distribution

mean particle radius of volume distribution

standard deviation of ln(r)

column volume of particles per column cross section

single scattering albedo

Biomass Burning

rg(µm)

rv(µm)

s

Vo(106cm3/cm2)

wo (0.67 µm)

Accumulation 0.61 0.130 0.500 -2.4 + 45t 0.90
Stratospheric 0.380 0.510 0.310 0.984 0.98
Coarse 1.0 - 1.3t 6.0 - 11.3t + 61t2 0.69 + 0.81t 2.4 - 6.3t + 37t2 0.84
Industrial/urban

rg(µm)

rv(µm)

s

Vo(106cm3/cm2)

wo (0.67 µm)

Accumulation 1 0.036 0.106 0.60 -2.0+70t-196t2+150t3 0.96
Accumulation 2 0.114 0.210 0.45 0.34-7.6t+80t2-63t3 0.97
Stratospheric 0.430 0.550 0.29 0.73 0.98
Sea Salt 0.990 1.300 0.30 -0.16+4.12t 0.92
Coarse 0.670 9.500 0.94 1.92 0.88

Dust

rg(µm)

rv(µm)

s

Vo(106cm3/cm2)

wo (0.67 µm)

1st mode 0.0010 0.0055 0.755 6.0x10-6 0.015
2nd mode 0.0218 1.230 1.160 1.0 0.95
3rd mode 6.2400 21.50 0.638 0.6 0.62

So, we have: t0.47new, t0.66new, t0.47cont and t0.66cont.

Now what?

The last step: the "correction"

t0.47 and t0.66 are corrected for the difference between the model chosen from the steps above and the preliminary t values obtained with the continental model. The following equation is used:

tlnew = (tlcont)(Plcontwlcont)/(Plnewwlnew)

Daily optical depth and mass concentration are stored on a resolution of 10x10 pixels (of 1 km nadir resolution). An interpolation is performed to obtain optical depth at 0.55 µm (as opposed to 0.47 or 0.66 µm).