Determination of the remote-sensing reflectance from above-water measurements with the "3C model": a further assessment.

The Three-Component Reflectance Model (3C) was primarily developed to improve the determination of the remote-sensing reflectance (Rrs) from above-water radiometric hyperspectral measurements performed during sub-optimal conditions (i.e., cloudy sky, variable viewing geometry, high glint perturbations, low illumination conditions). In view of further validating the model and showing its broad range of uses, this work presents the application of 3C to above-water radiometry data collected in oceanic and coastal waters with a variety of measurement conditions. Rrs derived from measurements performed during optimal and slightly sub-optimal conditions exhibit equivalence with Rrs obtained with an established above-water method that is commonly used to support ocean color validation activities. Additionally, the study shows that 3C can still provide relevant Rrs retrievals from field data characterized by low-light illumination, high glint perturbations and variable measurement geometries, for which the established method cannot be confidently applied. Finally, it is shown that the optimization residual returned by the 3C full-spectrum inversion procedure can be a potential relative indicator to assess the quality of derived Rrs.


Introduction
In-situ remote-sensing reflectance (R rs ) data from a variety of marine regions are required for the validation of satellite radiometric products. These R rs data can be confidently determined from above-water measurements performed during ideal (i.e., optimal) conditions, which require low-to-slight sea state, clear sky, strict observation and illumination geometries to minimize glint perturbations, and additionally that the footprint of the sea-viewing radiance sensor is well away from the superstructure [1].
Sub-optimal conditions are those that do not comply with the latter requirements. Still, measurements performed during sub-optimal conditions have relevance for applications such as bio-optical modelling and water quality monitoring, which exhibit less strict uncertainty requirements when compared to the validation of satellite data products. Because of this, the need for exploiting above-water radiometric data collected during sub-optimal conditions has recently got impetus [2]. R rs determinations from above-water measurements are commonly made from simultaneous spectral measurements of i. the total radiance from the sea L T , ii. sky radiance L i and iii. total downward irradiance E s . L T and L i are measured with the same relative azimuth φ from the sun, but supplementary zenith angles (i.e., the viewing angle θ v for L T and the angle 180-θ v for L i ). By making use of these parameters, the directional-to-hemispherical reflectance of the total radiance from the sea L T /E s can be modelled as the sum of R rs and a glint contribution R G : The quantification of R G is the main challenge of above-water radiometry. In the case of a flat sea surface, R G is determined by the sky radiance reflectance L i /E s multiplied by the Fresnel reflectance factor ρ f at the viewing angle θ v : However, natural water surfaces are not flat. This implies that R G can vary in a complex way with the sky radiance distribution (i.e., not only with L i from a given direction) as a function of wave heights and slopes [3]. Thus, ρ f L i may not accurately represent the glint contribution to L T for any measurement condition.
Mobley [4] (M99 hereafter) theoretically showed that for known wave slopes and illumination and observation geometries, Eq. (2) can still be applied with ρ f replaced by a coefficient ρ that, when excluding cloud perturbations, is expressed as a function of the wind speed and the illumination and observation geometries. Markedly, M99 implies the use of measurement geometries minimizing the effects of sun-glint and non-uniform sky radiance [4]. Recent investigations showed that ρ also depends on the aerosol load and type [5], which should then be measured as well for the glint correction. Nevertheless, during ideal measurement conditions, M99 was demonstrated to deliver R rs with uncertainties meeting requirements for satellite ocean color validation [6]. Because of this, despite of decades of investigations on alternative radiometric solutions for the quantification of R rs , M99 remains the basic method for above-water radiometry when relying on ideal measurement conditions [1].
In the presence of clouds, the determination of R rs faces increased difficulties. For instance, while the direct sunlight is unpolarized, the air molecules and clouds may induce partial polarization to sky radiance that would undermine the assumption of randomly polarized light [7]. In such sub-optimal conditions, M99 could lead to an incorrect estimate of glint contributions. This could be adjusted by determining and removing the residual glint. Over clear to moderately turbid waters (tentatively defined as those waters presenting concentrations of total suspended matter lower than 5 g m −3 , and chlorophyll-a lower than 5 mg m −3 ) this residual or offset is commonly estimated from R rs in the near-infrared (NIR) by assuming R rs (NIR)≈0. Conversely, in the case of highly turbid waters, where R rs (NIR) > 0, specific correction schemes should be applied [e.g., 8]. Both, the basic scheme suitable for clear or moderately turbid waters, and others for very turbid waters, provide wavelength-independent estimates of such an offset. However, irrespective of the water type, the offset affecting R rs determinations is unlikely to be spectrally constant and requires a spectrally-resolved correction [9].
In view of overcoming the above limitations and to support the exploitation of above-water hyperspectral radiometric measurements performed during sub-optimal conditions, the Three-Component Reflectance Model (3C) accounts for a first order term representing the Fresnel reflection of L i and an additional term accounting for the diffuse and direct downward irradiance contributions [2].
Regardless of the promising performance of 3C, both example applications and comparison to independent data were limited in the original publication [2]. Therefore, this work encompasses some model revisions and aims at further assessing 3C applied to field data from a wider range of marine environments and measurement conditions. Specifically, R rs determined with 3C and M99 from above-water measurements are assessed against R rs derived from in-water radiometric profile data for both ideal and slightly sub-optimal illumination conditions (where "slightly sub-optimal" here implies ideal conditions except for the presence of clouds, still not affecting either the direct sun beam or the field of view of the L i sensor). Additionally, the performance of 3C is investigated applying measurements heavily affected by a non-favorable observation geometry, low illumination, and high sea state. Finally, the sensitivity of 3C to glint perturbations is analyzed using sequences of measurements extending over several minutes.

