Estimation of near-infrared water-leaving reflectance for satellite ocean color data processing

The atmospheric correction algorithm employed by the NASA Ocean Biology Processing Group requires an assumption of negligible water-leaving reflectance in the near-infrared region of the spectrum. For waters where this assumption is not valid, an optical model is used to estimate near-infrared water-leaving reflectance. We describe this optical model as implemented for the sixth reprocessing of the SeaWiFS missionlong time-series (September 2009). Application of the optical model resulted in significant reductions in the number of negative water-leaving reflectance retrievals in turbid and optically complex waters, and improved agreement with in situ chlorophyll-a observations. The incidence of negative water-leaving reflectance retrievals at 412 nm was reduced by 40%, while negative reflectance at 490 nm was nearly eliminated. ©2010 Optical Society of America OCIS codes: (010.4450) Oceanic optics; (280.4991) Passive remote sensing References and Links 1. H. R. Gordon, and M. Wang, “Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm,” Appl. Opt. 33(3), 443–452 (1994). 2. D. A. Siegel, M. Wang, S. Maritorena, and W. Robinson, “Atmospheric correction of satellite ocean color imagery: the black pixel assumption,” Appl. Opt. 39(21), 3582–3591 (2000). 3. C. R. McClain, E. Ainsworth, R. Barnes, R. Eplee, Jr., F. Patt, W. Robinson, M. Wang, and S. W. Bailey, “SeaWiFS Postlaunch Calibration and Validation Analyses, Part 1,” NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2000). 4. A. Morel, “Optical modeling of the upper ocean in relation to its biogenous matter content (case 1 waters),” J. Geophys. Res. 93(C9), 10749–10768 (1988). 5. R. P. Stumpf, R. A. Arnone, J. R. W. Gould, P. M. Martinolich, and V. Ransibrahmanakul, “A partially coupled ocean-atmosphere model for retrieval of water-leaving radiance from SeaWiFS in coastal waters,” in Patt, F.S., et al., 2003: Algorithm Updates for the Fourth SeaWiFS Data Reprocessing. NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2003). 6. H. R. Gordon, O. B. Brown, R. H. Evans, J. W. Brown, R. C. Smith, K. S. Baker, and D. K. Clark, “A semianalytic radiance model of ocean color,” J. Geophys. Res. 93(D9), 10909–10924 (1988). 7. R. W. Gould, Jr., R. A. Arnone, and P. M. Martinolich, “Spectral dependence of the scattering coefficient in case 1 and case 2 waters,” Appl. Opt. 38(12), 2377–2383 (1999). 8. NIR, html, NIR Correction, http://oceancolor.gsfc.nasa.gov/REPROCESSING/SeaWiFS/R4/NIR.html 9. A. Morel, D. Antoine, and B. Gentili, “Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function,” Appl. Opt. 41(30), 6289–6306 (2002). 10. L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65-2.5 m spectral range,” Appl. Opt. 32(19), 3531–3540 (1993). 11. R. M. Pope, and E. S. Fry, “Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,” Appl. Opt. 36(33), 8710–8723 (1997). 12. Z. Lee, R. A. Arnone, C. Hu, P. J. Werdell, and B. Lubac, “Uncertainties of optical parameters and their propagations in an analytical ocean color inversion algorithm,” Appl. Opt. 49(3), 369–381 (2010). 13. P. J. Werdell, and S. W. Bailey, “An improved in-situ bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. 98(1), 122–140 (2005). 14. A. Bricaud, A. Morel, M. Babin, K. Allali, and H. Claustre, “Variations of light absorption by suspended particles with the chlorophyll a concentration in oceanic (case 1) waters: analysis and implications for bio-optical models,” J. Geophys. Res. 103(C13), 31033–31044 (1998). #122261 $15.00 USD Received 4 Jan 2010; revised 10 Mar 2010; accepted 11 Mar 2010; published 26 Mar 2010 (C) 2010 OSA 29 March 2010 / Vol. 18, No. 7 / OPTICS EXPRESS 7521 15. P. J. Werdell, S. W. Bailey, B. A. Franz, L. W. Harding, Jr., G. C. Feldman, and C. R. McClain, “Regional and seasonal variability of chlorophyll-a in Chesapeake Bay as observed by SeaWiFS and MODISAqua,” Remote Sens. Environ. 113(6), 1319–1330 (2009). 16. S. W. Bailey, and P. J. Werdell, “A multi-sensor approach for the onorbit validation of ocean color satellite data products,” Remote Sens. Environ. 102(1-2), 12–23 (2006).


Introduction
The retrieval of accurate geophysical data products (e.g., spectral remote sensing reflectance R rs (), and concentrations of the phytoplankton pigment chlorophyll-a, C a ) from satelliteborne radiometers requires the use of an 'atmospheric correction' algorithm [1].A critical component of this atmospheric correction algorithm is estimation of spectral aerosol reflectance,  a (), which are subtracted from the total measured signal as part of the process of determining R rs ().The algorithm implemented by the NASA Ocean Biology Processing Group (OBPG) [1] requires an assumption of negligible water-leaving reflectance in the near infrared (NIR) region of the spectrum (i.e., R rs (NIR) = 0 sr 1 ).The NIR bands of the NASA Sea-viewing Wide Field-of-view Sensor (SeaWiFS) are centered on 765 and 865 nm.With this 'black-pixel' assumption, the measured top-of-atmosphere (TOA) reflectance in two NIR bands can be used to estimate both the magnitude and spectral dependence of  a (NIR).Unfortunately, this assumption is rarely valid for waters with significant particulate (e.g., algal and mineral) backscattering [2].If R rs (NIR) can be modeled, however, it can be removed from the TOA signal, allowing the atmospheric correction to proceed in estimating  a (NIR) as if R rs (NIR) were negligible.
The OBPG implemented a bio-optical model to estimate R rs (NIR) as part of the third reprocessing of the SeaWiFS mission in May 2000 (R2000) [2,3].This approach required an initial run through the atmospheric correction process (ignoring R rs (NIR) > 0 sr 1 ) to estimate C a .This preliminary C a was used to estimate spectral particulate backscattering, b bp (NIR), which was used in turn to reconstruct R rs (NIR).With this modeled R rs (NIR) removed from the TOA signal, the atmospheric correction process was repeated, C a was recalculated, the process was iterated upon until convergence in C a was achieved [2,4].With the fourth SeaWiFS reprocessing in July 2002 (R2002), this C a -driven algorithm was replaced with a reflectancebased algorithm [5].The OBPG switched from [2] to [5] after observing that [2] depressed ratios of  a (NIR) in optically complex waters.This often resulted in the selection of spectrally flat aerosol models within the atmospheric correction process, which artificially depressed the final C a retrievals.The initial estimates of C a were unrealistically high in these complex water masses, which led to over-estimates of b bp (NIR) and, subsequently, R rs (NIR).The reflectancebased approach [5] was less prone to such over-correction.In this revised model, R rs () from an initial run through the atmospheric correction process and empirical estimates of the absorption coefficients for particles and dissolved materials at 670-nm were used to estimate b bp (670) via an ocean surface reflectance model [6].A spectral scattering function [7] was then used to derive b bp (NIR) from b bp (670).Finally, R rs (NIR) was reconstructed from b bp (NIR), the atmospheric correction process was repeated, and the process was iterated upon until convergence in R rs (NIR) was achieved.
Despite such efforts to account for non-negligible R rs (NIR), the atmospheric correction algorithm continued to demonstrate sub-optimal performance in highly scattering waters, as manifest by negative R rs (), particularly in the blue region of the spectrum, and erroneous spectral dependence in visible (VIS) part of the spectrum, with the latter leading to inaccurate C a retrievals.To address these issues, a re-evaluation of the implementation of [5] was initiated prior to the sixth SeaWiFS reprocessing in September 2009 (R2009) and several modifications were made to the R rs (NIR) model, which are presented below.First, an inconsistency in the relationship between R rs () and the inherent optical properties (IOP) was eliminated.Second, alternative estimates of the spectral dependence of b bp () and of the absorption coefficients for particles and dissolved materials at 670 nm were adopted.Finally, a revised iteration scheme was implemented.

