1. Introduction and background

Ocean color remote sensing requires in situ reference values of the spectral water-leaving radiance ${L_w}$ for the development of bio-optical algorithms, the validation of radiometric data products, or the indirect calibration (i.e., system vicarious calibration) of space sensors [1–3]. The in situ measurements necessary to compute ${L_w}$ can be acquired with in- or above-water methods. In-water methods determine ${L_w}$ from the regression of the upwelling radiance measured at different depths in the water column [4–6]. Above-water methods determine ${L_w}$ from sea-viewing measurements of the total radiance ${L_t}$ and estimates of the radiance reflected by the sea surface ${L_r}$ according to

The reflected radiance ${L_r}$ in any direction $(\theta ,\phi )$ is the result of the reflection by the sea surface of the sky-radiance ${L_{sky}}$ originating from the upper hemisphere including both the direct and the diffuse components (leading to sun- and sky-glint, respectively). By indicating with $r({{\theta^ \ast },{\phi^ \ast } \to \theta ,\phi } )$ the sea-surface radiance reflectance function, i.e., the fraction of incident radiance reflected by the sea surface from direction $({\theta ^ \ast },{\phi ^ \ast })$ to direction $(\theta ,\phi ),\,{L_r}$ is given by

The theoretical values of the reflected radiance L_r and of the incident radiance L_i, needed to compute the sea-surface reflectance factor through Eq. (2) and Eq. (4), are derived from sea-surface and the sky-radiance models [9–11]. It is emphasized that Monte Carlo (MC) ray-tracing simulations permit to compute $\rho $ accounting for multiple reflections at the sea surface [12–17].

Investigating the replicability of $\rho $ as a function of various sky-radiance and sea-surface models is the primary objective of this study. The Monte Carlo code for Ocean Color Simulations (so-called MOX) is utilized [18,19] restricting the analysis to clear sky conditions and accounting for the spatially averaged variability of sea-surface statistics.

The considered measurement geometry is that recommended by protocols for the validation of satellite ocean color data [7] and applied by the Ocean Color component of the Aerosol Robotic Network (AERONET-OC) [20]. As illustrated in Fig. 1, the measurement angles are $\phi = 90^\circ $ with respect to the sun azimuth, $\theta ^{\prime} = 40^\circ $ with respect to the zenith and $\theta = 180^\circ{-} 40^\circ $. Alternative geometries could certainly be considered and might provide relative advantages. For instance, choosing $\phi = 135^\circ $ instead of $\phi = 90^\circ $ might further minimize glint perturbations, but it would concurrently increase measurement perturbations induced by the deployment structure, which increase with the sun zenith.

 

Fig. 1. (a) Perspective-, (b) side- and (c) top-view of typical above-water measurement geometry. Note that the sea- and sky-viewing angles are defined with respect to the zenith.

Download Full Size | PDF

Look-up-tables of $\rho $ (e.g., [12,15]) allow deriving ${L_\textrm{w}}({\theta ,\phi } )$ from operational measurements of ${L_\textrm{t}}({\theta ,\phi } )$ and ${L_i}({\theta^{\prime},\phi } )$. The determination of ${L_\textrm{w}}({\theta ,\phi } )$ from in situ data presents, however, some caveats. In particular, L_t and L_i measurements exhibit dependence on sensor features such as field-of-view and sampling time. Additionally, L_i measurements from the direction $({\theta^{\prime},\phi } )$ might not accurately represent the ${L_{sky}}$ contributions leading to ${L_r}$. Consequently, ${L_\textrm{w}}$ determined from field measurements through Eq. (1) and Eq. (3) implies dependence on instrument features, measurement geometry, and obviously, sky-radiance distribution ${L_{sky}}$ and sea-surface radiance reflectance function r.

A core component of this study is the direct comparison between $\rho $ values computed with MOX and those included in the widely used $\rho $-tables generated with Hydrolight [12,21]. This benchmark analysis is performed assuming identical sky-radiance distributions and sea-surface statistics, consistently determined at the 550 nm center-wavelength neglecting polarization. Cross-comparisons based on independent data sources [9,10] are also presented.

It is emphasized that this work strictly focuses on: 1) the quantification of the replicability of ρ-factors, and 2) the assessment of processing methods and results through benchmarks. More general issues of relevance for above-water radiometry, such as the impact of measurement variability on the determination of ${L_t}$ values representing the actual observation conditions, the dependence of the sea-surface reflectance factor on wavelength and polarization [22], the evaluation of alternative methods for the reduction of above-water radiometric data [23,24], or even the provision of practical recommendations on measurements protocols, are all beyond the study objectives.

The structure of the work is as follows. Section 2 describes the simulation framework, the sky-radiance and sea-surface models, and the processing scheme to compute ρ. Section 3 documents the ρ-factor replicability and benchmark results. Findings from benchmark analyses and the impact of the study on the reduction of field measurements for the determination of ${L_w}$ are discussed in Section 4. Summary and remarks in Section 5 complete the work.

2. Methods

2.1 Simulation framework

A discretization of the upper hemisphere into quadrilaterals, called quads, is used in Hydrolight to perform MC ray tracing, collect simulation results and compute $\rho $. To simplify comparisons, MOX adopts the same scheme as Hydrolight. The partition accounts for ${N_\theta } \times {N_\phi }$ intervals in the zenith and azimuth direction, respectively (see Fig. 2). The i-th increment in the zenith direction is $\Delta {\theta _i}$, for $i = 1, \ldots ,{N_\theta }$ with $\mathop \sum \nolimits_i^{{N_\theta }} \Delta {\theta _i} = \pi /2$ (i.e., the $\mathrm{\Delta }{\theta _i}$ values can be different). The reference discretization is ${N_\theta } = 10$ with $\mathrm{\Delta }{\theta _i} = 5^\circ $ if i=1 or i=10, and $\mathrm{\Delta }{\theta _i} = 10^\circ $ otherwise. The j-th increment in the azimuth direction is $\mathrm{\Delta }{\phi _j} = 2\pi /{N_\phi }$ (i.e., the $\mathrm{\Delta }{\phi _j}$ segments are equal) for $N_\phi = 24$ and $\mathrm{\Delta }\phi = 15^\circ $. Grid points have zenith coordinates ${\theta _i} = \mathop \sum \nolimits_{i^{\prime} = 1}^{i - 1} \Delta {\theta _{i^{\prime}}}$, with ${\theta _1} = 0^\circ $ and azimuth coordinates ${\phi _j} = \Delta \phi \cdot ({j - 1/2} )$ for $j = 0, \ldots ,{N_\phi }$. Hydrolight unifies in a single polar cap the 24 slices converging to $\theta = 0^\circ $, while MOX keeps them separated (which is equivalent in terms of $\rho $ computation).

 

