Underwater DoLP was measured with the Multi-Angle Scattering Optical Tool (MASCOT) (Twardowski et al., 2012; Sullivan et al., 2013). It has an array of 17 silicon diode detectors covering 10°–170° scattering angles (source-sample-detector angle) in 10° increments. Incident light source is a 30 mW 658 nm laser diode passing through a wedge depolarizer to provide unpolarized incident radiation. The sampling rate is 20 Hz. Full angles of the Detector field of views (FOVs) are 0.8°, 2°, 3°, and 4° for the 10°, 20°, 30°, and 40° detectors, respectively. FOV is 5° for the rest of the detectors. MASCOT covers a large scattering angle range, enabling accurate measurement of VSF at mid- and back-scatter angles. The MASCOT was designed with minimal form factors and structural elements (moving parts, sample holder, etc.) to minimize stray light contamination. It is an in situ device designed for direct measurement in nominally unperturbed waters. A filter wheel mounted in front of the source window can generate unpolarized (empty space on wheel), horizontally or vertically polarized incident light, or a dark blank (opaque location on wheel). The wheel continually spins at a rate that allows each location on the wheel to be sampled for 1 s. Particulate DoLP for unpolarized light was computed with Eq. 8. Accurate measurements of polarization elements have been verified with Lorenz-Mie theory for microspherical beads. For details on the instrument calibration and correction procedures in VSF measurement of MASCOT see (Twardowski et al., 2012). The MASCOT has been deployed extensively since 2006 (Sullivan and Twardowski 2009; You et al., 2011b; Gleason et al., 2012; Twardowski et al., 2012; Sullivan et al., 2013; Moore et al., 2017). Table 1 lists the field sites and measurements relevant in this study. Locations cover inland lakes, coastal oceans and open oceans. Relevant particulate inherent optical properties (IOPs) were also measured. Particulate absorption coefficient a p [ m − 1 ] and scattering coefficient b p [ m − 1 ] were derived from measurements of non-water absorption a n w , non-water attenuation c n w , and absorption in the dissolved fraction a C D O M with a 25-cm flow cell WET Labs ACS or AC9 device following the protocol in Twardowski et al. (1999). Backscattering coefficient b b p [ m − 1 ] was derived from MASCOT unpolarized VSF measurements. The particulate backscattering coefficient is, b b p = 2 π ∫ 90 o 180 o β u n p o l , p ( ψ ) sin ( ψ ) d ψ , (10) where β u n p o l , p is the interpolated unpolarized particulate VSF measured with MASCOT. β u n p o l , p ( 180 ˚ ) =   β u n p o l , p ( 170 ˚ ) is assumed in Eq.10. Since the sin ( ψ )   weighting in the integral approaches 0 at 180°, precise accuracy near 180° is not critical. The particulate scattering coefficient b p measured with ACS or AC9 at 657 nm was used to compute the backscattering ratio b b p / b p . At each field site, measurements were collected at over a dozen stations spanning the region. For each station, all measurements were averaged to 1-m depth bins ranging from the surface down to the maximum measurement depth.

To prepare the DoLP angular curve for parameterization, the MATLAB Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) function was used to interpolate the measured 17-angle DoLP angular curve (10°–170° in 10° increment) to a 161-angle DoLP curve (10°–170° in 1° increment). Three parameters were then used to quantify the DoLP shape: the asymmetry parameter, the peak magnitude and the peak angle. The asymmetry parameter g is defined as, g = − ∫ 10 o 170 o P 0 ( ψ ) cos ( ψ ) sin ( ψ ) d ψ , (11) where P 0 ( ψ ) is the interpolated DoLP. When g = 0, the DoLP shape is symmetrical with respect to 90° and g >0 and g < 0 indicate that the DoLP shape is shifted towards the backward or forward direction, respectively. Parameter g can be computed for DoLP generated either by measurements or numerical simulations.