Glint characterization
The first step in the determination of R rs through 3C is the modelling of R G . The basic glint term ρ f L i E s corresponds to the glint contribution from a perfectly flat sea surface exhibiting Fresnel reflectance ρ f . This term is then expanded by adding a new term ∆ accounting for any glint contribution not quantifiable through L i : It is emphasized that ∆ in Eq. (3) is a spectrally dependent parameter. It is modeled as a function of the direct and diffuse ratios to the total downward irradiance, E sd E s and E ss E s , respectively: where E sd E s and E ss E s are scaled by their surface reflectance factors, ρ sd and ρ ss , which vary with the sun zenith angle. The factors f sd and f ss indicate the a-priori unknown fractions of the reflected direct and diffuse irradiances contributing to L T , which depend on the sky illumination and sea surface roughness. The term δ aims at facilitating retrievals with measurements affected by broken clouds or by high sea states causing bubbles and whitecaps. Specifically, it could help to model the spectral shape of glint contributions for overcast and scattered cloud conditions [10].
The irradiance ratios E sd E s and E ss E s are semi-analytically modeled according to Gregg and Carder [11]. While summary descriptions were formally provided [2,12], here it is only recalled that the model was developed for maritime clear skies and requires a number of input parameters such as the sun zenith angle, atmospheric pressure, relative humidity and the aerosol optical thickness τ a . This latter is expressed by the power law τ a = β(λ a ) −α , where α is the Ångström exponent and β the turbidity coefficient at the reference wavelength λ a = 550 nm. It is mentioned that the implementation of 3C only requires the diffuse and direct to total irradiance ratios, and not any absolute irradiance value. This strongly simplifies the model requirements because the input data related to extraterrestrial irradiance, aerosol, ozone, oxygen and water vapor transmissions cancel out.
Assuming a randomly polarized light, ρ sd is set equal to the Fresnel reflectance of a collimated and unpolarized sun ray with zenith angle θ s : The reflectance factor ρ ss for the diffuse irradiance, results from complex interactions between the sea surface and the light originating from any direction in the sky dome. Gege [12], relying on radiative transfer simulations for a clear sky, proposed its parametrization through a second-order polynomial function of 1 − cosθ s : sky light reflectance during clear sky conditions. Still, it was applied to measurements affected by haze, scattered clouds and overcast sky conditions [2]. Contrary to the original formulation [2], the separation of ρ sd and ρ ss , and of their unknown fractions f sd and f ss , reflects the deterministic nature of the former while the latter accounts for deviations from the analytically-derived model factors.
It is highlighted that the term ∆ includes the direct and diffuse fractions. Consequently, it can account for any glint contribution, which may make the term ρ f is a good first guess for the glint contribution. Still, when L i is poorly representing the glint affecting L T (e.g., with very high sun glint), it is possible to almost exclusively rely on the term ∆.

R rs modeling
The in-water part of 3C relies on a semi-analytical bio-optical model that yields R rs as a function of the optical properties of the significant water constituents and a number of boundary conditions. In view of relying on a comprehensive R rs model representation, a complex Case 2 water approach with several free parameters is applied to any water type (including Case 1 waters). This solution may lead to overfitting, which does not undermine the 3C performance and still ensures a spectral match regardless of the water type. In fact, the primary purpose of 3C is to separate the water and sky radiance contributions to L T . Consequently, any uncertainty affecting the concentration of the retrieved aquatic constituents is here of limited concern.
The semi-analytical bio-optical model by Albert and Mobley [13] (AM03 hereafter) is applied in 3C. Briefly, AM03 models R rs as a fourth-degree polynomial of b b /(a + b b ) multiplied by terms function of the sun and viewing geometry. Sole changes applied to AM03 are the parametrization of the absorption coefficient a g by colored dissolved organic matter (CDOM) and of the backscattering coefficient by particles b bp . In particular, a hyperbolic function a g = a g0 λ λ 0 −n g with λ 0 = 440 nm is applied for CDOM absorption due to its superior capability of fitting a g in the ultraviolet when compared to common exponential functions [14]. This solution minimizes potential cross-compensations involving CDOM absorption and Rayleigh sky glint whose spectral impact on L T are similar, and consequently allows for a more accurate separation of the two.
The particle backscattering coefficient b bp is also modeled by a power law function b bp = proportional to the concentration of the total suspended matter, C TSM , with λ 0 = 500 nm, hyperbolic slope η, and b * bp = 0.0042 m 2 g −1 [2,12]. The absorption coefficient by phytoplankton a ph = C a a * ph is determined by an average spectrum a * ph specific of Lake Constance with C a indicating the chlorophyll concentration [15].
The model choice is intended to embrace a wide range of water types. Yet, bio-optical parameterizations can readily be regionalized in view of further improving spectral matches.