Fig. 2. Partitioning of the upper hemisphere into quads. The thick solid and dashed lines indicate the measurement and the solar plane, respectively. The sensor-quad marked in black represents the virtual radiometer.

Download Full Size | PDF

To define a correspondence between field measurements and the simulation framework, the sun azimuth ${\phi _0}$ is set to $0^\circ $ with the solar plane crossing the hemisphere from $\phi = 0^\circ $ to $\phi = 180^\circ $. The measurement plane defined by the radiometer azimuth is determined by $\phi = 270^\circ $ and $\phi = 90^\circ $ (note that the assumption of ${\phi _0}$ = $0^\circ $ implies the equivalence of the relative azimuth $\phi $ with the azimuth identified by the measurement plane). The quad including $\theta = 40^\circ $ and $\phi = 270^\circ $ represents the virtual radiometer and is denoted as the sensor-quad.

Following the scheme implemented in Hydrolight [21], the sky radiance assumed homogeneous within each quad, is computed with the models presented in Sec. 2.2 for directions $({\theta_i^\ast ,\phi_j^\ast } )$, where $\theta _i^\ast\,{=}\, {\cos ^{ - 1}}({0.5 \cdot ({\cos {\theta_{i - 1}} + \cos {\theta_i}} )} ),\,\,\phi _j^\ast\,{=}\, 0.5 \cdot ({{\phi_j} + {\phi_{j + 1}}} )$, and (${\theta _i},\,{\phi _j}$) indicate the hemisphere grid coordinates. The quad including the sun zenith ${\theta _0}$ and sun azimuth ${\phi _0}$ angles is referred to as the sun-quad. The sea surface below the sky hemisphere permits simulating the light reflection according to the model representing the statistical properties of waves at the air-water interface (Sec. 2.4). Specific schemes adopted by Hydrolight and MOX for MC ray tracing are presented in Sec. 2.5 while detailing methods for the computation of $\rho $.

2.2 Sky-radiance models

The angular distribution of the sky radiance ${L_{sky}}$ is simulated with: 1) the Harrison-Coombes parametric model (HC) [25–27] embedded in the Hydrolight code, 2) the analytical model applied by Zibordi and Voss (ZV) providing realistic distributions in clear sky conditions [28], and 3) the highly accurate Finite Element Method radiative transfer numerical code (FEM) matter of several benchmarks and applications [29–32]. These methods are briefly introduced in the following sub-sections together with details on their implementation. The irradiance model from Gregg and Carder (GC) [33] is additionally applied to rescale the diffuse sky irradiance in HC and determine the direct radiance component for both HC and ZV. The main quantities necessary for the implementation of the various modeling solutions are summarized in Table 1 specifying their relevance for each scheme (i.e., HC, ZV and FEM).

Table 1. Main quantities applied for the simulation of sky-radiance distribution.

View Table | View all tables in this article

In view of providing absolute terms of comparison, simulated sky-radiance data are presented in radiance units (i.e., W m^-2 nm^-1 sr^-1). It is however recognized that, in agreement with Mobley [12], relative units could be alternatively considered for the computation of $\rho $ with the sole requirement to respect the ratio between direct and diffuse irradiance components.

2.2.1 Harrison-Coombes model (HC)

The HC model parameterizes the diffuse component of the sky radiance through a non-linear least-squares regression of a set of 300 measurements according to the scheme proposed by Pokrowski [39] and Kittler [40]. These measurements are integrated radiance values in the 300–3000 nm spectral range. The HC wavelength-independent parametric formulation of the diffuse sky-radiance component $L_{diff}^{HC}({\theta ,\phi } )$ for any direction $({\theta ,\phi } )$ of the sky dome and sun zenith ${\theta _0}$ is given by

The determination of spectral values for the direct and diffuse sky-radiance components follows the scheme used in Hydrolight [21]. Specifically, GC is applied to compute: 1) the direct sun irradiance $E_{dir}^{CG}({\theta _0})$ necessary to determine the direct sun radiance $L_{dir}^{HC}({\theta _0})$ = $E_{dir}^{CG}({\theta _0})/{\Omega _0}$, where ${\Omega _0}$ is the solid angle of the sun-quad, and 2) the diffuse sun irradiance $E_{diff}^{CG}({\theta _0})$ necessary to compute the scaling factor ${S^{HC}}({\theta _0}) = {{E_{diff}^{CG}({\theta _0})} / {E_{diff}^{HC}}}({\theta _0})$, with $E_{diff}^{HC}({\theta _0}) = \int_0^{2\pi } {\int_0^{\frac{\pi }{2}} {L_{diff}^{HC}({\theta ,\phi } )\sin \theta \cos \theta d\theta } } d\phi$. The HC total sky radiance is then determined as $L_{sky}^{HC}({\theta ,\phi } )= L_{diff}^{HC}({\theta ,\phi } )\cdot {S^{HC}}({\theta _0})$ for any $({\theta ,\phi } )$ different from $({{\theta_0},0} )$, and as $L_{sky}^{HC}({{\theta_0},0} )= L_{diff}^{HC}({{\theta_0},0} )\cdot {S^{HC}}({\theta _0}) + L_{dir}^{HC}({{\theta_0}} )$ for $({{\theta_0},0} )$.

2.2.2 Zibordi-Voss model (ZV)

The ZV model for the computation of the diffuse sky radiance relies on the analytical solution of the radiative transfer equation (RTE) proposed by Sobolev for a plane-parallel, homogeneous and cloudless atmosphere, bounded by a Lambertian surface [41].

The expression of the diffuse sky radiance for any direction $({\theta ,\phi } )$ of the sky dome and sun zenith ${\theta _0}$ is given by

The term $\sigma(\tau_d(\lambda))$ is the scattering transmission of the atmosphere given by

The molecular scattering is described by the Rayleigh phase function, while the aerosol scattering phase function is approximated by the two-term Henyey-Greenstein (TTHG) analytical function [44] with asymmetry parameters ${g_1} = 0.713,\,{g_2} = 0.759$, and $a = 0.962$ proposed for maritime aerosols [42].

