Effects of light polarization and waves slope statistics on the reflectance factor of the sea surface

Above-water radiometry depends on estimates of the reflectance factor ρ of the sea surface to compute the in situ water-leaving radiance. The Monte Carlo code for ocean color simulations MOX is used in this study to analyze the effect of different environmental components on ρ values. A first aspect is examining the reflectance factor without and by accounting for the sky-radiance polarization. The influence of the sea-surface statistics at discrete grid points is then considered by presenting a new scheme to define the variance of the waves slope. Results at different sun elevations and sensor orientations indicate that the light polarization effect on ρ simulations reduces from ∼17 to ∼10% when the wind speed increases from 0 to 14m s−1. An opposite tendency characterizes the modeling of the sea-surface slope variance, with ρ differences up to ∼12% at a wind speed of 10m s−1. The joint effect of polarization and the the sea-surface statistics displays a less systematic dependence on the wind speed, with differences in the range ∼13 to ∼18%. The ρ changes due to the light polarization and the variance of the waves slope become more relevant at sky-viewing geometries respectively lower and higher than 40◦ with respect to the zenith. An overall compensation of positive and negative offsets due to light polarization is finally documented when considering different sun elevations. These results address additional investigations which, by combining the modeling and experimental components of marine optics, better evaluate specific measurement protocols for collecting above-water radiometric data in the field. © 2016 Optical Society of America OCIS codes: (010.5620) Radiative transfer; (010.4450) Oceanic optics; (280.0280) Remote sensing and sensors; (010.5630) Radiometry. References and links 1. P. J. Werdell and S. W. Bailey, “An improved in-situ bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. 98, 122–140 (2005). 2. G. Zibordi, B. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J.-F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. E. Fabbri, S. Kaitala, and J. Seppälä, “AERONET-OC: A network for the validation of ccean color primary radiometric products,” J. Atmos. Oceanic Tech. 26, 1634–1651 (2009). 3. G. Zibordi, “Experimental evaluation of theoretical sea surface reflectance factors relevant to above-water radiometry,” Opt. Express 24, A446–A459 (2016). 4. D. D’Alimonte, G. Zibordi, J.-F. Berthon, E. Canuti, and T. Kajiyama, “Performance and applicability of biooptical algorithms in different European seas,” Remote Sens. Environ. 124, 402–412 (2012). 5. D. D’Alimonte, G. Zibordi, T. Kajiyama, and J.-F. Berthon, “Comparison between MERIS and regional highlevel products in European seas,” Remote Sens. Environ. 140, 378–395 (2014). #256116 Received 4 Jan 2016; revised 12 Feb 2016; accepted 21 Feb 2016; published 5 Apr 2016 © 2016 OSA 18 Apr 2016 | Vol. 24, No. 8 | DOI:10.1364/OE.24.007922 | OPTICS EXPRESS 7922 6. T. Kajiyama, D. D’Alimonte, and G. Zibordi, “Regional algorithms for European seas: a case study based on MERIS data,” IEEE Geosci. Remote Sens. Lett. 10, 283–287 (2013). 7. J. Chowdhary, B. Cairns, F. Waquet, K. Knobelspiesse, M. Ottaviani, J. Redemann, L. Travis, and M. Mishchenko, “Sensitivity of multiangle, multispectral polarimetric remote sensing over open oceans to waterleaving radiance: Analyses of RSP data acquired during the MILAGRO campaign,” Remote Sens. Environ. 118, 284 – 308 (2012). 8. C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements,” Appl. Optics 38, 7442–7455 (1999). 9. B. Bulgarelli and D. D’Alimonte, Optical radiometry for oceans climate measurements (Elsevier, 2014), chap. Simulation of in situ visible radiometric measurements, p. 43, Experimental Methods in Physical Sciences. ISBN10::0124170110. 10. J. Ramella-Roman, S. Prahl, and S. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13, 4420–4438 (2005). 11. T. Harmel, A. Gilerson, A. Tonizzo, J. Chowdhary, A. Weidemann, R. Arnone, and S. Ahmed, “Polarization impacts on the water-leaving radiance retrieval from above-water radiometric measurements,” Appl. Optics 51, 8324–8340 (2012). 12. C. D. Mobley, “Polarized reflectance and transmittance properties of windblown sea surfaces,” Appl. Opt. 54, 4828–4849 (2015). 13. M. Chami, B. Lafrance, B. Fougnie, J. Chowdhary, T. Harmel, and F. Waquet, “OSOAA: a vector radiative transfer model of coupled atmosphere-ocean system for a rough sea surface application to the estimates of the directional variations of the water leaving reflectance to better process multi-angular satellite sensors data over the ocean,” Opt. Express 23, 27829–27852 (2015). 14. T. Kajiyama, D. D’Alimonte, J. Cunha, and G. Zibordi, “High-performance ocean color Monte Carlo simulation in the Geo-info project,” in Parallel Processing and Applied Mathematics, Lecture notes in computer science, vol. 6068, R. Wyrzykowski, J. Dongarra, K. Karczewski, and J. Wasniewski, eds. (Springer, 2010), vol. 6068, pp. 370–379. 15. T. Kajiyama, D. D’Alimonte, and J. C. Cunha, “Performance prediction of ocean color Monte Carlo simulations using multi-layer perceptron neural networks,” in “Procedia Computer Science,” (Singapore, 2011), vol. 4, pp. 2186–2195. 16. T. Kajiyama, D. D’Alimonte, and J. C. Cunha, “Statistical performance tuning of parallel Monte Carlo ocean color simulations,” in “Parallel and Distributed Computing, Applications and Technologies 2012,” (Beijing, China, 2012), pp. 761 – 766. 17. D. D’Alimonte, G. Zibordi, T. Kajiyama, and J. C. Cunha, “Monte Carlo code for high spatial resolution ocean color simulations,” Appl. Optics 49, 4936–4950 (2010). 18. D. D’Alimonte, E. B. Shybanov, G. Zibordi, and T. Kajiyama, “Regression of in-water radiometric profile data,” Opt. Express 21, 27 (2013). 19. S. Chandrasekhar, Radiative transfer (Dover Publications, Inc., 1960). 20. G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphereocean with scattering according to a Rayleigh phase matrix: Effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989). 21. S. Bianchi, A. Ferrara, and C. Giovanardi, “Monte Carlo simulations of dusty spiral galaxies: Extinction and polarization properties,” Astrophys. J. 465, 127–144 (1996). 22. M. J. Raković, G. W. Kattawar, M. Mehru’́ubeoğlu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: Theory and experiment,” Appl. Optics 38, 3399– 3408 (1999). 23. B. G. Henderson, J. Theiler, and P. Villeneuve, “The polarized emissivity of a wind-roughened sea surface: A Monte Carlo model,” Remote Sens. Environ. 88, 453 – 467 (2003). 24. D. Goldstein, Polarized Light, Revised and Expanded (CRC, 2003). 25. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part II,” Opt. Express 13, 10392–10405 (2005). 26. R. W. Goosmann and C. M. Gaskell, “Modeling optical and UV polarization of AGNs. I. Imprints of individual scattering regions,” Astron. Astrophys. 465, 129–145 (2007). 27. C. Cornet, L. C-Labonnote, and F. Szczap, “Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud,” J. Quant. Spectrosc. Radiat. Transf. 111, 174 – 186 (2010). 28. J. E. Hansen and L. D. Travis, “Light scattering in planetary atmospheres,” Space Sci. Rev. 16, 527–610 (1974). 29. V. V. Sobolev, Light scattering in planetary atmospheres, vol. 76 of International series of monographs in natural philosophy (Pergamon, 1975). 30. R. Penndorf, “Tables of the refractive index for standard air and the Rayleigh scattering coefficient for the spectral region between 0.2 and 20.0 μ and their application to atmospheric optics,” J. Opt. Soc. Am. 47, 176–182 (1957). 31. H. C. Van De Hulst, Light Scattering by Small Particles (Structure of Matter Series) (Dover Pubns, 1957). 32. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (John Wiley and Sons, #256116 Received 4 Jan 2016; revised 12 Feb 2016; accepted 21 Feb 2016; published 5 Apr 2016 © 2016 OSA 18 Apr 2016 | Vol. 24, No. 8 | DOI:10.1364/OE.24.007922 | OPTICS EXPRESS 7923 1983). 33. C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. Discuss. 9, 17753–17791 (2009). 34. J. von Neumann, “Various techniques used in connection with random digits. Monte Carlo methods,” Nat. Bureau Standards 12, 36–38 (1951). 35. T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. Oc. 102, 15781 (1997). 36. M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a floating buoy,” in “Ocean Wave Spectra, proceedings of a conference, Easton, Maryland,” National Academy of Sciences (Prentice-Hall, 1963), pp. 111–136. 37. M. Hieronymi, “Monte carlo code for the study of the dynamic light field at the wavy atmosphere-ocean interface,” J. Eur. Opt. Soc, Rapid Publ. 8, 11 (2013). 38. J. Piskozub and W. Freda, “Signal of single scattering albedo in water leaving polarization,” J. Eur. Opt. Soc. Rapid Publ. 8, 13056 (2013). 39. G. Zibordi and K. J. Voss, “Geometrical and spectral distribution of sky radiance—comparison between simulations and field measurements,” Remote Sens. Environ. 27, 343–358 (1989). 40. B. Mayer and A. Kylling, “Technical note: The libRadtran software package for radiative transfer calculations— description and examples of use,” Atmos. Chem. Phys. 5, 1855–1877 (2005). 41. C. Cox


