Detailed validation of the bidirectional effect in various Case I and Case II waters

Simulated bidirectional reflectance distribution functions (BRDF) were compared with measurements made just beneath the water’s surface. In Case I water, the set of simulations that varied the particle scattering phase function depending on chlorophyll concentration agreed more closely with the data than other models. In Case II water, however, the simulations using fixed phase functions agreed well with the data and were nearly indistinguishable from each other, on average. The results suggest that BRDF corrections in Case II water are feasible using single, average, particle scattering phase functions, but that the existing approach using variable particle scattering phase functions is still warranted in Case I water. ©2012 Optical Society of America OCIS codes: (010.4450) Oceanic optics; (010.4588) Oceanic scattering; (010,5620) Radiative transfer; (280.4991) Passive remote sensing. References and links 1. H. R. Gordon and A. Y. Morel, Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery: A Review (Springer-Verlag, New York, 1983). 2. E. J. Kwiatkowska, B. A. Franz, G. Meister, C. R. McClain, and X. X. Xiong, “Cross calibration of ocean-color bands from moderate resolution imaging spectroradiometer on Terra platform,” Appl. Opt. 47(36), 6796–6810 (2008). 3. H. R. Gordon and D. K. Clark, “Clear water radiances for atmospheric correction of Coastal Zone Color Scanner imagery,” Appl. Opt. 20(24), 4175–4180 (1981). 4. T. Hirata, N. Hardman-Mountford, J. Aiken, and J. Fishwick, “Relationship between the distribution function of ocean nadir radiance and inherent optical properties for oceanic waters,” Appl. Opt. 48(17), 3129–3138 (2009). 5. H. Loisel and A. Morel, “Non-isotropy of the upward radiance field in typical coastal (Case 2) waters,” Int. J. Remote Sens. 22(2-3), 275–295 (2001). 6. A. Morel and B. Gentili, “Diffuse reflectance of oceanic waters: its dependence on sun angle as influenced by the molecular scattering contribution,” Appl. Opt. 30(30), 4427–4438 (1991). 7. A. Morel and B. Gentili, “Diffuse reflectance of oceanic waters. II Bidirectional Aspects,” Appl. Opt. 32(33), 6864–6879 (1993). 8. A. Morel and B. Gentili, “Diffuse reflectance of oceanic waters. III. Implication of bidirectionality for the remote-sensing problem,” Appl. Opt. 35(24), 4850–4862 (1996). 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. Z. P. Lee, K. Du, K. J. Voss, G. Zibordi, B. Lubac, R. Arnone, and A. Weidemann, “An inherent-opticalproperty-centered approach to correct the angular effects in water-leaving radiance,” Appl. Opt. 50(19), 3155– 3167 (2011). 11. K. J. Voss, A. Morel, and D. Antoine, “Detailed validation of the bidirectional effect in various Case 1 waters for application to ocean color imagery,” Biogeosciences 4(5), 781–789 (2007). 12. X. Zhang, M. Twardowski, and M. Lewis, “Retrieving composition and sizes of oceanic particle subpopulations from the volume scattering function,” Appl. Opt. 50(9), 1240–1259 (2011). 13. J. M. Sullivan and M. S. Twardowski, “Angular shape of the oceanic particulate volume scattering function in the backward direction,” Appl. Opt. 48(35), 6811–6819 (2009). #160488 $15.00 USD Received 23 Dec 2011; revised 5 Mar 2012; accepted 5 Mar 2012; published 20 Mar 2012 (C) 2012 OSA 26 March 2012 / Vol. 20, No. 7 / OPTICS EXPRESS 7630 14. IOCCG, “Remote sensing of inherent optical properties: fundamentals, tests of algorithms, and applications,” No. 5 (Reports of the International Ocean-Colour Coordinating Group, Dartmouth, Canada, 2006). 15. T. J. Petzold, “Volume scattering functions for selected ocean waters,” SIO Ref. 72–78 (University of California, San Diego, Scripps Institution of Oceanography Visibility Laboratory, 1972). 16. K. J. Voss and A. L. Chapin, “Upwelling radiance distribution camera system, NURADS,” Opt. Express 13(11), 4250–4262 (2005). 17. H. Claustre, A. Sciandra, and D. Vaulot, “Introduction to the special section bio-optical and biogeochemical conditions in the South East Pacific in late 2004: the BIOSOPE program,” Biogeosciences 5(3), 679–691 (2008). 18. S. Pegau, J. R. V. Zaneveld, B. G. Mitchell, J. L. Mueller, M. Kahru, J. Wieland, and M. Stramska, “Ocean optics protocols for satellite ocean color sensor validation, Revision 4, Volume IV: Inherent optical properties: instruments, characterizations, field measurements, and data analysis protocols,” NASA/TM-2003–21621/Rev4Vol IV (Goddard Space Flight Center, Greenbelt, MD, 2002). 19. M. S. Twardowski, E. Boss, J. B. Macdonald, W. S. Pegau, A. H. Barnard, and J. R. V. Zaneveld, “A model for estimating bulk refractive index from the optical backscattering ratio and the implications for understanding particle composition in case I and case II waters,” J. Geophys. Res. 106(C7), 14129–14142 (2001). 20. J. M. Sullivan, M. S. Twardowski, J. R. V. Zaneveld, and C. Moore, “Measuring optical backscattering in water,” in Light Scattering Reviews 7, A. Kokhanovsky, ed. (Springer Praxis Books, to be published). 21. A. L. Whitmire, E. Boss, T. J. Cowles, and W. S. Pegau, “Spectral variability of the particulate backscattering ratio,” Opt. Express 15(11), 7019–7031 (2007). 22. J. Ras, H. Claustre, and J. Uitz, “Spatial variability of phytoplankton pigment distributions in the Subtropical South Pacific Ocean: comparison between in situ and predicted data,” Biogeosciences 5(2), 353–369 (2008). 23. L. Van Heukelem and C. S. Thomas, “Computer-assisted high-performance liquid chromatography method development with applications to the isolation and analysis of phytoplankton pigments,” J. Chromatogr. A 910(1), 31–49 (2001). 24. J. L. Mueller, R. R. Bidigare, C. Trees, W. M. Balch, J. Dore, D. T. Drapeau, D. Karl, L. Van Heukelem, and J. Perl, “Ocean optics protocols for satellite ocean color sensor validation, Revision 5, Volume V: Biogeochemical and bio-optical measurements and data analysis protocols,” NASA/TM-2003(Goddard Space Flight Center, Greenbelt, MD, 2003). 25. H. R. Gordon, “Normalized water-leaving radiance: revisiting the influence of surface roughness,” Appl. Opt. 44(2), 241–248 (2005). 26. J. M. Sullivan, M. S. Twardowski, P. L. Donaghay, and S. A. Freeman, “Use of optical scattering to discriminate particle types in coastal waters,” Appl. Opt. 44(9), 1667–1680 (2005). 27. P. M. Teillet, “Rayleigh optical depth comparisons from various sources,” Appl. Opt. 29(13), 1897–1900 (1990). 28. C. D. Mobley, B. Gentili, H. R. Gordon, Z. H. Jin, G. W. Kattawar, A. Morel, P. Reinersman, K. Stamnes, and R. H. Stavn, “Comparison of numerical models for computing underwater light fields,” Appl. Opt. 32(36), 7484– 7504 (1993). 29. C. D. Mobley, L. K. Sundman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. 41(6), 1035–1050 (2002). 30. H. R. Gordon, “Sensitivity of radiative transfer to small-angle scattering in the ocean: Quantitative assessment,” Appl. Opt. 32(36), 7505–7511 (1993). 31. H. R. Gordon and K. Y. Ding, “Self-shading of in-water optical instruments,” Limnol. Oceanogr. 37(3), 491–500 (1992).