Model implementation and inversion
Following Eqs. (2)-(3) and the modelling of the atmospheric and aquatic components described in the previous sections, a theoretical estimate of L T /E s is obtained from: where the subscript "meas" is introduced to specify the measured term, opposite to the modelled ones identified by "mod," depending on a number of free parameters. Several ancillary parameters, such as the aerosol model (set to 4), relative humidity (set to 80%) and atmospheric pressure (set to 1013.25 mbar), are assumed constant due to their negligible impact on the direct and diffuse to total irradiance ratios exploited in 3C.
The goal function to be minimized (also known as objective, fitness, residual, or cost function), is equal to the squared Euclidean norm of the difference between modelled and measured L T /E s , weighted by a spectrally dependent function W: with W in units of sr. The process returns a minimized value of and a vector x opt , containing the optimized values of the various free parameters: x opt = (C a , C TSM , η, a g0 , n g , α, β, f sd , f ss )| opt (9) The optimized values of the atmospheric parameters (α, β, f sd , f ss )| opt are then applied to calculate the optimized value of the glint offset ∆| opt , added to the Fresnel reflection term ρ f L i E s meas and subtracted from L T E s meas , to determine the remote-sensing reflectance R rs from 3C: The goal function , which indicates the goodness of the model-to-data L T /E s fit, appears a relevant quality indicator for R rs retrievals. This is supported by the fact that the spectral shapes of R rs and R G are markedly different. Consequently, an accurate L T /E s fit generally suggests the accurate separation of the R rs and R G components.
Given the generic bio-optical model used in 3C, the derived aquatic parameters (C a , C TSM , η, a g0 , n g ) have the sole role of reconstructing a realistic R rs . As already stated, their values might be affected by overfitting and should therefore be considered with caution. Nevertheless, if their accurate determination is of concern, regionally-tuned parameters for the specific inherent optical properties should be applied, still recognizing that in most cases a generic parameterization is sufficient to reach a confident separation between the atmospheric and the aquatic contributions.
The 3C model applied in this work was implemented in MATLAB (The Mathworks, Inc.) using genetic algorithms (GA) for the minimization process [16]. Trial and error analysis led to the definition of the main tuning parameters of the GA algorithm: i. the crossover rate applied to generate new solutions from the existing ones by gene combination and ii. mutation controlling the rate of random gene changes during iterations expressed as a function of the shrink parameter. This heuristic analysis showed that a crossover fraction of 0.4 with respect to the default value of 0.8 and the reduction of the shrink parameter to 0.6 from the default value 1, lead to both satisfactory convergence and algorithm efficiency.

