Evolution of Ocean Color Atmospheric Correction: 1970-2005

Retrieval of water properties from satellite-borne imagers viewing oceans and coastal areas in the visible region of the spectrum requires removing the effect of the atmosphere, which contributes approximately 80–90% of the measured radiance over the open ocean in the blue spectral region. The Gordon and Wang algorithm originally developed for SeaWiFS (and used with other NASA sensors, e.g., MODIS) forms the basis for many atmospheric removal (correction) procedures. It was developed for application to imagery obtained over the open ocean (Case 1 waters), where the aerosol is usually non-absorbing, and is used operationally to process global data from SeaWiFS, MODIS and VIIRS. Here, I trace the evolution of this algorithm from early NASA aircraft experiments through the CZCS, OCTS, SeaWiFs, MERIS, and finally the MODIS sensors. Strategies to extend the algorithm to situations where the aerosol is strongly absorbing are examined. Its application to sensors with additional and unique capabilities is sketched. Problems associated with atmospheric correction in coastal waters are described.


Introduction
The literature on the optics of natural waters before 1970 suggested that the diffuse reflectance of the water contained information regarding the water itself and its constituents, i.e., mineral sediments as well as dissolved and suspended organic material (e.g., see Jerlov [1,2]). However, it was the pioneering work of Clarke et al. [3] that demonstrated that the upward spectral radiance measured above the sea contained information related to the concentration of chlorophyll a in the water. Chlorophyll a is the photosynthetic pigment contained in phytoplankton, which constitute the first link in the marine food chain. This direct observation led to the development of several aircraft imaging radiometers for viewing coastal (and ocean) waters and then to the first space-borne imager specifically designed for viewing the coasts (and ocean): the Coastal Zone Color Scanner (CZCS). Clarke et al. [3] also demonstrated that the information was degraded by the backscattering of the atmosphere between the surface and the sensor which they referred to as "air light". Regarding this "air light", at the end of their paper, they commented [3] "If such interference can be eliminated or identified and allowed for, spectroscopic procedures from aircraft (and perhaps from satellites) will be of great value in the rapid investigation of oceanic conditions, including conditions important for biological productivity." The effect of air light was further examined by Hovis and Leung [4] from altitudes as high as 16 km, where they observed that as much as 90% of the measured radiance actually resulted from scattering of sunlight between the surface and the sensor. This scattered light contained no information regarding the water. The removal of the air light radiance and the estimation of the radiance backscattered from within the water-the "water leaving radiance"-from that measured by the sensor is usually referred to as "atmospheric correction".
In this historical review I will describe the events leading up to the basic algorithm used for atmospheric correction of the modern U.S. sensors that followed the prototype CZCS, namely SeaWiFS, MODIS, and VIIRS (Gordon and Wang [5]). I devote much of the paper to this particular algorithm because of its widespread operational use, and refer to it as "GW94". I compare GW94 to the OCTS and MERIS algorithms, showing that they are essentially the same. I also show how (and why) the modifications to the algorithm were effected as a new understanding of the nature of aerosols over the oceans was obtained. A constraint on the original algorithm was that the aerosol exhibited only weak absorption compared to its scattering, and it often fails when strongly absorbing aerosols are present in the atmosphere, e.g., Saharan dust. Methods of dealing with strongly absorbing aerosols are described as well. Noting that the algorithm was developed for Case 1 waters, i.e., waters for which the optical properties all more-or-less covary with the concentration of chlorophyll a [6,7], and it often fails when applied to Case 2 waters (e.g., coastal waters), I will examine some efforts to extend the algorithm into the coastal regime; however, as this special issue of Remote Sensing is centered on atmospheric correction over such waters, the examination will be brief.
The review is divided into time periods that are quite naturally specific to individual sensors: 1970-1990CZCS operation;1990, SeaW-iFS/MODIS preparation (and operation); and 2000 to roughly 2005 (near the end of my active involvement in atmospheric correction development), algorithm enhancements and improvements. However, some work started before, but published after 2005, is described for completeness.
It is natural that in outlining the evolution of events in which one has participated, the result will reflect the author's perception of their import at the time. As such, there are likely contributions to the subject that, while important in their own right, had little impact on the development of the correction algorithm. I apologize for their omission.

1970-1980: Algorithm Conception and CZCS Development
It is important to note that up to and through the time period in which the CZCS was designed and fabricated there were no actual algorithms relating the color of the water (i.e., the spectrum of upwelling radiance just beneath the surface) to water constituents. The first such algorithm was that of Morel and Prieur [7] in which the ratio of the subsurface diffuse reflectance at 440 to that at 550 nm was related experimentally and theoretically to the "pigment concentration" (C P , the concentration of chlorophyll a and its degradation product phaeophytin a) in the water. This was closely followed by similar algorithms, based on a data set from coastal and ocean waters around the U.S. [8,9]. In addition, it is clear that in the design of CZCS little or no attention was paid to the problem of removing, or in some way accounting for, the masking effects of the air light. After its formation, the CZCS Experiment Team (Table 1) was in effect handed the specifications for the instrument and asked to "do what you can with this to make it work". surface and the sensor. Thus, one of the principal tasks set forth by the CZCS Experiment Team was to develop a method to account for the atmospheric effects. The other was to establish a quantitative relationship between the color of the water and the concentration of chlorophyll a. I volunteered to take the lead in the former.
Believing that aerosols presented the most significant problem for atmospheric correction, and applying the adage "one man's signal is another man's noise", I started by looking at efforts to retrieve aerosol concentrations from satellite-measured radiance, principally from the LANDSAT MSS sensor. Griggs [10] compared LANDSAT (then ERTS-1) radiances to the aerosol optical depth τ a (at 550 nm) over water surfaces and found a near-linear relationship between the two for the small range of τ a observed in the study. This suggested that given the aerosol optical depth (and possibly other parameters associated with the aerosol) and the known scattering by the air itself, it might be a simple matter to estimate the radiance at the sensor. Certainly this should be possible if a full multiple scattering solution to the radiative transfer equation were developed; however, with the enormous amount of data from the CZCS (one CZCS two-minute "scene" contained approximately 1.6 × 10 6 pixels with 4 usable wavelengths), full radiative transfer calculations at each pixel were clearly impossible. In addition, the aerosol properties required for such calculations are not known a priori. Thus, a simpler method was required to determine the atmospheric effect. The fact that the radiances for LANDSAT seemed to relate linearly to aerosol optical depth suggested that single scattering might be a usable approximation for the radiative transfer. Briefly, the (scalar) radiative transfer equation for the radiance L in a plane parallel medium illuminated from the top (τ = 0) by a beam of irradiance F 0 propagating in a directionξ 0 iŝ ξ •ê z ∂ ∂τ L(τ,ξ, λ) = −L(τ,ξ, λ) Allξ P(τ,ξ →ξ, λ)L(τ,ξ , λ) dΩ(ξ ) whereê z is a unit vector directed into the medium at the top,ξ is a unit vector in the direction of propagation of the radiance,"•" indicates the scalar product, τ is the optical depth, λ is the wavelength, ω 0 is the single-scattering albedo (scattering coefficient ÷ extinction coefficient), P is the scattering phase function, and dΩ(ξ ) is a small solid angle around the directionξ . The usual single-scattering solution to this equation involves simply leaving out the integral term and solving the remaining inhomogeneous differential equation. If the atmosphere is modeled as a plane parallel homogeneous slab of thickness τ 1 , the radiance exiting the top of the atmosphere (when no radiance enters either the top or bottom of the atmosphere) is where u =ξ •ê z and u 0 =ξ 0 •ê z (note that u < 0 for radiance exiting the top of the slab), and φ and φ 0 are the azimuth angles forξ andξ 0 , respectively. This is often referred to as the single-scattering radiance reflected from the slab; however, as it contains terms of order higher than the first in τ 1 , i.e., the exponential term, it is not really single scattering. True single scattering is obtained by expanding the exponential and retaining the lowest order term in τ 1 . The result is Often it is desirable to work with a normalized radiance ρ (referred to as "reflectance" here) rather than radiance. The reflectance is defined through ρ = πL/F 0 u 0 . (Note that ρ is dimensionless; however, the defining equation yields units of inverse solid angle, sr −1 . Since a solid angle is just the ratio of two areas, as a plane angle is the ratio of two lengths, and hence dimensionless, I prefer to leave ρ unit less; however, if the reader is uncomfortable with my convention, then the unit sr −1 should be attached to all reflectances.) Then, If we modify the model to account for the presence of a flat specularly reflecting water surface at τ 1 (but ignore the radiance exiting the water from below), it is easy to show that Equation (3) becomes where cos Θ ± = ∓uu 0 + 1 − u 2 )(1 − u 2 0 ) cos(φ − φ 0 ), and r F (x) is the Fresnel reflectance of the air-water interface for an incident angle (in the air) of arccos(x). How accurate is this approximation? Figure 1 shows the reflectance of a slab with Rayleigh scattering characteristic of Earth's atmosphere at 412 nm (Rayleigh optical thickness of 0.312) as a function of both the solar zenith angle θ 0 and the exiting angle θ ex in the perpendicular plane (φ − φ 0 = 90 deg). The angle the radiance exits the top of the atmosphere is θ ex , where θ ex = π − θ = cos −1 (|u|). The solar zenith angles are θ 0 = 20 deg (filed circles), 40 deg (filled squares), and 60 deg (filled triangles). (Right Panel) Error (%) in ρ r computed using Equation (4). Positive error indicates that the single-scattering ρ r is too high. (Taken from Gordon [11]).
The accuracy of Equation (4) in reproducing the results of exact (scalar) theory is truly remarkable. It suggested that this equation could be used to estimate the Rayleighscattering contribution to the radiance at the top of the atmosphere with an error of no more than about ±2% (note that the blue band on CZCS has λ = 443 nm, where the optical thickness is about 30% lower than that in Figure 1). Figure 2 shows a similar computation when aerosols are included in the slab. In this case, the aerosols are characteristic of a maritime aerosol with optical thickness of 0.2 at 412 nm. The computation was effected by computing the contributions from Rayleigh scattering (ρ r ) and aerosol scattering (ρ a ) separately and then adding them (ρ Total = ρ r + ρ a ).  (4). Positive error indicates that the single-scattering reflectance is too high. (Taken from [11]).
The error is now larger than in Figure 1, but so is the total optical thickness (>0.5). The larger error for small θ 0 and small θ ex owes to the large forward scattering by the aerosol, which is reflected from the lower surface. This suffers too little attenuation upon propagating to the upper boundary (zero attenuation in fact). In reality, this particular error in radiance is not very important as the Sun's glitter pattern for a roughened (by the wind) surface would likely overwhelm the radiance reflected from the atmosphere. Thus, it appeared that, within an error of only a few percent, one could write ρ T = ρ r + ρ a for the normalized radiance exiting the atmosphere, where ρ r and ρ a are computed individually using Equation (4). Note, this agrees with the observation [10] that LANDSAT radiances varied linearly with aerosol concentration (i.e., aerosol optical depth). It was taken as the working equation for the analysis of CZCS data [12], i.e., where ρ T is the total reflectance measured at the sensor, ρ w is the radiance exiting the water converted to reflectance and t is its diffuse transmittance [13] (see below). It is easy to underestimate the importance of the observation that single-scattering could be used to provide a reasonably accurate estimate of the radiance at the sensor that was due to the atmosphere. One need only consider the computing power that was available at that time. When CZCS was launched (November 1978) 1 MHz workstations with 1MB of memory were about a decade in the future. In fact some of the track line comparisons of ship measurements and satellite estimates of the concentration of chlorophyll we published (e.g., Gordon et al. [14]) were carried out by the author on a desktop computer with a single line display and 64 kB of memory (Hewlett Packard 9825). Thus, any complex computation of the atmosphere's influence on the sensor's radiance was out of the question. Mainframe computers at that time were too expensive to be dedicated to image processing. Actually, it was not until the early to mid 80s that even the single-scattering computation of ρ r could be effected on a pixel-by-pixel basis in an image processing environment (this application employed supermini computers such as the VAX 11/780).
There are many parameters in Equations (4) and (5) that were unknown in the application: ω 0 , P, and τ for the aerosol in Equation (4); and ρ a , t, and ρ w in Equation (5). These had to be estimated/determined from the imagery itself. In order to make up for the unknown quantities, several assumptions were required. The CZCS had four bands with sufficient sensitivity to be useful for ocean color. Their band centers were at 443, 520, 550, and 670 nm. The first assumption was that ρ w (670), the "water-leaving reflectance" in the red band, was negligible, i.e., ρ w (670) = 0. This is a good assumption for Case 1 waters with chlorophyll concentrations <0.5-1 mg/m 3 (the oceanic chlorophyll concentration has a broad distribution from 0.06 to about 0.3 mg/m 3 [15] with a global median around 0.13 mg/m 3 and higher values at higher latitudes). Then Equation (5) provided ρ a (670), as ρ r (670) could be readily computed from Equation (4). Note that in Equation (4) the angles are independent of wavelength, so that only ω 0 , P, and τ vary with λ. It was then assumed that ω 0 and P for the aerosol are independent of wavelength, so only the aerosol optical thickness τ a was spectrally dependent. Since ρ a is directly proportional to τ a , the aerosol optical depth, one has which shows that the water contribution to the total reflectance can be retrieved from Note that a term t(670)ρ w (670) is included on the right-hand side above even though it was assumed ρ w (670) = 0. This is for later use in situations in which ρ w (670) = 0.
As ρ r can be computed directly and ρ w (670) = 0, the only unknown on the right hand side of Equation (7) is the spectral variation of τ a . For many environmental aerosols where the (usually) positive number n is called the Ängstrom exponent. For a maritime aerosol, n ≈ 0. Thus, at this point all that was required to find tρ w in the blue and green bands was n. In more modern notation, the multiplier of ρ T (670) − ρ r (670) is called (λ, 670), and Equation (7) is usually written and for completeness, in the single-scattering approximation (λ 1 , λ 2 ) is given by where ω 0 , P, and τ refer to the associated aerosol properties. Once tρ w is determined, it remained to estimate t. In the absence of scattering in the atmosphere, the radiance that exited the water from a pixel being viewed by the sensor can arrive at the sensor by a single path: precisely in the direction of the sensor. This radiance is attenuated by a factor exp[−τ/|u|]. However, in the presence of atmospheric scattering radiance leaving the water in a different direction (i.e., from a different pixel) can also reach the sensor by being scattered into it. The later radiance depends on the angular distribution of radiance leaving the surface. If this radiance is totally diffuse (independent of the exiting direction) and the water is horizontally homogeneous, then for a pure Rayleigh-scattering atmosphere, and in the single-scattering approximation, the diffuse transmittance is given by (see [11]) where the "*" is to remind us that this formula is restricted to horizontal homogeneity of the water and to a totally diffuse ρ w at the surface. (The water-leaving radiance distribution is far from diffuse; however, the upward radiance just below the surface of the water is only weakly dependent on direction. My preference then was to compute t * assuming the upwelling radiance beneath the surface is totally diffuse, as in [13,14].) The dependence of t * r on r F owes to photons exiting the surface, backscattering toward the surface, and then reflecting from the surface in the direction of the sensor. The aerosol contribution to t is much smaller than the Rayleigh contribution because of the strong forward scattering exhibited by aerosols. It usually results in <1% change in t and was ignored. Thus, letting t = t * r , all that remained for determining ρ w (λ) from ρ T (λ) was the determination of (λ, 670).
In the initial validation of the CZCS-derived pigment concentration [16], (λ, 670) was determined by direct ship-based measurements of ρ w (λ) at the surface (i.e., at one pixel in the image). The derived -values were then used throughout the entire image. The uniformity of was justified by examination of the image provided in Figure 3. The image of (520, 670) was derived by using an estimate for L w (520) that was valid for the low-chlorophyll concentrations found in the interior of the Gulf of Mexico at the time of the image, Fall, 1978 (see later discussions of ρ w (λ) as a function of the chlorophyll concentration).   Figure 3. Note that varies by ±5% while L a (670) varies by about ±20%, and some of the variation in is due to variation in L w (520), albeit small. Thus, it was taken as a working hypothesis at that time that (λ, 670) was constant over an image (or sub-image), even in the presence of variations in the aerosol concentration. This is nearly equivalent to assuming that the aerosol size frequency distribution and refractive index are uniform over an image, even when the aerosol concentration varies (nearly because will undergo some variation because of its dependence on the Sun-viewing geometry). In sum, for atmospheric correction of the initial imagery from CZCS [16,17], the principle assumptions were: (1) the single-scattering approximation was sufficient to estimate ρ r ; (2) L w (670) = 0; (3) the parameter (λ, 670) was approximately constant over a CZCS scene; and (4) the upward radiance just beneath the water surface was totally diffuse.
Ship measurement of L w (λ) contemporaneous with the satellite overpass was carried out to determine (λ, 670), and these values were used for the entire image. The results suggested that the pigment concentration (C P ) calculated using the radiance ratios retrieved from the CZCS, i.e., was valid to within a factor of two. It should be noted that these relationships for C P [8,9,16] were developed from the sparce in situ measurements without regard to the Case 1-Case 2 distinction due to the paucity of surface measurements at that time. It is also important to note that within the literature, relationships such as these were dynamic in the sense that the coefficients, i.e., 0.5 and 1.3 in the first, were subject to variation as new in situ data were obtained, and as stratification of in-situ data between Case 1 and Case 2 waters became possible.

