Measuring and Modeling the Polarized Upwelling Radiance Distribution in Clear and Coastal Waters

: The upwelling spectral radiance distribution is polarized, and this polarization varies with the optical properties of the water body. Knowledge of the polarized, upwelling, bidirectional radiance distribution function (BRDF) is important for generating consistent, long-term data records for ocean color because the satellite sensors from which the data are derived are sensitive to polarization. In addition, various studies have indicated that measurement of the polarization of the radiance leaving the ocean can used to determine particle characteristics (Tonizzo et al., 2007; Ibrahim et al., 2016; Chami et al., 2001). Models for the unpolarized BRDF (Morel et al., 2002; Lee et al., 2011) have been validated (Voss et al., 2007; Gleason et al., 2012), but variations in the polarization of the upwelling radiance due to the sun angle, viewing geometry, dissolved material, and suspended particles have not been systematically documented. In this work, we simulated the upwelling radiance distribution using a Monte Carlo-based radiative transfer code and measured it using a set of ﬁsh-eye cameras with linear polarizing ﬁlters. The results of model-data comparisons from three ﬁeld experiments in clear and turbid coastal conditions showed that the degree of linear polarization ( DOLP ) of the upwelling light ﬁeld could be determined by the model with an absolute error of ± 0.05 (or 5% when the DOLP was expressed in %). This agreement was achieved even with a ﬁxed scattering Mueller matrix, but did require in situ measurements of the other inherent optical properties, e.g., scattering coefﬁcient, absorption coefﬁcient, etc. This underscores the difﬁculty that is likely to be encountered using the particle scattering Mueller matrix (as indicated through the remote measurement of the polarized radiance) to provide a signature relating to the properties of marine particles beyond the attenuation/absorption coefﬁcient.


Introduction
The in-water and water-leaving radiance in the ocean is partially polarized and this has implications for both biological activity [1] and for viewing the ocean from a satellite [2]. The biological implications of the polarization are a current topic of research [3] but the implications for ocean color remote sensing are just being exploited [2,4]. Various studies have indicated that measurement of the polarization of the radiance leaving the ocean can be used to determine particle characteristics [5][6][7]. In addition to using the polarization of the water leaving radiance for studies of the water properties, polarization may also be important for processing data from ocean color sensors that have unintended polarization sensitivities [8,9].
Measurements of the scalar (without regard to polarization) spectral upwelling radiance distribution have occurred more frequently since the development of the radiance distribution camera, RADS [10], and then the upwelling radiance distribution camera, NuRADS [11]. These instruments use electro-optic camera systems combined with filter changers and fisheye cameras to image the complete upwelling radiance distribution, for a specific wavelength, in one image. The CCD (charge coupled device) resolution and optics allow measurement of the radiance distribution with a 1 • angular resolution. These systems have been used in studies of the in-water light field [12] and studies of the angular radiance distribution variations in ocean color algorithms [13,14]. Recently, two new cameras, polarized radiance camera (PolRADS) [15] and downwelling polarized camera (DPOL) [16], have been used to make measurements of the polarized spectral upwelling radiance distribution, as described below.
Measurements of the upwelling, in situ, polarized radiance have typically been done with variations of a Gershun-tube radiometer, obtaining the angular distribution by changing the viewing direction of the radiometer [17,18]. These instruments have advantages such as the simplicity of calibration and the ability to do hyperspectral measurements. Hyperspectral measurements of the polarized light field have been suggested as a method to measure in situ solar stimulated fluorescence [19]. Unfortunately, to obtain the full radiance distribution requires many measurements during which the illumination conditions can change, and spatially varying angular features resulting from the sea surface may not be captured. By making measurements using fisheye lenses, while limited to multi-spectral measurements rather than hyperspectral measurements, higher angular spatial resolution can be obtained nearly simultaneously (typical exposure times are less than 1 s).
Models for the unpolarized bi-directional reflectance distribution function (BRDF) [20,21] have been validated in Case I [13] and Case II [14] waters, but variations in the polarization of upwelling radiance due to sun angle, viewing geometry, dissolved material, and suspended particles have not been systematically documented.

Methods
The light field polarization is easily described by use of the four element Stokes vector, I, as in, for example, Bohren and Huffman [22]. The four components of I are prescribed relative to some reference plane, which here we choose as the plane defined by the viewing direction and the nadir direction. Consider an electromagnetic wave of angular frequency ω. If the unit vectorˆ is parallel to, and the unit vectorr is perpendicular to, the reference plane (such thatˆ ×r is in the direction of propagation), the electric field can be written as → E = E ˆ + E rr , with: where E 0 and E r0 are real. The four components of the Stokes vector are then defined as: where δ ≡ δ r − δ . A derived parameter that we will use in our comparison is the degree of linear polarization, DOLP, defined as: Another useful parameter is the angle of the plane of polarization, χ, which is defined as: We measured the Stokes vector over the upwelling hemisphere using the PolRADS and DPOL cameras and simulated it using a Monte Carlo-based radiative transfer (RT) code. Data were available from three field campaigns: Hawaii, USA, in December 2005; the Ligurian Sea in March 2009; and the New York Bight in May 2009.

Measurements of the Polarized Upwelling Radiance Distribution
In Hawaii, measurements of the polarized upwelling radiance distribution were accomplished using the PolRADS camera system [15]. PolRADS is based on the NuRADS camera system [12], which is a compact (30 cm diameter, 30 cm length), multispectral camera that images the upwelling radiance distribution in six narrow (≈10 nm full-width at half-maximum, FWHM) spectral bands centered at 412, 436, 486, 526, 548, and 615 nm. In the PolRADS instrument, three synchronized NuRADS cameras are used, each with a linear polarizer, to simultaneously acquire images. Combining the images from three NuRADS cameras, when in PolRADS configuration, allows for retrieval of three elements of the Stokes vector (I, Q, and U, but not V), as well as DOLP.
In the Ligurian Sea and New York Bight, measurements of the polarized upwelling radiance distribution were accomplished with the DPOL camera system [16]. DPOL is similar to PolRADS in the sense that simultaneous measurements are made using multiple fisheye lenses with different polarizing filters in the optical path. DPOL differs from PolRADS by using optical-fiber bundles to project images from all lenses through a single spectral filter and onto a single CCD array. Eliminating the redundant filters and cameras makes DPOL much smaller than PolRADS, thereby reducing the instrument shadow. Furthermore, DPOL has four lenses, rather than the three on PolRADS, enabling retrieval of all four Stokes vector elements. The spectral filters used on DPOL were also ≈10 nm FWHM and were centered at 442, 488, 520, 550, 589, and 650 nm. For both PolRADS and DPOL, the channels above 600 nm will not be used because of instrument self-shading. In addition, because of its absence on DPOL, the 412 nm channel on PolRads will also not be discussed.
Data acquisition for both the PolRADS and DPOL systems used filter changers to rotate interference filters and thus sequentially acquire images in each wavelength band. Typical exposure times were less than one second. However, acquiring a set of images from all wavelengths took about two minutes due to the time required to read and store the data from the CCD. With both systems, typical deployments lasted from one to several hours, enabling multiple acquisitions over a range of solar zenith angles.
Reduction of the raw images consisted of applying calibration factors [15,16] and averaging images in both space and time to reduce environmental noise. After calibration, but 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. Calibrated images of the Stokes vectors were then averaged in 10-min 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 of 3 × 3-pixel windows was performed to produce final average images at a 1 • × 1 • resolution. Each pixel in the final, reduced image, therefore, could have been an average of up to 90 raw pixels (5 images × 2 image halves × 9-pixel window). The mean and coefficient of variation (CV = standard deviation divided by the mean) of the Stokes vector components were computed for each pixel in the reduced image using the up-to-90 raw pixels in the original images. The degree of linear polarization DOLP was then computed for each pixel from the mean I, Q, and U Stokes vector components (Equation (3)). Because V was not available from the PolRADS data, DOLP was calculated rather than the degree of polarization (DOP); however, in the upwelling lightfield, they are equivalent because the magnitude of V is negligible [23].
For in situ radiometric measurements, instrument self-shading must be considered [24]. In our case, instrument self-shading may also affect the DOLP measurements due to camera geometry. The physical displacement of the multiple lenses used means that the images through different polarizers will look into the instrument shadow in slightly different directions. When one lens views a portion of the radiance distribution with less shadow, it will appear brighter in that image. In the algorithm from which the polarization information is derived, for the scene to be unpolarized requires that region of the three camera images to be closely matched in intensity. If the region in one camera is less shaded, it will cause the algorithm to assume that the light field is polarized along the direction of the linear polarizer in that camera. Thus shadowing, which is not symmetric in all the images, will appear as an increase in the DOLP, and therefore, a negative model-data DOLP difference. This effect should be greater in turbid water and for the PolRADS instrument (because of its larger size) than in clear water or with the DPOL instrument. In our analysis, the area of direct shadow (the anti-solar point) was excluded from the data set. More subtle shadow areas may exist and could cause the measured DOLP to be larger than the modelled DOLP (as will be seen, this was not generally the case in our data set).

Modeling the Polarized Upwelling Radiance Distribution
A Monte Carlo model was used to solve the vector radiative transfer equation (Equation (5)) for the propagation of the Stokes vector, I, at wavelength λ, The four components of the column vector I = [I, Q, U, V] T were defined in Equation (2). The single-scattering albedo, ω 0 , the optical depth of the medium, τ, the scattering phase matrix, P, and the rotation matrix, R, are defined below. The right-handed coordinate system in Equation (5) has its origin at the top of the medium, z-axis pointed downward, and x-axis directed away from the sun. The nadir angle, θ, relative to the z-axis, and the azimuth angle, φ, in the x-y plane, define the direction of photon propagation. Photons originate in the solar beam, i.e., with θ = θ 0 , the solar zenith angle, and φ 0 = 0.
The model simulated a horizontally homogenous, two-layer system with a non-absorbing, Rayleigh-scattering atmosphere over an ocean comprised of Rayleigh-scattering water molecules plus scattering and absorbing hydrosols ("particles"). Fresnel reflectance at the air-water interface used the boundary condition defined in Gordon et al. [8]. The single-scattering albedo, ω 0 , and the optical depth of the medium, τ, were defined as per usual (e.g., Mobley [25]): and where the parameters a t (λ) + b t (λ) = c t (λ) are the total absorption, scattering, and attenuation coefficients (in m −1 ), respectively. For the atmospheric layer, a t (λ) = 0, therefore ω 0 = 1, and τ was taken from Teillet [26] assuming standard atmospheric conditions. For the oceanic layer, absorption (a w ) and scattering (b w ) coefficients of seawater were interpolated to PolRADS and DPOL spectral bands from Table 1.1 in Pegau et al. [27]. In the Ligurian Sea and New York Bight experiments, the absorption coefficient of dissolved and particulate constituents (a pg ) and the particle scattering (b p ) coefficients were measured in situ, depth-weighted, and interpolated as necessary to the PolRADS and DPOL spectral bands as described by Gleason et al. [14]. No in situ measurements were available from Hawaii, therefore for that dataset, a pg and b p were derived from Morel and Gentili [28] (Equations (7) and (8) in reference [28]) using an estimated total chlorophyll concentration of 0.1 g/m 3 , the long-term average value at the measurement site. Summing the water and particle contributions gave the total coefficients required for Equations (6) and (7): a t = a ω (λ) + a pg (λ) and b t = b ω (λ) + b p (λ). For each image data set, the model was run at the corresponding solar zenith angle. Scattering events in Equation (5) are represented by a 4 × 4 matrix called the Mueller matrix, M. For example, a photon travelling in direction ξ with Stokes vector I' would be scattered in a new direction ξ with Stokes vector I by the linear transformation I(ξ) = MI'(ξ'). Note that M is defined relative to the scattering plane (as is traditional), the plane defined by vectors ξ' and ξ, but we have defined I(θ, φ) relative to the viewing direction and the nadir direction. Pre-and post-multiplication by the rotation matrix R(α) is required to account for changes in reference frames: The rotation angle α = cos −1 (l l ·l r ) is measured clockwise from the vectorl l in the initial reference frame to the vectorl r in the rotated reference frame [8].
The scattering phase matrix, P, is a Mueller matrix normalized to the integral of the M 11 element over all solid angles: M ij represents the ith row and jth column of M, and Θ is the scattering angle. It is advantageous, when investigating the polarization effects of the Muller matrix, to normalize in a different way, specifically to the M 11 element at each angle Θ (Equation (11)).
Then, P = βS (12) where β = M 11 /b is the scattering phase function typically used in scalar radiance transfer models. For our modeling, P in the atmosphere and P w , the seawater component of the ocean, were both set to P r , given by Rayleigh scattering: (We examined a few cases using a code that included the (small) depolarization of Rayleigh scattering in the atmosphere and water by molecular anisotropy and concluded that its omission would not significantly affect the results.) Particle scattering in the ocean used β [25,29] and two alternative parameterizations of S. One set of model runs used S from Voss and Fry [30], which was experimentally determined from samples of seawater (referred to here after as V-F). Note that Voss and Fry found all off-diagonal elements of S equal to zero, within experimental error, except for S 12 and S 21 . Furthermore, Voss and Fry found S 12 ≈ S 21 . Subsequent laboratory measurements of the Mueller matrix elements for phytoplankton [31][32][33] have been similar to V-F. Specifically, the S 12 (Θ) element of phytoplankton samples tends to have a minimum value of between −0.6 to −0.8 at an angle between 90 • and 100 • . In contrast, experiments with solutions of suspended marine sediment have revealed smaller absolute values for S 12 (Θ) element near Θ = 90 • [32]. Because of this, we decided to carry out a second set of runs with a particle phase matrix having a smaller minimum in S 12 near 90 • (Figure 1). This was accomplished by replacing the V-F S 12 (Θ) with that for Rayleigh scattering with a depolarization factor ρ. ρ is the ratio of light scattered with polarization parallel to the plane of incidence to that scattered perpendicular to the plane of incidence, when the incident radiation is polarized perpendicular to the plane of incidence. The element S 12 (Θ) as a function of ρ is given by Equation (14) [34]. The other elements of S remained unchanged for this Mueller matrix (referred to as Mod-V-F). Setting ρ = 0.3 resulted in a minimum value of S 12 (90 • ) = −0.37, which is within the range of −0.38 to −0.25 observed by Volten et al. [32].
Kuik et al. [35] presented four inequalities, originally derived by Fry and Kattawar [36] and Hovenier et al. [37], that must be satisfied by the scattering matrix elements of randomly oriented particles having a plane of symmetry. These inequalities are all satisfied by the modified Mueller matrix. .
Kuik et al. [35] presented four inequalities, originally derived by Fry and Kattawar [36] and Hovenier et al. [37], that must be satisfied by the scattering matrix elements of randomly oriented particles having a plane of symmetry. These inequalities are all satisfied by the modified Mueller matrix.

Model-Data Comparison
Visual comparisons provided a qualitative check that the Stokes vector components from our RT model were in agreement with the camera data. Quantitative comparisons between the model and data were performed by computing the model-data difference every 10° in nadir from 0° to 80° and every 30° in azimuth from 0° to 180°. Thus, the difference in DOLP, DOLPdiff, between each model output, DOLPmodel, and the corresponding average PolRADS or DPOL image, DOLPdata, was computed at 63 angles using Equation (15).
where ( , ) and ( , ) are the DOLP predicted by the RT model and measured by the corresponding average PolRADS or DPOL image, respectively, at the same 63 nadir and azimuth angles.

Results
In total, approximately 528 individual images were selected and averaged in 10 min intervals to produce 219 data sets, which were compared with RT model runs (Table 1). Table 1. Number, N, of reduced images used for model-data comparisons in each of the three datasets in each of four spectral bands. Note, the PolRADS and DPOL cameras used different filters with

Model-Data Comparison
Visual comparisons provided a qualitative check that the Stokes vector components from our RT model were in agreement with the camera data. Quantitative comparisons between the model and data were performed by computing the model-data difference every 10 • in nadir from 0 • to 80 • and every 30 • in azimuth from 0 • to 180 • . Thus, the difference in DOLP, DOLP diff , between each model output, DOLP model , and the corresponding average PolRADS or DPOL image, DOLP data , was computed at 63 angles using Equation (15).
where DOLP model (θ, φ) and DOLP data (θ, φ) are the DOLP predicted by the RT model and measured by the corresponding average PolRADS or DPOL image, respectively, at the same 63 nadir and azimuth angles.

Results
In total, approximately 528 individual images were selected and averaged in 10 min intervals to produce 219 data sets, which were compared with RT model runs (Table 1). Table 1. Number, N, of reduced images used for model-data comparisons in each of the three datasets in each of four spectral bands. Note, the PolRADS and DPOL cameras used different filters with slightly different band centers; the center of each 10 nm-wide filter is listed in parenthesis.  Figure 2 shows an example image for clear water. In the right column, the data parameters are shown (I, Q, U, DOLP, and χ), while in the center column, the results from the RT model are shown, where the model inputs are based on the measured parameters for the data (b p , a pg , solar zenith angle). As mentioned earlier, during processing of this data, the symmetry about the principal plane was used to average the left and right side of the data (with proper handling of the sign of U/I). Thus, to make these images, this symmetry was used to generate the left side of the data images. The left column is a single scattering calculation, where the Mueller matrix was assumed to be V-F, and the solar zenith angle was the value appropriate for the data. As can be seen, the same broad patterns were visible in the single scattering model, the RT model and data. The I component of the Stokes vector was largest near the horizon in the direction towards the sun. The Q/I and U/I patterns were very similar in all three cases. For the DOLP, the single scattering model was significantly different than the data and the RT model. With single scattering, the maximum DOLP was both larger than in the RT model and the data, and occurred for all nadir angles at the scattering angle matching the minimum in the S 12 and S 21 elements (90-100 • ). Interestingly, while it was a subtle effect in the clear water case, the maximum DOLP in the model and data occurred at the same scattering angle, but at an azimuth of 90 • relative to the principal plane (the plane containing the anti-solar position and nadir). This can be contrasted with the case of the downwelling sky radiance distribution, where the maximum DOLP occurred in the principal plane, with a decrease in DOLP towards the horizon (for example, Liu and Voss [38]). The same cause was in effect in both cases: for downwelling radiance, the horizon had more multiple scattering than the principal plane at the 90 • scattering angle, while for the upwelling radiance, in the water, there was more multiple scattering for the principal plane at a 90 • scattering angle than for nadir angles closer to the horizon. The other feature in the data for this clear water case was the separation of the neutral points, areas of zero DOLP, on either side of the principal plane, as has been discussed previously [39,40]. Finally, the angle of the polarization plane, χ, was very similar in all three cases, perhaps less well-defined for the data in the area of the low polarization, around the anti-solar position.

Overall View of the Upwelling Polarization Signal
An example with more turbid water is shown in Figure 3. Once again, the patterns were generally the same, however in this case the data are noisier. Note that the data shown in this case has not been averaged using the left/right symmetry across the principle plane, thus appears less symmetric than the example in Figure 2. As has been pointed out before [41], the most stable parameter was χ, with the exception of the area with very low DOLP, in which case the plane of polarization was not well defined. In general, the images appeared much noisier in the data from the two more turbid sites (New York Bight and Ligurian Sea) than in the clear water off Hawaii. The pattern of Q/I and U/I tended to shift more between each image that goes into these data averages. In the scale of Q/I and U/I ranging from −1 to 1, the standard deviation obtained during averaging the data was approximately 0.05 for the New York Bight and Ligurian Sea data, while it was 0.012 for the clearer Hawaii data. One can also see that in this more turbid case, the single scattering was a poor approximation of all parameters, with the exception of χ.

Maximum DOLP and Nadir DOLP
How large can the maximum DOLP be in the upwelling radiance? Figure 4 shows this as a function of ct/at and ct. We used the variation of the optical properties (ct, bt, bbt, and at) with the wavelength to fill in the data, thus this figure includes all wavelengths with their appropriate optical properties. The quantity ct/at can be shown to equal the mean number of scattering events in the medium taken as a whole [28]. As the water gets more turbid, there is more multiple scattering, which will decrease the DOLP [7,42]. The maximum DOLP is often outside the Snell's cone, thus not retrievable from above the surface, so we also show the maximum value of the DOLP inside the Snell's cone in Figure 4. As can be seen, this was less, sometimes significantly less, than the maximum DOLP in the total upwelling field. Another point was that the maximum DOLP in the total upwelling field occurred at an azimuth nearly perpendicular to the principal plane, as this position was always at a 90° scattering angle to the sun, and had the largest component of single scattered light due to a combination of light attenuation with depth and the relative smoothness of the scattering function in

Maximum DOLP and Nadir DOLP
How large can the maximum DOLP be in the upwelling radiance? Figure 4 shows this as a function of c t /a t and c t . We used the variation of the optical properties (c t , b t , b bt , and a t ) with the wavelength to fill in the data, thus this figure includes all wavelengths with their appropriate optical properties. The quantity c t /a t can be shown to equal the mean number of scattering events in the medium taken as a whole [28]. As the water gets more turbid, there is more multiple scattering, which will decrease the DOLP [7,42]. The maximum DOLP is often outside the Snell's cone, thus not retrievable from above the surface, so we also show the maximum value of the DOLP inside the Snell's cone in Figure 4. As can be seen, this was less, sometimes significantly less, than the maximum DOLP in the total upwelling field. Another point was that the maximum DOLP in the total upwelling field occurred at an azimuth nearly perpendicular to the principal plane, as this position was always at a 90 • scattering angle to the sun, and had the largest component of single scattered light due to a combination of light attenuation with depth and the relative smoothness of the scattering function in the backward direction. For the maximum DOLP inside the Snell's circle, as can be seen in Figures 2  and 3, this azimuthal angle will move towards the principal plane, and hence got closer to the position where, above the surface, it will be in the region of the solar glitter pattern.  Another interesting parameter is the DOLP at the nadir angle, as this was the direction that many in situ radiometers make their measurement ( Figure 5). Here the RT model results are presented rather than the data, as in our definition of the reference frame, the nadir direction is a point of singularity, and the data become excessively noisy at that point.  Another interesting parameter is the DOLP at the nadir angle, as this was the direction that many in situ radiometers make their measurement ( Figure 5). Here the RT model results are presented rather than the data, as in our definition of the reference frame, the nadir direction is a point of singularity, and the data become excessively noisy at that point.  Another interesting parameter is the DOLP at the nadir angle, as this was the direction that many in situ radiometers make their measurement ( Figure 5). Here the RT model results are presented rather than the data, as in our definition of the reference frame, the nadir direction is a point of singularity, and the data become excessively noisy at that point.  Note that the nadir light field could be significantly polarized and the DOLP could be greater than 20%. Also, as can be seen, the largest values (22%) occurred at large solar zenith angles (when the scattering angle decreased to values closer to 90 • ), but there was also a trend of decreasing DOLP as c t /a t increased at a constant solar zenith angle due to increasing multiple scattering. The plane of polarization in all cases was perpendicular to the principal plane.

Comparison of DOLP Differences
For the first quantitative comparison of the RT model and the data, we looked at scatter plots comparing the measured and modeled DOLP (Figure 6). For each cruise, a best linear fit line was calculated, with the y-intercept = 0. The fit was calculated for all of the data, as well as the data inside the Snell's circle; however, the calculated slopes were not significantly different, so we will discuss only the fit for the total image. The slopes were 0.94 (±0.002, r 2 = 0.81), 0.88 (±0.006, r 2 = 0.84), and 0.82 (±0.008, r 2 = 0.76) for Hawaii, New York Bight, and Ligurian Sea, respectively.
Because the DOLP was significantly overestimated in the Ligurian Sea data, the alternative normalized Mueller matrix Mod-V-F (ρ = 0.30) was tried, with all other parameters the same, resulting in a slope of 1.20 (± 0.009, r 2 = 0.656); we concluded that the depolarization of 0.3 was too strong. The true average Mueller matrix for this location could be somewhere between the two we used.
Note that the nadir light field could be significantly polarized and the DOLP could be greater than 20%. Also, as can be seen, the largest values (22%) occurred at large solar zenith angles (when the scattering angle decreased to values closer to 90°), but there was also a trend of decreasing DOLP as ct/at increased at a constant solar zenith angle due to increasing multiple scattering. The plane of polarization in all cases was perpendicular to the principal plane.

Comparison of DOLP Differences
For the first quantitative comparison of the RT model and the data, we looked at scatter plots comparing the measured and modeled DOLP (Figure 6). For each cruise, a best linear fit line was calculated, with the y-intercept = 0. The fit was calculated for all of the data, as well as the data inside the Snell's circle; however, the calculated slopes were not significantly different, so we will discuss only the fit for the total image. The slopes were 0.94 (±0.002, r 2 = 0.81), 0.88 (±0.006, r 2 = 0.84), and 0.82 (±0.008, r 2 = 0.76) for Hawaii, New York Bight, and Ligurian Sea, respectively.
Because the DOLP was significantly overestimated in the Ligurian Sea data, the alternative normalized Mueller matrix Mod-V-F (ρ = 0.30) was tried, with all other parameters the same, resulting in a slope of 1.20 (± 0.009, r 2 = 0.656); we concluded that the depolarization of 0.3 was too strong. The true average Mueller matrix for this location could be somewhere between the two we used. Figure 6. Scatter plots of DOLP for model versus data, separated by cruise. Also shown is the best fit line (with y-intercept = 0) (black line) and 1:1 line (red line). The crosses correspond to data outside the Snell's circle (nadir angle larger than 48°), and the dots to points inside the Snell circle. V-F was used in the model for (a-c), Mod-V-F was used for (d). The correlation between DOLP diff and several environmental parameters was investigated. In Figure 7, we show the DOLP diff variation with solar zenith angle. In this figure, the mean DOLP diff was calculated for each data image, along with the standard deviation (shown as error bars). There was no trend in the data; however, there does appear to have been a slight increase in the magnitude of the deviation at solar zenith angles greater than 70 • . In general though, in almost all cases, the mean DOLP diff was within one standard deviation of zero.
The correlation between DOLPdiff and several environmental parameters was investigated. In Figure 7, we show the DOLPdiff variation with solar zenith angle. In this figure, the mean DOLPdiff was calculated for each data image, along with the standard deviation (shown as error bars). There was no trend in the data; however, there does appear to have been a slight increase in the magnitude of the deviation at solar zenith angles greater than 70°. In general though, in almost all cases, the mean DOLPdiff was within one standard deviation of zero. In Figure 8, we look at the DOLPdiff dependence on ct and ct/at. While there was no clear trend in DOLPdiff with ct, a slight trend occurred as a function of ct/at. Low values of ct/at (more absorption, less multiple scattering) had a slightly larger value of DOLPdiff, while DOLPdiff tended towards 0 with larger values of ct/at (more multiple scattering). However, the mean DOLPdiff was less than 5% for almost all cases, and this was within the measurement uncertainty of these instruments in the field.
Finally, in Figure 9, we look at the DOLPdiff as a function of bbt/bt (bbt/bt is the fraction of total scattering which was in the backwards direction, i.e., scattering angles from 90° to 180°). As bbt/bt increased, DOLPdiff appeared slightly higher, indicating less agreement at higher values of bbt/bt. Low values of bbt/bt were indicative of a low refractive index, phytoplankton-dominated environment, while higher values of bbt/bt indicated an environment with other higher refractive index particles, perhaps sediment. V-F was based on measurements in a variety of environments (clear ocean water to more turbid coastal water); however, measurements of the Mueller matrix of phytoplankton cultures also agreed with V-F. It has often been stated that the polarization resulting from higher index particles is less than V-F (i.e., larger depolarization at a 90° scattering angle) [32]. It is probably the case in our data set, that as the bbt/bt ratio increases, the S12 and S21 elements of the Mueller matrix should also decrease in an absolute value sense, but not as much of a decrease as used in Mod-V-F, where the decrease in this matrix would have caused a significant negative DOLPdiff. In Figure 8, we look at the DOLP diff dependence on c t and c t /a t . While there was no clear trend in DOLP diff with c t , a slight trend occurred as a function of c t /a t . Low values of c t /a t (more absorption, less multiple scattering) had a slightly larger value of DOLP diff , while DOLP diff tended towards 0 with larger values of c t /a t (more multiple scattering). However, the mean DOLP diff was less than 5% for almost all cases, and this was within the measurement uncertainty of these instruments in the field.
Finally, in Figure 9, we look at the DOLP diff as a function of b bt /b t (b bt /b t is the fraction of total scattering which was in the backwards direction, i.e., scattering angles from 90 • to 180 • ). As b bt /b t increased, DOLP diff appeared slightly higher, indicating less agreement at higher values of b bt /b t . Low values of b bt /b t were indicative of a low refractive index, phytoplankton-dominated environment, while higher values of b bt /b t indicated an environment with other higher refractive index particles, perhaps sediment. V-F was based on measurements in a variety of environments (clear ocean water to more turbid coastal water); however, measurements of the Mueller matrix of phytoplankton cultures also agreed with V-F. It has often been stated that the polarization resulting from higher index particles is less than V-F (i.e., larger depolarization at a 90 • scattering angle) [32]. It is probably the case in our data set, that as the b bt /b t ratio increases, the S 12 and S 21 elements of the Mueller matrix should also decrease in an absolute value sense, but not as much of a decrease as used in Mod-V-F, where the decrease in this matrix would have caused a significant negative DOLP diff .

Comparison of Q/I and U/I Differences
Rather than calculating DOLP, we could also show the difference between the model result for Q/I and U/I versus the data, shown in Figure 10 for the Ligurian Sea case. This case is shown because it was neither the best result, which was for Hawaii, nor the worse, which was New York Bight. The advantage of DOLP, and the reason most of the results are presented for this parameter, is that it is independent of the frame of reference, thus small errors in locating the solar plane in each image cause smaller differences in the comparison. Q/I and U/I comparisons are with respect to the frame of reference and depend on correctly locating the solar plane in the image, hence the comparison is not quite as good between the model and data result as with the DOLP. The line fits had a slope of 0.66 (±0.010, r 2 = 0.70) for Q/I and 0.71 (±0.013, r 2 = 0.67) for U/I. The fact that the scatter is qualitatively the same in these graphs as in the DOLP figure can be simply understood since Q/I and U/I can be derived from a combination of the DOLP and the plane of polarization, . As mentioned earlier, is a very stable parameter in the upwelling light field. The spread in the data is slightly larger for U/I versus Q/I because U/I has stronger spatial gradients, and small rotations affect this parameter more than Q/I.

Comparison of Q/I and U/I Differences
Rather than calculating DOLP, we could also show the difference between the model result for Q/I and U/I versus the data, shown in Figure 10 for the Ligurian Sea case. This case is shown because it was neither the best result, which was for Hawaii, nor the worse, which was New York Bight. The advantage of DOLP, and the reason most of the results are presented for this parameter, is that it is independent of the frame of reference, thus small errors in locating the solar plane in each image cause smaller differences in the comparison. Q/I and U/I comparisons are with respect to the frame of reference and depend on correctly locating the solar plane in the image, hence the comparison is not quite as good between the model and data result as with the DOLP. The line fits had a slope of 0.66 (±0.010, r 2 = 0.70) for Q/I and 0.71 (±0.013, r 2 = 0.67) for U/I. The fact that the scatter is qualitatively the same in these graphs as in the DOLP figure can be simply understood since Q/I and U/I can be derived from a combination of the DOLP and the plane of polarization, . As mentioned earlier, is a very stable parameter in the upwelling light field. The spread in the data is slightly larger for U/I versus Q/I because U/I has stronger spatial gradients, and small rotations affect this parameter more than Q/I.

Comparison of Q/I and U/I Differences
Rather than calculating DOLP, we could also show the difference between the model result for Q/I and U/I versus the data, shown in Figure 10 for the Ligurian Sea case. This case is shown because it was neither the best result, which was for Hawaii, nor the worse, which was New York Bight. The advantage of DOLP, and the reason most of the results are presented for this parameter, is that it is independent of the frame of reference, thus small errors in locating the solar plane in each image cause smaller differences in the comparison. Q/I and U/I comparisons are with respect to the frame of reference and depend on correctly locating the solar plane in the image, hence the comparison is not quite as good between the model and data result as with the DOLP. The line fits had a slope of 0.66 (±0.010, r 2 = 0.70) for Q/I and 0.71 (±0.013, r 2 = 0.67) for U/I. The fact that the scatter is qualitatively the same in these graphs as in the DOLP figure can be simply understood since Q/I and U/I can be derived from a combination of the DOLP and the plane of polarization, χ. As mentioned earlier, χ is a very stable parameter in the upwelling light field. The spread in the data is slightly larger for U/I versus Q/I because U/I has stronger spatial gradients, and small rotations affect this parameter more than Q/I.

Conclusions
We have compared measurements of the polarization properties of the upwelling light field in the marine environment to radiative transfer models based on the measured inherent optical properties (IOPs) and an average normalized Mueller matrix for the scattering. For the Ligurian Sea and New York Bight data sets, the base RT model used the same volume scattering phase function (Petzold) and Mueller matrix (V-F) for the standard model, with the scattering and absorption coefficients as measured coincidently on location. For Hawaii, the IOPs were not measured, but were estimated based on an average Chl value of 0.1 mg/m 3 at the measurement location. Even with the fixed scattering matrix parameters, on average, the DOLPdiff was less than 5%, and in almost every case was within one standard deviation of zero. Thus, to get more information, or discrimination, from the DOLP measurements, the uncertainty of the measured DOLP must be smaller than 5%, and the other parameters (scattering phase function, Mueller matrix, and absorption coefficient) need to be known well enough to constrain the DOLP model to this accuracy.
Some of the variability in our data was due to instrument noise (on the order of 2%); however, other noise was due to the variability in the environment resulting from the incident light field interacting with the wavy air-sea interface. This was more of a problem for the downwelling light field but was still evident in the upwelling light field. For remote sensing, this radiance (Stokes vector) would have been propagated through the (wavy) air-sea interface, which would introduce more noise.
The 5% agreement of the model with the data shows that the model will provide sufficient accuracy for correcting residual polarization effects in ocean color instruments. However, more work must be done at improving the measurement accuracy of the instrumentation if the goal is to be able to use polarization effects to characterize (and differentiate with respect to) the physical properties of the suspended particles other than c/a.

Conclusions
We have compared measurements of the polarization properties of the upwelling light field in the marine environment to radiative transfer models based on the measured inherent optical properties (IOPs) and an average normalized Mueller matrix for the scattering. For the Ligurian Sea and New York Bight data sets, the base RT model used the same volume scattering phase function (Petzold) and Mueller matrix (V-F) for the standard model, with the scattering and absorption coefficients as measured coincidently on location. For Hawaii, the IOPs were not measured, but were estimated based on an average Chl value of 0.1 mg/m 3 at the measurement location. Even with the fixed scattering matrix parameters, on average, the DOLP diff was less than 5%, and in almost every case was within one standard deviation of zero. Thus, to get more information, or discrimination, from the DOLP measurements, the uncertainty of the measured DOLP must be smaller than 5%, and the other parameters (scattering phase function, Mueller matrix, and absorption coefficient) need to be known well enough to constrain the DOLP model to this accuracy.
Some of the variability in our data was due to instrument noise (on the order of 2%); however, other noise was due to the variability in the environment resulting from the incident light field interacting with the wavy air-sea interface. This was more of a problem for the downwelling light field but was still evident in the upwelling light field. For remote sensing, this radiance (Stokes vector) would have been propagated through the (wavy) air-sea interface, which would introduce more noise.
The 5% agreement of the model with the data shows that the model will provide sufficient accuracy for correcting residual polarization effects in ocean color instruments. However, more work must be done at improving the measurement accuracy of the instrumentation if the goal is to be able to use polarization effects to characterize (and differentiate with respect to) the physical properties of the suspended particles other than c/a.