Inversion of remote-sensing reflectance
Remote sensing reflectance is proportional to the IOPs, specifically the marine absorption and backscattering coefficients, through: where The NIR reflectance model [5] implemented for R2002 and used through SeaWiFS reprocessing 5.2 in 2007 (R2007) suffered from an inconsistent application of Eq. ( 1).First, an approximation of this relationship, where G() was fixed at 0.51 and X() was simplified to b b () / a(), was used to obtain near-infrared IOPs from R rs (VIS).The quadratic form of [6], where G() = [g 1 + g 2 X()], g 1 = 0.0949, and g 2 = 0.0794, was later used to obtain R rs (NIR) from the modeled near-infrared IOPs.While a minor inconsistency, using dissimilar definitions for G() and X() was found to have significant effects in highly scattering waters.
Furthermore, the approach was inconsistent with the methods of [9] that are used elsewhere (and independently) in the atmospheric correction process to normalize R rs () for bi-directional effects.In [9], G() = f() / Q(), where f() is a function of the illumination conditions, water constituents, and sea-state, Q() is the ratio of spectral upwelling irradiance to spectral upwelling radiance, and ratios of f() to Q() are pre-tabulated from extensive radiative transfer simulations.For the R rs (NIR) model implemented in R2009, G() = f() / Q(), using the look up tables (LUTs) of [9], for both the conversion of retrieved R rs (VIS) to IOPs and the conversion of modeled IOPs to R rs (NIR).While these LUTs were originally generated for waters whose light field is predominantly influenced by phytoplankton, they proved suitable for use in this model, as will be demonstrated in Section 4. A full discussion of the applicability of [9] to all water types is beyond the scope of this paper.

Spectral backscattering
Given the relationship in Eq. ( 1), a measured R rs ( 0 ), and knowledge of a( 0 ), an estimate of b b ( 0 ) can be obtained, where  0 is any reference wavelength (e.g., 670 nm).This retrieved b b ( 0 ) can then be extrapolated to the NIR through a spectral model.The spectral dependence of b b () can be expressed as: where Y can take a number of forms, but is commonly parameterized as a power law function: In Eq. (2)b, the  term is typically in the range of [0, 2], with values around 2 observed in open ocean waters and turbid waters approaching zero [12].In contrast, the model employed in R2002-R2007 [3] used a spectral scattering relationship for Eq.(2)a developed by [7]: where the empirically derived coefficients  and  are 0.00113 and 1.62517, respectively.For use in the R rs (NIR) model, the latter relationship required an assumption of equivalent spectral dependency for backscattering and scattering.In the R rs (NIR) model employed in R2009, the more common power-law form of Eq. ( 2)b was adopted.Rather than relying on a fixed  value, however, the empirical relationship developed in [12] was adopted: The NASA bio-Optical Marine Algorithm Data set (NOMAD) [13] includes measurements of b b (VIS) at a number of wavelengths spanning a range of values from 0.0007 to 0.0133 m 1 .Unfortunately, this data set does not include measured b b (NIR).Nevertheless, this data set was used to evaluate the various forms of Eq. ( 2) and 3 using a reference wavelength of 670 nm and a retrieval wavelength of 555 nm.The formulation of [12] had the minimal bias and median percent difference relative to in situ measurements (Table 1, Fig. 1).and performed better than the linear, spectral scattering function [7].