Introduction
Reflectance factor ρ = L r /L i quantifies the fraction of reflected over incident radiance at the seaair interface (L r and L i , respectively).Above-water radiometric measurements require accurate ρ estimates to compute in situ water-leaving radiance L w = L t −ρ •L i , where the total radiance L t is collected in the sea-viewing mode (Fig. 1).In situ above-water systems have been used for the validation of radiometric data produced by space-born sensors like the Sea-viewing Wide Field-of-view Sensor SeaWiFS, the Moderate Resolution Imaging Spectroradiometer MODIS, the Medium Resolution Imaging Spectrometer MERIS and the Visible Infrared Imagery Radiometer Suite VIIRS [1].In the near future, field measurements acquired by the Ocean Color component of the Aerosol Robotic Network AERONET-OC [2,3] will also be applied for the analysis of the Ocean Land Colour Instrument OLCI deliverables.Additionally, above-water systems have relevance for the development of ocean color inversion schemes [4][5][6] and to address environmental monitoring tasks [7].
Radiative transfer studies have shown how the reflectance factor can vary with the sea-state and the viewing geometry [8,9].Look-up tables have then been compiled based on Monte Carlo MC simulations to report ρ values as a function of the wind speed, the radiometers orientation and the sun position [8].Recent extensions also account for the polarization of the incident light [10][11][12][13].The accuracy target for in situ radiometric measurements underpins additional analyses to comprehensively account for specific environmental conditions and extensively assess the performance of measurement protocols.
A new set of Monte Carlo experiments is performed in this study to detail the dependence of ρ values on light polarization and sea-surface statistics keeping the same benchmark for the inter-comparison of results.The first aspect is addressed without and by accounting for the sky-radiance polarization.The generation of the sea-surface at discrete grid points for MC ray tracing is then considered by presenting a new iterative scheme to improve the convergence between the computed and the target slope variance [12].Specific ρ variations due to the light polarization and changes in the sea-surface statistics are then detailed.Example of above-water radiometric measurements.The sensors orientation is that adopted for AERONET-OC [2].
This work relies on the MC code for Ocean Color Simulations MOX to model the reflectance factor in a three-dimensional domain using state-of-the-art High Performance Computing HPC solutions [14][15][16].Above-water simulation results integrate former MOX investigations on the precision of data products derived with in-water radiometric systems [17,18].The study is organized as follows.Methodological aspects are detailed in Section 2. Results are presented in Section 3 and discussed in Section 4. The study summary and conclusions of Section 5 end this work.