Assessment of 3C performance for ideal and slightly sub-optimal measurement conditions
The performance of 3C was evaluated against an alternative in-water radiometry method by comparing R rs determined from measurements performed during ideal or slightly sub-optimal conditions, relevant for the validation of satellite radiometric products. Specifically, measurements were collected during oceanographic campaigns ( Fig. 1) carried out in the Central Mediterranean Sea during Summer 2017, in the West Mediterranean Sea in spring 2014, and at the Acqua Alta Oceanographic Tower (AOOT) in the northern Adriatic Sea during June 2014. While the ship data only refer to ideal measurement conditions, the AAOT data are characterized by sub-optimal conditions. Above-and in-water measurements were generally performed within 15 min with b) Remote-sensing reflectances R rs derived from above-water measurements applying the 3C model. The low sensor sensitivity in the violet and the low signal in the near infrared, make R rs noisier in these extreme spectral regions. solar zenith angles between 20°and 60°. Wind speed was always lower than 5 m s −1 and sea state was low-to-slight, always below 4 of the Beaufort scale. L T and L i radiances (with viewing angle θ v of 40°and relative azimuth φ of 90°) were simultaneously collected using two RAMSES-ARC (TriOS Mess-und Datentechnik GmbH) hyperspectral radiometers. The use of φ = 90°, opposite to φ = 135°, is suggested by the need to minimize the impact of superstructure perturbations, which increase with the sun zenith angle [1]. The above-water downward irradiance E s was collected with a RAMSES-ACC radiometer operated above any superstructure element to avoid shading perturbations. In agreement with current protocols, data derived from sequences of successive L T , L i and E s measurements (i.e., 18 in this case) were filtered for excessive variability [1,17]. This was achieved by removing those measurement sequences i. characterized by tilts exceeding 2°for the E s sensor or ii. exhibiting variation coefficient exceeding 4% for L T , or 2% for L i or E S , in the 450-650 nm spectral interval (i.e., the region commonly exhibiting the highest signal-to-noise ratio for L T measurements). The former thresholds were intended to eliminate those measurement sequences affected by high glint (through L T filtering), significant sky-radiance inhomogeneity (through L i filtering) and ship motion or cloud contamination (through E s filtering). This filtering, integral part of the quality assurance process, naturally leads to the elimination of data collected with moderate to high wind speed (i.e., tentatively exceeding 5 ms −1 [6,18]). For each measurement sequence, L w was then computed from L T and L i data normalized with respect to E s by selecting, averaging and applying the three L T measurements exhibiting the lowest mean values in the 450-650 nm spectral range. It is emphasized that the application of the mean of relative minima of L T values is just an additional processing element aiming at further reducing the impact of sea surface perturbations in L w computed from L T already filtered for excessive variability and consequently not significantly affected by sun-glint perturbations [1,18].
In-water data were collected using a µPRO free-fall profiler equipped with OCR-507 multispectral radiometers (Satlantic Inc., Halifax) simultaneously measuring the spectral upwelling radiance Lu and the downward irradiance E d as a function of depth, and the above-water downward irradiance E s . After removing any in-water data affected by tilts exceeding 5°and E s with tilt higher than 2°, R rs was determined in agreement with current protocols [1,19].
Both, above-water M99-derived and in-water R rs were corrected for the anisotropy of the in-water light field through the scheme proposed for Case 1 waters [1,20].
When considering M99, the major sources of uncertainty affecting the retrieval of R rs from RAMSES are: i. the uncertainty of sensor calibration; ii. the uncertainty due to stray light effects; iii. the maximum sensitivity to polarization for L T and L i sensors; iv. the uncertainty in E s ; v. the uncertainty in sky-light correction as a function of the uncertainty in the surface reflectance factor ρ; and vi. the uncertainty in the determination of correction factors for minimizing the effects of anisotropy of the in-water water light field (i.e., the bi-directional effects as a function of the sensors viewing angle and sun zenith). Former investigations [21] indicated that the combined values of these uncertainties vary between 5 and 6% in the spectral range of interest.
Major contributions to combined uncertainties affecting R rs derived from in-water multispectral radiometric data are: i. the uncertainty in the absolute in-air calibration and immersion factors; ii. the uncertainty in the correction factors applied for removing self-shading perturbations; iii. the uncertainty in the determination of correction factors for minimizing the effects of anisotropy of the in-water light field (i.e., the bi-directional effects as a function of sun zenith); iv. the uncertainty in E s ; v. the uncertainty in the extrapolation of sub-surface data; and vi. the environmental variability resulting from the combination of wave-induced perturbations with seawater and illumination changes during profiling. A tentative estimate of these uncertainties leads to combined values approaching 5% across the blue-green-red spectral regions [22].
The quantification of the uncertainties associated to the determination of R rs with 3C has not been addressed to date and exceeds the aims of this work. Such a task is likely challenged by the complex combination of deterministic calculations, semi-analytical modelling and optimizations leading to R rs retrievals through 3C. Still, the following analyses will rely on the application of tests based on 3C model fit residuals as a first attempt to address the quality of derived R rs .
The 3C-derived spectra from the measurements performed in the Central and Western Mediterranean Seas show features typical of Case 1 waters, with the first indicating oligotrophic conditions and the second more representative of mesotrophic waters [ Fig. 1(b)]. Conversely, the spectra from the AAOT exhibit features more typical of Case 2 waters significantly affected by the presence of sediments and CDOM. For the comparison exercise, also 3C-derived data were corrected for bidirectional effects through the scheme proposed for Case 1 waters [20], still recognizing the limits of this approach when applied to optically complex waters such as those characterizing the AAOT site.
Following section 2.2, the generic Case 2 water model AM03 was applied to theoretically determine the R rs spectra. The parameter search was customized as detailed in Table 1. Retrievals showed improvements when setting a high spectral weight W for the goal function in the NIR and violet spectral regions. To optimize the convergence, the chlorophyll-a fluorescence peak in the red and the strong dip due to oxygen absorption in the NIR were excluded from the optimization process by setting the relevant spectral weights to zero. Also, wavelengths lower than 350 nm or higher than 920 nm were excluded as being heavily affected by noise. Overall, W was determined as: 3C optimizations were performed with 200 generations and a population of 100 elements. Using a laptop with an Intel Core i7-770HQ CPU at 2.8 GHz and 16 GB RAM, the run time was ∼12 s for each spectrum resolved with a 3.3 nm spectral sampling.
The computed M99 R rs were posteriorly corrected for any residual offset in the NIR by subtracting R rs (NIR), defined here as the minimum R rs in the 775-900 nm spectral interval (these products are hereafter indicated as M99 ). Such background correction is fully justified for the marine waters matter of this analysis, all characterized by low to moderate sediment concentrations. Figure 2 shows details on 3C-derived R rs for selected stations in the Central and Western Mediterranean Seas, and at the AAOT. Notably, the measured L t /E s spectra closely match the modelled ones. R G is the sum of the Fresnel reflectance ρ f L i /E s and ∆, with ρ f L i /E s dominant contribution and also first guess for the determination of R G itself. For ideal measurement conditions (Central and Western Mediterranean Seas), the Fresnel term nearly equals R G and the term ∆ is only a small correction. However, at the AAOT, where the measurement conditions were slightly sub-optimal due to the presence of clouds, the Fresnel term underestimates the glint contribution. In this case, ∆ significantly contributes to R G . In all cases, the M99 glint term (i.e., ρ M99 L i /E s + R rs (NIR)), which includes the bias adjustment R rs (NIR) required to force to zero R rs in the NIR, is plotted for comparison. 3C returns a glint spectrum R G very similar to the M99 glint term, but without the need to impose the condition of zero R rs in the NIR. The comparison between above-water and in-water-derived R rs is presented by grouping data according to the measurement conditions: ideal for the Central and Western Mediterranean Sea data (N=12), and slightly sub-optimal for the AAOT data (N=12). The R rs data from above-water measurements were re-sampled to match the spectral bands of the in-water multispectral sensors. Both TriOS and Satlantic sensors exhibit a close spectral resolution of 10 nm. Considering the 3 nm spectral sampling of TriOS sensors, the difference in center-wavelength characterizing the products from the different sensors is generally lower than 3 nm. Results in Fig. 3 show excellent agreement of 3C-and M99 -with µPRO-derived R rs for ideal measurement conditions, but also some degree of mismatch for the data collected during slightly sub-optimal conditions. Statistical comparisons are quantified through the average of unsigned percent differences (UPD) as a bias indicator, and the average of signed percent differences (SPD) as a dispersion indicator: where the subscripts "in" and "ab" indicate in-and above-water, respectively. Fig. 3. Comparison of 3C-derived R rs (top row) and M99 -derived R rs (bottom row) with R rs determined from in-water multispectral µPRO data at specific center-wavelengths. The columns refer to data grouped as a function of the measurement conditions: ideal (left column, Central and Western Mediterranean Seas) and slightly sub-optimal (right column, AAOT).
The statistical results in Table 2 indicate a good comparability of M99 -and 3C-derived R rs for the ideal measurements conditions characterizing the data collected in the Central and Western Mediterranean Seas. Results are slightly less favorable for the AAOT data likely due to a lower accuracy of both in-water and above-water measurements performed with sub-optimal conditions. For the AAOT inter-comparison, 3C slightly outperforms M99 .
The impact of determining R rs from above-water radiometry by removing a percentage of the data in measurement sequences was assessed by applying the median of all the L T and L i measured spectra within each sequence already quality checked for high glint perturbations. Results evidenced somewhat deteriorated M99 retrievals for ideal measurement conditions (with spectrally averaged SPD=7.08%), whereas 3C-derived R rs degrade only slightly (with spectrally averaged SPD=7.58%). Conversely, for the slightly sub-optimal measurement conditions, the same process indicates some improvement for both M99 (with spectrally averaged SPD=-1.26%) and 3C (with spectrally averaged SPD=0.46%), except at 412 nm. These results appear to confirm the appropriateness of the processing solutions applied for data collected with ideal measurement conditions. However, results from data collected with sub-optimal measurement conditions, for which the applicability of M99 (and consequently M99 ) is questionable, further show how a non-ideal sky-radiance distribution may affect R rs retrievals. An aspect relevant to 3C is the impact of the bio-optical parameterizations in R rs retrievals. In fact, 3C makes use of bio-optical models that may not fit specific regional aquatic features. It is thus of interest to evaluate the impact of different parametrizations of the water biooptical characteristics. Whereas the spectral shapes of backscattering and non-algal absorption coefficients embedded in 3C can be assumed almost universally applicable, the phytoplankton absorption coefficient may exhibit unique features during particular seasons or in some regions. As already anticipated, 3C relies on a specific phytoplankton absorption spectrum a * ph typical for Lake Constance. This default a * ph was replaced with the very atypical spectrum of a Cryptophyta type L species [15]. 3C retrievals of R rs performed with this alternative a * ph showed a degradation of comparison results with the in-water derived data leading to a spectrally averaged SPD = 19.87% for ideal measurement conditions and a spectrally averaged SPD = -10.57% for the slightly sub-optimal measurement cases. This result highlights the importance of relying on a bio-optical modelling tailored to the regional bio-optical features.
The 3C performance was further investigated by excluding wavelengths beyond 800 nm to assess the impact of NIR spectral bands on R rs retrievals. Results showed a significantly degraded performance for ideal measurement conditions with a spectrally averaged SPD=10.4% and higher mismatch in blue wavelengths. Conversely, the effects were less noticeable for slightly sub-optimal measurement conditions, except for blue wavelengths displaying a clear overestimation. These results underline the relevance of NIR data to improve separation between aquatic and atmospheric contributions.
Lastly, 3C data re-processing was carried out by performing optimizations with 2000 generations and 200 individuals. This reprocessing, carried out on the same computer applied for the basic processing, required ∼5-6 min for each spectrum. However, results showed very small to no improvement with respect to the former optimization process. This finding supports the application of 3C and specifically of the related convergence processes, with fairly low computing resources.