1980-1990: Algorithm Improvement and CZCS Operation
The basic concept of the atmospheric correction algorithm was to use a spectral band for which the radiance exiting the water is negligible to estimate the aerosol contribution and, by some means, scale it to provide an estimate in the other spectral bands. After the launch of CZCS and the processing of the initial data, it was clear that this algorithm concept was viable, and the effort then turned to improvements to the algorithm and to its implementation.

Ozone Absorption
The first improvement was to account for the absorption by Ozone in the stratosphere. Note that in Equations (5) and (9), all of the ρ's (or L's if we revert to radiances) are defined in the absence of Ozone. The Ozone correction was effected by assuming the Ozone was all in a non-scattering layer just above the slab atmosphere. Then every photon that arrives at the sensor had undergone two trips through the layer and therefore L T was reduced by a factor where τ Oz is the Ozone optical thickness, and M is the air mass The values of τ Oz for a concentration of 343 DU (the U.S. standard atmosphere) are: τ Oz (443) ≈ 0.; τ Oz (520) ≈ 0.02; τ Oz (550) ≈ 0.03.; and τ Oz (670) ≈ 0.02. Note however that the Ozone concentration varies from about 200-400 DU, so even though the absorption is small, it is highly variable, i.e., at 550 nm 0.932 < f Oz < 0.965 for 200 < DU < 400 and M = 2. Accounting for Ozone absorption, the measured radiance exiting the atmosphere is and in the presence of Ozone one must operate the algorithm with the "corrected" L T value above, to derive the value of L w extant without Ozone. The actual L w is related to that measured by L with Oz as the light that produced the water-leaving radiance traversed the Ozone layer only once. This is the value that should be compared with surface measurements of the waterleaving radiance.

Clear Water Radiance Concept
After deriving L with Oz w it is obviously useful to be able to compare measurements from one day to another or between different locations. Generally, in such comparisons, the solar zenith angles and the viewing angles will differ for each measurement. Gordon and Clark [18] attempted to normalize the recovered radiance for these differences by defining the normalized water-leaving radiance, [L w ] N as where a ⊕ is the Earth-Sun distance in AU (Gordon and Clark assumed a ⊕ = 1). Note if θ 0 = 0, τ r = τ Oz = 0, and a ⊕ = 1, then so if the Sun is at the zenith, and there is no atmosphere, [L w (λ)] N would be the radiance observed at the sensor. This is actually not quite correct, as without the atmosphere the downward radiance distribution at the water surface would simply be that of a Sun in a black sky, rather than attenuated direct sunlight plus diffuse sky radiance. In addition, [L w (λ)] N as defined in Equation (13) is based on the assumptions that (1) the water surface is flat and (2) the upwelling radiance just beneath the surface (i.e., within the water) is totally diffuse. Neither of these are extant for a real waterbody, so the statement is only approximate. However, Equation (13) defines [L w (λ)] N , and it is useful for comparing measurements made at different solar zenith angles. The normalized water-leaving radiance can easily be converted to normalized waterleaving reflectance [ρ w ] N through where F 0 is the extraterrestrial solar irradiance for a ⊕ = 1, i.e., F 0 = a 2 ⊕ F 0 . Then Equation (13) for L with Oz As an aside, it should be noted that modern algorithms use "remote sensing reflectance" (R rs = L w /E + d , where E + d is the downward irradiance just above the sea surface) rather than [ρ w ] N [19,20]. Fortunately, there is a simple relationship between the two: The importance of [L w ] N and [ρ w ] N , was the observation by Gordon and Clark [18] that for pigment concentrations (C P ) less than about 0.25 mg/m 3 [L w ] N and [ρ w ] N were approximately constant, i.e., independent of C P in the green spectral region (e.g., 520-550 nm). This is shown in Figure 5, which compares the variation in [ρ w ] N at 443 nm with that at 555 nm, derived from water-leaving radiances measured in Case 1 waters using in-situ radiometers [18]. These results suggested that for C P < 0.  This observation, referred to as the "clear water radiance concept," could be used to effect atmospheric correction in some situations. For example, if a region of an image could be classified as "clear water," the radiances at 550 and 670 nm could be used to estimate (550, 670) in that region, since L w at these wavelengths would be known. From this estimate, (λ, 670) could be inferred, e.g., by assuming that (λ, 670) ∝ (670/λ) n . The values of (λ, 670), as estimated, could then be used to correct the entire scene. This was the procedure carried out by Gordon et al. [14] in the region of the Gulf Stream in the Atlantic Ocean off the U.S. East Coast. In that validation study a warm core ring was observed in several images and the low pigment concentration within the ring constituted the "clear water" region.

Difficulties with, and Alternatives to, the Clear Water Requirement
There are obvious difficulties applying this procedure routinely, e.g., (1) the image of interest may contain no "clear water," or C P may not be small enough to take ρ w = 0 at 670 nm which can happen even in Case 1 waters at high C P , and (2) the ε's may vary over the image because of variations in aerosol type (or as a result of multiple scattering). Nevertheless, the algorithm, crude as it was, did demonstrate the feasibility of measuring pigment concentrations from space.
Attempts to remedy the obvious problems focussed on two different approaches: first, the use of empirical relationships between the L w 's to provide a forth equation to go with the three represented in Equation (7), allowing retrieval of the four ρ w 's; and second, the development of bio-optical models that relate the variation of L w (λ) to the pigment concentration intended to be valid for Case 1 waters. Briefly, in the first method Smith and Wilson [22] tried to overcome the requirement that L w (670) = 0 by employing an empirical relationship (which included both Case 1 and Case 2 waters): where L u is the upwelling radiance just beneath the water surface L w (λ) = t F L u (λ)/m 2 w , where m w is the refractive index of water and t F is the Fresnel transmittance of the air-water interface from the water side at the appropriate incidence angle . Equation (15) was used in an iterative manner in the following way. It was assumed that (λ, 670) = (670/λ) n and that n could be obtained in some alternate way, e.g., using the "clear water concept". In the first iteration, it was assumed that L w (670) = 0 and the algorithm returns L w (λ) for the other three bands. These were then used in Equation (15) to estimate a corrected value for L u (670). The corrected value of L u (670) provided a new value for L w (670), and operating the algorithm with this yielded new values for the three other L w (λ)'s. In this manner, the basic algorithm described in Section 2 could be operated in waters for which L u (670) = 0 as long as n could be estimated. Schematically, after converting all radiances to reflectances and with the superscript representing the iteration number starting with i = 1 (ρ (1) w (670) = 0), the procedure (with n fixed) was ending when ρ (i) w (670) became stable. Bricaud and Morel [23] developed a similar algorithm, however, they used measurements restricted to Case 1 waters alone as opposed to Equation (15): This relationship was used in a manner similar to that of Smith and Wilson to derive L w (670) given n, and then C P given empirical algorithms relating radiance ratios similar to those in Equation (12). In addition they also developed a bio-optical model relating ρ w (λ) at each of the CZCS spectral bands to the pigment concentration (C P ). This model was used to estimate n for image pixels that might not qualify as "clear water" by virtue of C P > 0.25 mg/m 3 . The estimate was obtained schematically as follows: Briefly, one starts with n = n (1) = 0 and enters the scheme in Equation (16), which returns values for ρ (1) w (λ) from which C (1) P is estimated. This pigment concentration is then used in Equation (18), where "Model" refers to the bio-optical model, to obtain a new estimate of the water-leaving reflectances ρ (2) w (λ) and hence ρ (2) a (λ) from which n (2) is obtained. The procedure given in Equation (16) is then operated with n (2) resulting in ρ (2) w (λ) and C (2) P , etc. In the last step in Equation (18), where Equation (12), is utilized, the equation actually used was a revised version of Equation (12) derived with more data and restricted to Case 1 waters.
Application of this algorithm for determining both n and C P showed that there was an upper limit in Case 1 waters, C P ∼1-2 mg/m 3 , beyond which there was incomplete decoupling between the ocean and atmospheric properties, i.e., there was a strong correlation between variations in n and C P , that was not seen at lower values of C P . At low values of C P there appeared spatial variations (noise) in n even in regions that satisfied the requirement for "clear water," suggesting that the assumption of a constant n for a relatively large region of an image should be used with caution. Such spatial variations in (550, 670) are also seen in Figure 3. In a later paper Andre and Morel [24] suggested that relationships provided by the bio-optical model should be used in place of Equations (12), (15) and (17), to assure internal consistency. They also suggested that with such a change, Equations (16) and (18) could be used to estimate n and C P on a pixel-by-pixel basis for C P 2 mg/m 3 .
These algorithms appear to have been rarely used in CZCS processing in the 1980's except for various demonstration studies. However, the basic approach was incorporated into the routine processing of SeaWiFS data (see Section 5.1), and in various reprocessing of CZCS data in the mid 2000's by the OBPG at NASA/GSFC.