Methods
The Stokes vector and the Muller matrix defining light polarization are introduced in Section 2.1.The sky-radiance and the sea-surface MOX components are presented in Section 2.2 and Section 2.3, respectively.

Polarized light
The electric field E expressing light polarization [10,[19][20][21][22][23][24][25][26][27][28] is defined in an orthonormal basis V = {ι ι ι , ι ι ι ⊥ } as with where A and A ⊥ are the amplitudes of the electric field in the ι ι ι and ι ι ι ⊥ direction, respectively; l is the distance along the light propagation direction; λ is the wavelength; ε and ε ⊥ are the phase components (the time dependence has been omitted).Following the schematics of Fig. 2, the ι ι ι vector is oriented in the meridian plane formed by the ray propagation direction ξ ξ ξ ι and the z-axis, while ι ι ι ⊥ satisfies the condition ι ι ι × ι ι ι ⊥ = ξ ξ ξ ι .In this reference system, the light polarization is described by the Stokes vector S S S ι as where δ = ε − ε ⊥ is the phase shift between E and E ⊥ .The I parameter represents the intensity; Q is the excess of polarization in the rather than in the ⊥ direction; U is the excess of polarization in the direction π/4 rather than in the directions 3π/4 with respect to ; and V indicates the excess of right-handed rather than left-handed polarized light.
The ray direction change from ξ ξ ξ ι to ξ ξ ξ χ , upon scattering or reflection, is expressed by the zenith θ * and azimuth angle φ * .Light polarization is then determined in a three-step process: 1) The first step defines the initial Stoke vector S S S ι in the change-of-direction plane determined by the old ξ ξ ξ ι and new ξ ξ ξ χ propagation directions.This is performed rotating counterclockwise (looking along ξ ξ ξ ι ) the basis V ι = {ι ι ι , ι ι ι ⊥ } by the angle φ * to reference S S S ι with respect to 2) The second step is to apply the Muller matrix-i.e., M M M κ (θ * )R R R ι,κ (−φ * )S S S ι -for light scattering or reflection, as detailed in the next Sections.
3) The last step specifies the transformed Stoke vector in the final meridian plane given by the z-axis and the new ray propagation direction ξ ξ ξ χ , which corresponds to the rotation R R R κ,χ (−ψ * ).The new polarization state S S S χ is hence computed as 2.2.Sky-radiance distribution MOX simulations of the sky-radiance distribution are performed at λ = 490nm by considering a plane-parallel, homogeneous and cloud-free atmosphere (e.g., [29]).The atmosphere bottom is idealized as a Lambertian surface with 5% albedo.The working hypothesis is that an accurate sea-surface modeling is required to analyze ρ values, but ρ values have a limited influence on the sky-radiance distribution.This allows for reducing computing time using the same sky-radiance distribution to initialize the photon trajectories when performing ρ simulations at different wind speed values and viewing geometries, while considering the same sun elevation.
The scattering optical thickness is a and τ G a specifying the contribution of the ozone O, water vapor W and permanent gases G, respectively.The absorption a atm and scattering b atm coefficients are defined scaling the corresponding optical thickness values by the atmosphere height and assuming that the scattering rather than absorption determines the optical contribution of the aerosol.Table 1 summarizes the parameters used in the MOX sky-radiance component.
The path length l of the photon trajectory is l = − log(u)/c atm (7) where c atm = a atm +b atm is the attenuation coefficient and u is a random number sampled from a uniform distribution between 0 and 1; i.e., u ∈ U(0, 1).The scheme to update the photon weight is w new = w old • ω atm (8) where ω atm = b atm /c atm is the single scattering albedo.This is equivalent to consider an initial ensemble of virtual photons and reduce its size based on the probability of each photon to be scattered.The ray-tracing of photons ending their trajectory in the atmosphere follows the rules of Rayleigh scattering for gas-molecules if u > µ, where u ∈ U(0, 1) and µ = τ M b /τ b (Section 2.2.1).The Mie theory instead applies for the aerosol scattering if u ≤ µ (Section 2.2.2).
The intensity and polarization of each virtual photon that arrives at the atmospheric bottom are recorded as a function of the propagation direction to define the sky-radiance distribution (Fig. 3).The photon is then redirected towards the sky.Photons reaching the atmosphere top are instead terminated.Simulation results include the diffuse-to-total r tot and the diffuse-to-direct r dir irradiance ratio and where E tot = E sun + E sky is the sum of sun and sky irradiance (E sun and E sky , respectively).Note that r dir can be expressed as a function of r tot [Eq.( 9)] as

.1. Muller matrix of Rayleigh scattering
The Muller matrix for Rayleigh scattering RLG due to anisotropic particles in random orientation is expressed as where being δ the depolarization factor [28,30].

Muller matrix of Mie scattering
The Mie theory models the scattering by an homogeneous spherical particle.The Muller matrix is formulated as where the amplitudes S 1 (θ * , n c , d) and S 2 (θ * , n c , d) are complex functions of the zenith scattering angle θ * , the complex refractive index n c = n r − in i , and the size parameter d = 2πr/λ defined by the particle radius r and the light wavelength λ .A description of the S 1 (θ * , n c , d) and S 2 (θ * , n c , d) functions is presented in the literature [19,28,[31][32][33].

Reflectance factor of the sea-surface
Upon generating the sky-radiance distribution L sky , the simulation of the reflectance factor ρ starts by sampling a random number q from an uniform distribution between 0 and 1: the photon origin is the sky if q < ρ sky , otherwise is the sun.The ρ sky threshold is given by with where the sun and the sky weighting factors W sun and W sky are and The direction cosines of photons that represent the direct light are set based on the sun zenith θ s and azimuth φ s , and also considering the sun angular size.The direction and the Stokes parameters of the photons that constitute the diffuse light are otherwise determined from skyradiance simulation results.The acceptance-rejection method [34] is in this case adopted for sampling from the two-dimensional distribution of tabulated L sky values.Photons leaving the simulation domain laterally are restarted from the opposite side [Fig.4(a)] assuming periodic boundary conditions.
The L r and L i fields are determined by tracking the photons trajectories.Each radiometric quantity is associated with an l-by-m-by-n matrix whose (i, j, k) entry accumulates values detected by the (i, j)-th photon collecting bin of the k-th horizontal layer.By indicating with ξ ξ ξ L i the direction of the radiometric sensor in the sky-viewing mode, the photon weight is added to L i (i, j, k) when the ray intersects the (i, j)-th bin of the k-th layer satisfying the constraint −ξ ξ ξ L i • ξ ξ ξ < cos(FOV), and an analogous scheme applies for L r .A single layer of collecting bins positioned at 10m above the sea level and a FOV of 5 • are considered in this work.

Sea-surface modeling
The sea surface is modeled following the scheme presented by Mobley [12], which provides a detailed description of the applied methods.The sea surface is discretized by N x and N y grid points along the xand y-axis respectively (Table 2, Fig. 4), with resolutions ∆ x = L x /N x and Table 2. MOX parameters for simulating the reflectance factor ρ.
where the harmonic components ẑ(u, v) are indexed by u and v in the wavenumber k domain.
Computing ẑ(u, v) requires two main quantities.One is the Elfouhaily omnidirectional spectral density S(k) [12,35], which determines the elevation variance based on the module of the wind speed and the wave age.The other is the "cosine-2s" spreading function [12,36], employed to distribute the sea-surface elevation to waves traveling in different directions.The wave age and spreading function parameters adopted for all ρ simulations considered in this work are respectively Ω c = 0.84 (i.e.,fully developed sea) and S= 2.8.In summary, the steps to generate the sea surface though Eq. ( 18) are: 1) compute the twosided discrete values of the elevation variance; 2) sample the amplitudes of the harmonic components with normal distributions; and 3) define the Hermitian coefficients of the IFFT to obtain real-valued z(r, s) elevations.An example of sea-surface generation is shown in Fig. 4(b) for a wind speed of 10 ms −1 .

Muller matrix of the reflected light
The incidence and transmission angles of a collimated light beam with respect to normal at the sea-surface intersection point are respectively denoted θ i and θ t , whereas their sum and difference are θ + = θ i + θ t and θ − = θ i − θ t .The Muller matrix expressing the polarization of the reflected light RFL is then [20,24]

Statistical figures of the sea surface
The main effect of the elevation variance on ρ values is shading portions of the sea surface to the light that propagates almost horizontally, and its relevance increases with the sun zenith [37].
The sea-surface slope variance has instead a direct influence on ρ because the angle of incidence θ i contributes to define the coefficients of the Muller matrix for light reflection [Eq.(19)].
In the ideal case, the target TRG values of the elevation ELE variance [σ ELE TRG ] 2 and the slope SLP variance where k FND is the fundamental wavenumber, k H denotes a wavenumber sufficiently high to account for the contribution of capillary waves, and assuming that the simulation domain includes the sea-surface wave with the largest wavelength of interest [12].Note that S(k) has different signification in Eq. ( 20) and (21).Specifically, in Eq. ( 20), it represents the spectral elevation variance in m 3 rad −1 , and hence [σ ELE TRG ] is in m 2 .Instead, the S(k) term of Eq. ( 21) is obtained upon spatial differentiation to evaluate the sea-surface slope statistics.Its units are m 3 rad −1 rad −2 (the rationale to add the extra rad −2 term is detailed below).This leads to a spectral slope variance k 2 S(k) in mrad −1 (not mrad as reported in the literature; e.g., [35] pag.15,795).And indeed the [σ SLP TRG ] term quantifies the variance of a slope rather then of an angle.Without loss of generalization, the sea-surface is defined for this dimensional analysis as z(x) = ∑ ∞ n=0 z n (x), where z n (x) = A n cos(k n x + φ n ).The term cos(k n x + φ n ) is dimensionless.The slope of this n-th harmonic component, computed with the chain rule for differentiation, is , which is the variance of a slope (i.e., not an angle) although the right term seems expressed in rad 2 at a first sight.The reasoning to formulate Eq. ( 21) in the correct units is analogous once the spreading function is accounted for relating S(k) to the Fourier coefficients [12].
The sea-surface statistics at high wavenumbers The IFFT [Eq.(18)] accounts for harmonic components with wavenumber in the range where k NYQ is the Nyquist value [12].Most of the elevation variance can be included in the discretized sea surface by choosing a proper grid resolution.It is however difficult to obtain the same result for the slope variance, which can still be relevant at high wavenumbers due to the k 2 term of Eq. ( 21).An efficient method has been presented in the literature [12] to modify the S(k) function and preserve the sea-surface slope variance at wavenumbers larger than k NYQ .
The scheme to model the High-Wavenumbers Statistics HWS of the sea surface defines a modified spectral density function Ŝ(k) as where and with k P indicating the peak frequency where the S(k) spectral density reaches its highest value.This approach ensures that The slope variance of the sea surface at the (r, s) grid points upon applying the standard HWS scheme of Eqs.( 22)-( 24) is where and ∆ yz (r, s) = z(r, s + 1) − z(r, s) y(s + 1) − y(s) .
The sea-surface statistics at the grid points of the simulation domain The HWS guarantees the validity of Eq. ( 25) in the continuous case.Numerical results however indicate that an underestimation tendency 26) and Eq. ( 21)] characterizes the IFFT seasurface realization [Eq.( 18)] due to truncation errors induced by the finite difference method.A novel algorithm is proposed in this work to model the Grid-Points Statistics GPS of sea surface by iteratively adjusting the spectral density function Ŝ(k) until converging to the target value [Eq.(21)].
2) Define the amplitude of the harmonic components directly from ŜGPS (k)-i.e., without sampling-and generate the sea surface with the IFFT scheme.

5) Repeat items 1 to 4 or terminate if ν ≥ 1
The η parameter is employed in step 4 as a trade-off for balancing the number of iterations and the accuracy of results, and an example is reported later when evaluating the algorithm performance.Once the δ GPS value has been set, the sea-surface can then be generated by sampling normally distributed elevation components as for the HWS scheme.

Variance of the sea-surface elevation and slope
The acronyms designating the main quantities considered in this analysis are highlighted below to facilitate the discussion of the results:   5 adopting the simulation setting of Table 2 and considering wind speed values v w in the range 2-14 m s −1 .In agreement with [12], Fig. 5(a) shows that both the HWS and the GPS schemes well approximate the target elevation variance (the results variability at v w = 14 m s −1 is the effect of statistical fluctuation).Instead, Fig. 5 2 even at low wind-speed values.The proof-ofconcept GPS algorithm presented in this study complements the HWS scheme demonstrating the possibility to reach a statistical convergence between the slope variance computed at the sea-surface grid points and the target value.The value of the η parameter employed in this case study is 0.02.Less than 20 iterations and a processing time of a few tens of seconds are required for generating each sea surface.The GPS convergence can be improved by updating the η value at run time, although this optimization is out of the scope of the present study.