NIOZ jetty station
The previous section provided an assessment of 3C-derived R rs from measurements performed during ideal or slightly sub-optimal conditions relevant for the validation of satellite data products [1]. However, often measurements are collected in regions heavily affected by cloudiness or from platforms such as towers, jetties or ships that cannot ensure optimal measurement geometries all the time. Still, derived R rs could be relevant for timely water quality monitoring purposes or bio-optical modelling. The suitability of 3C to produce R rs from above-water radiometry data collected with most illumination states and measurement geometries is hereafter addressed based on measurements performed at the NIOZ Jetty Station (NJS) (53°00 12 N, 4°47 16.5 E).
NJS, owned and managed by the Royal Netherlands Institute for Sea Research (NIOZ), is located nearby the Marsdiep inlet in the southern Wadden Sea. The nearby water is very turbid and affected by high resuspension, with Secchi disk depths generally smaller than 2 m. Spring blooms occur between April and May, followed by low production and turbidity in summer, and a late summer bloom. In autumn and winter, a high sediment resuspension due to strong winds leads to a higher turbidity. The tidal cycle largely drives the diurnal variability [23]. An optical radiometer system, located 5 meters above the water surface, has been operated almost continuously at NJS since 2001. This system [ Fig. 4(a)] includes six RAMSES hyperspectral radiometers: one ACC irradiance sensor measuring the downward irradiance (E s ); one ACC-UV irradiance sensor measuring the downward irradiance in the ultraviolet region whose data are not relevant for the following analysis; and two pairs of ARC radiance sensors operated in the south-east and south-west azimuthal directions, hereafter indicated as east and west for brevity. For both azimuth directions the total upwelling radiance L T (L T,e and L T,w ) with viewing angle θ v =35°and the sky radiance L i (L i,e and L i,w ) are measured accordingly. Fig. 4. a) Picture of the optical radiometer system operated at the NJS: 1. downward irradiance sensor with extended sensitivity in the ultraviolet (not relevant for this study); 2. downward irradiance sensor (E s ); 3. sea-viewing radiance sensor looking at south-east (L T,e ); 4. sea-viewing radiance sensor looking at south-west (L T,w ); 5. sky-viewing sensor looking at south-east (L i,e ); 6. sky-viewing sensor looking at south-west (L i,w ). b) Top-view schematic of the radiometers with respect to the sun position.
The measurement geometry leads to large sun glint perturbations in L T data during the course of the day [ Fig. 4(b)], and thus offers a distinctive case to assess the performance of 3C with data collected during extremely unfavorable conditions. More in general, the NJS data provide the opportunity to verify the consistency between R rs retrieved from measurements that were simultaneously collected over the same stretch of water with different measurement geometries for a variety of illumination conditions, and thus various levels of glint perturbation.
The hyperspectral above-water data collected at the NJS were processed with 3C applying the parameter bounds and initial guesses listed in Table 3. The spectral weights W of the goal function are determined by Eq. (13). Also in this case, a high weight at the NIR wavelengths has supported a better separation between marine and atmospheric contributions. Additionally, the spectral data affected by the oxygen absorption in the 750-775 nm interval as well as the highly noisy data beyond 920 nm were discarded by setting W = 0:

Analysis of one week of NJS data
This section focusses on 3C retrievals of NJS data collected at 5 minutes intervals from October 7 th until 13 th , 2012, and characterized by highly distinct illumination conditions. From October 7 th to 10 th the sky was patchy, with clouds intermittently covering the sun as indicated by the fluctuations affecting E s at 550 nm [Es(550)] shown in Fig. 5(a). On October 11 th , the smooth evolution of E s (550) indicates clear sky conditions, whereas on October 12 th and 13 th there were thick clouds as suggested by the permanently low E s (550) values. Intermittent rainfall occurred during these two last days as well.
Results from 3C retrievals of data from both the east and west NJS sensors are summarized in Fig. 5. The time series of the values displayed in Fig. 5(b) shows the cases for which 3C reliably reproduces the measured L T /E s ratio (a goal function threshold <0.02 was heuristically determined to identify successful retrievals). Results are considered satisfactory except for two specific conditions characterized by: i. highly variable sky illuminations, which may include the effects of cloud shadows; and ii. extreme sun glint leading to spectral values whose shapes cannot be reproduced by the models embedded in 3C. Occasional perturbations caused by foam, floating debris or rain drops, may not be either excluded. These are prone to a large increase of because 3C is unable to model their contributions to L T /E s . Similarly, sun glint may lead to saturated values or lens flare effects in the radiometer data, which cannot be handled with 3C. Overall, results appear to indicate that the goal function is a suitable quality indicator for these R rs retrievals. Figure 5(c) shows the weekly averages of the successfully retrieved R rs for the east and west directions. The shape of the spectrum indicates sediment-dominated waters and a strong presence of absorbing substances. High sediment concentrations are responsible for the non-zero R rs (NIR). The weekly averaged 3C retrievals for both directions are remarkably close, indicating that 3C Fig. 5. Results from 3C retrievals from NJS data. a) Time series of E s (550) nm evidencing different illumination conditions. The black symbols indicate all retrievals while the green symbols highlight the successful ones for at least one orientation (i.e., east or west). b) Time series of the goal function after each retrieval for the east and west sensors. c) Weekly averaged R rs for successful retrievals for both east and west orientations. d) Percent differences between the east and west individual successful R rs retrievals with respect to the mean of both, and the related inter-quartiles of the distribution. e) Time series of R rs percent differences with respect to their mean at 550 nm for the successful retrievals, before (black symbols) and after (red symbols) applying corrections for bidirectional effects. appears to successfully remove glint perturbations varying with the measurement geometry and illumination conditions, but the retrieved R rs values from both the east and west directions show instantaneous differences [ Fig. 5(d)] and daily patterns [ Fig. 5(e)] on sunny and partially sunny days. In general, R rs displays higher values for low relative azimuths φ between the sun and sensor, suggesting anisotropy effects. Indeed, the NJS time series shows a pattern with R rs,e > R rs,w in the morning, and vice-versa in the afternoon. Bidirectional effects were therefore compensated by applying a Case 2 correction scheme [24], restricted to θ s <70 • . Results indicate that the applied corrections reduce the differences during sunny (i.e., October 11 th ) and partially sunny (i.e., October 7 th to 9 th ) days. In particular, for October 11 th spectrally averaged SPD=-7.8% and UPD=11.9% between R rs,e (550) and R rs,w (550) determined without corrections, turned into SPD=-5.3% and UPD=9.0% when the corrections were applied. Conversely, as expected, the corrections induce artefacts in the data for overcast sky conditions (i.e., October 12 th and 13 th ). These are likely explained by the isotropic sky radiance distribution of a cloudy sky differing from that of a clear sky for which the correction scheme was developed.

Low light conditions
Based on observations collected during a clear sky day (October 11 th ), low light conditions defined here by the maximum spectral value of E s < 50 mW m −2 nm −1 , were also investigated. These conditions were met between 16:15 and 16:45 UTC, with extremely high sun-zenith angles θ s varying between 85°and 90°. For the processing of these data with M99, the values of ρ for θ s =80°were applied (i.e., 80°is the maximum value of θ s in the M99 look-up table). The relative azimuth φ varied between 116°and 122°for the east sensors and between 26°and 30°f or the west sensors. This led to the application of ρ = 0.0233 (ρ = 0.0282) for the east (west) sensor data. The wind data required for the determination of ρ were gathered from the De Kooy meteorological station of the Royal Netherlands Meteorological Institute located some 10 km away from NJS. Figure 6 shows both the corresponding 3C-and M99-derived R rs . Despite of the low light conditions, M99 (without NIR-bias removal) seems to produce relevant results for the east sensors due to a measurement geometry more adhering to consolidated protocols. For the west sensors, the small φ values indicate high glint affecting L T . Actually, in addition to major sun-glint contributions, this near-sun looking geometry is also challenged by the sharp changes of the sky radiance with φ, which make extremely inaccurate the quantification of the sky glint from L i measurements. Fig. 6. 3C-(left panels) and M99-derived (right panels) R rs from measurements performed during low light conditions. The upper panels refer to retrievals for the east sensors while the lower panels refer to retrieval for the west ones. The black lines indicate unsuccessful 3C retrievals (identified by >0.02) whereas the colored ones indicate the successful retrievals (determined by <0.02). The M99 retrievals in the right panels corresponding to the successful 3C ones, are displayed applying the same colors.
3C retrievals (left panels in Fig. 6) indicate that, among the measurements collected with challenging low light illumination, six out of eleven from the east direction and four out of eleven from the west one could be considered successful when applying <0.02 as the threshold for acceptable results. Overall, despite the measurement and illumination geometries being more challenging for the west sensors, 3C provides likely relevant R rs spectra from both viewing directions. This indicates that the air-water interface modelling applied in 3C [25,26] still performs with extreme sun zenith angles. Nevertheless, no conclusion can be drawn on the accuracy of these spectra and further investigations in this direction are required, ideally benefitting from reference measurements performed with alternative methods. However, it is important to note that at present, there are no consolidated protocols with instrument design, deployment method and uncertainty budget, defined for sub-optimal measurement conditions. Therefore, any R rs comparison relying on measurements performed in sub-optimal conditions with alternative methods would remain speculative.

