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
The spectral radiance exiting a natural water body (the water-leaving radiance, Lw) 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.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 .
Gordon and Clark  first defined the normalized water-leaving radiance (nLw) as the value Lw 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.e., they assumed a completely diffuse BRDF. Previous efforts to model / understand the actual BRDF [4–10] have produced the general insight that the variation can be modeled, at least in Case I waters (those for which the optical properties covary with the chlorophyll concentration , ).
The Morel et al.  model (referred to below as MAG2002) 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 , but is not suitable for Case II water for at least three reasons. First, the tables end at a maximum chlorophyll concentration of 10 μg l−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 (β(θ), where θ is the scattering angle), or equivalently, the absorption coefficient, the scattering coefficient (b), and the scattering phase function, (). Of these, measurements of β(θ) or the phase function have been reported in only a few studies, most of which were reviewed by  and references therein. Recently, using novel instrumentation, Sullivan and Twardowski  found remarkable consistency in the shape of the backward portion of β (i.e. for 90° < = θ < = 180°) measured in situ around the world. This is a promising development, because a “universal” particulate scattering phase function would simplify radiative transfer modeling in coastal waters. In particular, combining this phase function with a and bb (the backscattering coefficient), which are accessible from the remote sensing signal , would enable the BRDF to be determined.
One objective of this study was to assess the BRDF corrections from a recently published algorithm , which may apply in either Case I or Case II water. Lee et al.  (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 Lee2011 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 Lee2011 BRDF correction with those provided by MAG2002 as well as our own corrections generated with a radiative transfer model and the Sullivan and Twardowski  and Petzold  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 .
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, 526, and 548 nm with the NuRADS system , 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 Lee2011. Three additional sets of model output were generated using our own RT code (H.R. Gordon, unpublished) with different β(θ) parameterization. Data were available from six field campaigns: Chesapeake Bay, USA in 2004, the BIOSOPE cruise in 2004 , 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.
2.1. NuRADS image processing
The NuRADS system  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 ) 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 3 x 3 pixel windows was performed to produce final average images at 1 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.
2.2. IOP data processing
Inherent optical properties (IOPs) 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 (θs), 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 (aw) and scattering (bw) coefficients of seawater were interpolated to NuRADS spectral bands from Table 1.1 in . The absorption coefficient of dissolved and particulate constituents (apg) and the particle total scattering (bp) and backscattering (bbp) coefficients were measured in situ, depth weighted, and interpolated as necessary to the NuRADS spectral bands as follows.
Vertical profiles of the particulate and dissolved absorption, apg(λ, z), and particle scattering bp(λ, z) coefficients were measured at 9 wavelengths (λ) using WET Labs ac-9 instruments during the Chesapeake, BIOSOPE, LSCV’08, BP’09, and New York Bight experiments and at 83 wavelengths using a WET Labs ac-s in Monterey Bay. Particle backscattering coefficients bbp(λ, 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 bbp 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 apg(λ, z), bp(λ, z), and bbp(λ, z) followed the same four steps for each cast in each of the data sets. First, we computed depth-weighted values bbp-weighted(λmin) and bp-weighted(λclose) using Eq. (1) for the shortest wavelength data (λmin) from the backscattering sensor and the closest corresponding ac-9 or ac-s band (λclose) in each experiment. The λmin / λclose bands were 470 / 488 nm, 450 / 451 nm, 462 / 488 nm, 442 / 440 nm, and 462 / 488 nm in the BIOSOPE, Monterey Bay, LSCV’08, BP’09, and New York Bight data sets, respectively.
In Eq. (1), v(λ, z) is the depth-dependent variable being averaged (either apg(λ, z), bp(λ, z), or bbp(λ, z)), vweighted(λ) is the depth-weighted value of that variable, at(λ, z) and bbt(λ, z) are the total depth-dependent absorption (at(λ, z) = apg(λ, z) + aw(λ)) and backscattering (bbt(λ, z) = bbp(λ, z) + bbw(λ)) coefficients, respectively, and zmax is the maximum depth of the cast.
Second, we used the backscattering ratio, Bp = bbp-weighted(λmin) / bp-weighted(λclose) and the assumption that Bp did not vary with wavelength [19, 21] to compute bbp(λ, z) = Bpbp(λ, z) for all of the ac-9 / ac-s spectral bands. Third, we computed depth-weighted absorption apg-weighted(λ), scattering bp-weighted(λ), and backscattering bbp-weighted(λ) coefficients using Eq. (1) for all of the ac-9 / ac-s spectral bands. Finally, the depth-weighted values were interpolated to the NuRADS wavelengths. Linear interpolation was used for apg-weighted. A power law of the form αλ-γ was used to interpolate bp-weighted and bbp-weighted. Note, all further references to apg, bp, or bbp imply the depth-weighted value of those variables interpolated to a particular NuRADS spectral band.
2.3. Chlorophyll processing
Total chlorophyll concentration [Chl] (μg 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 [Chl] was determined via High-Performance Liquid Chromatography (HPLC).
Calculation of [Chl] from the BIOSOPE cruise used a slightly modified version (see ) of the method initiated by . Data for the other cruises were processed following NASA protocols . 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. (1) to generate a single [Chl] value.
2.4. MAG2002 table interpolation
MAG2002 provided look-up tables for f / 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 f is proportional to the irradiance reflectance; f = (at / bbt) * (Eu / Ed), where Eu and Ed 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 = (Eu / Lu), where Lu 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 f / Q at regularly spaced in-water view (θv) 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 f / Q at each view angle by the value at nadir (θv = 0, ϕ = 0), as described in Section 2.7, produced a normalized radiance distribution, Lu(θv, ϕ) / Lu(0, 0), for comparison with the normalized NuRADS images.
2.5. Lee2011 table interpolation
Lee2011 provided look-up tables for dimensionless parameters G0 and G1, which relate remote sensing reflectance to the water and particle plus dissolved matter absorption and backscattering coefficients via:
where κ = at + bbt and G0 and G1 are empirical coefficients determined separately for water, Gw, and particles, Gp, with dependence on solar zenith angle, above-water view angle (θva), and azimuth. Tables of coefficients G were defined for θs < = 75°. Remote sensing reflectance, Rrs, is defined as the ratio of the water leaving radiance, Lw, to the downwelling irradiance, Ed, measured just above the surface.
Normalizing Rrs at each view angle by the value at nadir (θva = 0, ϕ = 0) produced a normalized above-water radiance distribution, Lw(θs, θva, ϕ) / Lw(θs, 0, 0). Note that Ed is canceled when Rrs is normalized in this way. Correction for the reflection-transmission properties of the air-sea interface was necessary to compare normalized Lw with normalized Lu, which is the quantity measured by NuRADS:
The dimensionless factor in principle depends on both θs and θva, but, as shown in , can be approximated by:Eq. (4), Tf(θva) is the Fresnel transmittance for a flat air-water interface and m is the refractive index for water, taken to be 1.341.
Solar zenith angle, aw, apg, bbw, and bbp measured and processed as described in Section 2.2 were used to retrieve Rrs at regularly spaced angles (θva, ϕ) via interpolation of the Lee2011 look-up tables (i.e. Rrs was derived from IOPs and the tables, not measured directly in the field). Rrs(θva, ϕ) was converted to normalized Lu(θv, ϕ) for comparison with the NuRADS images. As mentioned in Section 2.2, no bbp measurements were made during the Chesapeake cruise, yet bbp was required for Eq. (2). To address this limitation, we assumed Bp = 2%, a common value for coastal waters [21, 26], and computed bbp = Bp*bp. To check the validity of this assumption, we ran the same analysis using Bp = 1% and 3% and computed normalized Lu from Eqs. (2-4) using the measured bp and Bp = 1% and 3%. The normalized Lu values for Bp = 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.
2.6. Our RT model parameterization
Our 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 , and a particle volume scattering function, βp(θ). Varying the choice and / or normalization of βp(θ) generated three sets of model output for comparison with each NuRADS image.
To recreate the proper βp(θ) for input to the RT code under method 1 we multiplied by the depth-weighted, spectrally-interpolated bp. Note that, due to the acceptance angle of the ac-9 (0.93 degrees), our depth-weighted bp underestimates the value that would be measured by a perfect instrument for the Petzold phase function by 20-30% .
Method 2 used the Sullivan and Twardowski  average phase function. The MASCOT instrument used in  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 . 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 apg, and bp from 0.01 to 1.0 m−1 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 either 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. 1 ).
The observation that truncating the phase function does not affect the calculated upwelling radiance significantly is consistent with Mobley et al. , who concluded that the exact shape of the phase function in the forward direction was not critical as long as bbp was correct. Furthermore, Gordon  showed that β(θ) truncated at scattering angles < 15 degrees (such as the one used here) did not alter calculations of irradiances by more than a few percent.
Note that Sullivan and Twardowski  used a normalization in the backward direction only, i.e. their measured βp(θ) were normalized by bbp and not by bp as had been done in . Therefore, the for method 2 satisfies the normalization:
To recreate the proper βp(θ) for input to the RT code under method 2 we multiplied by the depth-weighted, spectrally-interpolated bbp.
Method 3 also used the Petzold  turbid-water β(θ), but rather than scaling it by bp, as in method 1, it was normalized in the backward direction, as in method 2. Thus, for method 3 satisfied Eq. (6) and was multiplied by the depth-weighted, spectrally-interpolated bbp to produce the β(θ) input to the RT calculations. The β(θ) used for method 3 has the same shape as β(θ) used for method 1, but, like method 2, forces bbp in the model to agree with the measurements.
2.7. Evaluation of the model-data difference computation
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 5 to 45 degrees and every 15 degrees in azimuth from 0 to 180 degrees. Thus, the difference, D, between each model output, LuM, and the corresponding average NuRADS image, LuD, was computed at 117 angles based on Eq. (7).
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  model modified to incorporate the three-dimensional 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 (ds) was estimated and used with the total absorption (at) to compute the transmission (T) through the shadow: T = exp(-atds). 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 m−1) 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(θv, ϕ) were plotted as grouped by chlorophyll, ϕ, and θs. 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 apg or bp and chlorophyll, as described below.
3.1. Case I and Case II divisions within the data
The optical properties of Case I waters are determined by the chlorophyll concentration . Our data set was divided into Case I and Case II subsets, therefore, by visual inspection of plots of inherent optical properties as a function of chlorophyll concentration (Fig. 2 ). The IOP model of  (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) [Chl] > 10 μg l−1, or (b) bp > 1 m−1, or (c) apg > 0.2 m−1.
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 apg and bp were plotted in Fig. 2.
3.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 MAG2002 tables because the remainder, all of which were from either the Chesapeake or Monterey Bay data sets, had total chlorophyll values greater than 10 μg / l (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 Lu(θv, ϕ) / Lu(0, 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 least-squares 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, D(θv, ϕ), 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 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(θv, ϕ) were highly skewed or contained extreme outliers; box plots were useful for visualizing such features.
The population (N) of each distribution is given along the top of each plot. Due to the data acquired with chlorophyll values greater than 10 μg l−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 D(θv, ϕ) = 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 D(θv, ϕ). 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.
For chlorophyll concentrations < 1.0 μg l−1, MAG2002 had both the smallest overall range of D(θv, ϕ), 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(θv, ϕ) for Lee2011 was as close to zero as MAG2002 for chlorophyll concentrations < 0.18 μg l−1 (Fig. 4).
At chlorophyll concentrations > 1.0 μg l−1 in the Case I waters we sampled, all of the models had similar distributions of D(θv, ϕ). It was not obvious that one model performed consistently better than the others for chlorophyll concentrations > 1.0 μg l−1 (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 μg l−1, all of the models had close to zero median and similar extreme values of D(θv, ϕ). In the chlorophyll bins with smaller sample sizes, 1.00 - 1.78, 17.8 - 31.6, 56.2 - 178, and > 178 μg l−1, 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 μg l−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 μg l−1) and the bin above (10.0 - 17.8 μg l−1) the one in question.
In Case I conditions, all of the models had positive biases in and largest extreme values of D(θv, ϕ) 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 bbp 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 bp values had larger biases and extreme values than MAG2002 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 MAG2002 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 Lee2011 and our RT model using the Sullivan-Twardowski phase function produced median values of D(θv, ϕ) 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 D(θv, ϕ) 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 θs > = 20° and for our RT model for θs > = 25°. All models most closely agreed with the data at solar zenith angles between 25° to 35° (Fig. 8). At θs > = 40° a positive bias in D(θv, ϕ) and highly skewed distributions were observed for all models (Fig. 8). The median value of D(θv, ϕ) computed with MAG2002 exhibited the least bias for θs > = 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 D(θv, ϕ) that were within the average NuRADS coefficient of variation for solar zenith angles as large as 65°.
The Case II results, plotted against solar zenith angle (Fig. 9), show much greater model-data agreement as well as consistency among the models than the data from Case I conditions (Fig. 8). The only examples of median D(θv, ϕ) falling outside the range of the average NuRADS coefficient of variation were for situations with small sample sizes: all models with 35° < = θs < 40° (N = 29), MAG2002 for 45° < = θs < 50° (N = 9), and MAG2002 for 50° < = θs < 55° (N = 24). No systematic patterns of bias varying with solar zenith angle were observed (Fig. 9). Note, however, that there were no observations available for Case II water with θs < 35°. Therefore if a shadow effect such as the one observed in Case I water for θs < 25° occurs in Case II water, the data would not have captured it.
Across all of the observed Case I conditions, MAG2002 agreed with much of the available BRDF data more closely than any of the other models. In particular, three portions of the Case I data space were best fit by MAG2002: (a) chlorophyll concentrations in the range 0.18 to 1.0 μg l−1 (Fig. 4), (b) azimuth angles near the principal plane, i.e. outside of the range 90° - 135° (Fig. 6), and (c) solar zenith angles greater than 35° (Fig. 8). At both low and high chlorophyll concentrations, and for solar zenith angles from ~20-40 degrees, however, BRDF correction factors from Lee2011 and our method 2, using the Sullivan and Twardowski  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 . 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 μg l−1 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 f 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 backscattering (i.e. the value of bbp) and the detailed shape of the backward portion of the particulate phase function affected the BRDF correction in Case I water. At small azimuth angles, methods 2 and 3, in which the value of bbp matches measurements, both agree more closely with the data than method 1, in which the value of bp matches measurements (Fig. 6). This result indicates that the magnitude of bbp is more important for BRDF corrections at small azimuth angles than the detailed shape of the backward portion of the particulate 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 bbp correct was necessary but not sufficient to match the data at large azimuth angles, and that the shape of the particulate 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 Sea-viewing 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.
Two questions addressed were: How well does the recent Lee2011 BRDF correction agree with a wide suite of in situ measurements? Is the choice of phase function used in Lee2011 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  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 different particulate phase functions did not always produce similar output suggests that in Case I water the shape of the particulate 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 particulate phase function was important at 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 Lee2011 in Case II conditions.
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  was superior to MAG2002  in Case II water at small azimuth angles and for chlorophyll concentrations greater than the 10 μg l−1 limit of MAG2002. MAG2002 , however, remained superior to Lee2011  for Case I water.
This work was supported by NASA under grants NNX08AH93A (Voss) and NNX06AH14G (Trees). The acquisition of BIOSOPE data were funded through Centre National de Recherche Scientifique - Institut National des Sciences de l'Univers grants. The Monterey Bay data were collected as part of the COAST experiment, supported by NOAA. The LSCV’08 cruise was funded by the European Space Agency (21916/08/I-OL, Trees) and the NATO Undersea Research Centre. BP’09 was funded by the NATO Undersea Research Centre. Clark supported by NIST contract SB134110SE0990. Scott Freeman is thanked for his technical assistance in the field.
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). [CrossRef] [PubMed]
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). [CrossRef] [PubMed]
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). [CrossRef]
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). [CrossRef] [PubMed]
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). [CrossRef] [PubMed]
10. Z. P. Lee, K. Du, K. J. Voss, G. Zibordi, B. Lubac, R. Arnone, and A. Weidemann, “An inherent-optical-property-centered approach to correct the angular effects in water-leaving radiance,” Appl. Opt. 50(19), 3155–3167 (2011). [CrossRef] [PubMed]
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). [CrossRef]
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). [CrossRef] [PubMed]
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).
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). [CrossRef]
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/Rev4-Vol 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). [CrossRef]
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).
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). [CrossRef]
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). [CrossRef] [PubMed]
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).
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). [CrossRef] [PubMed]
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). [CrossRef] [PubMed]
31. H. R. Gordon and K. Y. Ding, “Self-shading of in-water optical instruments,” Limnol. Oceanogr. 37(3), 491–500 (1992). [CrossRef]