Introduction
The spectral radiance exiting a natural water body (the water-leaving radiance, L,) is dependent on the solar zenith angle and the viewing direction, as well as the optical properties of the water and its constituents.When this radiance is imaged with sensors on earth-orbiting satellites (ocean color remote sensing), every pixel in the image is associated with a different illumination-viewing geometry.Most algorithms that are used to associate water-leaving radiance with geophysical parameters, such as the concentration of chlorophyll a, have been developed based on in situ measurements of upwelling spectral radiance propagating toward the zenith [I], i.e., at a single viewing angle.Thus, it is necessary to be able to relate the water-leaving radiance at any sun-viewing direction to the nadir view.In addition, today there are many earth-orbiting sensors in operation, and each views a given location under different illumination-viewing geometries.Therefore, in order to compare or merge the imagery from different sensors, it is necessary to relate each to a common geometry [2].
Gordon and Clark [3] first defined the normalized water-leaving radiance (nl., i as the value L.. would have if the sun were at the zenith, the atmosphere absent, and the sensor looking at nadir.This estimate was obtained by assuming that the angular distribution of upwelling radiance just beneath the water's surface (the bidirectional reflectance distribution function, BRDF) was independent of the viewing direction, i. ) is currently the standard model used to account for the angular variation of the upwelling radiance in satellite processing.It works well in Case I water [II], but is not suitable for Case II water for at least three reasons.First, the tables end at a maximum chlorophyll concentration of 10 ug I' 1 , but chlorophyll can exceed this in places.Second, the MAG2002 BRDF tables were developed from a radiative transfer (RT) model that used scattering particle phase functions that assumed that all the particles are biogenic, and thereby are applicable only to Case I waters.Third, it requires that the absorption contribution from other optically active agents (e.g.yellow substance) co-vary with concentration of chlorophyll a, which is often invalid for Case II waters.
Radiative transfer computations require the absorption coefficient (a) and the volume scattering function (ß(0).where 0 is the scattering angle), or equivalently, the absorption coefficient, the scattering coefficient (b).and the scattering phase function, (ß(9) = ß(9)lb).Of these, measurements of ß(6) or the phase function have been reported in only a few studies, most of which were reviewed by [12] and references therein.Recently, using novel instrumentation, Sullivan and Twardowski [13] found remarkable consistency in the shape of the backward portion of/? (i.e. for 90° < = 8 < = 180°) measured in situ around the world.This is a promising development, because a "universal" paniculate scattering phase function would simplify radiative transfer modeling in coastal waters.In particular, combining this phase function with a and b b (the back scattering coefficient), which are accessible from the remote sensing signal [14], would enable the BRDF to be determined.
One objective of this study was to assess the BRDF corrections from a recently published algorithm [10], which may apply in either Case I or Case II water.Lee et al.
[10] (referred to below as Lee2011) used in situ measurements from just 3 locations to validate their model; here we used a much larger data set across a wide variety of inherent optical properties.
A second objective of this study was to assess the effects of using different phase functions when modeling the BRDF.MAG2002 and Lee2011 used different phase functions.MAG2002 used a phase function that varied as a function of chlorophyll concentration.In contrast, the phase function used by Lee20l 1 did not vary systematically as a function of the BRDF look-up-table input parameters.However, the synthetic data set used to generate the Lee2011 BRDF look-up-table had been computed with a randomly varying linear mixture of two phase functions (one for mineral particles, one for phytoplankton).One question posed in Lee2011 but not addressed with experimental data was how critical the choice of phase function is for the purposes of a BRDF correction.We investigated that question here by comparison of the Lee20l 1 BRDF correction with those provided by MAG2002 as well as our own corrections generated with a radiative transfer model and the Sullivan and Twardowski [ 13] and Petzold [ 15] turbid-water phase functions.
The overall goal of this work was to assess potential BRDF corrections applicable in Case II water.As mentioned above, two specific objectives were (a) to validate Lee2011 using a large data set of in situ hemispherical radiance measurements and (b) to test the sensitivity of results to the shape of the phase function.For comparison, we also included data from Case I sites and corrections from MAG2002 in the evaluation.Although the focus was on Case II water and the new model Lee2011, the older model MAG2002 provided a useful benchmark, since it had been previously validated in Case I conditions using methods similar to those presented here [II].

Methods
The approach was to compare measured hemispherical upwelling radiance distributions with model simulations.Each comparison consisted of a set of average images measured at 412, 436, 486, 326, and 548 nm with the NuRADS system [16], and sets of results from each of five numerical models.The first set of model output was interpolated from the tables described by MAG2002.The second set of model output was interpolated from the tables described by Lee20l I. Three additional sets of model output were generated using our own RT code (H.R. Gordon, unpublished) with different ß(0) parameterization.Data were available from six field campaigns: Chesapeake Bay, USA in 2004, the BIOSOPE cruise in 2004 [17]

/. NuRADS image processing
The NuKADS system [16] is a compact (30 cm diameter, 30 cm length), multispectral camera that images the upwelling light field in six narrow (-10 nm FWHM) spectral bands centered at 412, 436, 486, 526, 548 and 615 nm.This study did not use data from the 615 nm channel due to the large influence of instrument self-shading at this wavelength.
NuRADS was deployed on the surface with a nadir view direction thereby collecting upwelling radiance -30 cm below the surface.The camera sequentially acquired images in each band with the use of a rotating filter changer.Typical exposure times were less than one second.Acquiring a set of images from all six wavelengths took about two minutes, however, due to the delay reading the data from the CCD in between exposures for each spectral band.A typical deployment lasted from one to several hours, enabling multiple acquisitions over a range of solar zenith angles.
Reduction of the raw NuRADS images consisted of applying calibration factors (see [16]) and averaging images in both space and time to reduce environmental noise.Before averaging, every image was inspected to find the anti-solar point, to correct the geometry of the image, and to check for obstructions in the field of view such as fish, the power/data cable, the side of the ship, or other instruments.Images were then averaged in 10-minute bins, excluding those that had been flagged as unacceptable in the inspection stage.The symmetry of the images about the principal plane was exploited to further average both halves of each image.In addition, spatial binning over 3x3 pixel windows was performed to produce final average images at I x 1 degree resolution.Therefore, each pixel in the processed image over a 10-minute window for a given band could be an average of up to 90 raw pixels (5 images x 2 image halves x 9 pixel window).The mean and coefficient of variation (CV = standard deviation divided by the mean) were computed for each pixel in the processed image from the up to 90 raw pixels in the original images.

IOP data processing
Inherent optical properties (lOPs) measured during the field experiments along with solar zenith angle, calculated for the time of each NuRADS image, were used to index the look up tables (Sections 2.3, 2.5) and parameterize the RT model (Section 2.6).The solar zenith angle i a i. defined as the angle between the vector to the sun and local vertical, was computed based on the time, latitude, and longitude of each NuRADS image.Absorption in,.) and scattering (/)") coefficients of seawaler were interpolated to NuRADS spectral bands from Table 1.1 in [18].The absorption coefficient of dissolved and paniculate constituents (a n ) and the particle total scattering </>,,) and backscattering </>,.,,)coefficients were measured in situ, depth weighted, and interpolated as necessary to the NuRADS spectral bands as follows.
Vertical  ,(Ä, ;) = Bft^X, :) for all of the ac-9 I ac-s spectral bands.Third, we computed depth-weighted absorption a".»riskuAX).scattering ^»«"wU), and backscattering V«^) coefficients using Eq.(1) for all of the ac-9 I ac-s spectral bands.Finally, the depth-weighted values were interpolated to the NuRADS wavelengths.Linear interpolation was used for a n »«"tan/.A power law of the form aX' was used to interpolate bp.^^ and iv«,,,*,«/.Note, all further references to a n , bp, or biy imply the depth-weighted value of those variables interpolated to a particular NuRADS spectral band.