High sun glint perturbations
NJS radiometry data from a clear sky day (i.e., June 29, 2019) are presented to showcase 3C-derived R rs from measurements affected by high sun glint. At the NJS latitude, during this period of the year, the sun zenith and the L i -viewing angle are close (i.e., θ v ≈ θ s when φ ≈ 0). These geometries lead to very high glint perturbations in L T measurements, and may saturate the L i sensors when these look close to the sun disk. These extreme measurement geometries are not compliant with the requirements set by consolidated protocols that aim at minimizing sun glint contributions to L T . For this reason, only the 3C retrievals are attempted.
The input data for 3C are shown in Fig. 7. In agreement with expectations for a clear sky day, E s (550) exhibits a bell shape [ Fig. 7(a)]. The L i (550)/E s (550) values [ Fig. 7(b)] increase for the east or west observations as φ approaches zero, reaching the saturation between 9 and 10 UTC for the first and between 13 and 14 UTC for the latter. Additional evidence of glint perturbations is given by the temporal development of the L T (850)/E s (850) ratio in Fig. 7(c). 3C-derived R rs determined with fixed ρ=ρ f are displayed in Fig. 8 (left panels). The glint correction performed using the saturated L i values causes negative R rs . When L i is not saturated, but still closely oriented toward the sun disk, the retrieved R rs (insets) display negative values in the blue spectral region, indicating an overestimated sun glint contribution. Figure 9 (left panels) shows the time series of and of the various terms composing the computed glint contributions. The peaks of normally indicate a degraded performance explained by the high glint (Fig. 9, upper left), especially when 3C is affected by the saturated L i values.    (Table 3) while the right panels refer to ρ values varying between 0 and ρ f ( Table 4).
The lower left panel in Fig. 9 indicates that L i /E s is dominant in the glint removal with ρ fixed to the Fresnel reflection coefficient ρ f . Therefore, to improve the retrieval of R rs from data affected by extreme sun glint perturbations, the 3C settings were modified as detailed in Table 4. Specifically, the reflection coefficient ρ applied to L i /E s was allowed to vary from zero to the value of the Fresnel coefficient. Additionally, the f sd and f ss coefficients, defining the contribution to the glint signal, were allowed to vary over a larger range to account for higher glint fractions.
Results obtained applying these less restrictive bounds show a significant improvement of the R rs retrievals from measurements heavily affected by sun glint (e.g., compare the values of in the left and right panels of Fig. 9). Notably, the values of still show an increase for conditions affected by high glint, but they are always lower than 0.04. This significant improvement of 3C-derived R rs displayed in the right panels of Fig. 9, was achieved in an unsupervised manner by determining near-zero ρ values through optimizations, thus de-facto discarding the L i data while alternatively determining the glint contribution from the direct irradiance fraction E sd /E s . Markedly, the east and west R rs values exhibit a different magnitude. The east and west cases were recorded at nearly equal time differences from the local noon, thus making their respective sun zenith angles almost equal. Additionally, the values of φ were close to 0 for both cases, thus leading to nearly equivalent measurement geometries. Therefore, the difference in R rs between east and west based measurements, is confidently explained by changes in the turbidity due to tidal effects [23].
Noteworthy is the noise affecting some of the spectra. This is likely due to retrievals performed with extremely perturbed measurements (i.e., the glint contributions is approximately 10 times higher than R rs in the green region of the spectrum). Yet, the fact that apparently relevant R rs spectra can be determined from near sun-looking measurement geometries is remarkable. Definitively, these data are unlikely to meet requirements for the validation of satellite data products. However, fully recognizing the impossibility of assigning any uncertainty estimate, when accompanied by suitable quality control procedures they may still support applications with less demanding accuracy needs.
Comparison between the left and right upper panels in Fig. 9 reveals that allowing for a bounded determination of ρ worsens the retrievals for conditions not affected by high glint perturbations, opposite to the application of a fixed value equal to the Fresnel reflection coefficient. This is evidenced by a slight increase, and a high variability, of the goal function values indicating a more problematic convergence. This finding suggests that the Fresnel reflection term is generally a good first estimate in an operational context. However, the implementation of alternative configurations, one with fixed ρ = ρ f for low to moderate levels of glint and the other with variable ρ for very large glint perturbations, appears a viable solution. Since high glint contributions can be anticipated from measurement geometries and are easily detectable (e.g., by L i /E s >>1/π), the handling of the two configurations can certainly be automated.