Compensation for Radiometric Calibration Error (Vicarious Calibration)
The top-of-atmosphere determination of ρ T in Equation (5) involves measurement of L T and division by F 0 , both of which quantities have associated uncertainties (errors). It is easy to show that if these quantities have uncorrelated errors of the order of 5%, then the resulting L w 's are useless (they may even be negative in some cases). In fact, the accuracy required for the CZCS in the blue region of the spectrum far exceeded the capabilities of even carefully controlled laboratory calibration at that time. The approach taken to remedy this problem was to adjust the calibration constants of the CZCS (i.e., the constants relating digital counts to radiance in radiometric units) to bring the CZCS-estimated water-leaving radiances into confluence with the corresponding surface-measured values [25][26][27]. This procedure was referred to as "vicarious calibration" (VC) [28]. Briefly, and ignoring Ozone corrections, using measurements of ρ w (λ) the algorithm was used to estimate (λ, 670): It was also assumed that there was no error in ρ T (670). Then ρ T (λ) is varied for each band until a "smooth" variation in (λ, 670) is obtained. At that time a "smooth" variation was subjective and deemed to be something like (λ, 670) ≈ (670/λ) n , where 0 < n < 0.5 for a maritime atmosphere. Although it might not seem so, this is a very sensitive method of finding a good "vicarious calibration" (see [29,30] for a complete analysis of the accuracy of this procedure as applied to SeaWiFS). Using this procedure with CZCS resulted in good agreement between surface and retrieved values of ρ w for the times of the year near the time of the vicarious calibration, but for earlier times errors were increasing with increasing time difference between calibration and application. Applying the vicarious calibration procedure multiple times from launch showed that the radiometric sensitivity of the CZCS blue band exhibited a rapid decrease in sensitivity with time-a 20% decrease in a little over 4 years [27]. An account of this decrease was included in all CZCS processing. Although varying in detail from that here, vicarious calibration procedures have been employed with all ocean color scanners since CZCS, usually after the long-term drift of the radiometric responses has been assessed in other ways, e.g., see [31].

Computation of ρ r including Polarization and Multiple Scattering
Although the single-scattering approximation provides a good estimate to ρ r , the error can still reach as much as 4%. To reduce the error an exact computation of ρ r using vector radiative transfer theory ("vector" implies that it includes polarization of the radiance, as opposed to "scalar" theory in which polarization is ignored) was carried out [32] providing "exact" values for this quantity as a function of both the solar position and the viewing direction, i.e., ρ r (λ; u 0 , φ 0 ; u, φ). The results were in the form of look-up-tables from which the reflectance for any geometry could be calculated. In addition, as ρ r is a function of the atmospheric pressure (τ r ∝ P 0 , where P 0 is the surface pressure), a simple method was developed (by trial and error) to allow the accurate computation of ρ r given small changes in pressure, through where τ r is the optical thickness at mean atmospheric surface pressure P 0 , and τ r = (P /P 0 )τ r , with P being the actual surface pressure.

Effect of Aerosol Multiple Scattering and Rayleigh-Aerosol Interaction
Although the correction algorithm based on the single-scattering approximation was quite effective, Figure 2 clearly shows that there can be considerable error in the computation of ρ r + ρ a . In fact, a significant error remains even when ρ r is computed as in Section 3.5. Deschamps et al. [33] studied the multiple scattering effects through the following computations for an atmosphere with a totally absorbing lower boundary. First, (1) compute ρ r (including multiple scattering, but in the absence of aerosols) with an optical thickness τ r ; (2) compute ρ a (including multiple scattering, but in the absence of the air, i.e., Rayleigh scattering molecules) with an optical thickness τ a ; and (3) compute ρ T for an atmosphere composed of both molecules (τ r ) and aerosols (τ a ). Then, form (3) − [(1) + (2)], which in the single scattering approximation is zero. Finally, define the remainder ρ ra ≡ (3) − [(1) + (2)]. It is the error in the single-scattering approximation due to photons that scatter one or more times from both aerosol particles and molecules, and depends on all of the parameters of the computation: Sun-viewing geometry; τ r and τ a ; the aerosol phase function, etc. Deschamps et al. [33] showed that ρ ra was close to the limit of detectability of CZCS, but would have to be addressed for the sensors with higher radiometric sensitivity that were being proposed at that time. Thus, it was clear that a more correct formulation of Equation (5) would have to be used when more sensitive instruments were developed: where it is understood that the term ρ r in this equation is computed from an exact solution of the radiative transfer equation (including polarization) for an aerosol-free atmosphere, and ρ A = ρ a + ρ ra . Gordon and Castaño [34] approached the problem of assessing the effects of multiple scattering in a different manner. They used the second form of Equation (20) and studied the variation of over CZCS scan lines for a horizontally homogeneous aerosol. They also assumed that the aerosol phase function was independent of wavelength. Their results suggested that the error in atmospheric correction could be reduced close to the quantization interval of CZCS if the algorithm as originally developed, i.e., Equations (5)-(9), was used with ρ A replacing ρ a and replacing , and at each pixel the resulting value of (443, 670), i.e., that provided using Equation (8), was reduced by 5%. Note that these changes could be easily accommodated within the iterative algorithms in Section 3.3.

Algorithm Modifications for Processing the CZCS Global Data Set
After it became clear that the CZCS retrievals of pigments agreed well with in situ measurements in the open ocean, acquisition of a global data set became a priority of the Experiment Team. This required that a significant fraction of the CZCS operation time (2 h/day) be dedicated to global data collection.
In 1986, after the CZCS ceased operation, NASA begin a project to process the entire data set [35], a total of approximately 68,000 two-minute scenes (1968 pixels/line × 970 lines/scene). The spatial resolution of the CZCS was approximately 0.825 km at nadir. These data were subsampled, i.e., samples were taken from every 4th line and every 4th pixel for a pseudo spatial resolution of about 4 km. The data were atmospherically corrected using Equations (5)-(9) with the Gordon and Castaño [34] modifications described in Section 3.6 and with ρ r computed as described in Section 3.5. The Ozone correction (Section 3.1) was effected as well, using Ozone concentration data from the TOMS Ozone sensor that was carried on the NIMBUS-7 satellite along with the CZCS. The water-leaving radiance at 670 nm was taken to be null, as the alternatives to this assumption (Section 3.3) were too costly in computational resources for such an undertaking. Owing to that fact that maritime aerosols were thought to exhibit little spectral variation in extinction (Shettle and Fenn [36], referred to as "SF79"), the exponent n in Equation (8)  The effect of the aerosol on the diffuse transmittance was ignored and t was replaced by t * r in Equation (9). Calibration variations over the CZCS mission were based on the assessment of Evans and Gordon [37]. The results from processing the global data set from radiance to pigment concentration helped demonstrate the value of ocean color data to the broader oceanographic community and paved the way for the new sensors, e.g., SeaWiFS and MODIS [38].
Summarizing, the significant atmospheric correction activities of the period 1980 to 1990 consisted of the following: • accounting for absorbing gases, e.g., Ozone; • development of clear water radiance (CWR) concept; • using CWR for determination; • using surface measurements to effect vicarious calibration; • including exact computations of ρ r (including polarization); • assessing the error in neglecting multiple scattering; and • applying the concomitant improvements/understanding to the global data set.

1990-2000: Preparation for Improved Instruments
Even before the launch of CZCS it was understood that, although the sensor was more-or-less optimized for observing the principal variations in ocean color, i.e., those resulting from variations in the concentration of chlorophyll-like phytoplankton pigments, it was not well designed for atmospheric correction [39]. That is, there were no spectral bands for which the sensor recorded only radiance resulting from atmospheric and surface interactions. In fact, the radiance exiting in the red band is negligible only in low to moderate chlorophyll waters, and hardly ever near the coasts. However, for bands in the near infrared (NIR, 700-1000 nm), it is usually the case that L w ≈ 0, except in regions with high suspended particle load (non-living or living). The CZCS had a wide (100 nm) spectral band at 750 nm that had a radiometric sensitivity similar to a corresponding LANDSAT band, but that radiometric sensitivity was too low to be of value for atmospheric correction. The opportunity of increasing the sensitivity of this band electronically was presented to the author after CZCS had been integrated with the NIMBUS satellite; however, as this involved considerably risk, the decision was made not to effect it. A working group [39] proposed in the summer of 1978 that ocean color sensors following CZCS have at least two spectral bands in the NIR with which the effect of the atmosphere on the radiance at the sensor could be assessed and extrapolated to the visible bands. Such spectral bands have been placed on all post-CZCS sensors, with the MODIS and VIIRS (and some proposed) sensors having additional bands that are still further into the infrared (the short-wave infrared, SWIR), although their presence was justified for other purposes.

Required Accuracy for Atmospheric Correction
In the case of CZCS, where there were no bio-optical algorithms to use as a guide, the goal was to perform atmospheric correction as best as one could (given the meager computational resources). The goal of the second-generation SeaWiFS sensor was stated (SeaWiFS Prelaunch Report V1 [40]) "To achieve ... water-leaving radiances to within 5% absolute, and chlorophyll a concentration to within 35% over the range of 0.05-50.0 mg m −3 ". However, the description of this goal was incomplete. For example, at 670 nm L w ≈ 0, and retrieving a very small radiance to within 5% is clearly impossible. The correctly stated goal was to retrieve the water-leaving radiance at 443 nm with an accuracy of 5% in oligotrophic waters. But even this more complete description of the goal still depends on the definition of "oligotrophic." The author took that to mean a pigment concentration of about 0.05 mg/m 3 . Figure 5 shows that [ρ w (443)] N ≈ 0.04 for this pigment concentration, so I took the goal of atmospheric correction to be ∆[ρ w (443)] N < 0.002. Linear regression of the log-transformed data in Figure 5 provides and when C P ≤ 0.6 mg/m 3 the standard error of estimate in C P , using the second form for regression, is about 27%. Using this we find so for the above-stated goal of atmospheric correction the atmospherically-induced error would correspond to an error in C P of approximately 18% at C P = 0.1 mg/m 3 . Combining this with the error in the bio-optical algorithm, in a root-sum-of-squares (RSS) addition of the two, predicts a total error in C P of approximately 32% for C P ≈ 0.1 mg/m 3 ; however, for C P = 0.6 mg/m 3 the atmospheric correction error approximately doubles and the combined RSS error increases to about 45%. In order to achieve an uncertainty of less than 35% in an RSS sense over the range 0.05 ≤ C P ≤ 0.6 mg/m 3 , one would need ∆[ρ w (443)] N ≤ 0.0012. With the final SeaWiFS algorithm, this goal was achievable sometimes, but not always; however the goal of ∆[ρ w (443)] N ≤ 0.002 appeared to be achievable under most circumstances as long as the aerosol is non-absorbing or weakly absorbing. Antoine and Morel [41] obtained similar criteria using the requirement that a sensor (MERIS in their case) be capable of distinguishing 30 equally spaced classes of log 10 C (in their paper they used the chlorophyll a concentration, C, rather than C P ) in the range 0.03 ≤ C ≤ 30 mg/m 3 .

An Extended Single-Scattering Algorithm
It was important to implement the equations for atmospheric correction, in the single scattering approximation, including these additional bands in the NIR as an initial substitute for the yet-to-be-developed multiple scattering algorithm, which was expected to be orders of magnitude more computer intensive.
Consider two spectral bands in the NIR positioned at λ s and λ , where "s" and " " stand for "short" and "long," respectively. For SeaWiFS, λ s = 765 nm and λ = 865 nm. Thus, since ρ w (λ) = 0 for λ = λ s and λ = λ , the equivalent of Equation (9) is where ω a , P a , and τ a refer to the associated aerosol optical properties. Writing Equation (23) for λ = λ s , we have which would be a known quantity. It should be noted that in Equations (23) and (25), ρ r (λ) is provided in the form of look-up-tables (LUTS) prepared from exact computations (including polarization) of the radiance exiting a purely Rayleigh scattering atmosphere bounded by a totally absorbing ocean with a flat Fresnel reflecting surface, i.e., the multiplescattering analog to Equation (4). In the case of SeaWiFS, an estimate of the atmospheric pressure was to be available and would modify ρ r in a manner prescribed by Equation (19). The basic idea of the algorithm was to estimate (λ, λ ) given (λ s , λ ). In the case of a non-absorbing aerosol for which the phase function is independent of wavelength, this is particularly simple (in this single-scattering approximation): (λ s , λ ) = τ a (λ s )/τ a (λ ) or (λ, λ ) = τ a (λ)/τ a (λ ), and the spectral variation of the aerosol optical depth alone deter-mines the extrapolation into the visible. An example we considered earlier was when the aerosol satisfies Ängstrom's law (τ a (λ) ∝ λ −n ), in which case (λ s , λ ) = (λ /λ s ) n determines n, and then (λ, λ ) = (λ /λ) n . Unfortunately, the phase function of most aerosols depends on λ, so this simple example (used here only as an illustration) is not realized in practice, although it was universally assumed in the earlier algorithms (e.g., [14,22,23]). Wang and Gordon [42] studied the efficacy of Equations (23)-(25) when the realistic aerosol models developed by Shettle and Fenn [36] were used to characterize the aerosol. They found that for these aerosol models where k is a constant that depends on the particular model and Sun-viewing geometry. It can be determined via where (λ s , λ ) is found from Equation (25). Using simulations, they concluded that (1) if the aerosol was non-absorbing (or weakly absorbing, e.g., maritime), (4) if the resulting value of (443, λ ) is reduced by 4.6%, then ∆[ρ w (443)] N is usually less than 0.002. The required reduction in (443, λ ) by 4.6% is consistent with the approximate correction of Gordon and Castaño (1987) [34] for the effects of multiple scattering. The extended single-scattering algorithm was to be a "fall-back" algorithm if its multiple-scattering counterpart proved too costly in terms of computational resources.