Total chlorophyll concentration [Chi] (ug l 1 ) represents the sum of the concentrations of the following suite of pigments: chlorophyll a, divinyl chlorophyll a, chlorophyllid a, and chlorophyll a allomers and epimers. Calculation of [Chi] was determined via High-Performance Liquid Chromatography (HPLC).
Calculation of [Chi] from the BIOSOPE cruise used a slightly modified version (see [22]) of the method initiated by [23].Data for the other cruises were processed following NASA protocols [24].All of the Chesapeake water samples were acquired at 0.5 m depth.Water samples for the other cruises were acquired within 2 m of the surface and, for some stations, at depths as great as 20 m.Stations for which a depth cast was available were optically weighted using Eq.(I) to generate a single [Chi] value.

Lee20I I table interpolation
Lee2011 provided look-up tables for dimensionless parameters G" and G t , which relate remote sensing reflectance to the water and particle plus dissolved matter absorption and backscattering coefficients via: GS(0"6 m ,4)+G>(0"0 m ,t)2!-p- where * = a, + b* and G" and G, arc empirical coefficients determined separately for water, C. and particles, G'\ with dependence on solar zenith angle, above-water view angle (9"), and azimuth.Tables of coefficients G were defined for 9, < = 75°.Remote sensing reflectance, ft,,. is defined as the ratio of the water leaving radiance, /.".u> ihe downwelling irradiance, /.',,.measured just above the surface.
( 4 ) 0.985 »m 2 to within a few percent for 0 m < 60° and to better than 10% for 9" = 70°.In Eq. (4), T/,9 m ) is the Fresnel transmittance for a flat air-water interface and m is the refractive index for water, taken tobe 1.341.Solar zenith angle, a", a n , b hw , and b^ measured and processed as described in Section 2.2 were used to retrieve R n at regularly spaced angles (#",, ^) via interpolation of the Lee20l 1 look-up tables (i.e.R" was derived from IOPs and the tables, not measured directly in the field).R n (0 w f) was converted to normalized L^(0 V , f) for comparison with the NuRADS images.As mentioned in Section 2.2, no b^ measurements were made during the Chesapeake cruise, yet b^ was required for Eq.(2).To address this limitation, we assumed B p = 2%, a common value for coastal waters [21,26], and computed b^ = B^b^ To check the validity of this assumption, we ran the same analysis using B p = \% and 3% and computed normalized L* from Eqs. (2-4) using the measured b r and B p = 1% and 3%.The normalized Z« values for B p = 1% and 3% agreed within 3.8% for 95% of the Chesapeake data set and within 5.6% for 99% of the Chesapeake data set.The observation that truncating the phase function does not affect the calculated upwelling radiance significantly is consistent with Mobley et al. [29], who concluded that the exact shape of the phase function in the forward direction was not critical as long as b hp was correct.Furthermore, Gordon [30] showed thai ß(0) truncated at scattering angles < 15 degrees (such as the one used here) did not alter calculations of irradiances by more than a few percent.

Our
Note that Sullivan and Twardowski [13] used a normalization in the backward direction only, i.e. their measured ß^O) were normalized by b^ and not by b p as had been done in [28].
Therefore, the ß p (0) for method 2 satisfies the normalization: Method 3 also used the Petzold [ 15] turbid-water ß{ 0), but rather than scaling it by />". as in method I, it was normalized in the backward direction, as in method 2. Thus, ß (0) for method 3 satisfied Eq. ( 6) and was multiplied by the depth-weighted, spectrally-interpolated bt,,, to produce the ß{0) input to the RT calculations.The ß(0) used for method 3 has the same shape as ß(6) used for method 1, but, like method 2, forces b^ in the model to agree with the measurements.

Each of the radiance distributions was normalized by its value at nadir. Comparison between the modeled and measured normalized radiance distributions was performed by computing the model-data difference at every 5 degrees in nadir from
^,(0.0) J The NuRADS shadow affected some of the view angles, which were subsequently discarded based on an estimate of the shadow contamination.Instrument self-shading was estimated using the Gordon and Ding [31] model modified to incorporate the threedimensional shape of NuRADS rather than the flat circular disk used in the original calculations.For each view angle, the path length through the instrument shadow (</,) was estimated and used with the total absorption (a,) to compute the transmission < T) through the shadow: T = e\p(-a,d,).For most of the data, angles for which the transmission was less than 95% were discarded from the comparison.The total absorption coefficient in the Chesapeake Bay data set, however, was so high (mean +/-std.dev.■ 3.1 +/-1.3 rrf') that almost the entire data set was flagged as shadow.Therefore, the shadow threshold for the Chesapeake data set was reduced to 90% transmission.
After omitting points thought to be contaminated by instrument shadow, the distributions of the remaining D(0,, </» were plotted as grouped by chlorophyll, *t, and 0,.These plots were constructed using all images in the data set and also after dividing the images into "Case I" and "Case II" subsets based on the visual relationship between a" or b p and chlorophyll, as described below.The entire data set contained 1474 averaged images, 80% of which were acquired in Case I conditions.The example shown (Fig. 2) is for 526 nm, which had 319 averaged NuRADS images.Typically, many images were acquired at each IOP sampling station, which is why only 112 combinations of u pf and b p were plotted in Fig. 2.

Model-data comparison
The number of average NuRADS images processed and matched with the corresponding runs from our RT model was 265, 284, 290, 319, and 316 in each of the five NuRADS bands (412, 436, 486, 526, and 548 nm, respectively).Of these, only 233, 254, 257, 284, and 279 in the respective bands were matched with the M AG2002 tables because the remainder, all of which were from either the Chesapeake or Monterey Bay data sets, had total chlorophyll values greater than 10 ug /1 (the upper limit available from MAG2002).Also, only 263, 279.284, 306, and 305 images in the respective bands were matched with the Lee2011 tables because the remainder were acquired with solar zenith angles > 75°.For Case I conditions, the correlation between the MAG2002 modeled L u (0 r , «*) / /."(().0) and the NuRADS data was similar to that between the Lee2011 model and the NuRADS data (Fig. 3A, B).The Lee2011 model more closely agreed with the NuRADS data than MAG2002 in Case II, however (Fig. 3C, D).For these Case I data (Fig. 3A, B), the leastsquares fits fell closer to the 1:1 line than the eye might predict because there were so many overlapping points.At 412 nm there were 22,185 matched points and at 526 nm there were 28,491; other bands fell within this range.Aggregate plots such as Fig. 3 gave an overview of the results, but details of the models' performance were apparent when looking at subsets of the data, which we did using box plots.
Box plots show the distribution of model-data differences, />(#,, #), across the observed chlorophyll levels (Figs. 4, 5), azimuths (Figs. 6, 7), and solar zenith angles (Figs. 8, 9).The plots were prepared for all bands but are only shown here for the 526 nm channel because the The population (AO of each distribution is given along the top of each plot.Due to the data acquired with chlorophyll values greater than 10 ug 1"', the number of model-data differences was not always the same for the different models.When two values of N are reported, the top one refers to the number of model-data differences for the MAG2002 look-up table and the bottom one refers to the number of model-data differences for the Lee2011 tables and our three RT model calculations.If only one value of N is reported, the number of model-data differences was the same for all five models. An area around Dtfh. ^ = 0 in the box plots has been shaded to show +/-the mean coefficient of variation of the NuRADS images at the points used to compute />(#,, <j>).Section 2.1 describes how the CV of the NuRADS images was calculated.The CV gives a measure of the limit to which the data could agree with even the best model, given environmental variability, due to effects such as wave focusing within the images.At chlorophyll concentrations > 1.0 ug I" 1 in the Case I waters we sampled, all of the models had similar distributions of £>(0" f)-It was n °l obvious that one model performed consistently better than the others for chlorophyll concentrations > 1.0 ug I"' (Fig. 4).
Likewise, in the samples taken from Case II waters, it was not obvious that one model performed consistently better than the others (Fig. 5).In the chlorophyll bins for which relatively large sample sizes were available, namely 1.78 -3.15, 5.62 -10.0, 10.0 -17.8, and 31.6 -56.2 ug I ', all of the models had close to zero median and similar extreme values of D(0," 4»).In the chlorophyll bins with smaller sample sizes, 1.00 -1.78, 17.8 -31.6, 56.2 -178, and > 178 ug I ', all of the models exhibited similar biases and extreme values of approximately equal magnitude.The only exceptions to these generalizations were MAG2002 and, arguably, Lee2011, in the 5.62 -10.0 ug I" 1 bin.These differences, however, were not part of any trend.MAG2002 and Lee2011 performed just as well as the others in the bin below (1.78 -3.16 ug I" 1 ) and the bin above (10.0-17.8 ug I" 1 ) the one in question.
In Case I conditions, all of the models had positive biases in and largest extreme values of 0(0,, 4) at small azimuth angles and minimum model-data differences between 105° -120° azimuth (Fig. 6).At larger azimuth angles, from 135° to 180°, the model-data agreement for MAG2002 and Lee2011 did not decrease much relative to their optimum angles between 105° -120° degrees, but the biases and largest model-data deviations for the other models steadily increased (Fig. 6).
MAG2002 had the smallest bias and smallest extreme values over the widest range of azimuth angles of any of the models tested (black boxes in Fig. 6).Our RT model using either the Sullivan-Twardowski (red boxes) or Petzold (blue boxes) phase functions scaled to give correct b^ values had approximately the same bias, though slightly larger extreme values, as MAG2002 for small azimuth angles, but they had notably larger biases than MAG2002 at large azimuth angles (Fig. 6).Conversely, Lee2011 had a larger bias than MAG2002 at small azimuth angles, but the two models had very similar performance at azimuth angles > = 90°.Our RT model using the Petzold phase function scaled to give correct />,, values had larger biases and extreme values than MAG2O02 for both small and large azimuth angles (green boxes in Fig. 6).
One way in which the Case II data differed from the Case I results, however, was that MAG2O02 exhibited the largest biases of all the models, which were found from 0° -45° azimuth.A second way in which the Case II data differed from the Case I results was that the range of azimuth angles over which there was good model-data agreement was much larger for the Case II data.For example, both Lee20ll and our RT model using the Sullivan-Twardowski phase function produced median values of D(0,.^) that were within the average NuRADS coefficient of variation across all azimuth angles (Fig. 7).Third, for the Case II data, the range of extreme values of O(0" f) fell between -10% to + 20% of the upwelling radiance at nadir for azimuth angles 45° to 180°, whereas in Case I the smallest extreme values were about +/-20% only between azimuth angles of 90° to 135°.Finally, between 0 and 45° the largest differences were positive in Case I and negative in Case II.
Large negative biases were apparent in all models at small solar zenith angles because the NuRADS data were normalized to nadir (Fig. 8).At small solar zenith angles, near-nadir measurements are reduced by the instrument shadow thereby artificially increasing values at other angles when the data were normalized to nadir.For the Case I data, the near-nadir shadow bias was eliminated for MAG2002 and Lee2011 for 6, > = 20° and for our RT model for 0, > = 25°.All models most closely agreed with the data at solar zenith angles between 25° to 35° (Fig. 8).At 6, > = 40° a positive bias in D(0,, 4) and highly skewed distributions were observed for all models (Fig. 8).The median value of D(0,, *>) computed with MAG2002 exhibited the least bias for 0, > = 40°, falling within the average NuRADS coefficient of variation for all solar zenith angles between 15° to 70°.Lee2011 and our RT model using the Sullivan-Twardowski phase function produced median values of O(0" ^) that were within the average NuRADS coefficient of variation for solar zenith angles as large as 65°.

Chi (ug/l)
Hg. 4. Summary of IKO., 4) as a function of chlorophyll (Chi) for the 526 nm NuRADS images in Case I water.See text for a description of box plots and of the population sizes < V i Prom left lo right within each group of boxes: Black boxes correspond to values from MAG2002 minus the data.Light blue boxes correspond to values from Lec20l I minus the data.Green, red.and dark blue boxes correspond to the RT model using VSF parameterization methods 1-3.respectively, minus the data.The shaded grey area around 0 difference is the mean coefficient of variation of the NuRADS images at the points used to compute IX8" #).The Case II results, plotted against solar zenith angle (Fig. 9), show much greater modeldata agreement as well as consistency among the models than the data from Case I conditions (Ftg.8).The only examples of median D(6\, of) falling outside the range of the average correction factors from Lee2011 and our method 2, using the Sullivan and Twardowski [13] phase function, each matched the NuRADS data as well as MAG2002.All of the models matched the NuRADS data well for azimuth angles near 90°.It is not surprising that MAG2002 worked well in Case I waters, as that had been shown previously [II].The new results here, for Case I water, are that the other, models with simpler parameterization for their particle phase functions worked nearly as well as MAG2002 under some conditions.
For Case II conditions, in contrast, the other models tested matched the NuRADS data at least as well as the corrections from MAG2002.Lee2011 and all three versions of our RT code were hardly distinguishable in Case II conditions across the entire range of chlorophyll concentrations (Fig. 5), azimuth angles (Fig. 7), and solar zenith angles (Fig. 9) sampled by the NuRADS data.The BRDF correction generated from MAG2002 worked well in Case II water, up to the 10 ug I ' limit of the tables, except for small azimuth angles (Fig. 7).Detailed results have been presented for only the 526 nm data (Figs.4-9), but the patterns observed were consistent across all spectral bands sampled (Fig. 3).
Use of the MAG2002 tables without Raman scattering would not likely have changed these results significantly.MAG2002 results with and without Raman scattering differed for/ only for wavelengths > 600 nm, which were not considered in this study, and only by a few percent for Q The agreement of our RT methods 1-3 with the data revealed that both the total backscattcring (i.e. the value of b^) and the detailed shape of the backward portion of the paniculate phase function affected the BRDF correction in Case I water.At small azimuth angles, methods 2 and 3, in which the value of />,,,.matches measurements, both agree more closely with the data than method 1, in which the value of />,, matches measurements (Fig. 6).This result indicates that the magnitude of b^ is more important for BRDF corrections at small azimuth angles than the detailed shape of the backward portion of the paniculate phase function.
At large azimuth angles, however, method 2, using the Sullivan-Twardowski phase function, agreed better with the data than either method using the Petzold phase function (Fig. 6).This result indicates that having the magnitude of b^ correct was necessary but not sufficient to match the data at large azimuth angles, and that the shape of the paniculate phase function in the backward direction was also important.The Sullivan-Twardowski and Petzold phase functions differ by about 5% near 145° and 10% at 170° (Fig. 1

inset).
In Case II waters, the differences between our RT methods 1-3 were much smaller than in Case I (Figs.6, 7) suggesting that multiple scattering dampened out the most of the effects of the phase function shape for determining the BRDF.
The observation that the largest model-data discrepancies were close to the principal plane (Figs. 6, 7) is relevant to correcting ocean color satellite imagery for bidirectional reflectance effects because the area of largest errors is also the portion of satellite images most affected by sunglint.Ocean color processing typically discards data collected within and close to sunglint.Therefore, Moderate Resolution Imaging Spectroradiometer (MODIS) and Seaviewing Wide Field-of-view Sensor (SeaWiFS) data that have passed quality control checks for sunglint will have been acquired over azimuth angles with the smallest model-data BRDF discrepancies.The results of this study suggest that operational correction of remotely sensed data for bidirectional reflectance effects could be performed with two models, depending on the water type.Lee2011 [10] was superior to MAG2002 [9] in Case II water at small azimuth angles and for chlorophyll concentrations greater than the 10 ug I" 1 limit of MAG2002.MAG2002 [9], however, remained superior to Lee20l 1 [10] for Case I water.

Conclusion
e., they assumed a completely diffuse BRDF.Previous efforts to model / understand the actual BRDF [4-10] have produced #160488 -$15.00USD Received 23 Dec 2011; revised 5 Mar 2012; accepted 5 Mar 2012; published 20 Mar 2012 (C) 2012 OSA the general insight that the variation can be modeled, at least in Case 1 waters (those for which the optical properties covary with the chlorophyll concentration [ I),).The Morel et al. [9] model (referred to below as MAG2002 , Monterey Bay, CA, USA, in October 2006, the Ligurian Sea in both October 2008 (LSCV'08) and March 2009 (BP'09), and the New York Bight in 2009.
profiles of the paniculate and dissolved absorption, a"(k, z), and panicle scattering b p (k, z) coefficients were measured at 9 wavelengths (JL) using WET Labs ac-9 instruments during the Chesapeake, BIOSOPE, LSCV08, BP'09, and New York Bight experiments and at 83 wavelengths using a WET Labs ar-s in Monterey Bay.Panicle backscattering coefficients b, r U. z) were measured with ECO-BB3 instruments during BIOSOPE, Monterey Bay, LSCV'08, and New York Bight cruises and with a HOBI Labs Hydroscat-6 instrument for BP'09.No b^ measurements were made during the Chesapeake cruise.Quality control and calibrations were applied to these data sets by the investigators who acquired them [19, 20].Processing the calibrated a Pf (X, z), b p (X, z), and b^X, z) followed the same four steps for each cast in each of the data sets.First, we computed depth-weighted values b hp ." it i ttrt i(X mi ,) #160488 -$15.00USD Received 23 Dec 2011; revised 5 Mar 2012;accepted 5 Mar 2012; published 20 Mar 2012 (C) 2012 OSA and bp.^fiuJiKiitu) using Eq.(1) for the shortest wavelength data (Ä«*) from the backscattering sensor and the closest corresponding ac-9 or ac-s band (X rkl ") in each experiment The X^ I X iio " bands were 470 / 488 nm, 450 / 451 nm, 462 / 488 nm, 442 / 440 nm, and 462 / 488 nm in the BIOSOPE, Monterey Bay, LSCV08, BP'09, and New York Bight data sets, respectively.r~ vO.z)exp(-(a, W. z)+K (X, z))z)dz J o exp(-(a,(/l.z)+ * h (/l,z))2Vfe In Eq. (1), v{\, z) is the depth-dependent variable being averaged (either a n {X, z), b p (X, z), or b^X, z)), Vwtighmlik) is the depth-weighted value of that variable, a,(X., z) and b^X, z) are the total depth-dependent absorption («,(/, ;) ■ a"(X, z) + aJ,X)) and backscattering {b^iX.z) = bhp(X.z) + /»/»</.))coefficients, respectively, and Zmax >s the maximum depth of the cast.Second, we used the backscattering ratio, B p ■ b^^^X*,,,) I b p .wrif[klt ji,X, l , Klt ) and the assumption that 1< did not vary with wavelength [19, 21] to compute K, provided look-up tables for// Q as a function of wavelength, solar zenith angle, chlorophyll concentration, and view angle.One set of MAG2002 look-up tables included the effects of Raman scattering whereas a second set of tables did not.The results here used the tables that did include Raman scattering.The dimensionless parameter/is proportional to the irradiance reflectance; / = (a, I b^) * (£, / Ej), where E" and E d are the upwelling and downwelling irradiances just below the water surface, respectively.The bidirectional function Q is defined as the ratio of the upwelling irradiance to radiance; Q = (E m l LJ, where L" is the upwelling radiance just below the water surface.Total chlorophyll concentration and solar zenith angle, calculated for the time of each NuRADS image, were used to retrieve / / Q at regularly spaced in-water view (0,,) and azimuth ($ angles via interpolation of the MAG2002 look-up tables.Azimuth is defined relative to the principal plane, the plane containing the anti-solar point, local vertical, and the sun.The tables were not interpolated spectrally because they had been computed at wavelengths close to the NuRADS bands (412.5,442.5,490,510, 560 nm).Normalizing// Q at each view angle by the value at nadir (&,.= 0, <p = 0), as described in Section 2.7, produced a normalized radiance distribution.L U (0 V .<f>) I /,"((), 0), for comparison with the normalized NuRADS images.#160488 -$15.00USD Received 23 Dec 2011; revised S Mir 2012; accepted S Mar 2012; published 20 Mar 2012 (C) 2012 OSA

0
RT code required the following input parameters: the IOPs and solar zenith angle described in Section 2.2, the Rayleigh optical depth of the atmosphere, computed from [27], and a particle volume scattering function, /!,,<<!).Varying the choice and / or normalization of ßp(0) generated three sets of model output for comparison with each NuRADS image.Method I used the Petzold [15] turbid-water ß p (0) as normalized in [28] to produce the Petzold turbid-water phase function, ß p (9).Note that this ß p (9) satisfies the normalization: #160488 -$15.00USD Received 23 Dec 2011; revised S Mar 2012; accepted S Mar 2012; published 20 Mar 2012 (To recreate the proper ßJO) for input to the RT code under method 1 we multiplied ß p (0) by the depth-weighted, spectrally-interpolated />,,.Note that, due to the acceptance angle of the ac-9 (0.93 degrees), our depth-weighted b p underestimates the value that would be measured by a perfect instrument for the Petzold phase function by 20-30% [15].Method 2 used the Sullivan and Twardowski [ 13] average phase function.The MASCOT instrument used in [13] collected data over the range from 10 to 170 degrees, but Sullivan and Twardowski extrapolated the average phase function in the backwards direction from 170 to 180 degrees [13].Extrapolation of the phase function in the forward direction is more difficult; therefore we performed some initial experiments to test the importance of the exact shape of the forward direction on the BRDF.These tests consisted of a set of model runs at 412 nm, over a range of«,,,, and b p from 0.01 to 1.0 m ' and solar zenith angles 0, 30, and 60 degrees, in which we truncated the Petzold phase function at 2, 5, and 10 degrees.The normalized upwelling radiance distributions created from the phase function truncated at cither 5 or 10 degrees were always within +/-3% of the corresponding result using the phase function truncated at 2 degrees.Therefore, we did not extrapolate the Sullivan and Twardowski phase function in the forward direction, instead using it truncated at 10 degrees (Fig. I).

-
properßJO) for input to the RT code under method 2 we multiplied ß p (0) by the depth-weighted, spectrally-interpolated b hp .
Sullivan-Twardowski and Petzold turbid water phase functions used in methods 1-3.The phase function for method 1 is normalized over all angles (Eq.(5)) and multiplied by ft, to generate ßjO) in the RT model.The phase functions for methods 2 and 3 are normalized in the backward direction (Eq.(6)) and multiplied by Kp to generate ß/8) in the RT model.Inset shows details of methods 2 and 3 over the backscattering directions.