The peak magnitude D o L P m a x is simply the magnitude of the DoLP peak. For the measured 17-angle DoLP, the peak angle falls on the detector angles such as 80°, 90°, 100°, etc. The MATLAB PCHIP interpolation retains the local maximum and minimum of the 17-angle DoLP so the peak angle of the 161-angle DoLP remains the same. To obtain higher angle precision than 10°, ψ m a x was obtained by fitting a skewed Rayleigh DoLP to the 161-angle DoLP. Following (Kokhanovsky 2003), the skewed Rayleigh DoLP is given in this study as: P ( ψ ) = D o L P m a x ⋅ s i n 2 ( ψ − ψ m a x + π 2 ) 1 + D o L P m a x ⋅ c o s 2 ( ψ − ψ m a x + π 2 ) . (12) With D o L P m a x already obtained, ψ m a x is determined by minimizing the RMSE between P ( ψ ) and the measured 161-angle DoLP ( P 0 ( ψ ) ). Graphically, this means horizontally shifting P ( ψ ) to fit the measured DoLP. The RMSE ( ε ) is defined as, ε = ∑ θ = 10 θ = 170 ( P ( ψ ) − P 0 ( ψ ) ) 2 161 . (13) The use of this method is justified by the overall small RMSE between measurements and Eq. 12, and small deviation of ψ m a x from 90° (see Measured Particulate DoLP Angular Shape second paragraph). For DoLPs with high enough angular resolution such as those generated by numerical simulations, there is no need to apply this method.

Light scattering simulations with homogeneous spheres were performed to help interpret observations. Also, we attempt to retrieve the particulate refractive index by fitting simulated DoLP to measurements, i.e., through inversion. In modelling, The bulk particulate DoLP is given by, D o L P p , m o d e l ( ψ , n p , ξ ) = − M 12 ( ψ , n p , ξ ) M 11 ( ψ , n p , ξ ) , (14) where the bulk M 12 and M 11 are computed for Junge type particle size distribution as, M 12 ( ψ , n p , ξ ) = ∫ D m i n D m a x S 12 ( ψ , n p , D ) D − ξ d D , (15a) M 11 ( ψ , n p , ξ ) = ∫ D m i n D m a x S 11 ( ψ , n p , D ) D − ξ d D . (15b) where n p is the particulate refractive index relative to water, n p = n + i n ′ , n and n ′ are the real and imaginary parts. D is the particle diameter, ξ is the slope of the Junge-type particle size distribution. The size averaging was computed in a diameter range from 0.01 μ m to 163 μ m with 195 logarithmically spaced abscissae. Range of ξ was [2,4] in 0.1 increments. Scattering angle ψ range was [0°,180°] in 0.25° increments. In this study, the incident wavelength is 658 nm in vacuum (MASCOT detection wavelength (Twardowski et al., 2012)), that is around 495 nm in water. S 12 and S 11 are outputs from the single scattering models. For the homogeneous spheres, the Lorenz Mie computational program adapted from the Bohren and Huffman formulation (Bohren and Huffman 1998; Twardowski et al., 2001) was used. A Lookup Table (LUT) was constructed with n in the range [1.001,1.3] in 0.002 increments, and n ′ fixed at 0.005 (relative to water).

Figure 1 shows the simulated DoLP angular functions of the homogeneous sphere model for different ξ and n . Figure 2 shows variation of the particle refractive index n with the three shape parameters (peak magnitude D o L P m a x , peak angle ψ m a x , asymmetry parameter g ) of the simulated DoLPs. The overall trend in ξ between 2.8 and 4 is that the DoLP peak decreases and the DoLP shape shifts to large scattering angles (towards ψ ∈ [90°,180°]) with increasing particulate refractive index. Negative branches in DoLP can be spotted near 180° in Figure 1.

Another light scattering model, the asymmetric hexahedral particle scattering model was examined in this study. Originally, the Volume Scattering Function (VSF) from the asymmetric hexahedral model dataset was used in another study to retrieve particle size distributions (Twardowski et al., 2012). In this study, we examine its DoLP element. In the asymmetric hexahedral model, at each particle size, the scattering properties of an ensemble of randomly distorted hexahedra were computed with the Discrete Dipole Approximation (DDA) method (Yurkin and Hoekstra 2007) and the Improved Geometric Optics Method (IGOM) (Yang and Liou 1996). After ensemble averaging at each particle size and a subsequent size averaging, the scattering properties for a polydispersion of particles were obtained. The particulate refractive index range was [1.02,1.2] in 0.02 increments, with the imaginary part fixed at 0.002. The equivalent sphere diameter increased exponentially from 0.01 to 163 μm. Incident wavelength was 658 nm in vacuum ( ≈ 495 nm in water).

Figure 3 panels show the simulated DoLP angular functions of the asymmetric hexahedron model for different ξ and n p . Figure 4 panels show variation of the particle refractive index n with the three shape parameters (peak magnitude D o L P m a x , peak angle ψ m a x , asymmetry parameter g ) of the simulated DoLPs. In contrast to the sphere model, there is little variation of the DoLP shape with the refractive index for ξ values of 3.6 and 4 (Figure 3). At ξ = 4, the DoLP shape does not vary with refractive index, the peak magnitude stays at 1, and the DoLP shape stays symmetrical with respect to 90° throughout the [1.02,1.2] refractive index range. Figure 4 panels show that for ξ ≥ 3.6, there is little variation of D o L P m a x , ψ m a x and g with refractive index. For ξ ≥ 3.6, the asymmetric hexahedra model cannot be used effectively for retrieval of the refractive index. For ξ < 3.6, the hexahedral DoLPs is problematic as well, as there are multiple solutions in some cases. This is in contrast to Figures 1, 2 of the sphere model, where variations in DoLP shape with refractive index are more monotonic and unique. Due to these properties with the hexahedral model, we did not use it for the retrieval part of our study. Nonetheless, the comparison between the DoLPs of the homogeneous sphere and the asymmetric hexahedron model showed the DoLP shape is very sensitive to particle shapes. Thus, for simulation purposes, this brings the question of which particle shape may be the most suitable for reproducing the measured DoLPs in the field. We will discuss this further below.

Figure 5A panels show measured particulate DoLP at each field site. A considerable number of DoLP curves exhibited large, highly intermittent spikes at 30°, 40°, 140°, 150°, and 160°. These isolated spikes at specific angles are generated by particulates such as large organisms and aggregates drifting into the sample volume, and isolated light paths for individual detectors. At those angles, spikes with magnitudes larger than 0.2 were flagged and these DoLPs were removed from the analysis. After that screening, at each field site, outliers with magnitude falling out of the [15%, 85%] inter-percentile range at any one of the 17 scattering angles were removed. This percentile range eliminated DoLPs with potentially high uncertainty and maintained enough DoLP samples for each field sites. The HI data also showed more variability in the backward direction due to the very low measured scattering signals in this clear water; this was from a higher contribution from MASCOT instrument noise. Remaining fluctuations in the DoLP angular curves are from small scale environmental variability in particle fields; note measurements of vertical and horizontal polarized scattering in the sample volume were not made exactly simultaneously.

Figures 5B–G panels show vertical profiles of the three shape parameters, percent RMSE and backscattering parameters at each field site. Linear least-squares fit to the vertical profiles in Figures 5B–G panels are displayed if the coefficient of determination R 2 >0.5. Note that Eq. 12 is introduced to deal with the lower angular resolution of the instrument and to estimate a proper ψ m a x for each measured DoLP as accurately as possible. Horizontally shifting the standard Rayleigh DoLP can bring non-zero values at 0° and 180°. Nonetheless, the RMSEs are generally small (<10%) between measurement and Eq. 12 (see Figure 5E panels), and the obtained ψ m a x typically show small deviation from 90° (within 5°), justifying the use of Eqs 12, 13 to treat this special situation. The peak angle ψ m a x and the asymmetry parameter g will complement each other in the task of quantifying the relative asymmetry of measured DoLP with respect to 90°.

In Figure 5B panels, a relatively high proportion of g are positive among the six field sites. Correspondingly, Figure 5C panels also show a high proportion of ψ m a x that are greater than 90°. This means the measured particulate DoLPs of the six field sites are not symmetric to 90° and generally shift to angular range greater than 90°. This observation agrees with previous measurements (Voss and Fry 1984; Koestner et al., 2020). In Koestner et al. (2020), from the analysis of ocean water samples in the San Diego area, ψ m a x ranged from 91° to 98° with a mean value of 94°. The right-shifted tendency in DoLP is also observed in the simulated DoLP in Modeling, and in other single scattering simulations of the DoLP of marine-like particulates (Lotsberg and Stamnes 2010; Xu et al., 2017; Stegmann et al., 2019). At field sites ES and SBC, D o L P m a x decreases with increasing depth.

Figure 6 shows the box-whisker plots of D o L P m a x , ψ m a x and g of the six field sites. For the six field sites, g generally fell in the range [0,0.04], ψ m a x in the range [90°,95°], with a ensemble-averaged mean of around 92°. D o L P m a x were in the range [0.6,0.9] with a ensemble-averaged mean around 0.7. LE, LS and NYB show very low D o L P m a x , with the minimum at the three sites approaching 0.5. Simulations have shown low DoLP peaks are usually associated with particles with high bulk refractive index, as shown in Figures 1, 2. However, in Zhai et al. (2020), we also showed that a densely-packed aggregate of cyanobacteria cells (with low refractive index) can also have a very low DoLP peak (<0.5) due to the colony structure. The Lake Erie measurements were made when the lake was in a state of cyanobacteria (mainly Microcystis) bloom (Moore et al., 2017). The origin of the low D o L P m a x values for LS and NYB are more difficult to interpret relative to the simulations, as the backscattering ratios at LS and NYB were relatively low (see Figure 5G). A possible explanation could be the complication of having a mixture of multiple particle types contributing to polarized scattering; further investigation is needed to better understand the relative contributions of both high refractive index suspended sediments and low refractive index biological material in the manifestation of DoLP for complex particle mixtures, and some contribution from bubbles may possibly play a role.

DoLPs from all field sites in Figure 5A were averaged into 1 m depth bins and sorted according to their measurement depth to obtain a ensemble averaged DoLP in Figure 7. Figure 7 shows the ensemble averaged particulate DoLP along with the vertical profiles of g , ψ m a x and D o L P m a x . Ensemble average DoLPs mostly exhibited positive asymmetry parameters and peak angles greater than 90°, meaning the DoLPs were generally right shifted. The vertical profile of P m a x shows a discontinuity at around 30 m depth; this jump is an artefact since the water depth at several field sites end at 30 m. Within the depth range [0 m,30 m], D o L P m a x showed decreasing trends with increasing depth.

Figure 8 compares the ensemble averaged DoLP in this study with the Voss and Fry global mean DoLP (Voss and Fry 1984). Figure 8 left panel shows the ensemble averaged particulate DoLP of this study and the seawater DoLP computed with Eqs. 9a–b at a temperature of 20°C and a salinity of 35 ppt. Figure 8 right panel shows the ensemble averaged total (seawater plus particulates) DoLP, Voss and Fry DoLP and seawater DoLP. The red error bars and the grey area indicate standard deviation of our DoLP and the Voss and Fry DoLP, respectively. The pure seawater DoLP has a peak of around 0.92, the total DoLP has a higher peak (0.74) than the particulate DoLP (0.72), but they are both higher than the peak of the Voss and Fry DoLP (0.66). The Voss and Fry DoLP and our total and particulate DoLPs all have highest variability around 90°. Note Voss and Fry DoLP did not converge to zero in far backward angles (beyond about 140°); this may be due to reflection errors related to the cuvette in the bench top apparatus that was used.

As mentioned in Modeling, there is potential with using a DoLP LUT to retrieve the particulate refractive index through inversion. For a measured DoLP with a corresponding PSD slope value, a bulk particulate refractive index can be retrieved by finding a best-fit to the measurement from the simulated DoLP LUT. We attempt to retrieve particulate refractive index with this approach. With the co-located backscattering ratio measurements, another bulk particulate refractive index retrieval approach (Twardowski et al., 2001) (referred to as “the Twardowski model” in texts and figures below) was used for comparison.

At each measured DoLP depth, the Junge-type PSD slope ( ξ ) is estimated from the particulate attenuation ( c p ) spectral curve following (Boss et al., 2001). The measured c p spectrum in the range 440–712 nm can be parameterized by a power-law relationship with a beam attenuation slope γ , c p = c p ( 440 ) ( λ 440 ) − γ . (16) The PSD slope ξ is related to γ by (Boss et al., 2001), ξ = γ + 3 − 0.5 exp ( − 6 γ ) . (17) The Twardowski et al. model relates ξ and backscattering ratio (denoted as b b p ˜ here) to the bulk particulate refractive index n (Twardowski et al., 2001), n = 1 + b b p ˜ 0.5377 + 0.4867 ( ξ − 3 ) 2 [ 1.4676 + 2.2950 ( ξ − 3 ) 2 + 2.3113 ( ξ − 3 ) 4 ] . (18) Figure 9 shows the results of refractive index retrieval at field site East Sound (ES). In Figure 9B, The second panel shows relative errors of the decimal part of the refractive indices retrieved with the two models, ε n p = 100 ( n D o L P − n T w a r d o w s k i n T w a r d o w s k i − 1 ) . (19) Figures 10–14 show the same contents for the other field sites, HI (Hawaii), LE (Lake Erie), Ligurian Sea (LS), NYB (New York Bight) and SBC (Santa Barbara Channel). Overall, at ES, HI, LE and SBC, refractive indices retrieved with the DoLP model and the Twardowski model were roughly comparable. At ES (Figure 9), n retrieved with the DoLP model are mostly lower than the n retrieved with the Twardowski model. The combination of relatively high b b / b and relatively low PSD slope ξ resulted in lower n from the DoLP model compared to the Twardowski model. Moderate to low PSD slope means that larger particles take up relatively more proportion in the particle assemblage. In terms of the quality of the fitting between measured and simulated DoLP in Figure 9A, we see that a lot of the measured DoLPs possess negative branches near 10° and 170°, and the simulated DoLPs generally do not have this feature around the same scattering angles. Simulated DoLPs generally have higher peak angle ψ m a x and higher peaks. At the other field sites, we saw a similar pattern.

At HI, LE and SBC, refractive indices retrieved with the DoLP model and the Twardowski model were comparable. We capped the PSD slope values at 4 for these retrievals because oceanic PSD slope values do not typically exceed 4 in natural waters (Reynolds et al., 2010; their Figure 11), thus removing some potential uncertainty in the derivation of ξ from γ . Among the six field sites, HI and SBC exhibited strongest agreement between the refractive indices retrieved with the two models. One reason for this is the DoLP peaks were generally high (0.7–0.8) at the two sites; another reason is the b b / b and ξ values were in copacetic ranges.

At LE, the DoLP peaks are low, with extreme values approaching 0.5. Thus, the DoLP model returns high refractive indices. Values of b b / b at LE are the highest among the six field sites, so the Twardowski model also returns high refractive indices similar to the DoLP model. At LS and NYB, the situation is in contrast to the aforementioned sites. On one hand, b b / b at LS and NYB are generally lower compared to the other sites, resulting in low refractive indices (<1.1) retrieved from the Twardowski model. On the other hand, a considerable amount of DoLPs possess peaks around 0.6, leading to retrieved n values around 1.15–1.3 for the DoLP model. The extreme value of 1.3 returned by the DoLP model at LS and NYB is due to the combination of high ξ (approaching 4) and low DoLP peak (∼0.6). In Figure 2 first panel, for the ξ = 4 curve, if D o L P m a x is lower than 0.6, the corresponding n approaches 1.3.

Inconsistencies in the models for LS and NYB could be due to several factors not accounted for that may affect the two models differently, including complex particle mixtures (i.e., broad n _p distributions), complex particle shapes, aggregations and bubbles. Further work is needed to interpret the influence of these factors. The Twardowski model and the DoLP model here are both based on the Lorenz-Mie (homogeneous sphere) model, although the key scattering parameters of the two models, DoLP and b b / b , might have different sensitivity to changes in the particle size, shape and refractive index. Systematic study is needed to investigate the feasibility of using non-spherical and inhomogeneous particle shapes such as coated spheres, spheroids, coated spheroids, etc. in DoLP simulations. In Zhai et al. (2020), for cellular contents with very low particle refractive index ( n =0.75 for gas vacuole, n = 1.035 for cytoplasm, n = 1.05 for cell wall), we observed significant decreases in DoLP with increasing aggregate size for densely packed aggregates of layered spheres. At large aggregate sizes (diameter ∼24  μ m ), simulated DoLP peaks decrease to values below 0.5. This aggregate model produces very high backscattering ratios, consistent with the values observed for LE (see Figure 11), and is an example of how accounting for complex particle mixtures and aggregation can significantly affect scattering parameters. Besides modeling issues, the negative branches in measured DoLP curve at around 10° and 170° were the most prominent at NYB among all field sites (see Figure 13A first panel). This feature helped in forming very low DoLP peaks in the simulated DoLPs (see Figure 13A second panel). Negative branches are more prominent at NYB, LS, HI and SBC (see Figure 5A); these locations have clear ocean waters compared to ES and LE. This indicates that in clear waters, MASCOT instrument noise or calibration error could be magnified around 10° and 170°.

Although model discrepancies here require further work to elucidate impacts of these complicating factors, the DoLP model provides another tool to assess oceanic particle composition that may help constrain bulk refractive index estimates. It is also progress toward operational application of remote polarimetry data in determining bulk refractive index and closely related particle density. Particle density is essential for determinations of particle sinking rates (Briggs et al., 2020; Nayak and Twardowski 2020; Omand et al., 2020) and relationships between particulate organic carbon and chlorophyll concentrations (Loisel et al., 2007), and is currently not a parameter that can be derived remotely.