Estimating absorption
In the NIR, a() is dominated by the absorption due to water, a w () [10,11], such that absorption coefficients for particles, a p (NIR), and dissolved material, a g (NIR), can be ignored.However, a reasonable estimation of the absorption coefficient is crucial for the inversion of R rs () to retrieve b b (670), where a p (670) and a g (670) can be a significant portion of the total.
In [5], the estimate of a(670) involved two functions, one based on C a [14] and the other on a green-to-red reflectance ratio [5].In the modified R rs (NIR) model, this was updated to a single, C a -based relationship derived from NOMAD [13] (Type II linear regression using spectroscopy only, limited to C a > 0.2 mg m 3 ; Fig.The a p (670) and a g (670) values contained within the NOMAD data set span a range of 0.00001-0.03m 1 for a p (670) and 0.00057 0.8413 m 1 for a g (670).

Iteration scheme
Since the retrieval of R rs () in highly scattering waters requires knowledge of R rs (NIR), and estimation of R rs (NIR) requires knowledge of R rs (VIS), the R rs (NIR) algorithm must be implemented iteratively.The iteration begins by assuming R rs (NIR) = 0 sr 1 , so that initial R rs () can be derived and a first estimate of R rs (NIR) can be achieved.The iteration converges when R rs () changes by less than 2%.In the majority of cases, the method converges after 3-4 iterations, but the algorithm allows up to 10 iterations.If, however, the initial atmospheric correction results in negative or an otherwise nonphysical spectral distribution for R rs () (criteria based on in situ observations), the iteration is re-initialized assuming  a (NIR) = 0 (i.e., all reflectance except that from Rayleigh scattering and Sun glint originates from the water mass).This is effectively the opposite extreme from the initial condition of R rs (NIR) = 0 sr 1 .The iteration process is then allowed to proceed to convergence as previously described.For the rare cases where the iteration fails to converge after this reset, an additional iteration is forced, again with an assumption of  a (NIR) = 0.This pixel is then flagged with an 'atmospheric correction warning' to alert users to the questionable nature of the pixel, although a qualitatively useful retrieval may still result.
The controls on this iteration scheme differ slightly from previous implementations [8], especially with regards to the how the iteration is reset, which is currently more robust than that of past methods.It should also be noted that the R rs (NIR) is forced to zero if the initial iteration results in a C a < 0.3 mg m 3 .This is done to avoid the addition of noise into clearwater pixels through the estimation of otherwise negligible R rs (NIR) contributions.Furthermore, to avoid the introduction of transitional artifacts, the R rs (NIR) estimation is linearly weighted from 0 to 1 for 0.3  C a  0.7 mg m 3 (Fig. 3).This C a -based weighting of R rs (NIR) was first introduced in R2002 and is effectively unchanged for R2009.

Impact
A long-term satellite and in situ data set was recently compiled to study C a variability in Chesapeake Bay [15].These data were used to evaluate the modifications to the NIR correction described above, as implemented for the most recent SeaWiFS reprocessing.Specifically, the NIR correction implemented with R2002 (and used through R2007) was compared to the method described in this work and employed for R2009.The satellite and in situ data acquisition and processing methods are described in [15].
The revised correction reduced the percentage of negative R rs (412) for the time series of SeaWiFS retrievals over the Bay by nearly half (18.5% vs. 31.9%).The incidence of negative R rs (490) was reduced to negligible levels (Table 2).The satellite estimates of C a showed a substantial improvement in the agreement with ground truth measurements in the turbid and highly productive regions of the Middle and Upper Bay (Fig. 4).The revised correction also improved reflectance retrievals relative to coincident in situ measurements.Following the method of [16], and using a global validation data set compiled by the OBPG, we compared R rs () retrievals for SeaWiFS after limiting the data set to where the NIR correction was used in the processing of the satellite data (e.g., C a  0.3 mg m 3 ).Improvement in the satellite agreement with in situ R rs was evident for all bands, with the greatest improvement seen for the 412 and 443 nm bands (Table 3).

Conclusion
The revised R rs (NIR) model that is presented in this paper has been implemented in the operational processing of satellite ocean color sensor data by the NASA Ocean Biology Processing Group.This revised correction eliminated a number of inconsistencies present in the previous version.A variable estimate for the spectral slope of the backscattering coefficient was implemented in place of a fixed linear relationship.A simplified estimate of absorption in the red region of the visible spectrum was also implemented.This resulted in a significant reduction in the incidence of negative R rs for satellite retrievals in the 412 to 490 nm range, while improving the agreement of satellite-derived R rs () and C a with in situ measurements.
The results and discussion presented here have been focused on ocean color retrievals from the SeaWiFS sensor.It should be noted, however, that the same algorithm has been applied in the processing of MODIS-Aqua and other ocean color sensor data sets archived by NASA's Ocean Biology Processing Group.
The approach presented, while a significant improvement over the previous implementation, cannot correct all cases where the black pixel assumption is invalid.The use of empirical optical models (e.g.spectral backscattering, particulate and dissolved absorption) limits the applicability of this algorithm to waters that are similar to those over which these empirical models were developed.In addition, the practical implementation of the iteration scheme can result in conditions of non-convergence.Fortunately, the optical domains covered by the data sets used in the generation of these algorithms are quite diverse and the cases of convergence failure are small.Therefore, these limitations do not detract from the benefits realized by the use of this algorithm.

G
() is a function of the illumination conditions, water constituents, and sea-state, b b () is the total backscattering coefficient [ b bw () + b bp ()], and a() is the total absorption coefficient [ a w () + a p () + a g ()].The additional subscripts w, p, and g indicate specific contributions by water, particles, and dissolved material, respectively.

Fig. 1 .
Fig.1.Modeled vs. measured bb(555 nm) from the NOMAD data set.Gray symbols are the power law formulation of[12].Red symbols are the linear function of[7].

Fig. 2 .
Fig. 2. Particulate and dissolved absorption at 670 nm vs. Ca derived from the NOMAD data set.

Fig. 3 .
Fig. 3. Global map of likely application of the NIR correction.This was created from the SeaWiFS mission cumulative climatology.Black indicates land; grey Ca < = 0.3 mg m 3 ; white Ca > 0.3 mg m 3 .The white area indicates where the NIR correction is likely to be applied.

Fig. 4 .
Fig. 4. Histograms of in situ and satellite derived Ca for the Lower, Middle and Upper Chesapeake Bay