5 to 45 degrees and every 15 degrees in azimuth from 0
to 180 degrees.Thus, the difference, D, between each model output, / .,,. and the corresponding average NuRADS image, L^,, was computed at 117 angles based on Eq.
The optical properties of Case I waters are determined by the chlorophyll concentration [I].Our data set was divided into Case I and Case II subsets, therefore, by visual inspection of #160488-$15.00USD (O2012OSA Received 23 Dec 2011; revised 5 Mar 2012; accepted S Mar 2012; published 20 Mar 2012 plots of inherent optical properties as a function of chlorophyll concentration (Fig. 2).The IOP model of [6] (their Eq. (7-10) was used as a reference during this process.Although the division between Case I and II was done visually, in effect the criteria for classification as a Case II site were (a) [Chi] > 10 ug I" 1 , or (b) b p > 1 nf', or (c) a" > 0.2 m"'.

Fig. 2 .
Fig. 2. Depth-weighted b, (top) and a" (bottom) at 526 nm plotted against total chlorophyll concentration.Left plot« on both the top and bottom show details near the origin.Solid lines plot the IOP model used as a reference for Case I [6].Solid dots and open circles represent conditions designated as Case I and Case II.respectively.

#
160488-SI 5.00 USD (Q2012OSA Received 23 Dec 2011; revised 5 Mar 2012; accepted 5 Mar 2012; published 20 Mar 2012 general trends were the same at all wavelengths (see Fig. 3).Box plots depict the distribution of model-data differences as follows: the box shows the interquartile range (the middle 50% of the data); the circle with dot shows the median model-data difference; thin lines known as "whiskers" extending out of the box show the full range of data values.Some of the distributions of D{0" D were highly skewed or contained extreme outliers; box plots were useful for visualizing such features.
U9.. 4) I U0, 0) for MAG2002 vs. NuRADS data (A and C) and Lee20l I vs. NuRADS data (B and D).The plots for each of the 5 wavelengths have been offset vertically by 1.0 for clarity.Solid lines on each plot show the I: I slope and dashed lines show a linear fit to the data.Note that there is no obvious spectral dependence of these results.For chlorophyll concentrations < 1.0 ug 1 *', MAG2002 had both the smallest overall range of D(0 t .*f), i.e. the least extreme values, and the smallest interquartile range (Fig. 4).Furthermore, over this range MAG2002 generally had median values closest to zero, although the median D(0,, f) for Lee2011 was as close to zero as MAG2002 for chlorophyll concentrations < 0.18 ug 1"' (Fig. 4).#160488-$15.00USD (O2012OSA Received 23 Dec 2011; revised S Mar 2012; accepted S Mar 2012; published 20 Mar 2012 Two questions addressed were: How well does the recent Lee20l I BRDF correction agree with a wide suite of in situ measurements?Is the choice of phase function used in Lee20l I a limitation to its broad implementation?The answers to these questions differed in Case I and Case II water.In Case I conditions, MAG2002 produced the closest agreement with field measurements of the models considered in this study.Relative to earlier work [II] that also validated the accuracy of MAG2002 in Case I water, this study added additional data sets in other parts of the world and five competing models for comparison.The fact that MAG2002, using a variable particulate phase function, produced the best results and that the other models with «160488 -$15.00USD Received 23 Dec 2011; revised S Mar2012; accepted 5 Mar 2012; published 20 Mar 2012 (C) 2012 OSA different paniculate phase functions did not always produce similar output suggests that in Case I water the shape of the paniculate phase function is important for modeling BRDF of the ocean.Specifically, as indicated by the results of our RT methods 1-3, the detailed structure of the backward portion of the paniculate phase function was import.miat large azimuth angles, and the relative amount of forward to backwards scattering was important at small azimuth angles.In Case II conditions, several models produced similar results; of the ones considered, Lee2011 is probably the most convenient to implement, since published tables are available already.Lee2011 compared their model with two NuRADS images, but this study validated the model with more than 1000 NuRADS images across a wide range of absorption and scattering coefficients, solar zenith angles, and wavelengths.Models with different, but fixed, phase functions produced similar output, all in agreement with the data, suggesting that in Case II water the precise shape of the phase function is not important for modeling BRDF of the ocean and therefore is also not a limitation to the adoption of Lee20ll in Case II conditions.