Glint filtering in a sequence of measurements under high wind
The impact of glint filtering in above-water radiometry data over a highly roughened sea (determined by a wind speed w = 9 m s −1 ) was also investigated. Data were collected at a station in the Dutch Wadden Sea on October 12 th 2017 at 12:30 UTC. The water was turbid, dominated by sediments and with a total absorption coefficient of approximately 4 m −1 at 440 nm. E s , L i and L T measurements were made with TriOS RAMSES sensors operated on a portable frame, with observation geometry determined by relative azimuth φ = 135 • and viewing angle θ v = 35°. The sequence consisted of 20 consecutive measurements for each sensor and was completed in 40 seconds. Figure 10 shows the acquired data. The illumination conditions were ideal with negligible variability affecting E s and L i and therefore their medians were taken as representative of the measurement sequence. As expected, the L T spectra show a high variability as a result of the glint contributions generated by wave facets. 3C was applied to each individual L T measurement and to the medians of E s and L i . Each retrieved R rs was evaluated as a function of the returned value of the goal function . The parameter search bounds and spectral weights were the same applied for the NJS case analysis [Eq. (13), Table (3)].  Figure 11 displays the retrievals for each L T /E s in the measurement sequence. Results show that the model effectively matches the measured L T /E s over the dynamic range of the various measurements, except for one case for which the convergence was not achieved (see the outlier displayed in Fig. 12(a) exhibiting a value of approaching 0.05). However, when L T /E s increases in magnitude due to higher glint contributions, 3C encounters increasing difficulties in matching the fine spectral features characterizing L T /E s , as highlighted by the correlation between L T /E s and . This finding suggests that the general scheme for glint filtering, introduced in section 3.1, is advisable as a pre-processing step to 3C retrievals when multiple measurements are collected.
Indeed, a condition for glint filtering is the availability of a sufficient number of L T values exhibiting low variance within each measurement sequence: a condition satisfied for the L T data showing the lowest values [see Fig. 12(b)]. In particular, for the specific example, if the lowest 20% L T are applied to compute a median L T and determine R rs using 3C, the optimization yields = 6.5 · 10 −3 , close to = 6.1 · 10 −3 obtained for the individual measurement exhibiting the best optimization. Fig. 11. 3C retrievals for each L T of the measurement sequence. Measured L T /E s , modelled L T /E s , R G and R rs are identified by the blue, green, red and black lines. The scaling is the same for all graphs. b) 3C-derived R rs for each L T measurement of the sequence (black lines), and 3C-derived R rs from the median of the 20% lowest L T (green lines). The red spectrum identifies the R rs associated to the lowest .
The analysis indicates that 3C retrievals from each individual measurement of a sequence may allow the quality control of the R rs through statistical approaches. This naturally suggests the need for further investigations on the quality control of 3C-derived R rs by relying on time series of L T data.
Despite the likely favorable measurement geometry (yet, not that specifically recommended by consolidated protocols), high wind speed can produce significant sun-glint perturbations in L T . Therefore, the analysis of this measurement sequence was repeated, allowing ρ to vary from zero to ρ f as previously done in section 3.2.3. The optimized ρ are negatively correlated with f sd (results not shown), thus confirming the complexity of the glint contributions. Interestingly, the optimizations returning the lowest , for which ρ≈ρ f and f sd ≈0, are also those pertaining to the 20% lowest L T . For those retrievals, the values of are comparable to those obtained with a fixed ρ=ρ f . This suggests that an optimization scheme with fixed ρ=ρ f can be confidently applied to process those measurements exhibiting low sun-glint contamination as resulting with low wind speed. Additionally, the removal of some of the L T spectra exhibiting larger values prior to 3C processing appears a suitable option in the case of high glint contamination.
Finally, this latest analysis clearly indicates that cannot provide a universal answer on the quality of R rs retrievals. For instance, the threshold = 0.02, suggesting a satisfactory quality of R rs from the NJS data is definitively overestimated for the retrievals from the measurement sequence applied to investigate glint effects with relatively high wind speed. A goal function normalized to the reflectance intensity may mitigate such a dependence.

Conclusions
The 3C model exhibits remarkable flexibility in the determination of R rs from above-water radiometry obtained in a wide range of illumination conditions and measurement geometries. When compared to above-water methods recommended for satellite ocean color validation, which rely on ideal illumination conditions (i.e., clear sky), low wind, and specific measurement geometries minimizing glint perturbations, 3C exhibits similar accuracy. However, for less favorable illumination and measurement geometries, 3C still shows the ability to produce significant R rs retrievals. This qualifies 3C for the determination of R rs from autonomous measurement platforms under almost all weather conditions and geometries, still recognizing the difficulty to assign strict quality levels and more importantly uncertainties, to the retrieved R rs .
Thresholds applied to the goal function appear valuable to qualitatively rank the R rs retrievals. Appropriate thresholds for may be chosen based on the application at hand, by trading off between valid retrievals and required accuracy. However, the values of may largely vary with the quality and magnitude of the input data. Because of this, the use of universal thresholds relying on to quality assure 3C retrievals is discouraged and should be matter of future investigations.
Notably, 3C retrievals do not require ancillary data such as wind speed, cloud cover, or relative azimuths. The performance of 3C relies on its flexible design that leads to the determination of glint contributions by matching the total radiance reflectance L T /E s to semi-analytically modelled L T /E s data only. The inclusion of L i is beneficial to reach convergence, except in extreme cases (e.g., with near sun-looking geometries).
It is emphasized that 3C is not intended to be a full inversion method for water constituents. The retrieved model parameters (e.g., chlorophyll) are byproducts obtained by the reconstruction of glint and R rs spectra. In this context, possible overfitting does not have negative effect on R rs as long as the atmospheric and aquatic contributions can be reliably separated. Convergence is ensured as far as the bio-optical modelling can realistically reproduce the R rs spectrum, which is generally achieved with the aquatic optical properties included in the AM03 model described here. For unusual events such as cyanobacteria blooms, absorption spectra may fall outside the ranges specified for AM03. In this case, the bio-optical modelling included in 3C might then not be able to accurately reproduce observations and therefore return accurate R rs . This issue could be mitigated by a choice of aquatic optical properties, representative of a number of water types, combined with the application of statistical approaches ranking the likelihood of the retrieved R rs spectra [27]. Additional benefit is obtainable from realistic search bounds for each 3C parameter as well as weights for the goal function avoiding spectral areas affected by noise or processes unaccounted for in the bio-optical modelling (e.g., fluorescence). In future developments, experimental determinations of the diffuse and direct irradiance fractions could be considered to simplify the modelling of the atmospheric components.
In view of supporting a wider use of 3C, a Matlab implementation of the model is made accessible in the GitLab repository (https://gitlab.com/jaipipor/rrs_model_3c_matlab).

Funding
NIOZ: Ministry for Science and Culture of Lower Saxony, Germany: Project "Coastal Ocean Darkening" (VWZN3175); Joint Research Centre.