Overall tendencies of the reflectance factor of the sea-surface
The ρ simulations for the azimuth viewing geometry φ v = 90 • with respect to the sun plane are summarized in Fig. 6.This azimuth angle is selected in agreement with the measurement protocol of AERONET-OC [2,3] and MC simulation results [38].Each panel presents results for θ v = 0 • , 10 • , 20 • , . . ., 80 • with respect to zenith (the same values define the sea-viewing mode with respect to nadir).Wind speed values are in the range 0-14m s −1 .The column panels from left to right refer to θ s = 0 • , 30 • and 60 • , respectively.The row panels from top to bottom instead consider the HWSUNP, HWSPOL, GPSUNP and GPSPOL simulation cases (Table 3).Repeated ρ simulations display differences due to the MC intrinsic noise in the order of 1% (details not presented).Results have been additionally verified through replicated sea-surface generations.The ρ values analyzed in this work for different environmental conditions (i.e., selected sun zenith and wind speed values) meet ranges documented in the literature [12].Trends of Fig. 6 indicate a significant ρ variability, as well as a strong dependence on the wind speed values, when the L i radiometer is oriented close to zenith.This depends on the sun-glint contribution, as further discussed in Section 4. Larger viewing angles lessen the glint effect, and the overall ρ variability at different wind speed and sun elevations reaches a minimum at θ v ∼ 65 • .The ρ dependency on θ v indicates larger values when the radiometer orientation approaches the horizon due to the enhanced Fresnel reflectance.
Results indicate that increasing the wind-driven slope variance of the sea-surface corresponds to: • lower ρ values for θ v < 20 • and θ s ∼ 0 • .
These tendencies are summarized in Fig. 7.
Selected cases addressed in Fig. 6 are compared against each other in Fig. 8. Results at different wind speed values are highlighted adopting the same color scheme of Fig. 6.The triangle , the circle and the square symbols denote results for θ s = 0 • , 30 • and 60 • , respectively.The benchmark quantity in the abscissa axis is the ρ value obtained without accounting for the light polarization and using the HWS scheme to generate the sea surface.The ordinate axis refers instead to ρ simulations performed by: 1) including light polarization (HWSPOL, Fig. 8(a)); 2) applying the iterative algorithm to generate the sea surface (GPSUNP, Fig. 8(b)); and 3) jointly considering the effect of light polarization and the sea-surface modeling revision (GPSPOL, Fig. 8(c)).
Fig. 8. Scatter plots for different ρ simulation settings.The benchmark case is HWSUNP.The ρ variations induced by the light polarization (HWSPOL), the use of the novel iterative algorithm for the sea-surface generation (GPSUNP), and their combined effect (GPSPOL) are highlighted in panel (a), (b) and (c), respectively.Colors identify results at different wind speed values as in Fig. 6.
The scattering ε and bias δ between the compared quantities are where N = 216 is the number of data points (i.e., 8 wind speed values, 3 sun elevations and 9 viewing geometries).The determination coefficient r 2 is also reported.The polarization effect of on ρ simulations detailed in Fig. 8(a) indicates an average scattering ε = 10% and almost no bias.Accounting for light polarization produces lower ρ values at θ s = 0 • .This tendency is reversed for θ s = 60 • .The effect of light polarization reduces for θ s = 30 • .
The comparison of results obtained by generating the sea surface with the HWS and the GPS scheme is presented in Fig. 8(b).The overall scattering and bias are ε = 5% and δ = 3%, respectively.The tendency ρ HWSUNP > ρ GPSUNP can be observed for ρ values below 0.15, while ρ HWSUNP < ρ GPSUNP appears for ρ > 0.15.This dual feature is likely a consequence of the manifold effect of the sea-surface statistical figures on ρ when considering different viewing geometries and sun elevations, as summarized in Fig. 7.
The composite effect of light polarization and the use of the GPS scheme to generate the sea surface is addressed in Fig. 8(c).Results indicate an increase of both scattering and bias (i.e., ε = 13% and δ = 4%).
An additional analysis on the influence of the wind speed is performed for the θ v = 40 • and φ v =90 • viewing geometry by averaging differences between ρ simulations obtained at different sun elevations (Fig. 9).Results of Fig. 9(a) for the HWSPOL vs. HWSUNP case indicate that raising v w from 0 to 14m s −1 reduces the light polarization effect on the reflectance factor from ∼ 17 to ∼ 10%.This can be explained considering that, when the variance of sea-surface slope becomes larger, a bigger fraction of the sky dome with different levels of polarization is reflected in the field of view of the radiometer measuring L t , which translates in a depolarization tendency.
The GPSUNP vs. HWSUNP comparison of Fig. 9(b) shows that differences between ρ simulations increase up to δ = 12% when the wind speed reaches v w = 10m s −1 , and then reduce slightly.This tendency is likely due to the fact that ρ variations becomes lower once the wind Colors identify results at different wind speed values as in Fig. 6.
speed reaches a certain level, as it can be noticed by carefully observing results of Fig. 6.The comparison of GPSPOL and HWSUNP results shown in Fig. 9(c) jointly consider the effects of light polarization and the use of the GPS scheme to generate the sea-surface.The differences between ρ simulations display a less systematic dependence on the wind speed, with values in the range ∼13 to ∼18%.The overall trend addressed in Fig. 9(c) is approximately the composition of results presented in Fig. 9(a) and 9(b).Note however that ρ differences between the GPSPOL and the HWSUNP simulation results tend to be lower than the sum of the offsets documented for the HWSPOL vs. HWSUNP and the GPSUNP vs. HWSUNP cases.This likely depends on the depolarization effect due to the larger slope variance of the sea-surface generated with the GPS with respect to the HWS scheme.

Low wind speed case
A wind breeze of about 6 m s −1 starts the formation of white caps, limiting the accuracy of in situ measurements of the water-leaving radiance.The effect of white caps is also not taken into account in the present work.A case study for wind speed values between 2 and 6 m s −1 is then presented.The top panels of Fig. 10 highlight the comparison of the HWSUNP and GPSUNP results.The panels of the central row refer to the HWSUNP and HWSPOL cases.The comparison between the HWSUNP and GPSPOL simulations is finally addressed in the bottom panels.The column panels from left to right are addressed to θ s = 0 • , 30 • and 60 • , respectively.The ρ variations between the considered cases confirm the complexity of details that emerge when light polarization, or refinement of the sea-surface statistics, are taken into account even at relatively low wind speed.
To better identify overall tendencies, results of Fig. 10 are summarized in Fig. 11 by considering the mean ρ difference at wind speed values between 2 and 6 m s −1 .Figure 11(a) shows that light polarization induces larger variations at θ v ≥40 • .It is also confirmed that the polarization effect is reduced when considering negative and positive differences between ρ simulations performed at different sun elevations.Larger ρ differences are instead reported at θ v ≤ 40 • in Fig. 11(b) for the comparison of ρ values obtained using the HWS and the GPS scheme to generate the sea surface.Figure 11(c) shows how the results presented in Fig. 11 methods for the regression of optical profile values [18].In the present study, this code has been featured with new functionalities to model the sky-radiance distribution and undertake abovewater radiometric simulations.Atmospheric parameters have been set in agreement with an analytical scheme validated with experimental data [39].The simulated sky-radiance patterns are similar to those reported in the original study (details not shown).By the same token, the Mie scattering, included for a more realistic representation of light polarization effects, provides patterns analogous to published results [40].
As in former investigations [14][15][16], high-performance computing methods have been strategic to satisfy Monte Carlo simulation requirements by relying on a large-scale computer cluster for production runs.The Navigator supercomputer, University of Coimbra, Portugal, has been utilized in the present study.The number of virtual photons traced to simulate each reflectance factor is N ph = 10 9 , and a summary of execution times for different values of the wind speed and sun elevation is presented in Fig. 12 (i.e., HWSUNP case).
Ray-tracing has been executed excluding direct sun-light contribution in L i .Both the diffuse and the direct sky-radiance inputs to L r are instead accounted for.This leads to large ρ values for the θ s = 0 • , v w = 0m s −1 and θ v = 0 • (left column panels of Fig. 6).An example of ρ simulations adding the direct radiance to L i for θ s = 0 and θ v = 0 is presented in Fig. 13 (i.e., HWSUNP case).Lower values now characterize ρ simulations for θ v = 0 and v w > 0 m s −1 .This can be explained considering that waves at the sea surface significantly lessen the sun-glint contribution to L r even at low v w , while the L i sensor records the direct light.The conclusions of this work do not depend on the sensor response in the zenith orientation.The choice to exclude the direct light from L i records (which only matters for θ s = 0 • and θ v = 0 • ) is then justified.
The HWS scheme offers the advantage to update the slope variance analytically, but without accounting for the specific case where the sea surface is generated at discrete grid points.The present study has shown that complementing the HWS scheme with the GPS iterative algorithm enables a better convergence to the target value of the sea-surface slope variance.Considering for instance ρ simulations at v w = 6m s −1 , the target values based on the Elfouhaily omnidirectional spectral density variance S(k) [35]  GPS ] 2 =0.0576m 2 and [σ SLP GPS ] 2 = 0.0366, showing an improved convergence with the target slope variance.Note that the HWS and GPS values can be subject to some statistical variability due to the sampling adopted to generate the sea-surface with the IFFT scheme.This study concerns a fully developed sea (i.e., Ω c = 0.84) and follows the parameterization adopted in [12] for the implementation of the Elfouhaily model [35].The similarity between the sea-surface statistics considered here and Cox-Munk results [41], which probably refer to a mature sea state [12], would then improve upon setting Ω c =1. Besides, slightly varying agreements between the Cox-Munk and Elfouhaily slope variance are documented in the literature (e.g., [42,43]).Preliminary analyses performed in this work (details not presented) indicate that the convergence between the two models further enhances defining the friction velocity at the sea surface through Eqs. ( 60) and (66) of [35] (e.g., as in [44]), rather than Eq. ( 61) of [35] (as both in [12] and in the current study).These elements address future analyses of the sea-surface statistics under different modeling assumptions.
Additional accuracy improvements will be part of future MOX developments.It has been observed that the omnidirectional density function of the sea-surface elevation variance depends on the age of the wind-generated waves.The spreading function [12,36] adopted in the two-dimensional case can also influence ρ values.The definition of the sea-surface statistics only based on the wind speed is then an approximation of the environmental conditions, and this contributes to the uncertainty budget of L w values determined with above-water systems in the field.The synergy between theoretical and experimental marine optics is needed for understanding better these aspects, as well as to evaluate additional uncertainties depending on measurement protocols.Analyses shall also be conducted to account for the elevation at which the system is deployed above the sea level, as well as the field-of-view of the radiance sensors.Radiative transfer simulations at different wavelengths and atmospheric optical properties are relevant as well.On a longer term, MOX extensions are planned to address the closure problem by comparing L w results obtained simulating the above-and in-water measurements protocols.

Summary and conclusions
The look-up tables of ρ values to derive the in situ water leaving radiance from above-water radiometric measurements were originally implemented based on the uniform sky-radiance distribution and adopting the Cox-Munk parameterization of the sea-surface elevation [8,41].These tables have been recently re-formulated taking into account the polarization of an idealized Rayleigh sky, as well as using the IFFT method for modeling the sea-surface elevation [12,28,30].The look-up table revision document has presented ρ changes by considering the joint effects of these updates [12].The rationale of the present work was performing a comparative study using a common reference to investigate the specific influence that light polarization and the sea-surface statistics can have on the ρ factor.A dimensional analysis has  also allowed for verifying that the spectral slope variance is expressed in rad m −1 , and not in radm as reported in the literature (e.g., [35]).
Test cases indicate that differences between ρ simulations performed without or by taking light polarization into account tend to reduce when considering their overall effect for different sun zenith values.The relevance of light polarization would then increase in the case of a limited number of field measurements systematically collected at high or low sun elevations.Additionally, the effect of polarization on ρ values determined at the θ v = 40 • and φ v = 90 • viewing geometry by averaging results for different sun elevations is lessened when the wind speed increases.The opposite wind speed effect on ρ differences instead appears when including the GPS scheme to adjust the sea-surface slope variance at high wavenumbers.Underestimated ρ values correspond to overestimated L w data derived from in situ measurements.And this highlights the importance of the GPS algorithm from an operational perspective.
In conclusion, this work confirms the capability of the MOX code to perform state-of-thearts marine optics investigations and provide relevant information for the collection of abovewater radiometric measurements in the field.Additional analyses, joining the modeling and experimental study components, are then planned to better understand uncertainty budgets by considering different measurement protocols for the data collection in the field.
Fig.1.Example of above-water radiometric measurements.The sensors orientation is that adopted for AERONET-OC[2].

Fig. 2 .
Fig.2.The Stokes vector is at first referenced to the initial meridian plane specified by the base V ι = {ι ι ι , ι ι ι ⊥ }.Before applying the Muller transformation, light polarization needs to be expressed in the basis V κ = {κ κ κ , κ κ κ ⊥ }, where the κ κ κ vector is within the change-ofdirection plane.The new polarization state has afterwards to be formulated in the basis V χ = {χ χ χ , χ χ χ ⊥ }, where χ χ χ is in the final meridian plane.Each change of base is a rotation with a different Euler angle.The figure refers to the reflection of light at a facet of the sea surface (θ i indicate both the incident and the reflection angle with respect to the normal n n n).The schematic for the photon scattering is analogous.

Fig. 3 .
Fig.3.Sky-radiance simulations in the panels from left to right refer to θ s = 0 • , 30 • and 60 • , respectively.The I, Q, U Stokes parameters (in units of W m −2 nm −1 sr −1 ) and the polarization degree P (in percent) are presented from the top to the bottom row.Results are scaled so that the total irradiance is 1 Wm −2 nm −1 for θ s = 0 • .

Fig. 4 .
Fig. 4. Schematic of photon tracing and an example of sea-surface generated at the grid points of the simulation domain (v w = 10m s −1 ) in panel (a) and (b), respectively.

Fig. 5 .
Fig.5.Elevation and slope variance of sea surfaces generated for wind speed values in the range 2-14 m s −1 are in the left and right panel, respectively.The legend acronyms are defined as follows: TRG denotes target quantities; HWS and GPS refers to the highwavenumbers statistics and to the grid-points statistics, respectively; RAW indicates the case where no correction is applied.

Fig. 9 .
Fig. 9. Summary of ρ variations when considering the θ v = 40 • viewing geometry and averaging results for θ s = 0 • , 30 • and 60 • .Effects due to the light polarization, the seasurface statistics, and their combination are presented in panel(a), (b) and (c), respectively.Colors identify results at different wind speed values as in Fig.6.

Fig. 11 .
Fig. 11.Summary of ρ variations by averaging results for v w = 2, 4 and 6 ms −1 .The effect of light polarization, sea-surface statistics and their composition are considered in panel (a), (b) and (c), respectively.

Fig. 12 .
Fig. 12. Execution times for simulating the reflectance factor as a function of θ s and v w .Results refer to the tracing of 10 9 photons for production runs of the HWSUNP case using 24 CPU cores on the Navigator supercomputer, University of Coimbra, Portugal.

Fig. 13 .
Fig.13.Example of ρ simulations for θ s = 0 • , v w = 0 m s −1 and θ v = 0 • by including the direct light contribution to L i .The symbol and color scheme for results at different wind speed v w values are the same as in Fig.6.See text for details.
where τ M b and τ A b indicate the molecules M and aerosol A contribution, respectively.The absorption optical thickness is instead τ a

Table 3 .
Test cases to evaluate the variability of ρ simulation results.