Recalling that ZV only provides the diffuse component of the spectral sky radiance, the determination of the spectral direct radiance component also benefits from the diffuse and direct irradiances computed with GC. Specifically, assuming equality of the GC and ZV total irradiances, the ZV direct sun irradiance component is determined as $E_{dir}^{ZV}({{\theta_0}} )= E_{dir}^{GC}({{\theta_0}} )+ E_{diff}^{GC}({{\theta_0}} )- E_{diff}^{ZV}({{\theta_0}} )$ where the explicit dependence on $\lambda $ has been omitted for simplicity, and $E_{diff}^{ZV}({\theta _0}) = \int_0^{2\pi } {\int_0^{\frac{\pi }{2}} {L_{diff}^{ZV}({\theta ,\phi } )\sin \theta \cos \theta d\theta } } d\phi$ is the ZV diffuse irradiance for sun zenith angle ${\theta _0}$ implicit in $L_{diff}^{ZV}({\theta ,\phi } )$. The direct radiance is then given by $L_{dir}^{ZV}({{\theta_0}} )$ = $E_{dir}^{ZV}({\theta _0})/{\Omega _0}$, with ${\Omega _0}$ the solid angle of the sun-quad. Computationally, the sky radiance from ZV is given by $L_{sky}^{ZV}({\theta ,\phi } )= L_{diff}^{ZV}({\theta ,\phi } )$ for any $({\theta ,\phi } )$ different from $({{\theta_0},0} )$, and by $L_{sky}^{ZV}({{\theta_0},0} )\,= \, L_{diff}^{ZV}({{\theta_0},0} )+ L_{dir}^{ZV}({{\theta_0}} )$ for $({{\theta_0},0} )$. It is finally noted that $L_{sky}^{ZV}({{\theta_0},0} )$ cannot be computed for ${\theta _0}$ and consequently it can only be approximated by values determined nearby the sun zenith angle. This singularity, however, only affects the determination of the diffuse sky radiance in the sun direction and does not represent a limitation to the present work. In fact, by adopting the same scheme utilized in Hydrolight, the sun is located at the center of the quad, while the diffuse sky radiance is computed for the direction $({\theta_i^\ast ,\phi_j^\ast } )$, as illustrated in Section 2.1.

2.2.3 FEM numerical code

The Finite Element Method radiative transfer code (FEM) [29], applies the finite element deterministic approach to model the propagation of unpolarized monochromatic radiation in the ultraviolet, visible and near-infrared electromagnetic regions within a vertically inhomogeneous atmosphere-ocean system. It fully accounts for multiple scattering, reflection and refraction at the air-water interface, and has been extensively tested and validated through benchmark with other popular codes [29–32].

The atmosphere-ocean system is divided into ${\rm {\frak L}}$ optically homogeneous horizontal layers. For each layer, the phase function is expanded in Legendre polynomials with a finite number of terms M, and the radiance L is correspondingly represented by a Fourier cosine series

The remaining integrals are further approximated by projecting each radiance harmonic ${L^m}\; $onto a set of basis functions, each different from zero only in a finite interval (hence, the finite element designation of this method). Specifically, the solution of Eq. (10) is sought in the form

The resulting system of matricial equations is solved analytically, imposing the necessary boundary conditions, and the solutions are lastly improved by applying an interpolation procedure.

It is highlighted that FEM ensures flux conservation independently of the number of grid points (N = 2 is already sufficient) and of the shape of the phase function (i.e., of the order of polynomials M used to correctly represent it), opposite to alternative numerical methods (e.g., the Discrete Ordinate Method), which require that the number of Gaussian quadrature points N used to solve the RTE integral is at least equal to the number of Legendre polynomials M. Additionally, by selecting basis functions which are finite, FEM generates a piecewise approximate solution that may converge faster than other global approximations (e.g., those adopted in the Spherical Harmonic Method). A number of grid points N = 64 is usually sufficient to accurately model the radiance angular distribution even when characterized by extreme gradients. The approach implemented in FEM is thus particularly suitable for highly asymmetric phase functions, like those characterizing the scattering by aerosols and hydrosols.

For the present work, the atmosphere is divided into 14 plane-parallel layers (i.e., ${\rm {\frak L}}$=14) resolving the vertical distribution of aerosol, ozone, and other gas molecules. Each atmospheric layer contains a variable mixture of ozone (only absorbing), other gas molecules (only scattering), and aerosols (scattering and absorbing). The vertical profile of ozone and Rayleigh molecules are determined in agreement with Guzzi et al. [45], and Marggraf and Griggs [46], respectively. The vertical profile of aerosols is modeled according to Elterman [47] assuming a mean aerosol scale height of 1.2 km. The molecular and aerosols scattering phase functions are defined as in ZV.

2.3 Simulated sky-radiances

Examples of diffuse sky-radiance distributions obtained with HC, ZV and FEM, are displayed in Fig. 3 for different ${\theta _0}$ (i.e., $20^\circ ,\,40^\circ $ and $60^\circ $). These radiance distributions show differences among the various models that become more pronounced with increasing ${\theta _0}$. This is confirmed by the diffuse sky-radiance values shown in Fig. 4 for the solar and measurement planes (solid and dashed lines, respectively), with HC exhibiting the largest deviations with respect to both ZV and FEM, and the latter being regarded as the most accurate approach. These differences are explained by the diverse capability of the applied models to account for forward and multiple scattering, more challenged in the circumsolar region and near the horizon, respectively.

 

Fig. 3. Sky-radiance distribution L_sky at 550 nm (restricted to the diffuse component) obtained with HC, ZV and FEM displayed from the top to the bottom row, respectively. Panels from left to right correspond to ${\theta _0}$ = 20°, 40° and 60°.

Download Full Size | PDF

 

Fig. 4. Sky radiance L_sky at 550 nm (restricted to the diffuse component) determined for the solar and measurement planes (solid and dashed lines, respectively) with HC, ZV and FEM for different ${\theta _0}$ (i.e., $20^\circ ,\,40^\circ $ and $60^\circ $). Dots highlight the zenith angles at which the sky-radiance values are computed (see Sec. 2.1). The direct radiance components from the sun-quad at $20^\circ ,\,40^\circ $ and $60^\circ $ (not displayed) are respectively 77.03, 37.32 and 21.16 $\textrm{W }{\textrm{m}^{ - 2}}\textrm{n}{\textrm{m}^{ - 1}}\textrm{s}{\textrm{r}^{ - 1}}$ for both HC and ZV. Corresponding FEM values are 77.13, 37.36 and 21.15 $\textrm{W }{\textrm{m}^{ - 2}}\textrm{n}{\textrm{m}^{ - 1}}\textrm{s}{\textrm{r}^{ - 1}}$.

Download Full Size | PDF

Among the three sky-radiance models applied, FEM ensures highly accurate spectral and angular distributions of the diffuse radiance under any condition. This includes an accurate determination of the sky radiance from both the circumsolar region (mostly affected by a correct treatment of the aerosol forward scattering) and from the horizon (mostly affected by a correct treatment of multiple scattering effects). This further implies accurate estimates of both diffuse and direct irradiances. Conversely, HC exhibits intrinsic limitations due to its wavelength-independent parametrization of the diffuse sky-radiance distribution through Eq. (5). Finally, ZV may provide realistic spectral sky-radiance distributions whose accuracy is expected to be higher when the multiple scattering contribution is lower (i.e., for low values of the aerosol optical depth).

2.4 Sea-surface models

The effect of the sea-surface statistics on $\rho $ replicability is analyzed accounting for: 1) the Cox-Munk sea-surface model (CM) [48], and 2) the Elfouhaily, Chapron, Katsaros and Vandemark sea-surface model (ECKV) [49]. Additional models (e.g., [50,51]) were not considered by restricting the analysis to statistical solutions either already adopted in the determination of the ρ tables currently utilized for the reduction of field data (CM), or characterized by a better capability to represent the sea-surface slope and elevation distributions (ECKV). The values of wind-speed ${v_w}$ considered in the present study vary from 0 to 10 $\textrm{m}{\textrm{s}^{ - 1}}$ with 2 $\textrm{m}{\textrm{s}^{ - 1}}$ increments.

2.4.1 Cox-Munk Sea-surface model (CM)

The CM equations express the slope distribution of surface waves as a Gaussian function with variance components along ($\mathrm{\sigma }_{\textrm{up}}^2$) and across ($\mathrm{\sigma }_{\textrm{cr}}^2$) wind direction given by

 

Fig. 5. CM2 (top panels) and ECKV (bottom panels) representations of the sea surface. The left and right columns refer to ${v_\textrm{w}}$ of 3 and 7 $\textrm{m}{\textrm{s}^{ - 1}}$, respectively.

Download Full Size | PDF

2.4.2 ECKV sea-surface model

The ECKV sea-surface model is based on the distribution of waves spectral energy. The elevation $z({i,j} )$ is expressed at the grid-points of the sea surface, discretized with ${N_x}$ and ${N_y}$ intervals along the x- and y-axis, respectively, under the assumption of linear wave theory, applying [15]

Processing steps for generating the sea surface require: 1) computing the two-sided discrete values of the elevation variance, 2) sampling the amplitudes of the harmonic components with normal distributions, and 3) defining the Inverse Fast Fourier Transform Hermitian coefficients to obtain the sea-surface elevation. The spectral density function ${\cal \textrm{S}}(k )$ is iteratively adjusted until converging to the target value of the waves slope distribution [17]. In the present work, the propagation direction of waves is assumed isotropic. Sea-surface representations generated with the ECKV model are illustrated in panels (c) and (d) of Fig. 5.

2.5 Computation of the sea-surface reflectance factor

The partition of the upper hemisphere into quads presented in Sec. 2.1 is here used to define the values of a radiance reflectance matrix ${\cal {R}}$ towards the sensor-quad (i.e., the quad at $\theta = 40^\circ$ and $\phi = 270^\circ $—see Fig. 2), so that

2.5.1 Hydrolight scheme

Hydrolight adopts forward ray tracing from the source-quad to the sea surface. Rays are initiated with unitary statistical weight. One quad at a time is selected as a source until the entire hemisphere is considered. The sea surface is represented by a hexagonal domain divided into triangular facets, called triads, tilted according to the Cox-Munk waves slope distribution (CM2, Eq. (16); see also Fig. 6). The Snell and Fresnel equations are applied to update the weight and direction of rays hitting the sea surface. To enhance the computing efficiency (an essential feature when the algorithm was conceived [53]), rays are forced to converge to the center of the hexagon as illustrated in Fig. 6(b). Rays initiated from a quad close to the horizon might be reflected by a triad located between the starting and the target points (the latter corresponding to the central triad) without reaching the hexagon center. This occlusion results in a phenomenon usually termed as “shadowing effect”, which nonetheless affects only a small number of rays.

 

Fig. 6. Illustration of the Hydrolight ray-tracing scheme. (a) Rays initialization from a random position of each emitting quad. (b) Perspective view, where dashed and solid lines illustrate incident rays and their reflection towards the sensor-quad, respectively. Rays reflected out of the sensor-quad are not shown. Each ray corresponds to a different hexagon realization (only one surface is however displayed).

Download Full Size | PDF

For each quad, the accumulated weight of the rays reflected toward the sensor-quad is divided by the number of emitted rays to determine the ratio $\Phi ({{i^\ast },{j^\ast }} )$ between reflected and incoming radiant fluxes. The process is repeated for several hexagon realizations by resampling the inclination of all triads to comprehensively represent the overall statistical properties of the sea surface. The computed average flux ratios $\bar{\mathrm{\Phi }}({{i^\ast },{j^\ast }} )$ are then transformed into radiance ratios (and hence elements of the reflectance matrix ${\cal {R}}$) by accounting for the geometry and solid angles of the emitting and receiving quads (for details, see [21]). Note that surface reflectance and transmission properties are inherent optical properties of the sea surface, parameterized as a function of the wind speed. Their values are pre-computed with the ray-tracing scheme formerly detailed and are not generated for each Hydrolight simulation.

2.5.2 MOX scheme

The MOX code relies on backward MC starting the rays from the sensor-quad, see Fig. 7(a). The sea surface is defined over a square domain of side $L = 50$ m discretized in $N \approx 500$ intervals along the x- and y-axis [17]. The spatial resolution and coordinates along the x-axis are ${\Delta _x} = L/{N_x}$ and $x(i )= i \cdot {\Delta _x}$, with $i = 0, \ldots ,{N_x}$. Equivalent quantities for the y-axis are ${\Delta _y} = L/{N_y}$ and $y(j )= j \cdot {\Delta _y}$, for $j = 0, \ldots ,{N_y}$. Each $({{\Delta _x},{\Delta _y}} )$ element is split into two triads whose vertices are elevated based on the CM1, CM2 or ECKV sea-surface model (Sec. 2.4).

 

Fig. 7. As in Fig. 6, but for the MOX backward ray-tracing scheme. The solid and dashed lines denote rays emitted from the sensor-quad and reflected by the sea surface, respectively. Rays that exit the domain from one side continue their path at the opposite side based on the periodicity of the sea surface.

Download Full Size | PDF

Rays are initiated at random horizontal positions from the same height above the sea surface—see Fig. 7(b). Those exiting the domain from one side continue their paths at the opposite side based on the periodicity imposed at the sea surface boundaries. The weight of the reflected ray is determined with the Fresnel equation. Upon the last reflection, the ray weight is accumulated at the hemisphere quads following the same scheme described for Hydrolight. Namely, the result of the weight accumulation, scaled by the number of emitted rays, represents the radiant flux ratio for the considered sea-surface model and wind speed. The MC simulations are repeated by resampling ten times the triads slope to compute the mean radiant flux ratio $\bar{\mathrm{\Phi }}({{i^\ast },{j^\ast }} )$ to better represent the sea-surface statistical properties (see Fig. 8). Considering that ray-tracing is performed in backward mode, ${\cal {R}}({{i^\ast },{j^\ast }} )= \bar{\mathrm{\Phi }}({{i^\ast },{j^\ast }} )$ as a result of the reciprocity principle [54].

 

Fig. 8. Slope variance as a function of the wind speed for each sea-surface model addressed in the study (i.e., CM1, CM2 and ECKV labeled with a circle, square and triangle, respectively). The solid lines connecting the symbols represent the expected theoretical values. The mean and the standard deviation of the slope variance resulting from multiple sea-surface realizations are displayed as dashed lines and shaded regions, respectively. Overlapping lines confirm the satisfactory agreement between expected and computed trends.

Download Full Size | PDF

The ray-tracing executed in this work to analyze the $\rho $-factor replicability considers 3 different sea-surface models and 6 wind-speed values. As anticipated, each sea surface, made of approximatively $5 \cdot {10^5}$ triads, is resampled 10 times. Individual MC simulations account for ${10^7}$ rays, and for each it is necessary to find the target triad along the ray trajectory, and then determine the intersection point. The overall possible number of ray-triad intersections is then approximately ${10^{15}}$. MOX simulations were carried out using the Navigator supercomputer of the University of Coimbra, Portugal. The system features 164 identical computing nodes equipped with 2 twelve-core Intel Xeon E5-2697v2 processors at 2.70 GHz, 96 GB RAM and InfiniBand interconnect. Multiple simulation cases were executed simultaneously using 24 CPU cores in parallel. This allowed running all simulation cases in about 15 minutes. Examples of radiance reflectance matrix obtained with the MOX code are presented in Fig. 9.

 

Fig. 9. Examples of radiance reflectance matrix based on the CM1, CM2 and the ECKV sea-surface statistics in the top, middle and bottom row panels, respectively. Column panels from left to right correspond to v_w = 2, 6 and 10 ms^-1, respectively. The solar plane crossing the hemisphere from $\phi = 0^\circ $ to $\phi = 180^\circ $and the radiometer orientation from $\phi = 270^\circ $ to $\phi = 90^\circ $are indicated by the dashed line and the arrow, respectively.

Download Full Size | PDF

3. Results

3.1 Replicability of the sea-surface reflectance factor

The mean and standard deviation of $\rho $ values are analyzed in Fig. 10 for different combinations of sky-radiance (HC, ZV, and FEM) and sea-surface models (CM1, CM2 and ECKV). The radiance reflected at the sea surface accounts cumulatively for the contribution originating from the sun-quad (S) and that from all other quads (D). The total $\rho $ value (where $\rho = {\rho _D} + {\rho _S}$) is then examined in terms of ${\rho _S}$ and ${\rho _D}$, respectively. The two components of $\rho $ are modeled through Eq. (18) and (19) accounting for the entry of the radiance reflectance matrix corresponding to the sun-quad.

 

Fig. 10. Variability of $\rho $ values as a function of wind speed, as resulting from the application of various models for sky-radiance distribution at 550 nm (first row of panels) and sea-surface statistics (middle row of panels), and the combined effects of sky-radiance and sea surface (bottom row of panels). Column panels from left to right report results for ${\theta _0}$ = 20°, 40° and 60°, respectively. The mean and the standard deviation of $\rho $ values are represented with a trend line and a shaded region, respectively.

Download Full Size | PDF

The capability of separating $\rho $ into its direct- and diffuse-related components, and additionally discussing the dependence of the two with wind speed, is of interest for above-water applications. In this respect, it is noted that the sun-quad has dimensions larger than the sun disc. As such, the glint contribution of the sky radiance originating from the sun-quad with respect to the case where the exact size of the sun disc is considered, translates into a larger number of rays characterized by a lower weight. In terms of operational above-water radiometry, the impact of this approximation is indeed lessened by the application of measurement geometries minimizing the sun-glint contribution. In fact, the contribution of ${\rho _S}$ is the most difficult to be statistically represented in field measurements and, for this reason, its perturbing effects are minimized through the adoption of suitable viewing geometries.

The HC and CM2 are the sky-radiance and sea-surface reference models applied in Hydrolight. These models are here considered as baselines to investigate the replicability of $\rho $ values. Namely, panels in the top row of Fig. 10 refer to the same sea-surface model (CM2) to detail the variability of $\rho $ induced by different sky-radiance models. Conversely, panels in the middle row rely on the same sky-radiance distribution (HC) to illustrate the effects due to different sea-surface representations. These results hence permit to document the specific sensitivity of the $\rho $ factor to the selection of the sky-radiance model (top-row panels) or of the sea-surface model (middle-row panels). Panels in the bottom row finally show the variability obtained when considering all sky-radiance and sea-surface models. Panels in the left, central and right columns refer to ${\theta _0} = 20^\circ ,40^\circ \textrm{ and }60^\circ$, respectively (results are presented adopting different scales to better visualize trends).

The solid lines and the shaded regions in Fig. 10 represent the mean and the standard deviation of different factors. The scrutiny of results indicates that the mean value of $\rho $ is driven by two distinct trends of ${\rho _D}$ and ${\rho _S}$ as a function of ${v_w}$. An increase approximately proportional to ${v_w}$ characterizes ${\rho _D}$ for the considered wind-speed interval. This trend displays minor variations among different sky-radiance and sea-surface models. Instead, at low ${\theta _0}$ and for ${v_w}$ above 4 ms^-1, ${\rho _\textrm{S}}$ grows rapidly with ${v_w}$ showing an exponential-like tendency as a result of increasing contributions to L_r of the reflected sun-quad radiance. This effect reduces with large ${\theta _0}$ leading to less pronounced sun glint perturbations.

As displayed in Fig. 10(a), both ${\rho _D}$ and ${\rho _S}$ show dependence on the selected sky-radiance model at high sun elevations (e.g., ${\theta _0} = 20^\circ $). With reference to ${\rho _S}$, its variability is mainly ascribed to the different sky-radiance values in the circumsolar region, where the HC ones are significantly larger than those from ZV and FEM (see Fig. 4). For the same sun elevations, Fig. 10(d) shows that ${\rho _\textrm{D}}$ is almost unaffected by the choice of the sea-surface model. The variability of ${\rho _\textrm{S}}$ appears instead similar to that documented in Fig. 10(a), although the mean value of ${\rho _S}$ is larger in Fig. 10(d) than in Fig. 10(a) for ${v_w}$ = 10 $\textrm{m}{\textrm{s}^{ - 1}}$. This confirms the relevance of the sea-surface slope variance in defining the fraction of direct light reflected into the sensor-quad.

At low sun elevations (e.g., ${\theta _0} = 60^\circ $), the variability of ${\rho _D}$ due to different sky-radiance models is larger than that due to different sea-surface models—c.f., Fig. 10(c) and Fig. 10(f). This results from the different capabilities of the various methods to represent the atmospheric multiple scattering, which exhibits more pronounced effects with sun zenith (see Fig. 4).

The above considerations on the different impacts of sky-radiance distributions and modeled sea surfaces are confirmed by the analysis of Fig. 10(g-i) showing their combined effects. Variations of $\rho $ are similar to those induced by alternative sea-surface models with low ${\theta _0}$ (c.f., Fig. 10(g) and Fig. 10(d)), and by different sky-radiance models with high ${\theta _0}$ (c.f., Fig. 10(i) and Fig. 10(c)).

Table 2 primarily aims at reporting the values of $\rho $ and providing additional elements to better understand peculiarities embedded in Fig. 10. Note that any application relying on the $\rho $ values included in this table should consider that simulations have been performed at a specific wavelength (i.e., 550 nm). Data in the top, middle and bottom data-block refer to the impact of different sky-radiance models, sea-surface models, and their combined effects, respectively. Row- and column-wise data contain statistical figures for the considered wind-speed ${v_w}$ and sun-zenith ${\theta _0}$ values, respectively. For selected ${v_w}$ and ${\theta _0}$, data sub-sets report the mean $\mu $, the standard deviation $\sigma $ and the coefficient of variation CV in percent ($100 \cdot \sigma /\mu $) of $\rho $, in the first, second and third entry from top to bottom.

Table 2. Statistical figures summarizing the variability of ρ at 550 nm. For each selected ${\theta _0}$ and ${v_w}$, data sub-sets report the mean $\mu $, the standard deviation $\sigma $ and the coefficient of variation CV in percent ($100 \cdot \sigma /\mu $) in the first, second and third entry from top to bottom.

View Table | View all tables in this article

Data indicate that CVs characterizing $\rho $ due to the application of different sky-radiance models are generally below 2% for ${v_w}$ up to 4 $\textrm{m}{\textrm{s}^{ - 1}}$. Conversely, CVs vary between 3% and 6% when ${v_w}$ reaches 10 $\textrm{m}{\textrm{s}^{ - 1}}$. At low wind speed, the sky-radiance contributions mainly originate from regions of the sky dome where the models mostly agree. An increased wind speed, instead, conveys sky-radiance contributions from a larger number of quads into the sensor field of view. These contributions may include those from the circumsolar and the horizon regions where the accuracy in the simulation of forward and multiple scattering, respectively, is more challenged (see Fig. 4).

When considering the diverse models applied to represent the sea surface, CVs of approximately $1.7$% characterize the values determined with ${v_w}$ = 0 $\textrm{m}{\textrm{s}^{ - 1}}$. This non-null value is explained by the fact that, in absence of wind, the CM2 and ECKV models assume a completely flat sea surface, while CM1 considers the sea surface to be always perturbed by some smooth wave, even in the absence of local wind (as in reality, see Fig. 8).

Figure 8 is also relevant to explain another detail that emerges from the analysis of Table 2. Still referring to the case of different sea-surface models, statistical figures indicate a general CV increase with ${v_w}$. However, an opposite tendency can be seen varying ${v_w}$ from 4 to 6 $\textrm{m}{\textrm{s}^{ - 1}}$. For instance, when ${\theta _0} = 10^\circ $, CV varies from 7.4% for ${v_w}$ = 4 $\textrm{m}{\textrm{s}^{ - 1}}$ to 6.5% for ${v_w}$ = 6 $\textrm{m}{\textrm{s}^{ - 1}}$. This can be explained with the help of Fig. 8. The variance of the sea-surface slope distribution foreseen by CM1 and CM2 linearly varies with wind speed (the former is higher, and a constant bias separates them). The variance predicted by the ECKV model, instead, is non-linear with wind speed. Specifically, the ECKV value is above that of CM1 with ${v_w}$ = 4 $\textrm{m}{\textrm{s}^{ - 1}}$, while it is in between those of CM1 and CM2 with ${v_w}$ = 6 $\textrm{m}{\textrm{s}^{ - 1}}$. The slope variance hence varies more in the former than in the latter case, which also affects the $\rho $ replicability.

Finally, considering the combined effects of the various sky-radiance and sea-surface models on the replicability of $\rho $, results further confirm that cases with higher wind speed, and sun elevation close to the zenith or to the horizon, are those where the sky-radiance and the sea-surface modeling is more influential.

3.2 Benchmark results

The values of $\rho $ determined with Hydrolight based on the HC sky-radiance and CM2 sea-surface models [12] are largely applied by the scientific community and currently recommended into the validation protocols [7] to derive ${L_w}$ from measurements of ${L_i}$ and ${L_t}$. A comparison between these values and equivalent ones obtained with MOX using the HC sky-radiance and the CM2 sea-surface models is presented in Fig. 11.

 

Fig. 11. Values of $\rho $ determined with (a) Hydrolight and (b) MOX adopting the same sky-radiance and sea-surface statistics, with solar zenith angle varying from $10^\circ $ to $80^\circ $. (c) Direct comparison between Hydrolight and MOX derived $\rho $. (d) Percent difference of $\rho $ values determined with MOX with respect to Hydrolight.

Download Full Size | PDF

Results show that the values of $\rho $ computed with MOX are always appreciably higher than those from Hydrolight. This difference also suggests an impact on the experimental determination of the water-leaving radiance through above-water radiometry, which relies on the application of $\rho $. The source of the observed difference is addressed in the next discussion section, benefitting of simplified simulation conditions and additional benchmark cases.

4. Discussion

4.1 Analysis of benchmark results and explanation of differences

In view of addressing differences between Hydrolight and MOX results, additional $\rho $ simulations were performed for an isotropic sky-radiance distribution, still modeling the sea-surface statistics with CM2 for measurement geometries defined by $\theta = 140^\circ $ and additionally $\theta = 180^\circ $. This simulation set permits to consider, as an additional and fully independent term of comparison, the $\rho $ values determined by Austin [10] for an isotropic sky and the Cox-Munk sea-surface statistics applying the ergodic-cap scheme proposed by Gordon [9]. The ergodic-cap approach is based on the construction of a three-dimensional surface with a slope distribution equivalent to that of the surface waves. The center and the boundary of this surface are associated with the lower and the larger waves slope, hence the term cap. Note that the cap is stretched to account for the slope probability distribution of waves, and a correction term is applied to resolve its projection in the viewing direction when computing $\rho $.

For notation compactness, the sea-surface reflectance factors are here indicated as $\rho _{{v_w}}^C(\theta )$, where: 1) the symbol C is H for Hydrolight, M for MOX, and A for Austin, 2) the values of the zenith viewing angle $\theta $ are $180^\circ $ and $140^\circ $, and 3) the wind-speed values ${v_w}$ are $4$ and 10 $\textrm{m}{\textrm{s}^{ - 1}}$. As an example, $\rho _4^H(180) = 0.0201$ denotes the sea-surface reflectance factor computed with Hydrolight assuming an isotropic radiance for ${v_w}$=4 $\textrm{m}{\textrm{s}^{ - 1}}$and nadir view.

Table 3 confirms that the values of $\rho $ computed with Hydrolight are lower than those obtained with MOX also in the presence of isotropic sky-radiance distribution, e.g., $\rho _4^H(140) = 0.0251$ and $\rho _4^M(140) = 0.0267$. The latter MOX values show a significant agreement with those from Austin. It is further highlighted that the Hydrolight results $\rho _4^H(180) = 0.0201$ and $\rho _{10}^H(180) = 0.0194$ are below the minimum expected values of $\rho = 0.0211$ for a flat Fresnel surface in nadir view. The latter finding, which does not appear to have a physical explanation, suggests investigating ray-tracing details.

Table 3. Comparison of ρ values obtained from Hydrolight, MOX and by Austin for an isotropic sky-radiance distribution and the CM2 sea-surface statistics.

View Table | View all tables in this article

This is undertaken with the help of Fig. 12, where: 1) the sky-radiance is assumed isotropic, 2) the sea surface is one-dimensional and for simplicity chosen to coincide with the measurement plane (see Fig. 2), 3) multiple reflections are not accounted for, and finally 4) waves occluding effects are neglected (the last assumption is supported by the viewing direction of the virtual radiometer [9,11]).

 

Fig. 12. Schematic of ray tracing in the measurement plane.

Download Full Size | PDF

The discretized sea surface shown in Fig. 12 is composed by a set of segments with coordinates $({x_j},{z_j})$ with $j = 0, \ldots ,N$. These segments are determined by sampling angles ${\alpha _j}$ from the waves slope distribution, and by defining ${z_0} = 0$ and ${z_j} = {z_{j - 1}} + \Delta {z_j}$ with a fixed interval $\Delta x$ and elevation $\Delta {z_j} = \Delta x\tan {\alpha _j}$. The length of the j-th segment is ${A_j} = \Delta x\sec {\alpha _j}$. The projection of the j-th segment onto the plane perpendicular to the observation direction ${\mathbf v}$, which is the segment seen by the virtual radiometer, has length ${w_j}$

With $W = \sum\limits_{j = 1}^N {{w_j}} $ and recalling the assumption of isotropic sky-radiance, the value of ${L_r}$ is computed as

The ray initialization performed in MOX from a random point of the domain top permits to statistically consider the relative weight of each triad according to its projection in the ray direction (see Fig. 7). In fact, the “visible” area (i.e., projection) of each facet in the sensor direction varies according to the facet orientation. This feature is also considered by Austin with the ergodic-cap approach utilizing a projection correction, which supports the agreement between $\rho $ determined with these simulation schemes. Conversely, the Hydrolight scheme implies that about the same number of rays (i.e., excluding the minority that is affected by the shadowing effect) hits the central triad for repeated hexagon realizations, as illustrated in Fig. 6. When computing the mean of the radiant flux ratio (see Sec. 2.5.1), the former assumption hampers the capability to account for the relative weight of the reflecting triad at the center of each hexagon realization due to its projection towards the sensor-quad. In other words, directing the rays towards the central triad does not allow to properly account that the “visible” area of the reflecting facet varies among sea surface realizations, with clear implications on the computation of the average contribution.

Additional MC tests have confirmed the above findings by showing that: 1) backward and forward Monte Carlo simulations converge to the values of $\rho $ that are close to those computed with MOX, and 2) the forward scheme performed with the simple arithmetic mean instead of a weighted mean as in Eq. (22), leads to smaller values of $\rho $. In this case, the underestimate is analogous to that observed for Hydrolight—e.g., $\rho = 0.0250$ vs. $\rho _H^4(140) = 0.0251$ and $\rho = 0.0192$ vs. $\rho _H^{10}(180) = 0.0194$. Differences in the ray-tracing schemes are hence identified as the key factor explaining the lower values of $\rho $ computed by Hydrolight with respect to MOX.

4.2 Relevance of the statistical representativity of in situ measurements

The statistical representativity of ${L_t}$ measurements is a fundamental aspect in above-water radiometry. Remarking that ${L_t} = {L_r} + \textrm{ }{L_w}$, and recalling that ${L_r}$ applied for the determination of $\rho $ is quantified numerically, the actual ${L_r}$ characterizing the in situ measurements may exhibit differences from the theoretical value of ${L_r}$. The working hypothesis that underlies the theoretical determination of $\rho $, which should be also featured by in situ measurements, is that the reflective properties of the sea surface are averaged over time, space and propagation direction of waves. Individual ${L_t}$ measurements may instead vary over time and as a function of the sensor features (e.g., field-of-view and sampling time). Consequently, each measurement can only capture part of the variability induced by a wavy surface over time and space.

This aspect is illustrated through Fig. 13. Dots in the right side of the sky hemisphere exhibiting HC sky-radiance distribution with ${\theta _0} = 20^\circ $, represent ten independent rays trajectories randomly sampled in backward mode from a virtual sensor without any restriction applied to sky- and sun-glint contributions. Their dispersion is explained by the slope distribution of waves here generated with CM2 for ${v_w}$ = 4 $\textrm{m}{\textrm{s}^{ - 1}}$. Assuming that each trajectory corresponds to a different measurement performed by the sensor with infinitesimal integration time and field-of-view, the mean ${L_r}$ values of the two simulated data records are 0.0020 and 0.0022 $\textrm{W}\; {\textrm{m}^{ - 2}}\; \textrm{n}{\textrm{m}^{ - 1}}\; \textrm{s}{\textrm{r}^{ - 1}}$, thus indicating a $10$% difference. The former result firmly suggests the need to rely on the averaging of a statistically significant number of ${L_t}$ measurements to enhance consistency between simulated and actual ${L_r}$.

 

Fig. 13. Example of repeated determinations of ${L_r}$ through a limited set of virtual measurements. The sky radiance and the sea surface are represented by the HC (${\theta _0} = 20^\circ$) and the CM2 (${v_w} = 4\textrm{m}{\textrm{s}^{ - 1}}$) models, respectively. The ${L_r}$ values corresponding to panel (a) and (b) are 0.0020 and 0.0022 $\textrm{W }{\textrm{m}^{ - 2}}\textrm{n}{\textrm{m}^{ - 1}}\textrm{s}{\textrm{r}^{ - 1}}$, respectively.

Download Full Size | PDF

Figure 13 illustrates an idealized case that only accounts for the variability of ${L_r}$ due to the slope variance of surface waves. Differences in repeated ${L_t}$ measurements, however, also depend on the sky-radiance distribution (that can vary as a function of several parameters including wavelength). Still, acknowledging the impossibility to account for any specific parameter that may impact the spatio-temporal averaging of surface effects in ${L_t}$ such as the sensor sampling time and field-of-view or its height above the sea level, the former considerations indicate the need for $\rho $ factors relying on those sky-radiance (also accounting for spectral dependence) and sea-surface models exhibiting the highest accuracy.

4.3 Further considerations on operational above-water measurements

In addition to the statistical representativity of ${L_t}$ measurements addressed in the former section, uncertainties and biases affecting $\rho $ unavoidably impact the determination of ${L_w}$. Such an impact can be quantified in relative terms following Gergely and Zibordi [56], but excluding the effects of the uncertainties affecting ${L_t}$ and ${L_i}$, as

Specifically, by assuming 1) ${v_w}$ = 4 $\textrm{m}{\textrm{s}^{ - 1}}$, 2) ${\theta _0}$ larger than $20^\circ $, and 3) setting $u(\rho )$ equal to the standard deviation of $\rho $ given in Table 2, the impact on ${L_w}$ (i.e., $100\textrm{ } \cdot \textrm{ }u({L_w})/{L_w}$) at 550 nm due to the adoption of diverse sky-radiance and sea-surface models in the determination of ρ, is lower than 0.5% for most water types and clear sky conditions.

More importantly, in agreement with Eq. (1), any underestimate (overestimate) of $\rho $ leads to an overestimate (underestimate) of ${L_w}$. Regarding the data shown in Fig. 11, a negative bias of approximately 0.002 affecting $\rho $ at 550 nm (comparable to that shown by the values of $\rho $ determined with MOX and Hydrolight for ${v_w} = \textrm{ }4\,\textrm{m}{\textrm{s}^{ - 1}}$ and ${\theta _0}$ larger than $20^\circ $), induces an overestimate in ${L_w}$ approaching 1.5%, (see Eq.s (1) and (3)).

The above considerations have relevance for the current processing solutions applied to AERONET-OC data [57,20], which rely on a modified version of Eq. (1) and the application of $\rho $ determined using Hydrolight [12]. Anticipating that the AERONET-OC data largely refer to sea-state conditions driven by low-to-mild wind speeds generally not exceeding ${v_w} = 5$ $\textrm{m}{\textrm{s}^{ - 1}}$ as a result of strict quality assurance and quality control schemes, the value of ${L_t}$ applied for the determination of ${L_w}$ in the AERONET-OC processing is the mean of a percentage of the individual sea-radiance data exhibiting the lowest values in a measurement sequence. This solution, only supported by experimental evidence provided by the comparison of ${L_w}$ data determined from above-water and in-water methods [7,58], was tentatively justified by the need to exclude those sea-viewing measurements exhibiting excessively large glint contamination and expected to lead to an overestimate of ${L_t}$. The $\rho $ values displayed in Fig. 11 give novel elements to explain the better performance of Eq. (1) when relying on the averaging of relative minima from a series of sea-viewing measurements opposite to the averaging of all the measurements. Given the underestimate of the $\rho $ values applied in data processing, the related impact on the determination of ${L_w}$ is heuristically compensated by an underestimate of ${L_t}$ as resulting from the averaging of a fraction of the sea-viewing data exhibiting the lowest values in each measurement sequence.

5. Summary and remarks

Results from the study, restricted to a consolidated measurement geometry [7], indicate that the variability of $\rho $ derived from the application of alternative methods to determine the sky-radiance distribution and the sea-surface statistics is generally below 2% when considering wind speeds up to 4 ms⁻¹ and sun zenith angles larger than approximately $20^\circ $. The variability of $\rho $ increases with wind speed, especially when the sun is close to the zenith or the horizon. This dependence has been interpreted acknowledging that, when the variance of the sea-surface slope distribution increases due to wind speed, the sky region that contributes to the radiance reflected towards the sensor-quad widens to include areas where sky-radiance distributions can show significant differences (e.g., the circumsolar and horizon regions) both in intensity and spectrally.

The variability of $\rho $ resulting from the replicability exercise can exceed 5% for wind speeds higher than 4 $\textrm{m}{\textrm{s}^{ - 1}}$ and sun elevations close to the zenith (${\theta _0} = 20^\circ $ or less). This result indicates that at low sun zenith angles, the determination of ${L_w}$ through above-water radiometric methods is challenged by the capability to properly represent the relative contribution of the direct and diffuse sky-radiance in the computation of $\rho $.

The ρ replicability has been evaluated at 550 nm to specifically account for reference values as benchmark [12], whilst analyses at different wavelengths are beyond the scope of this work. The study has hence shown that the values of $\rho $ computed with Hydrolight are lower than those obtained with MOX regardless of the adoption of the same sky-radiance and sea-surface models. The reason for this difference has been identified in the ray-tracing schemes employed in the two radiative transfer codes. The solution applied in Hydrolight to improve simulation efficiency by fixing the point towards which rays are emitted, presents some caveats when combining results from different sea-surface realizations. Specifically, the relative weight of the reflecting sea-surface facets due to their projections towards the sensor-quad is not taken into consideration in the determination of $\rho $. Conversely, the MOX code implicitly accounts for this effect in statistical terms as a result of the randomness of the starting point from the top of the simulation domain for rays traced in backward mode. This allows improving the accuracy of the computed $\rho $ although additional processing resources are demanded.

The study also offered discussion elements on the operational determination of the water-leaving radiance with above-water systems. One aspect is that the uncertainties affecting ${L_w}$ naturally depend on the actual statistical representativity of the replicated sea-viewing measurements applied for the determination of ${L_t}$, besides the accuracy of $\rho $.

Results from the study also support the rationale for the AERONET-OC processing applying ${L_t}$ values determined from the averaging of relative minima, in alternative to the averaging, of a series of sea-viewing measurements. This solution is shown contributing to compensate for the impact of the underestimated $\rho $ values from Hydrolight simulations.