The SeaWiFS Multiple Scattering Atmospheric Correction Algorithm
In 1989, the author and M. Wang began working on an extension of the algorithm to include the effects of multiple scattering in a less ad hoc manner than Gordon and Castaño [34]. They started with the Deschamps et al. [33] observation that there was significant contribution to ρ T from photons that scatter from both molecules and aerosol particles, i.e., ρ ra in Equation (20) or the replacement of ρ a by ρ A . To avoid confusion regarding what refers to single scattering and what refers to multiple scattering, we now change the notation slightly. Let the single-scattering reflectance exiting the atmosphere after scattering by the aerosol henceforth be indicated by ρ as , where with the subscript "a" indicating "aerosol", and subscript "s" indicating "single scattering." Note that, in a given Sun-viewing geometry, ρ as is directly proportional to the aerosol optical thickness τ a . The success of the single-scattering algorithm at low values of τ a suggested that it should be mimicked as closely as possible in a multiple-scattering version. That was the path followed. The only link between the aerosol and its physical properties that can be determined from ρ T (λ) are the values of ρ A (λ s ) and ρ A (λ ), both of which contain the effects of multiple scattering. Because ρ as is proportional to τ a , it is clear that for a given wavelength, Sun-viewing geometry, and aerosol model (say the jth), one can write ρ A as a function of ρ as . An example is where a (j) , b (j) and c (j) are constants that depend on the model and the Sun-viewing geometry. In the single-scattering approximation, a (j) = 1 and b (j) = c (j) · · · = 0. For now we use the following notation for the relationship designated in Equation (27), where the big curly brackets mean ρ A (λ) for the jth model given as a function of ρ as for the jth model. (Later I will simplify this notation to ρ A = g(ρ as ) or ρ A = g (j) (ρ as ) for the jth model, where the g's are the polynomials in ρ as .) Equation (28) can be inverted to yield ρ as as a function of ρ A for the jth model, i.e., ρ (j) The relationships in Equations (28) and (29) must be established by way of rigorous solutions to the radiative transfer equation. These in turn require models for the aerosol and the structure of the atmosphere. Gordon and Wang [5] simplified the vertical structure of the atmosphere to just three layers: (1) next to the Fresnel-reflecting, flat, water surface, was a layer containing only aerosols of a specific type and characterized by the optical thickness τ a ; (2) a layer above the aerosol containing only Rayleigh-scattering molecules; and (3) a layer above the second displaying only absorption to account for Ozone (and which can be neglected in the radiative transfer computations). The bulk water itself was considered to be totally absorbing, so no radiance crosses the water surface from below. Calculations showed that, when the aerosol was non-or weakly absorbing, the radiance exiting the top of the atmosphere was insensitive to its vertical structure [43,44], so this atmospheric structure should yield radiances that are sufficently accurate to establish the relationships in Equations (28) and (29) as long as realistic aerosol models are employed.
The aerosol was modeled according to Shettle and Fenn [36]. Their aerosol models consisted of two components: a log-normal size distribution of small particles (modal diameter < 1 µm) and a similar distribution of large particles (modal diameter > 1 µm). Over the oceans the small fraction represented either the background aerosol or an aerosol transported from land. In contrast, the large-sized component represented locally generated aerosol (e.g., by breaking waves). Each component swells with increasing humidity and its refractive index decreases (but always remains greater than that of pure water). The background aerosol has a refractive index (dry) of approximately 1.53 with a modest absorption, while the large-sized component has an index close to water and virtually no absorption. Gordon and Wang [5] used the Maritime models in SF79 at four values of relative humidity (RH = 50, 70, 90, and 99%) and an additional model with RH = 80% was used to test the final algorithm. These are referred to as M50, M70, M90, M99, and M80. They also used models without the large-sized component and referred to these as terrestrial or tropospheric aerosols with the label "T" (i.e., T50, T70, etc.). Finally, they added a class of aerosols that had half of the large fraction compared to the maritime models and referred to these as coastal models (C50, C70, etc.). Summarizing, Gordon and Wang [5] used aerosol models with a background aerosol and a varying component representing the locally-generated marine contribution.
Full multiple scattering solutions to the (scalar) radiative transfer equation (using the successive order of scattering method) were carried out for each aerosol model (twelve in total) at eight aerosol optical depths (0.05 to 0.8), eight wavelengths (412 to 865 nm) and 33 values of the solar zenith angle (θ 0 = 0 to 80 deg in 2.5 deg increments), for a total of more than 25,000. Carrying out these simulations was a large undertaking, as each solution required about one hour on a 1 MHz workstation. It was effected by utilizing a large number of workstations at the laboratory (RSMAS) during hours when they would otherwise be idle. The result of each computation was twenty Fourier coefficients (in viewing azimuth relative to the Sun, φ v ) of ρ A at 100 viewing angles and the solar zenith angles above, i.e., A (λ) was then established as a power series in ρ (j) A via least-squares fitting. This provided look-up-tables from which . Examples of the relationship in Equation (27) are provided in Figure 6 for aerosols T50 and M99 which have the smallest and largest aerosol particles of the collection of models. The near linearity and similarity of the curves for λ = 765 and 865 nm suggests that the first term in Equation (27) is dominant and the multiple scattering effects are approximately the same at both NIR wavelengths. These multiple scattering results were used in the following manner. For a given pixel we determine the measured values of ρ A at λ s and λ l which we denote as ρ M A (λ s ) and ρ M A (λ l ), using the assumption that ρ w (λ) = 0 at both wavelengths. Then, for the jth model, we can find i.e., the values of ρ as at λ s and at λ l that would be valid if the actual aerosol were identical to that for the jth model. For the jth model then, we can find , the value of ε(λ s , λ l ) if the jth aerosol model were correct. As mentioned above, since multiple-scattering effects are similar in both NIR bands, e.g., Figure 6, the retrieved value of ε (j) (λ s , λ l ) will be close to the true value, independent of the aerosol model. That is, if we assume several (N) different models ("candidates") for the aerosol, but one model say the ith is correct, each model, e.g., the jth, will yield an ε (j) (λ s , λ l ) that will be close in value to the true ε (i) (λ s , λ l ). Since the ε's resulting from each of the trial aerosol models are close to one another, rather than pick a single model to estimate this quantity, as it is likely none of the aerosol models are correct, it seemed more reasonable to estimate ε(λ s , λ l ) through where ε (j) (λ s , λ l ) is the value of ε(λ s , λ l ) derived by assuming that the jth aerosol model is correct. To obtain a more robust estimate of ε(λ s , λ l ), the two models with ε(λ s , λ l ) values farthest (above and below) the mean are rejected, and a new average is computed omitting these two models. This is executed several times until only 4 models remain, the average of which gives the final value of ε(λ s , λ l ). Such an estimate will be most closely bracketed by two of the models, i.e., one with a slightly smaller value and one with a slightly larger value. Call these ε (−) (λ s , λ l ) and ε (+) (λ s , λ l ), respectively. Then as (λ l ) and ρ (+) as (λ i ) . Thus, we have two values for ρ A (λ i ). The procedure usually adopted was to assume that ρ A (λ i ) falls between ρ A (λ i ) in the same proportion as ε(λ s , λ l ) falls between ε (−) (λ s , λ l ) and ε (+) (λ s , λ l ). Some kind of assumption such as this is required, but this one is not universally true, in particular when the aerosols are strongly absorbing. The N aerosol models are called "candidate" models, and to the extent that the candidate models were similar to the actual aerosol, this particular assumption worked quite well.
This procedure yielded tρ w , where t is the diffuse transmittance from the surface to the sensor. Recall that earlier algorithms estimated t by assuming ρ w was totally diffuse and the atmosphere was aerosol free, i.e., t ≈ t * r . However, admitting that t depends on the aerosol concentration as well, and assuming that the upwelling radiance just beneath the water surface was totally diffuse, Gordon and Wang [5] used the principle of reciprocity (Yang and Gordon [45] In the Gordon and Wang algorithm, the diffuse transmittance was derived from t * (+) and t * (−) in the same manner as ρ a (λ i ) in the same proportion as ε(λ s , λ l ) falls between ε (−) (λ s , λ l ) and ε (+) (λ s , λ l ).
Tests of the algorithm with simulated data for θ 0 = 20, 40, and 60 deg. and θ ≈ 0 and 45 deg. suggested that the algorithm met the criterion that |∆ρ w | ≤ 0.002 for most cases when the true aerosol was SF97's M80, C80, or T80 (not among the candidate models). Using relationships similar to those in Equations (12) and (21), it was found that for pigment concentrations C P ≈ 0.1 or 1 mg/m 3 , for 95% of the cases tested |C P |/C P < 30%. Thus, the algorithm met the goals specified for the SeaWiFS sensor.

The OCTS and MERIS Multiple Scattering Atmospheric Correction Algorithms
During this time period one ocean color sensor was launched, OCTS in 1996 by NASDA (Japan), and another was in the planning phase, MERIS by ESA. Atmospheric correction algorithms that were similar in concept to the GW94 algorithm were developed for both sensors . The OCTS algorithm developed by Fukushima et al. [46] (F98) was virtually identical to GW94 in that it employed nine of the SF79 aerosol models (and one Asian dust model, which was strongly absorbing) with radiative transfer (RT) calculations to account for multiple scattering. They did not specify the vertical structure used in the computations, although it is largely irrelevant for the nine SF79 models, which are weakly absorbing.
The principle difference from GW94 was that F98 developed RT look-up-tables for the expansion of ρ A (λ) as a function of τ a (λ), i.e., ρ A = h(τ a ), where h is a polynomial function of τ a , as opposed to the tables of ρ A = g(ρ as ) utilized in GW94. Note that since τ a ∝ ρ as for a given Sun-viewing geometry, the g-tables can be converted to h-tables and vice versa. F98 chose λ s = 670 nm rather than 765 nm to avoid the problem of having to correct the 750 nm band for absorption by atmospheric Oxygen [47]. This choice restricts use of the algorithm to the open ocean because of the requirement that ρ w (λ s ) = 0 (see Section 5.1).
To operate the algorithm they first inverted the τ a -ρ A relationship in the red and NIR for each aerosol model: .
Then the γ (j) (λ s , λ )'s were averaged (or combined in a weighted average) to find the "derived" estimate of γ(λ s , λ ). Finally, the two "best" models among candidates were found based on comparison of the derived γ(λ s , λ ) with the individual γ (j) (λ s , λ )'s, and the value of γ(λ, 865) was determined at each λ by interpolation between the two "best" models based on their deviations from γ(670, 865). (Recall GW94 used a similar procedure replacing the γ's with 's.) The performance of this algorithm with simulated data was similar to that in GW94.
In the case of MERIS there were three bands that could be used for atmospheric correction: 705, 775, and 865 mn; however, the basic algorithm for correction utilized the 775 and 865 nm bands. The algorithm developed by Antoine and Morel [41] was also similar to that of GW94, in that it utilized aerosol models, but differed in more details than the F98 version.
Using radiative transfer theory they confirmed that the "path" reflectance (ρ path = ρ r + ρ a + ρ ra = ρ r + ρ A ) varied almost linearly with the aerosol optical thickness. In their algorithm they chose to develop ρ path /ρ r , rather than ρ path alone, as a function of the aerosol optical depth, i.e., where f (τ a ) is a polynomial in τ a . Comparing this to GW94, where ρ A was expanded as a function of ρ as , i.e., ρ A = g(ρ as ), it is easy to see that Thus, for any aerosol model and Sun-viewing geometry the two representations contain exactly the same information (as they must), i.e., for a given aerosol model and geometry f (τ a ) could be derived given g(ρ as ) and vice versa in much the same manner as h(τ a ) could be derived from g(ρ as ).
As in GW94, Antoine and Morel [41] (AM99) developed look-up-tables of ρ path /ρ r based on radiative transfer computations for an atmosphere bounded by a Fresnel reflecting surface, below which the water was totally absorbing. They used the Monte Carlo method of solving the transfer equation (radiance uncertainty about 2%). The vertical structure of the atmosphere was more elaborate than that of GW94, consisting of fifty 1 km thick layers with the molecular scattering properties assigned to each layer based on the air density. Ozone absorption was assigned to each layer based on the models of Elterman [48]. The aerosol was assumed to reside in three homogeneous layers (0-2 km, 2-12 km and 12-50 km), but unlike GW94, they assigned (different) aerosols to all three layers.
The layer nearest the surface (the marine boundary layer) contained locally generated aerosols with properties as modeled by SF79. The middle layer (the free troposphere) contained the "background" aerosol that had weak absorption (similar to the SF79 Tropospheric models) and a fixed optical thickness of 0.025 at 550 nm. The upper layer (the stratosphere) contained the non-absorbing stratospheric aerosol composed mostly of sulfu-ric acid and water and had an optical thickness of 0.005 at 550 nm. These combined with the variable marine boundary layer constituted what AM99 referred to as an "assemblage." The marine boundary layer could be composed of the SF97's M70, M80, M90, or M99 aerosols, and the choice of one of these models and the other two (fixed) layers led to one of four "standard assemblages." Other "assemblages" contain different components, usually strongly absorbing, in the marine boundary layer and/or the free troposphere, e.g., dust or the SF79 Urban models.
The standard assemblages were first used to perform an atmospheric correction as follows. First, the measured ρ path (865) is used to determine the set τ path (775), which is then compared with the measured counterpart ρ path (775). Finally, in a manner similar to GW94, the two best models, i.e., those that bracket most closely the observed ρ path (775), are used in the correction. The values of ρ path (λ) are then estimated by mixing the prediction of the two models in the same proportion as the excursions of the two "best" ρ (j) path (775)'s from ρ path (775). If the correction is successful i.e., the resulting waterleaving reflectances pass certain tests (requirements) for Case 1 waters, then the procedure is finished. If the retrieved water-leaving reflectances fail the tests, then the procedure is repeated with other (absorbing) assemblages until success is obtained. Simulations with pseudo-data for cases with weakly-absorbing aerosols suggest that this algorithm has approximately the same accuracy as that of GW94. An interesting difference between AM99 and the earlier two algorithms is the GW94 and F98 had to average the results of the calculations of (j) (λ s , λ )' and γ (j) (λ s , λ ) to obtain (λ s , λ ) and γ(λ s , λ ), respectively, for comparison with the individual models, while in AM99, the values of ρ path /ρ r , to be compared to the model results, came directly from the measured reflectance.
The principal differences between AM99 and GW94 are that (1) AM99 placed aerosols of different types as well as Rayleigh-scattering molecules in all layers with the variable aerosol component in the lowest layer, while GW94 had a single (variable) aerosol type only in the lowest layer, with all of the Rayleigh-scattering molecules in an (aerosol-free) layer above the aerosols; (2) AM99 modeled Ozone absorption using a continuous profile in the generation of the LUTs (discretized into 50 layers, as was atmospheric density for Rayleigh scattering), while GW94 did not include Ozone in their LUTs, but added it during the atmospheric correction process as a single layer at the top of the atmosphere; and (3) AM99 generated their LUTs using Monte Carlo methods with estimated (statistical) errors of the order of 2% while GW94 generated theirs with the successive order of scattering method with much smaller error (order of 0.1%). (With regard to (1) above, it should be noted that many of the GW94 models included both a pure maritime component (locally generated by wave breaking) and a terrestrial component; however, in contrast to AM99, these two components were in the same (lower) layer.) The procedures used in the three atmospheric correction algorithms are summarized in Table 2.
One should note that in the cases of weakly-or non-absorbing aerosols, e.g., the standard assemblages, the vertical structure of the aerosol is irrelevant at the accuracy required for atmospheric correction. While, for the strongly-absorbing assemblages, e.g., those with dust or the SF97 Urban models, the correction will fail in their presence unless the assumed vertical structure closely approximates the actual structure of the atmosphere (see Section 4.5). Table 2. Comparison of multiple scattering algorithms for SeaWiFS, OCTS, and MERIS. The superscript "j" refers to the jth aerosol model.

RT Tables
the N models to the N models to for each model find (λ s , λ ) find γ(λ s , λ )

Final
Interpolate between Interpolate between Interpolate between Selection the 2 "best" models the 2 "best" models the 2 "best" models * in nm.
There were two visible-NIR sensors developed mostly for land and aerosol applications during this era that could also be used for ocean color observations. These instruments were based on the new concept of measuring the angular distribution of radiance exiting the Earth's atmosphere and as such had the potential to determine the water-leaving radiance in several directions, albeit at slightly different times (up to 7 min). MISR [49] provided measurements of the radiance exiting the atmosphere at nine viewing directions distributed fore and aft of the satellite's orbit track, plus the nadir, and at four wavelengths (nominally 440, 550, 670, and 865 nm). POLDER [50] was capable of measuring the exiting radiance over a continuous range of viewing directions up to about 50 deg. from nadir, at wavelengths similar to MISR, but with the additional feature that the polarization of the radiance could be determined at 443, 670, and 865 nm as well. The atmospheric correction algorithms described in earlier sections could also be used with these sensors by applying them directly with λ s = 670 and λ = 865 nm. However, the availability of multi-directional reflectance provides further possibilities. In the case of MISR, Wang and Gordon [51,52] showed that the angular variation in ρ T measured by MISR, when compared to that produced by modeled aerosols, e.g., SF97, could be used to obtain the aerosol properties, and Gordon [44] built on this to show that atmospheric correction could be effected using a single band (865 nm, for which ρ w (λ ) = 0) by using the angular distribution of ρ T (λ ) to select SF79 aerosol models. Atmospheric correction based on a single band in the NIR is attractive because of its potential for application to Case 2 waters. Simulations showed that the efficacy of such an algorithm was slightly better than the SeaWiFS-type algorithm for Northern Hemisphere Summer geometries (small solar zenith angle), but much worse for Northern Hemisphere Winter geometries (large solar zenith angle). Note that such an algorithm requires that the phase function of at least one of the candidate aerosol models must closely match that of the actual aerosol extant in the atmosphere, i.e., a much stronger constraint than that of the GW94-type algorithm.
Frouin et al. [53] developed an approach with POLDER that was similar to that of GW94, but with some exceptions. They proposed to utilize the band at 670 nm along with λ s = 765 nm and λ = 865 nm. In essence, for each pixel they computed where i is an index over the set of viewing directions extant for the given pixel, and w is a weighting factor (set equal to 1 for the first instrument). This was then compared to its counterpart for RT simulations using the SF79 models to select the two "best" models, which are mixed by interpolation similar to GW94, with similar results [53]. To the author's knowledge, the polarization capabilities of POLDER have not been used to assist in the estimation of ρ A .
Comparison of the performance of the SeaWiFS/MODIS, OCTS, MERIS, and POLDER algorithms operating on a synthetic data set can be found in IOCCG Report No. 10 [54].

Absorbing Aerosols
Aerosols that have strong absorption, either neutral with wavelength, e.g., carboniferous aerosols from coal-fired power plants, or selective with wavelength such as desert dust with a significant Iron content (strong absorption in the blue), pose significant challenges for atmospheric correction. The first, and probably the most important, challenge is that all algorithms discussed so far utilized the observation that the radiance exiting a weakly-absorbing atmosphere can be considered to be independent of the aerosol vertical structure, at least at the accuracy required for atmospheric correction. This is no longer true when the aerosols are absorbing. An example of this is provided in Figure 7 for the Shettle and Fenn [36] Urban model at relative humidity 80% (U80). For this model ω a ≈ 0.78-0.74 for λ = 412-865 nm, respectively, and also, τ a ∝ λ −1.06 . The behavior seen in the figure, a progressive decrease in ρ A as a fixed amount of aerosol is mixed higher and higher into the atmosphere, is a characteristic behavior of an absorbing aerosol and demonstrates the strong dependence on vertical structure. One should also note from Figure 7 that measurement of ρ A in the NIR provides little or no information regarding the vertical structure even when the physical properties of the aerosol are known. However, if the aerosol is suspected to be absorbing, and candidate models are restricted to those with similar optical properties (e.g., U50, U70, U90, and U99, when the actual aerosol is U80), Figure 8 shows that an acceptable atmospheric correction can be obtained with the Gordon and Wang [5] algorithm, but only when the actual vertical structure is close to that of the candidate models, for which all of the aerosol is in a thin layer at the surface. Thus, it appears that atmospheric correction in the presence of an absorbing aerosol is possible if (1) the general nature of the aerosol is known, i.e., in the Figure 8 example the aerosol type is known to be Urban, and (2) the vertical distribution of the aerosol is known to be similar to that in the candidate models.  Figure 7. Influence of the physical thickness of the aerosol layer on the spectrum of ρ a + ρ ra = ρ A . For U80 the aerosol is confined to a thin layer near the surface, while for U180, U280, U480, and UU80, the aerosol is uniformly mixed with air to a height of 1 km, 2 km, 4 km, and the whole atmosphere, respectively. Viewing is near nadir and θ 0 = 60 deg and τ a (865) = 0.2. From Gordon, Du, Zhang [55].
How does one know that absorbing aerosols are present? The only way seems to be the appearance of unexpected results following the usual atmospheric correction, which assumes the aerosol is non-absorbing. For example, consider the U480 case in Figure 7. If one assumed that the aerosol was non-absorbing, then ρ A (443) would be approximately the extrapolation from that at 765 and 865 nm in Figure 7, or approximately 0.026; however, the true value would be about 0.014, so the estimated value is 0.012 too high, or tρ w (443) would be 0.012 too low. Figure 5 shows that if the true value of C P were 0.1 mg/m 3 , then an error of this magnitude would lead to a retrieved C P ≈ 1 mg/m 3 , so one manifestation of the presence of the absorbing aerosol might be a much higher-than-expected C P . Another manifestation would be a retrieved tρ w (443) < 0, that could occur in waters with even moderate C P .  Curves from top to bottom refer to situations in which U80 is confined to a layer just above the surface, between the surface and 1, 2, 4, and 6 km, and uniformly mixed throughout the atmosphere. From Gordon [44].
Another method of detecting atmospheric correction problems such as the unknown presence of absorbing aerosols is based on the right panel of Figure 5, which shows that [ρ w (550)] N varies within a restricted range of values. The Antoine and Morel [41] algorithm actually used [ρ w (510)] N , which shows even less variation than [ρ w (550)] N , to detect failure of the atmospheric correction using their standard assemblages, i.e., retrieved values of [ρ w (510)] N falling outside the range of acceptable values. When failure was detected, atmospheric correction was carried out using another set of assemblages tailored to the nature of the particular situation, e.g., dust, urban pollution, etc. Thus, their procedure for the detection of absorbing aerosols was based on the fact that one has an a priori expectation of the range of water-leaving reflectance values in some cases. That is, one has a model, either mathematical or experimental, of what to expect for [ρ w (λ)] N for a given location, e.g., for Case 1 waters.
Carrying this expectation further, Gordon, Du, and Zhang [55] proposed a method of solution, called the "spectral matching algorithm" (SMA), that required [ρ w (λ)] N for all values of λ to fit a particular semi-analytic bio-optical model [21]. In their algorithm, first the bio-optical model provided [ρ w (λ)] N as a function of C P and a particle scattering parameter b 0 , i.e., it was used to find [ρ w (λ, C P , b 0 )] N . Using the same approximation as in Equation (30), these provided trial values of for the ith aerosol model and waters characterized by particular values of C P and b 0 . Next, using the usual assumption that ρ w (865) = 0, τ (i) a (865) for the ith aerosol model was determined from ρ T (865), and this yielded t (i) (λ)ρ (i) w (λ) for the other spectral bands, i.e., the value t(λ)ρ w (λ) would have if the ith aerosol model were correct. Finally, the residual where n is the number of visible wavelengths, was examined. The value of δ was computed for each aerosol model and discrete set of ocean parameters. Note that, if for a particular parameter combination C P and b 0 and a particular aerosol model, δ = 0, then an exact solution of the problem has been found, i.e., the top-of-atmosphere ρ T (λ) is consistent with the water parameters and the atmospheric model. This is most unlikely, so normally the set of parameters i, C P , and b 0 that yielded the smallest δ(i, C P , b 0 ) would be chosen as the best, i.e., the solution to the problem. However, knowing that it is also unlikely that the "correct" model is one of the set of candidates, averaging over some number of the best retrievals (i.e., retrievals with the lowest values of δ(i, C P , b 0 )) to obtain the retrieved ocean and aerosol parameters seemed more reasonable. Gordon, Du, and Zhang [55] averaged over the best ten.
To test this algorithm they generated pseudo data using the M80, C80, T80, and U80 aerosol models for the atmosphere and the Gordon et al. [21] water-leaving reflectance model. Note that the same water model is used to generate the pseudo data and to operate the SMA; so this represented the optimum situation, i.e., a correct bio-optical model. Sixteen aerosol models are used as candidates: M, C, T and U models with RH = 50, 70, 90, and 99%. Note the aerosol models used to generate the pseudo data were not part of the candidate set. For this test, the vertical structure of the aerosol in the generated pseudo data is the same as that in the candidate models -again, an optimum situation. The test demonstrated that with the above limitations -correct water model and correct aerosol vertical distribution -the SMA performed very well. Of particular importance was the observation that the algorithm had no difficulty determining the presence of strongly absorbing aerosols and performed as well as the Gordon and Wang [5] algorithm when the aerosol was non-absorbing. In a real-world operation of the SMA, the major stumbling block was that one would have to incorporate information regarding the aerosol vertical structure when it was determined that the aerosol was absorbing.
Chomko and Gordon [56] proposed a different approach, but similar in spirit, to effecting atmospheric correction both in the presence of absorbing and non-absorbing aerosols. Rather than using the realistic aerosol models of SF79, they assumed a simple power-law size distribution: with fixed values of D 0 , D 1 , and D 2 . The aerosol was assumed to be homogeneous in nature (all particles have the same composition), and was completely specified by K, ν, m r , and m i , where the refractive index of the (assumed-spherical) particles is m r − im i . Thus, ρ A = ρ A (λ, K, ν, m r , m i ). (Note, for a given set ν, m r , and m i , τ a ∝ K.) They used the same water model as the SMA, so ρ w = ρ w (λ, b 0 , C P ). For SeaWiFS this provided eight equations of the form the set of which were "solved" by standard optimization techniques, hence the name spectral optimization algorithm (SOA). Applying this to simulated "data" developed using the SF79 aerosol models T80, C80 and M80, they found good retrievals of C P (range, 0.1 ≤ C P ≤ 1.0 mg/m 3 .), showing that accurate aerosol size distributions are not required to effect atmospheric correction; however, on the contrary, as one would expect, the derived aerosol optical thickness showed large errors. Excellent retrievals were also obtained for strongly absorbing aerosols (U80) but, again, only when the correct vertical distribution of the aerosols was known and employed in the algorithm. The SOA was tested using SeaWiFS imagery from the Middle Atlantic Bight and the Sargasso Sea by Chomko and Gordon [57]. They examined imagery from two days (DOY 279 and 281, 1997), one of which had a turbid atmosphere and the other a clear atmosphere. They found good continuity between the derived pigment concentrations over the two days, and a complete decoupling between the atmosphere and the ocean patterns for a relatively turbid atmosphere, but for the clear day water patterns were clearly evident in the retrieved aerosol products. This atmosphere-ocean coupling was significantly reduced by using a more sophisticated water model (Chomko et al. [58]). In that study, the water was modeled using radiative transfer along with the absorption coefficients of phytoplankton a ph (λ) and colored detrital material CDM, a cdm (λ) , and the backscattering by phytoplankton b bp (λ) . These parameters were further related to wavelength through where C is the concentration of chlorophyll a, a * ph (λ) is the chlorophyll-specific absorption coefficient of phytoplankton, S characterized the spectral variation of CDM, η specified the spectral variation of phytoplankton backscatter, and λ 0 = 443 nm. The parameters a * ph (λ), S, and η were optimized for Case 1 waters by Maritorena et al. [59] using in-situ data. Thus, the unknown parameters of the water model were C, a cdm (λ 0 ), and b ph (λ 0 ), i.e., ρ w (λ) = ρ w λ, C, a cdm (λ 0 ), b ph (λ 0 ) . When this model was tested with SeaWiFS imagery, a more complete decoupling between the water and the atmosphere was obtained, and good agreement was found between a cdm (λ 0 ) and its proxy measured contemporaneously by airborne lidar, and between the derived C and that from standard SeaWiFS processing. To apply such an algorithm to an absorbing aerosol that was distributed higher in the atmosphere, rather than in a layer at the surface, would require look-up-tables tailored to each likely aerosol vertical distribution. The required resolution in the vertical structure for effecting atmospheric correction is yet to be determined.
It should be pointed out that the models in this section that relate ρ w to the absorption and backscattering coefficients of the water's constituents are approximate in that the resulting upward radiance distribution is taken to be uniform, an assumption that must be invoked because the scattering phase function of the constituents in the water is not used in the calculation. A more realistic model has been presented by Morel, Antoine, and Gentili [60], in which the angular distribution of ρ w is directly related to C. This model required a full solution to the radiative transfer equation (including the Raman scattering) and the results were given in the form of lookup tables. They are for Case 1 waters only, and in application to atmospheric correction, have been used in the SeaWiFS implementation to remedy failures of the "black pixel" approximation (see Section 5.1).

Second-Order Effects
There are many effects that were omitted in the development of atmospheric correction algorithms as their influence was thought to be much smaller than the effect of aerosols, e.g., they were second-order effects. Among these were whitecaps from breaking waves, Sun glint cause by surface waves, departure from the plane-parallel atmosphere assumption due to the curvature of the Earth, and polarization of light in the ocean-atmosphere system. There were of course issues that were due to idiosyncrasies of specific sensors, e.g., in the case of SeaWiFS interference due to the Oxygen "A" band absorption near 765 nm, non-ideal spectral response functions, or excessive instrument response to the polarization of ρ T as in SeaWiFS and MODIS. Although important, these sensor-specific issues will not be discussed here (e.g., see [47,61,62]).

Reflectance Augmentation from Whitecaps
Whitecaps from breaking waves introduce additional radiance reflected from the surface over and above the water-leaving radiance. Their effect is small but can exceed [ρ w (λ)] N in the red and NIR, and therefore interfere with atmospheric correction [63]. The fraction of the water surface covered by whitecaps is strongly dependent on wind speed (∼ W a , where a ∼ 2.5 to 3.5, and W is the surface, or ten-meter-height, windspeed). Thus, their reflectance is strongly dependent on W as well. (A complete discussion of early studies aimed at relating whitecap coverage to wind speed and other parameters, such as atmospheric stability, etc., can be found in Stramska and Petelski [64]). Most studies attempted to relate the whitecap coverage (in % of surface covered) to the wind speed, and then to the augmentation of reflectance by multiplying the coverage by their effective reflectance (e.g., Koepke [65]). Whitecap coverage was usually determined from shipboard photographs of the water surface [65,66], but some measurements were obtained radiometrically [67][68][69]. Both determinations of coverage are difficult and the "derived" increased reflectance has significant variance (standard deviation of coverage ∼ coverage). Frouin et al. [67] showed that sea foam displays a clear decrease in reflectance from the visible to the NIR. Moore et al. [69] gave the radiometrically-determined increase in reflectance due to whitecaps as 3 × 10 −6 × S(λ) × W 2.55 , where S is the spectral variation being unity for all bands in the visible and 0.95 and 0.85, respectively for 750 and 860 nm. In contrast, from the Stramska and Petelski [64] photographically-determined windspeed-coverage relationship and the average reflectance of whitecaps (22% [65]), the reflectance was given by a wc × 19.25 × 10 −5 × (U 10 − 6.33) 3 , where a wc = 0.889, 0.765, and 0.645 at 670, 765, and 865 nm, respectively, (and unity elsewhere) was taken from Frouin et al. [67], and U 10 is the wind speed 10 m above the surface. The Stramska and Petelski [64] formulation was used to estimate and remove the whitecap contribution from ρ T in the SeaWiFS/MODIS processing.

Sun Glint
Sun glint, the specular reflection of direct solar radiation from the water surface to the sensor, can either be a zeroth, first, or second order effect. Looking near the specular image of the Sun, i.e., the position of the Sun in an image for a flat water surface, results in a huge radiance that exceeds all other terms in the radiance budget and can even saturate the sensor. This region, the radiance and spatial size of which depends on the wind speed (Cox and Munk [70]) was always simply masked from processing (not processed). Further away from the specular image the radiance is still large and, although it could be estimated given the wind speed and direction using the Cox and Munk [70] surface slope probability density, the relationship between the probability density and wind speed was not well enough known to be able to predict the radiance accurately, so this region was masked as well. In fact, the masks were so conservative that Sun glint was generally ignored. There were many Sun glint studies in the early 2000's aimed mostly in assessing the glint and reducing the size of the glint-masked regions, e.g., Wang and Bailey [71] implemented a correction using the Cox-Munk slope distribution that operated when the glitter reflectance was predicted to be less than about 0.005-0.01.

Polarization of Radiance in Ocean-Atmosphere System
In the Gordon and Wang [5] algorithm, the Rayleigh scattering contribution ρ r was computed using full vector radiative transfer theory, while the aerosol look-up-tables for the relationships in Equations (28) and (29) were generated using scalar radiative transfer theory. Understanding that ignoring polarization in radiative transfer can result in errors of a few percent, Gordon [44], in a restricted study, showed that the algorithm could perform somewhat better (applied to pseudo data) when vector theory was employed, but did not implement it for SeaWiFS. In a later study, Wang [72], recomputed the entire set of tables using vector theory. Use of the new tables showed significant improvement when the algorithm was operated with synthetic data; however, after vicarious recalibration, a procedure that is required whenever an algorithm is changed, there was a negligible difference in the derived products, i.e., the ρ w 's. In part, this result was a demonstration of the power of vicarious calibration and in part a demonstration of the robustness of the algorithm concept. When a new set of candidate aerosol models were introduced [73] (Section 5.4 below) the computation of the aerosol look-up-tables was carried out using a vector radiative transfer procedure.

Curvature of Earth's Surface
The GW94 atmospheric correction algorithm (Section 4.3) was based on radiative transfer (single and multiple scattering) within an assumed plane-parallel geometry (PPG) for the atmosphere as opposed to the correct spherical-shell geometry (SSG). Although one would expect the PPG to be an excellent approximation, its applicability will clearly break down when the solar zenith angle becomes large, i.e., at high latitudes. Gordon and Ding [74] tested the applicability of the PPG by creating pseudo data in SSG and then using the GW94 algorithm (based on PPG) to effect atmospheric correction. They found that for θ 0 ≤ 70 deg the algorithm performed as well as with pseudo data created in PPG; however, for θ 0 > 70 deg, significant deviations were observed. In addition, they discovered that the additional error could be almost eliminated completely by computing ρ r (λ) in SSG rather than in PPG. They also presented a simplified computation scheme for ρ r in SSG based on the observation [75] that the ratio of the true ρ r to that computed assuming single scattering was approximately the same in PPG and SSG.

Activities in the Early 2000's
In the early years of the new millennium effort was focussed on addressing atmospheric correction issues that were problematic. These included relaxing the requirement of negligible water-leaving radiance in the NIR; demonstrating that absorbing aerosols could be dealt with, at least on a case-by-case basis; improvement of "candidate" aerosol models for use in the GW94 algorithm; and addressing the problem of atmospheric correction in sediment-dominated Case 2 waters.

Relaxation of the "Black Pixel" Assumption
An assumption usually made in discussing atmospheric correction was that the water leaving radiance was negligible in at least one spectral band: the band at 670 nm for CZCS and the two bands in the NIR (765 and 865 nm) for SeaWiFS. If this were not the case, then the retrieved values of ρ w (λ) would be in error. Consider the CZCS. If ρ w (670) = 0, then the retrieved value of ρ w (λ), which we will denote by {ρ w (λ)} Ret , is easily seen to be and is thus too small (and may even be negative). For CZCS, the reflectance change (at θ 0 = 0) for one digital count (DC) was as follows: 443 nm, 0.00075; 520 nm, 0.00053; 550 nm 0.00042; and 670 nm, 0.00024. Direct measurements and modeling (e.g., see Bricaud and Morel [23]) show that for a pigment concentration of 0.1 mg/m 3 in Case 1 waters the value of ρ w (670) ≈ 1 DC, while at C P = 1 mg/m 3 , ρ w (670) ≈ 3-4 DC. Thus, ignoring the t's in the above equation, and assuming that (λ, 670) can be approximated as (670/λ) n with n = 1, {ρ w (443)} Ret will be essentially unchanged for C P = 0.1 mg/m 3 , but too small by approximately 1-2 DC for C P = 1 mg/m 3 . This would lead to an overestimate for C P by ∼30% or more for C P = 1 mg/m 3 Equations (21) and (22) . For larger C P the error in {ρ w (443)} Ret is larger still. In comparison to CZCS, the SeaWiFS sensor had roughly four times the sensitivity: a factor of 2 2 through ten-bit digitization compared to eight-bit for CZCS. It also had spectral bands placed in the NIR specifically for atmospheric correction. Because of the much larger absorption of pure water in the NIR, the water-leaving reflectance in these bands is much lower than that at 670 nm. In fact, the water absorption coefficient at 765 nm is about five times that at 670 nm, while at 865 nm the corresponding factor is about 10. Since the reflectance is proportional to the backscattering-to-absorption ratio [6], and particle backscattering is slowly varying with respect to wavelength, one would expect that ρ w (765) ≈ ρ w (670)/5 and ρ w (865) ≈ ρ w (670)/10. However, as the SeaWiFS sensitivity is 4× that of CZCS, ρ w for the two sensors will be approximately the same at 670 nm (CZCS) and 765 nm (SeaWiFS) when measured in units of digital counts. Thus, the comments in the preceding paragraph apply to SeaWiFS as well, but with one important difference. For SeaWiFS, the relevant factors are determined using the NIR, and the fact that the water contribution is more significant at 765 than at 865 nm will result in an overestimation of causing in an additional error in atmospheric correction. Siegel et al. [76] showed that the overall error in C resulting from this effect could range from approximately −20 to −160% for C = 5 mg/m 3 , depending on the bio-optical algorithm used for the conversion ρ w (λ) → C. Thus, it was important to correct for the non-blackness of the water-leaving reflectance in the NIR.
Siegel et al. [76] referred to the assumption the ρ w = 0 in the NIR as the "black pixel" assumption, and developed a method for its correction. Briefly, they used empirical models of the dependence of the inherent optical properties of Case 1 waters and their suspended particles (the latter being plankton and their particulate detritus) on the concentration of chlorophyll a (C), and a simple radiative transfer model [21] relating these to the waterleaving reflectance to estimate the ρ w (NIR) as a function of C. Their scheme was iterative beginning with ρ w (NIR) = 0, which provided an initial estimate of (λ s , λ l ) and therefore ρ w (λ), then The iteration usually converged after 3-4 estimates of C. Using this correction, Siegel et al. [76] showed good agreement with the distribution of C (frequency of occurrence of C versus C) in the Chesapeake Bay, and the scheme was incorporated into the standard SeaWiFS processing. Later, Bailey, et al. [77] provided an updated, and more accurate, model that eliminated the need for C in the above iteration scheme by relating the particle backscattering's spectral variation to the values of ρ w at 443 and 555 nm, in a manner first formulated by Carder et al. [78].

Absorbing Aerosols Case Study: Saharan Dust
In order to employ the SMA in an atmosphere with absorbing aerosols, one needs a set of aerosol physical models that capture the absorption and scattering properties of the particular absorbing aerosol extant and have some knowledge of its vertical structure. These would need to be developed for each particular geographical region to be studied. Moulin et al. [79] performed a case study providing a procedure designed to meet these requirements. The case they chose was that of Saharan dust transported by the winds off the West Coast of Africa out over the Atlantic Ocean.
Such large dust plumes often flow off the African Continent over the Northeast Atlantic and are clearly visible in imagery from many space-borne sensors. In the regions of highest concentration, the reflectance of the dust is so high that atmospheric correction would be impossible; however, such regions do provide the opportunity to gain some information regarding the dust's optical properties. Moulin et al. [79] ignored the water contribution to the reflectance at the SeaWiFS sensor (in regions of high dust concentration) allowing the determination of ρ A (λ). The results displayed a decreasing ρ A progressing from NIR through the blue wavelengths similar to that shown in Figure 7. This is in contrast with the behavior of ρ A for non-absorbing aerosols which shows a steady or increasing ρ A with decreasing λ. The inhomogeneous nature of the plume allowed the establishment of a relationship between ρ A (λ) and ρ A (865) for each spectral band. Well-known physical models of Saharan dust, e.g., a size distribution from Shettle [80] and a refractive index from Patterson [81], were then combined with simple vertical distributions of the dust (also similar to those in Figure 7) to produce theoretical values of ρ A (λ) as a function of ρ A (865) for each model (i.e., each combination of the physical properties and the vertical distribution). (Note that τ a (λ) can be found from τ a (865) given the model of the aerosol, so such calculations do not require knowing the actual optical depth of the aerosol in the image.) For the vertical distribution in their computations they examined cases with the aerosol located just above the surface and the aerosol mixed uniformly with air to an altitude of 2, 4, 6, or 8 km with only the air above it (respectively, V02, V04, V06, or V08); and the aerosol uniformly mixed with air the entire atmosphere (VUU). Comparison of the predicted ρ A (λ) vs. ρ A (865) relationship with that from the imagery showed that for some regions of the image one of the extant models (referred to as BDS1, BDS2, and BDS3) provided relationships in good agreement with the experimental observations, while in other regions all of these models had too much absorption. Thus, Moulin et al. [79] developed a new set of models based on the original set, but with less absorption (BDW1, BDW2, BDW3).
To apply the SMA to imagery with Saharan dust, Moulin et al. [82] combined the six physical models (BDS1, . . . BDW3) with three vertical distributions (V02, V04, V06) for a total of 18 possible atmospheres. They applied the SMA to retrieve the pigment concentration off the African Coast from SeaWIFS imagery. Although there were no surface measurements, comparisons between one "clear" and two "dusty" images obtained over seven days showed excellent continuity between pigment retrievals over the three images, and good agreement with the monthly means of chlorophyll a. [Recall that post-CZCS bio-optical algorithms return chlorophyll a, while the Gordon et al. [21] model returns C P .] Further evidence of the SMA's efficacy was provided in Banzon, et al. [83], wherein it was used to study the evolution of the pigment concentration in the Northern Arabian Sea during the southwest monsoon. Again, good retrievals were obtained over dusty regions of images that had been masked (unprocessed) by the standard SeaWiFS processing. These studies supported the contention that atmospheric correction could be effected in regions with strongly-absorbing aerosols. In the future, when aerosol vertical structure is available from other sensors, such retrievals could be carried out routinely.

Improved Aerosol Models for the Gordon and Wang Algorithm
To achieve a better understanding of the aerosol in a "clean" maritime environment, Smirnov and coworkers [84,85] analyzed a large set of size distribution retrievals from three remote areas: Lanai in the Pacific (Hawaii), Bermuda in the Atlantic, and Kaashidoo in the Indian Ocean (Maldives). The size distributions were retrieved from measurements of the Sun and sky radiances made in the Aerosol Robotic Network (AERONET) [86]. These retrievals yield the columnar size distribution, i.e., the fictive distribution assuming the particle size frequency distribution of whole column of atmosphere is independent of altitude. Although a broad distribution of aerosol optical depth at 500 nm was observed (mostly less than 0.3) the values retrieved with the highest frequency were 0.06, 0.09, and 0.11 at Lanai, Bermuda, and Kaashidoo, respectively. Limiting their analysis to retrievals with τ a (500) ≤ 0.15, they retrieved the mean columnar volume size distribution shown in the left panel of Figure 9, and fit these to a volume size distribution where the individual dV /d nR's are log-normal. Here V is the volume of particles contained in a vertical column of atmosphere with a base of 1 µm 2 , i.e., the units of V are µm 3 /µm 2 , and R is the particle "radius." The "coarse" mode is primarily sea salt. The right panel of Figure 9 compares the average Lanai distribution with the analytical fit. Figure 10 compares the fine and coarse modes with the similar components from the Shettle and Fenn models with RH ≤ 50%. We see that the modal sizes are similar; however, the Shettle and Fenn distributions are considerably wider, especially the fine mode. It should be noted that the measurements are for ambient RH, and the particles are expected to increase in size with increasing RH. Figure 9 is for what is considered to be a pure maritime aerosol; however, when compared to the full data set from the three locations, most of the variability is in the relative contributions of the two modes, i.e., the relative values of the total particle volume (V T ) for the fine and coarse modes. An analysis of the correlation between the size distribution and the wind speed at Midway Island showed that the modal size of both modes was virtually independent of W , the wind speed averaged over 24 h, as is V T for the fine mode. In contrast, the total volume of particles of the coarse mode increased roughly linearly with W at a rate of 0.0073 µm 3 /µm 2 per m/s of wind. (Lanai) size distributions for the fine and coarse modes. Each mode has been normalized to unity at the size of maximum dV/d nD. Note that in contrast to Figure 9, the horizontal axis is particle diameter as opposed to radius. The values of V T for the fine and coarse modes at Lanai are 0.010 and 0.039 µm 3 /µm 2 , respectively. (Taken from [11]).
More recent data [87] from AERONET island stations provide mean values for the aerosol optical thickness τ a (500) and the Ängstrom power p (τ a (λ) ∝ λ −p ) for several of the world's oceans. The global average is τ a (500) = 0.108 and p = 0.573. The Southern Ocean has the clearest atmosphere (τ a (500) = 0.060 and p = 0.380) and lowest variability, while the Atlantic Ocean has the most turbid atmosphere (τ a (500) = 0.190 and p = 0.604) and the highest variability. An important conclusion to be drawn from these studies is that the for most of the World oceans, τ a (500) ≤ 0.2, and so with the Atlantic value of p, τ a (865) ≤ 0.14. In Gordon and Wang [5] the efficacy of the atmospheric correction algorithm was assessed taking τ a = 0.20. These results show that τ a (865) = 0.20 for the open ocean is quite conservative and explains the success of the single-scattering variant applied to CZCS in the open ocean.
The AERONET-derived models [84,85,87] lead to optical properties that are similar to those of Shettle and Fenn [36]. The scattering phase function and spectral variation of the extinction coefficient, computed assuming spherical particles with a complex refractive index of 1.37 − 0.001i and the Lanai size distribution are shown and compared with the same quantities for several SF79 models in Figure 11. The Smirnov et al. models clearly yield optical properties that fall within the range of variability of the Shettle and Fenn models. Building on the work of Smirnov, et al. [84,85,87], Ahmad et al. [73] developed a set of aerosol models for deriving more accurate values of aerosol optical thickness (τ a ) from ocean color sensor data. As with the SF79 models, in the Ahmad et al. [73] models the modal particle sizes varied with relative humidity (RH) as did the particle refractive index. In their models, the "coarse" mode was primarily sea salt and water and, as such, non-absorbing. All absorption was placed in the "fine" mode, which was taken to be continental in origin. The composition of the fine mode was adjusted to achieve the spectral dependence of absorption (actually ω 0 ) seen in the AERONET retrievals.
The Ahmad et al. [73] models consisted of eight values of RH from which the parameters of the size distributions of the fine and course modes were established. For each RH value, the relative concentration by mode volume of the fine mode was varied from 0 to 1, defining ten models. Therefore they employed a total of 80 distinct aerosol models compared to 12 in the original GW94 algorithm which were defined by varying the RH and the relative concentration of the two modes by the total number of particles (rather than volume) in each mode. In their pixel-by-pixel implementation, they first selected the values of RH in the model set bracketing the RH of the observation that was taken from NCEP analysis for the given pixel, and then selected the best two models based on the reflectance in the NIR in a manner similar to GW94. Before global processing, they carried out a complete vicarious calibration with the new models in place. Application of the algorithm with the new models showed significant lowering (compared to the original models) in the retrieved τ a (λ) in better agreement with AERONET observations. The generally lower τ a (λ) is partially explained by a fact that the newer aerosol models, based on the AERONET models over the open ocean, generally result in larger values for the scattering phase function at scattering angles greater than about 140 deg. An example is shown in Figure 11 for the Smirnov et al. model. Since the bulk of the aerosol reflectance is proportional to τ a P, a larger P requires a smaller τ a for the same aerosol reflectance in the NIR ρ T (NIR) − ρ r (NIR) .
It was found in utilizing the new models that the retrieved water-leaving reflectances were essentially unchanged from those using the SF79 models. This again demonstrated the robustness of GW94 and the power of vicarious calibration: as long as the aerosol models used in atmospheric correction reasonably approximate the range of possible variations of (λ, λ l ), they will provide acceptable values for ρ w (λ), after vicarious calibration is effected. This is a less stringent requirement than that for deriving τ a , for which a realistic phase function is required as well.

Atmospheric Correction for Turbid (Case 2) Waters
This Special Issue of Remote Sensing focusses in the problem of atmospheric correction in coastal regions. These regions are generally classified as Case 2 waters. As mentioned earlier, Case 2 waters are defined as those for which the optical properties cannot be estimated from the chlorophyll concentration alone [6,7]. Such waters include those with high concentrations of suspended inorganic particles (e.g., resuspended sediments), waters laden with riverine-sourced dissolved organic matter (as distinct from dissolved organic matter derived from phytoplankton degradation in Case 1 waters), or both of these constituents. They can also contain very high concentrations of phytoplankton. In the following I shall adopt the terminology of Kirk [88] and refer to the concentration of the retinue of phytoplankton, particulate detritus, and inorganic suspended particles as the seston concentration (S). Seston-dominated Case 2 waters provide a challenge for atmospheric correction because they can have significant water-leaving reflectance in the NIR. Furthermore, the dependence of ρ w (NIR) on the seston concentration depends on the nature of the suspended material and thus may be location specific. For waters dominated by dissolved organic material, the question of atmospheric correction is almost moot, as such waters tend to be strongly absorbing throughout the visible, and even under the most favorable circumstances, ρ w (λ) may be too small to be accurately derived.
In this section I will briefly introduce concepts relating to atmospheric correction in seston-dominated Case 2 waters, leaving the detailed description to other papers in this issue. Nevertheless, I will describe a correction method utilized in the MODIS and VIISR sensors, which have additional spectral bands in the short-wave infrared (SWIR, λ > 1000 µm) portion of the spectrum where the water is so absorbing that ρ w (SWIR) = 0 is almost always an excellent approximation. One should note, however, that the closer a band is to the NIR-SWIR boundary (on the long-wave side) the higher the probability that there may still be detectable ρ w with modern, more-sensitive, instruments.

Traditional Approaches (Homogeneous Aerosol Type: Invariant )
Consider the number of equations and the number of unknowns in Equation (23) (including nonzero values for ρ w in the NIR) for the 8-band SeaWiFS sensor. There are seven equations, eight unknown values of ρ w and seven unknown values of . Six values of are determined from (λ s , λ ) (as long as the aerosol is at most weakly absorbing), leaving a total of nine unknowns and seven equations. The GW94 algorithm closes the system with ρ w (765) = 0 and ρ w (865) = 0. However, the latter two relationships are not valid in seston-dominated Case 2 waters, so in these waters two more relationships among the various unknowns are required. (Note, the diffuse transmittance t(λ) is not regarded as an unknown here as it is only weakly dependent on the aerosol, and can easily be corrected/updated after an estimate of the aerosol model(s) becomes available.) One way of closing the Equation (23) system, and that uses only information from the sensor itself, has been developed by Ruddick et al. [89]. It was similar to the Gordon et al. [14] method for the CZCS that uses "clear water" in a scene to determine correction parameters and then uses these parameters for the whole scene. Their idea was to use very turbid water areas in a scene to estimate α = ρ w (765)/ρ w (865) and to use areas with low turbidity water to estimate (765, 865) = ρ A (765)/ρ A (865). Then, letting ρ rc = ρ T − ρ r , where the subscript "rc" stands for "Rayleigh corrected," one has ρ rc (865) = ρ A (865) + t(865)ρ w (865) ρ rc (765) = (765, 865)ρ A (865) + αt(765)ρ w (865).
The two parameters α and are determined by plotting ρ rc (765) as a function of ρ rc (865). The behavior of this plot for large values of ρ rc (865) determines α and for low values determines . (Actually, α ≈ a w (865)/a w (765) = 1.72, since ρ w ∝ backscattering/absorption, backscattering is a weak function of wavelength, and the absorption is due almost entirely to water itself.) For the method to work, one must have a wide range of turbidities in a given scene. In particular, there must be areas for which ρ w (NIR) ≈ 0. In addition, (765, 865) was assumed to be approximately constant over the scene. Since the t's are easy to estimate, the above equations can then be solved for ρ A (865) and ρ w (865). Then the usual GW94 algorithm was used to estimate (765, 865) and finally two bracketing aerosol models as before, effecting an atmospheric correction in the visible bands. A similar technique was described in Hu et al. [90].
Lavender et al. [91] and Moore et al. [92] developed a Case 2 algorithm for MERIS that they tested with SeaWiFS. Through laboratory measurements, they developed relationships between the NIR (and red) ρ w 's and S (which in their case was mostly mineral particles) to close the system as described above. Using these, the atmospheric contribution to ρ A + tρ w in the red and NIR was deduced by solving the resulting nonlinear equations numerically. The ρ A 's corrected by this procedure were then entered into the GW94 algorithm to complete the atmospheric correction. Such algorithms can in general be summarized as follows. First, develop a relationship between the seston concentration and its inherent optical properties. Use these in radiative transfer models to estimate ρ w (λ) as a function of the seston concentration (S). Then through an iterative procedure that conceptually looks like · · · S→ ρ w (NIR) with the final ρ w (λ) being used to derive other products, e.g., chlorophyll a, as well. This is much the same as the modification of the black pixel assumption in Section 5.1. As this special issue of Remote Sensing deals extensively with such algorithms, they will not be discussed further here.

Use of More Detailed Water Models Coupled with Radiative Transfer
There were also attempts at using models of the water reflectance as in the spectral optimization algorithm. Kuchinke and coworkers [93,94] modified the SOA to handle cases in which ρ w (NIR) = 0. Briefly, they operated the SOA assuming that ρ w (NIR) = 0 deriving the three water parameters C, a cdm (443) and b bp (443), with which they computed estimates for ρ w (NIR) and subtracted these from ρ T (NIR) and reentered the "corrected" ρ T (NIR) into the SOA, etc. They found an excellent decoupling between the atmosphere and the water. Because of the sparse nature of the surface truth, the retrievals of C were compared to surface measurements on a frequency distribution-to-frequency distribution basis rather than location-to-location basis, and showed similar distributions. A similarly based water model for Case 2 waters was developed by Land and Haigh [95,96] for atmospheric correction and retrieval of water properties. The model was tested with pseudo-data and performed well.
These studies suggest that the combination of aerosol and water models with the "best" parameters derived through iteration or optimization (or in the case of the SMA "brute force") can produce reasonable results in Case 2 waters.

An Alternative to the NIR: The SWIR
In the case of MODIS (and VIIRS) there was an alternate approach to dealing with seston-dominated Case 2 waters by virtue of the spectral bands in the short-wave infrared (SWIR) region of the spectrum, present on this sensor [44]. In atmospheric windows MODIS has SWIR bands at 1240, 1640, and 2130 nm. At these wavelengths the absorption coefficient of pure water is approximately 88, 498, and 2228 m −1 (Hale and Querry [97]), respectively. Compared to the "atmospheric correction" band at 865 nm for which the water absorption coefficient is about 5 m −1 , one sees that the the assumption ρ w (SWIR) ≈ 0, will hold at much higher seston concentrations than in the NIR. Wang [98] studied the efficacy of the GW94 algorithm, when λ s and λ l were shifted from the NIR to the SWIR. He concluded that SWIR combinations λ s = 1240 and λ l = 1640 nm, as well as λ s = 1240 and λ l = 2130 nm performed as well as the combination of the two NIR bands; however, the signal-to-noise ratio required for such performance was an order of magnitude higher than that extant in the MODIS sensor: the sensitivity of the SWIR bands was specified following requirements for use in land remote sensing. Nevertheless, Wang and Shi [99] demonstrated the efficacy of the SWIR bands in retrieving ρ w (NIR) in turbid waters. Given ρ w (NIR), the GW94 algorithm could then be applied to effect atmospheric correction in the visible.

Other Atmospheric Correction Methods
In the early 2000's considerable efforts were undertaken to apply sophisticated methods of optimization (e.g., as used in the SOA), such as neural networks, to the problem of atmospheric correction and the general problem of water property retrievals. These will not be discussed here, as they deviate significantly from the operational methods used to process global data which were the principal focus of this historical review. The reader can deduce the "flavor" of such methods from the (now vast) literature, e.g., [100][101][102][103][104][105][106][107][108].

Concluding Thoughts
There is little doubt that the success of the CZCS led to the development of ocean color as a significant component of research related to living marine resources and to climate change on a global scale. At present count there have been 29 ocean color sensors following CZCS placed in space and 8 in the planning stage (IOCCG website). Its success led directly to the follow-on sensors SeaWiFS, OCTS, and MERIS. It seems reasonable to assert that the success of the CZCS followed directly from the ability to execute atmospheric correction to an accuracy that allowed application of the bio-optical algorithms that retrieved the pigment concentration in the global oceans to an accuracy well within a factor of two. This ability resulted partly from a good measure of serendipity.
The choice of the CZCS spectral bands, although well chosen for retrieving pigments given ρ w , was inadequate for atmospheric correction, a fact already known prior to launch [39]. Apparently there was no consideration given to atmospheric correction in the design of the CZCS, only to ways of possibly minimizing the atmospheric contamination, e.g., using band ratios ρ T at various wavelengths [109]. The red band, that ultimately was central to atmospheric correction in the open ocean, was actually included on CZCS for the detection of chlorophyll at high concentrations and had significant water-leaving radiance only at the very high chlorophyll concentrations that one might find in the coastal zone, estuaries, and lakes. (Recall that chlorophyll a has an absorption band at 670 nm, so any enhanced radiance there had to be due to the backscattering of plankton containing the pigment, their associated detritus, and whatever other kinds of particulates are present-a very indirect and insensitive method of measuring chlorophyll.) The NIR band on CZCS was intended for surface vegetation and was too insensitive for oceanic application. Interestingly, the excess radiance in the red band, although minimal, prevented atmospheric correction in high pigment areas, e.g., the coastal zones, so we had the conundrum that the Coastal Zone Color Scanner performed poorly near the coasts, but, in the end, very well in the open ocean. It is interesting to note that the CZCS deficiencies led most of the oceanographic community (except the experiment team, Table 1) in the prelaunch era to consider the CZCS project a waste of precious funding. Had the red band not been included, atmospheric correction on the open ocean may not have been possible, and the CZCS would likely have been considered a failure. As the time between the launches of CZCS and SeaWiFS was 19 years, one is led to wonder how much that time interval would have expanded if CZCS atmospheric correction had failed. But it did not fail, the CZCS mission proved that chlorophyll (actually pigments) could be measured on a global scale and, concomitantly, an indication of primary productivity on a global scale, and the CZCS truly became the "proof-of-concept" mission.
What are the optical characteristics of the atmosphere that led to the possibility of atmospheric correction of CZCS? First and foremost is the fact that the atmosphere is basically optically thin, meaning that single-scattering of sunlight provides a reasonable approximation to the atmosphere's contribution to ρ T (Figures 1 and 2). [It also means that one can image through the atmosphere, an obvious requirement for space-borne ocean color sensors.] In particular, molecular scattering, the largest contribution to ρ T in the blue (the most important spectral band for estimating C P in the open ocean, Figure 5, left panel), could be computed using single-scattering theory, as it was in the initial validation studies [14,16]. It enabled atmospheric correction on a pixel-by-pixel basis to be carried out with the computational resources available at that time. Second, because radiative transfer in the atmosphere was close to single scattering, Rayleigh and aerosol scattering were independent (ρ ra ≈ 0) so ρ A ∝ τ a ( [10], and Equation (4)). This meant that ρ A at different wavelengths were proportional to one another (Equation (5)), and therefore could be determined from one another. Third, although the proportionality "constant" (λ 1 , λ 2 ) (Equation (10)) depends on the aerosol phase function, the fact that the dependence is weak, especially for marine aerosols, which are the main aerosol in the open ocean, meant that (λ 1 , λ 2 ) ≈ τ a (λ 1 )/τ a (λ 2 ), which can be accurately modeled by Ängstrom's law. Thus, if the aerosol's Ängstrom exponent remained constant over a scene (a constant aerosol "type," i.e., uniform aerosol size frequency distribution and refractive index) and could be determined in one part of the scene, then ρ A (λ 1 ) could be determined from ρ A (λ 2 ).
[Although not a characteristic of the atmosphere, a characteristic of Case 1 waters that facilitated this process was the "clear water" radiance concept [18].] The presence of the red band (670 nm), where ρ w ≈ 0 for most of the world's oceans, facilitated the determination of ρ A (670), from which ρ A (λ) could be found leading to ρ w (λ) for use in the bio-optical algorithms (Equation (12)). These characteristics of the atmosphere (and "clear" water), enabled the arguably "rapid," correction of imagery obtained early in the mission and preparation and publication of initial results in a prestigious journal [16,17], fending off critics. In addition, the fact that (λ, 670) ≈ 1 for a marine aerosol (e.g., M99) enabled processing of the full CZCS archive without the cumbersome process of locating clear water to determine in several sub-scenes of each of the approximately 68,000 images [35,38].
What would happen if one of the links in the above chain were broken, e.g., the weak dependence of the aerosol phase function on wavelength. The manifestation of failure of this particular link would be that (λ 1 , λ 2 ) determined from one portion of the image, e.g., "clear water," could not be "carried over" to another portion of the image (without explicit knowledge of the phase functions, Equation (10)), and the resulting "corrected" image with the assumed constant would contain easily-identifiable atmospheric features (from spatial scale differences between oceanic and atmospheric features)-a failure of atmospheric correction in particular and the CZCS mission in general.
I cannot over stress the importance of the fact that atmospheric radiative transfer, at the level of accuracy that was required for CZCS, is well-approximated by single scattering. In early 1979, when CZCS imagery became available, the image processing system (the Atmospheric and Oceanic Image Processing System-AOIPS at NASA/GSFC) available for our use with CZCS imagery was powered by a DEC PDP 11/55 mini computer (circa 1976). Although state-of-the-art, it could only perform simple manipulation of the various spectral bands (image planes) of an image. Initially, I computed the Rayleigh reflectance ρ r in scan coordinates (line#-pixel#) using a Univac 1106 mainframe (circa 1965) computer, which would be a "toy" by today's standards. In this form ρ r could be simply subtracted from ρ T to form ρ A using AOIPS. Then using ρ A (670) determined by assuming ρ w (670) = 0 and ρ A (λ) = (λ, 670)ρ A (λ), the desired ρ w (λ) was found (Equation (9), with t = 1). The first CZCS image to be atmospherically corrected in this manner (or any manner) is provided in Figure 12. The correction was based on measured L w 's at a location in the upper portion of the image, enabling (λ, 670) to be determined there, and then used throughout the sub-scene. It was photographed from the AOIPS monitor screen by the author as atmospheric correction of the 443 nm band of CZCS took place, line-by-line. The time to compute this simple correction for a single scan line (512 pixels in the subscene) was of the order of ten seconds. As I watched the correction take place, it was as if a veil were slowly being pulled downward off the image -and it was exhilarating.
Were the algorithm to require multiple scattering computations, it is hard to envisage (1) how imagery could be processed in a fashion timely enough to forestall designation of the CZCS mission as a failure, and (2) how more than a few CZCS scenes could have been processed in the 1980s. "Fortunately, by the late 1980's, with the introduction of the DEC VAX line of super mini-computers, and software developed at the University of Miami and installed at NASA/GSFC, it was possible to process the entire CZCS archive of images (∼68,000 images, each with 1968 pixels/line and 970 lines) and demonstrate the potential of ocean color imagery [35,38]. Comparison of the radiometric performance of SeaWiFS with CZCS showed that a correction algorithm for SeaWiFS had to be significantly more accurate than the simple CZCS algorithm, and in fact required a full multiple-scattering algorithm. Although the spectral bands in the NIR facilitated atmospheric correction and indeed a single-scattering version using these bands was quite accurate [42], it was not accurate enough. M. Wang and I begin working on a multiple-scattering version (GW94) in 1989 when SeaWiFS was approved by NASA. Although multiple scattering was required, success with the singlescattering algorithm suggested that we should try and use the same structure. Early tests were carried out with aerosol models for which both the phase function and the extinction coefficients were independent of wavelength to approximate a marine aerosol (the aerosol was also non-absorbing). Results were encouraging, suggesting that a correction that was roughly an order of magnitude better than CZCS could be effected [110,111]. However, when spectral variations of the aerosol optical properties expected for realistic aerosol models, e.g., SF79, were incorporated into the algorithms and its pseudo data, the resulting accuracy was reduced by a factor of two to three-fortunately, still accurate enough for SeaWiFS.
Implementation of the algorithm for operational SeaWiFS processing posed problems. First, the generation of the LUTs containing the expansion coefficients a, b, and c in Equation (27) was a significant computational problem. The computations required solving the radiative transfer equation at eight wavelengths, 12 aerosol models, 33 values of the solar zenith angles, and eight values of the aerosol optical depth for a total of 25,344 simulations. Our radiative transfer code (successive order of scattering) required one hour per simulation on a DECstation 3100, so the total time required for a single work station was approximately 2.9 years. The actual computations were carried out in about one month using a large number of work stations at RSMAS during times when they would otherwise be idle (thanks to Robert Evans and James Brown for developing a scheduler that effected this task). Second, the LUTs had to be available for rapid access for pixel-by-pixel processing, i.e., stored in the computer's main memory. Because our successive-order code used Fourier decomposition in the azimuth, rather than store tables of a, b, and c for a set of azimuth angles, we stored the first 15 of their Fourier coefficients. The final LUTs contained over 100 Mb of data: 3 × 15 Fourier coefficients for each case of 8 λ's, 100 viewing angles, 33 solar zenith angles, and 12 aerosol models. Systems with such large memory were not readily available at that time. But, by launch (1997) the systems with sufficient memory were available and the GW94 algorithm could be operated on SeaWiFS data in a timely manner. Lastly, the prototype algorithm was only capable of processing three to four pixels per second on a DECstation 3100. This was too slow to even consider processing imagery in real time, a goal for the SeaWiFS mission; however, by launch, processors were fast enough to process the imagery efficiently, and by the end of the mission, the entire SeaWiFS data set could be reprocessed in a matter of weeks. Saved by Moore's law again.
A final development that became an indispensable development tool for SeaWiFS and for future sensors was the SeaDAS computer program developed and maintained at NASA/GSFC which was freely available to all, and allowed processing of ocean color imagery without the significant complications of navigation, calibration, etc., using a personal computer in the comfort of one's office. As I can personally attest, the availability of the source code for this program made it possible to try new approaches to atmospheric correction with a minimum of effort [56][57][58]79,82,93,94].
Here, I have tried to describe the ocean color atmospheric correction algorithm and its development within the context of the time, i.e., the available computational resources. Through a good measure of luck and some insight, the CZCS and follow-on sensors worked sufficiently well to enable ocean color remote sensing to achieve respectability, and become a significant player in global marine-biological and climate-related studies. Additional reflections on the journey from CZCS to SeaWiFS can be found in Gordon [112] and Acker [113]. lem, and Menghua Wang in particular with whom the GW94 algorithm was developed. Third, I acknowledge the key players at NASA/GSFC during the CZCS (post-launch), SeaWiFS and MODIS eras: Gene Feldman, Chuck McClain and Wayne Esaias. They assumed the responsibility of processing the data, maintaining instrument calibration, and shepherding the projects through the concept/design/approval process. Finally, NASA for supporting my research in this area for 35 years.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: