Validation of a spectral correction procedure for sun and sky reflections in above-water reflectance measurements

A three-component reflectance model (3C) is applied to above-water radiometric measurements to derive remote-sensing reflectance Rr s (λ). 3C provides a spectrally resolved offset ∆(λ) to correct for residual sun and sky radiance (Rayleighand aerosol-scattered) reflections on the water surface that were not represented by sky radiance measurements. 3C is validated with a data set of matching aboveand below-water radiometric measurements collected in the Baltic Sea, and compared against a scalar offset correction ∆. Correction with ∆(λ) instead of ∆ consistently reduced the (mean normalized root-mean-square) deviation between Rr s (λ) and reference reflectances to comparable levels for clear (∆: 14.3 ± 2.5 %, ∆(λ): 8.2 ± 1.7 %), partly clouded (∆: 15.4 ± 2.1 %, ∆(λ): 6.5 ± 1.4 %), and completely overcast (∆: 10.8 ± 1.7 %, ∆(λ): 6.3 ± 1.8 %) sky conditions. The improvement was most pronounced under inhomogeneous sky conditions when measurements of sky radiance tend to be less representative of surface-reflected radiance. Accounting for both sun glint and sky reflections also relaxes constraints on measurement geometry, which was demonstrated based on a semi-continuous daytime data set recorded in a eutrophic freshwater lake in the Netherlands. Rr s (λ) that were derived throughout the day varied spectrally by less than 2 % relative standard deviation. Implications on measurement protocols are discussed. An open source software library for processing reflectance measurements was developed and is made publicly available. © 2017 Optical Society of America OCIS codes: (280.0280) Remote sensing and sensors; (240.6690) Surface waves; (120.5700) Reflection; (010.7340) Water; (010.4450) Oceanic optics; (010.1110) Aerosols; (010.1615) Clouds. References and links 1. C. Mobley, Light and water: Radiative transfer in natural waters (Academic Press, 1994). 2. C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements.” Appl. Opt. 38, 7442–7455 (1999). 3. G. Zibordi, S. B. Hooker, J. F. Berthon, and D. D’Alimonte, “Autonomous above-water radiance measurements from an offshore platform: A field assessment experiments,” Journal of Atmospheric and Oceanic Technology 19, 808–819 (2002). 4. S. B. Hooker, G. Lazin, G. Zibordi, and S. Mclean, “An evaluation of aboveand in-water methods for determining water-leaving radiances,” Journal of Atmospheric and Oceanic Technology 19, 486–515 (2002). 5. C. Cox and W. Munk, “Measurement of the Roughness of the Sea Surface from Photographs of the Sun’s Glitter,” J. Opt. Soc. Am. 44, 838–850 (1954). 6. C. D. Mobley, “Polarized reflectance and transmittance properties of windblown sea surfaces,” Appl. Opt. 54, 4828–4849 (2015). 7. E. Aas, “Estimates of radiance reflected towards the zenith at the surface of the sea,” Ocean Science Discussions 7, 1059–1102 (2010). 8. K. Ruddick, V. De Cauwer, and B. Van Mol, “Use of the near infrared similarity reflectance spectrum for the quality control of remote sensing data,” in “Remote Sensing of the Coastal Oceanic Environment,” , vol. 5885-1 (2005), vol. 5885-1. 9. K. G. Ruddick, V. De Cauwer, Y.-J. Park, and G. Moore, “Seaborne measurements of near infrared water-leaving reflectance: The similarity spectrum for turbid waters,” Limnology and Oceanography 51, 1167–1179 (2006). 10. S. G. Simis and J. Olsson, “Unattended processing of shipborne hyperspectral reflectance measurements,” Remote Sensing of Environment 135, 202–212 (2013). 11. V. Martinez-Vicente, S. G. H. Simis, R. Alegre, P. E. Land, and S. B. Groom, “Above-water reflectance for the evaluation of adjacency effects in earth observation data: Initial results and methods comparison for near-coastal waters in the western channel, UK,” Journal of the European Optical Society 8 (2013). 12. Z. Lee, Y.-h. Ahn, C. Mobley, and R. Arnone, “Removal of surface-reflected light for the measurement of remotesensing reflectance from an above-surface platform,” Opt. Express 18, 171–182 (2010). 13. T.-W. Cui, Q.-J. Song, J.-W. Tang, and J. Zhang, “Spectral variability of sea surface skylight reflectance and its effect on ocean color.” Opt. Express 21, 24929–24941 (2013). 14. L. G. Sokoletsky and F. Shen, “Optical closure for remote-sensing reflectance based on accurate radiative transfer approximations: the case of the Changjiang (Yangtze) River Estuary and its adjacent coastal area, China,” International Journal of Remote Sensing 35, 4193–4224 (2014). 15. D. A. Toole, D. A. Siegel, D. W. Menzies, M. J. Neumann, and R. C. Smith, “Remote-sensing reflectance determinations in the coastal ocean environment: impact of instrumental characteristics and environmental variability.” Appl. Opt. 39, 456–469 (2000). 16. X. Zhang, S. He, A. Shabani, P.-W. Zhai, and K. Du, “Spectral sea surface reflectance of skylight,” Opt. Express 25, A1 (2017). 17. P. Gege, “The water color simulator WASI: an integrating software tool for analysis and simulation of optical in situ spectra,” Computers & Geosciences 30, 523–532 (2004). 18. W. Gregg and K. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limonology And Oceanography 35, 1657–1675 (1990). 19. P. Gege, “Analytic model for the direct and diffuse components of downwelling spectral irradiance in water.” Appl. Opt. 51, 1407–1419 (2012). 20. K. Dörnhöfer, A. Göritz, P. Gege, B. Pflug, and N. Oppelt, “Water Constituents and Water Depth Retrieval from Sentinel-2A A First Evaluation in an Oligotrophic Lake,” Remote Sensing 8, 941 (2016). 21. P. Gege, “WASI-2D: A software tool for regionally optimized analysis of imaging spectrometer data from deep and shallow waters,” Computers and Geosciences 62, 208–215 (2014). 22. P. Gege, “A Case Study at Starnberger See for Hyperspectral Bathymetry Mapping Using Inverse Modeling,” in “WHISPERS 2014,” (2014), pp. 1–4. 23. J. T. O. Kirk, Light and photosynthesis in aquatic ecosystems (Cambridge University, 1994). 24. C. Stedmon, S. Markager, and H. Kaas, “Optical Properties and Signatures of Chromophoric Dissolved Organic Matter (CDOM) in Danish Coastal Waters,” Estuarine, Coastal and Shelf Science 51, 267–278 (2000). 25. A. Albert and Mobley, “An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters.” Opt. Express 11, 2873–2890 (2003). 26. H. Buiteveld, J. H. M. Hakvoort, and M. Donze, “Optical properties of pure water,” SPIE 174, 2258 (1994). 27. A. Morel, “Optical properties of pure water and pure sea water,” in “Optical Aspects of Oceanography,” , N. G. Jerlov and E. Steemann Nielsen, eds. (Academic Press, 1974), chap. 1, pp. 1–24. 28. P. Gege, “Characterization of the phytoplankton in Lake Constance for classification by remote sensing,” ArchHydrobiolSpecIssues AdvLimnol 0, 179–193 (1998). 29. T. Heege, “Flugzeuggestutzte Fernerkundung von Wasserinhaltsstoffen im Bodensee,” Ph.D. thesis, Freie Universitat Berlin (2000). 30. Z. Lee, K. L. Carder, C. D. Mobley, R. G. Steward, and J. S. Patch, “Hyperspectral remote sensing for shallow waters. I. A semianalytical model.” Appl. Opt. 37, 6329–6338 (1998). 31. H. R. Gordon, J. W. Brown, O. B. Brown, R. H. Evans, and R. C. Smith, “A semianalytic radiance model of ocean color,” Journal of Coastal Research 93, 10909–10924 (1988). 32. P. Gege and P. Groetsch, “A spectral model for correcting sun glint and sky glint,” in “Proceedings of Ocean Optics XXIII,” (2016). 33. R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A Limited Memory Algorithm for Bound Constrained Optimization,” SIAM Journal on Scientific Computing 16, 1190–1208 (1995). 34. H. J. Gons, “Optical teledetection of chlorophyll a in turbid inland waters,” Environmental Science and Technology 33, 1127–1132 (1999). 35. D. Stramski and J. Dera, “On the mechanism for producing flashing light under a wind-disturbed water surface,” Oceanologia 25 (1988). 36. Y. You, D. Stramski, M. Darecki, and G. W. Kattawar, “Modeling of wave-induced irradiance fluctuations at near-surface depths in the ocean: a comparison with measurements.” Appl. Opt. 49, 1041–1053 (2010). 37. M. Hieronymi and A. Macke, “Spatiotemporal underwater light field fluctuations in the open ocean,” Journal of the European Optical Society: Rapid Publications 5, 1–8 (2010).


Introduction
The spectral signature of a natural water body contains information about its optically active components, such as phytoplankton and other suspended matter, and coloured dissolved organic matter.Quantitative knowledge of these constituents and their dynamics is crucial to environmental monitoring and understanding ecosystem dynamics, and can inform how these environments are affected by human activities and climate change.Spectral signatures can be obtained at high accuracy from water samples in a laboratory environment, and acquired remotely from satellites, aeroplanes, unmanned aerial vehicles, ships-of-opportunity, coastal laboratories, and handheld devices.Passive optical remote-sensing methods yield particularly fast and cost-effective observation coverage compared to relatively laborious water sample collection and laboratory analysis.However, changing illumination conditions and wave patterns induce temporally and spatially highly variable sun and sky reflections on the water surface, which pose a particular challenge to accurate retrieval of spectral signatures of natural water bodies from remotely sensed observations.
Remote-sensing reflectance R r s (λ) is a widely used parameter to express the spectral reflectance signature of a water body and is defined as with water-leaving radiance L w (λ), and downwelling irradiance E d (λ) just above the water surface.For an overview of symbols and units see Table 1.E d (λ) may be divided into two spectrally distinct components: solar irradiance directly transmitted through the atmosphere from the direction of the sun, and diffuse irradiance scattered in the upper hemisphere by e.g.air molecules, aerosols or clouds.The corresponding specular reflections on the water surface constitute water surface-reflected radiance L r (λ).Neither L w (λ) nor L r (λ) can be measured directly from above the water surface, but their sum is recorded simultaneously as upwelling radiance L u (λ) [1], such that Eq. (1) becomes: A widely used approach to correct for surface reflections is to express L r (λ) as the fraction ρ s of sky radiance L s (λ, θ view , φ), measured in specular direction to L u (λ, θ view , φ).Here, θ view is the viewing zenith angle and φ indicates the azimuth angle to the sun.For this widely used measurement protocol, R r s (λ) is then expressed as (omitting viewing geometry for brevity): For a perfectly flat water surface or a perfectly uniform sky, L s (λ) is spectrally representative of L r (λ) and the air-water interface reflectance factor ρ s equals the Fresnel reflectance factor ρ f (θ view , n w ), which is a function of θ view and the refractive index of water n w [1].Neither of these conditions is commonly encountered, and thus L r (λ) will need to be determined by other means in operational contexts.Moreover, L r (λ) can be orders of magnitude higher than L w (λ), especially when sun glint or cloud reflections contribute.A good estimate of L r (λ) is therefore essential to derive R r s (λ).
Mobley [2] reported that sun glint contributions for most wind-roughened water surfaces are at a minimum when observing with a viewing geometry of approximately θ view = 40 • and φ = 135 • (Fig. 1).Sun glint contributions can be further minimized by selecting the lowest of several consecutive measurements of L u (λ) [e.g.3,4].However, even with the most strict filtering thresholds it is not possible to rule out contamination by sun glint because typical integration times -from tens of milliseconds to seconds -are too long to resolve short exposure to sun Table 1.Free model parameters (upper section) in 3C and L10 model optimization.Parameters not allowed to vary in the optimization were set to their start value (indicated by *).Parameters that were variable, but were determined outside of the optimisation, are marked with **.Parameters that were not used in a model run are marked with -.In the following, chlorophylla is abbreviated as chl, suspended particulate matter as SPM, and coloured dissolved organic matter as CDOM.The lower section lists further relevant parameter names and acronyms that are used in this manuscript.Units for radiances (mWm −2 nm −1 sr −1 ) and irradiances (mWm −2 nm −1 ) are abbreviated with 1 glint, e.g.glitter patterns on waves.Sun glint contributions can not be minimized if redundant measurements are not viable or viewing geometry is fixed, e.g. for airborne and satellite remote sensing applications, or fixed position installations.
air water Parameter dependencies were omitted for brevity and for further information on parameters and units, please refer to Table 1.
For surface-reflected sky radiance, several approaches exist to determine the value of ρ s in Eq. (3) for non-flat water surfaces.Water surface slope distributions [5,6] can be used to statistically derive ρ s for various illumination conditions, wind speeds, and viewing geometries [2,7,6].Ruddick et al. [8,9] derived from theory and observations that near-infrared (700-900 nm) reflectance measurements are primarily shaped by pure water absorption, provided the water is sufficiently turbid to yield exploitable signal in this spectral domain.Sky radiance reflected at the water surface also changes the reflectance signal in the near-infrared.The reported shape of near-infrared reflectance can therefore be used to determine reflected skylight in reflectance measurements ('similarity spectrum approach') when assuming no other sources of error.Simis and Olsson [10] exploited the property of water reflectance spectra to exhibit only broad and spectrally slowly varying features, compared to more finely featured downwelling irradiance and radiance spectra.The authors proposed an optimization scheme that determines ρ s such that atmospheric gas absorption features are minimized in the reflectance spectrum R r s (λ) ('fingerprint approach').Martinez-Vicente et al. [11] applied the latter two approaches to reflectance observations from a moving research-vessel in coastal waters of the Western Channel, UK, and found large variations in the resulting ρ s values.Both solutions assume that ρ s is spectrally constant, whereas Lee et al. [12] point out that simulated L r (λ) can differ by a factor 8 between 400 and 800 nm from the approximated L r (λ) = ρ s • L s (λ), and conclude that ρ s should be treated as a function with strong wavelength dependency [e.g.12,13,14].
Despite the relative success of these approaches under a wide range of measurement conditions, the product of ρ s and L s (λ) is only truly representative of L r (λ) when sky radiance is uniformly distributed, e.g. for fully overcast conditions [2,15].This is because water surfaces are practically always wind-roughened to the extent that wave facets reflect different parts of the sky and even the sun disc.These elements contribute to observed L r (λ) and typically differ from L s (λ) in both spectral appearance and intensity [16].Lee et al. [12] acknowledged this [12,Eq. (8)], yet state that the spectral dependence of radiance reflected from a roughened water surface is generally unknown.The authors then introduced a spectrally independent (scalar) correction factor ∆ to Eq. ( 3) (further referred to as L10), which can be estimated through bio-optical modelling [12] or with the similarity spectrum approach [8]: It had been proposed earlier [e.g.17] that the spectral dependence of L r (λ) may be resolved by accounting for spectrally distinct sun glint and surface-reflected sky irradiance.This approach distinguishes three irradiance components: E dd (λ) for direct solar irradiance, E dsr (λ) for diffuse molecular-scattered irradiance (Rayleigh component), and E dsa (λ) for diffuse aerosol-scattered irradiance.These components can be calculated based on the Eqs. of Gregg and Carder [18] (further referred to as the GC90 model) and atmospheric gas absorption spectra [19].The GC90 model has been derived for cloudless maritime atmospheres and requires aerosol optical properties as input.These properties can be measured (using e.g. a sun photometer), or adopted from atmospheric correction results in case of airborne or satellite remote-sensing observations [20].Alternatively, aerosol properties can be derived from from L u (λ) or from L u (λ) E d (λ) ratios in combination with bio-optical modelling [21].Application of this approach to hyperspectral airborne [22] and Sentinel-2A satellite imagery [20] has been demonstrated for cloudless atmospheres.This three-component approach is further referred to as the 3C model.
In the present study, we apply 3C to in situ remote-sensing reflectance data sets obtained under a broad range of illumination and wave conditions.The model replaces the scalar offset ∆ in Eq. ( 4) with a spectrally resolved correction factor ∆(λ). Validation of this spectral correction scheme is approached in two steps.First, radiometric measurements of Baltic Sea waters are processed to R r s (λ) with scalar and spectrally resolved offset correction factors, ∆ (L10) and ∆(λ) (3C), respectively.Reference reflectances derived from sub-surface observations are used to compare model performance under clear, partially clouded, and overcast sky conditions.Second, 3C and L10 are applied to semi-continuous daytime radiometric observations recorded by a fixed-position hyperspectral reflectance sensor in a eutrophic lake, to assess whether a spectrally resolved correction factor ∆(λ) relaxes measurement geometry requirements.Implications on measurement protocols are discussed.

Study sites
In situ data were obtained prior to this study during spring and summer research cruises of RV Aranda at 44 stations in the Baltic Sea.Their locations are shown in Fig. 2.Only stations where above-and below-water radiometric measurements were simultaneously recorded were selected.Several accompanying parameters, such as absorption of light by coloured dissolved organic matter (CDOM, see section 2.1.4)and wind speed, were recorded for each station.
A separate data set was used to assess the effect of non-optimal viewing geometries on sun glint and surface-reflected sky radiance under otherwise stable conditions.Semi-continuous radiometric measurements for two azimuth angles were recorded during daylight hours by a fixed-position instrument in a fresh water lake in the Netherlands.Variability in derived reflectances and model parameters were analysed with respect to viewing geometry.Fig. 2. Stations that were sampled in the Baltic Sea between 2010 and 2012 and for which both above and below surface radiometric measurements were carried out.Four cases that represent a range of sky conditions are marked with an x, and are presented in detail in Fig. 4.

Sub-surface reflectance
Up-and downwelling spectral irradiance E u (z, λ) and E d (z, λ) were recorded as a function of depth z from 0 to 15 m with TriOS Ramses-ACC-VIS sensors (λ = 320 − 950 nm).An additional sensor located at a high point on the ship recorded E d (λ) to account for fluctuations in downwelling irradiance.During acquisition, the ship was positioned such that no shadow was cast onto the sensors.Instrument self-shading correction was applied to measurements of E u (λ), based on spectral absorption measurements that were simultaneously recorded by a Wetlabs AC-S instrument.The self-shading correction factor for the wavelength range 400-700 nm averaged to 1.04 ± 0.02.Sampling depth was determined from pressure readings at sensor level (Sea Bird Electronics SBE-50).Irradiance reflectance just below the water surface R − (λ) was calculated by extrapolating exponential fits to depth profiles of E u (λ) and E d (λ) to zero depth.The CDOM-rich water of the Baltic Sea efficiently absorbs light in the ultra-violet and near-infrared spectral regions and reflectance is typically weak, so that measurements in these bands suffered from a low signal to noise ratio preventing accurate fits to the data.Evaluation of reflectance profiles was therefore restricted to 400-700 nm.

Above-surface reflectance
Above-surface measurements of spectral upwelling radiance L u (λ, θ view , φ), spectral downwelling radiance in the specular direction L s (λ, θ view , φ), and E d (λ) were recorded continuously from a prototype Rflex system.Rflex keeps viewing azimuth angles as close as possible to 135 • and always > 90 • with respect to the solar azimuth [10], and was installed on the bow of RV Aranda.Observations were considered matching sub-surface reflectance measurements if they were recorded within ± 30 minutes and 500 meters of the cast point.The system was equipped with inter-calibrated TriOS Ramses-ACC-VIS irradiance and two ARC-VIS radiance sensors (λ = 320 − 950 nm).Viewing angles θ view , θ view of down-and upwelling radiance sensors were fixed to 40 • from zenith and nadir (Fig. 1), respectively, and were not corrected for ship roll or tilt.From this data set, observations were discarded if either of L s (λ), L u (λ), or E d (λ) deviated more than 30% from the respective station mean spectrum, indicating the measurement was likely affected by sea spray, white caps, sensor tilt, or obstacles within the sensor field-of-view.This quality control filter on the data is expressed as: with ì x(λ) indicating a series of consecutive spectra.The term z( ì represents a series of z-score normalized spectra ì x(λ) with corresponding means ì µ and standard deviations ì σ. z( ì x(λ)) is the mean spectrum of such a series of z-score normalized spectra.If any of the radiance or irradiance spectra belonging to a single observation exceeded the selection criterion, the whole observation was omitted.In addition, observations were considered suspect (e.g.affected by foam-caps, sea-spray, or surface scum) if the ratio L u (λ) E d (λ) exceeded an empirical threshold of 0.025 sr −1 at any wavelength between 800 and 950 nm.This quality control procedure generally preserved in the order of 70 observations per station.
Semi-continuous measurements were recorded by a fixed-position Ecowatch instrument in the eutrophic Lake Paterswoldsemeer (Paterswoldsemeer), the Netherlands (Latitude: 53.1631 • ; Longitude 6.5773 • , mean depth in summer 0.8 m).Chlorophyll-a concentrations in Paterswoldsemeer average to 30-50 µgl −1 in summer and cyanobacterial blooms are frequently observed.The lake bed consists of peat and mud, which frequently become re-suspended and contribute to low average Secchi-depths of 0.5 to 0.6 m.The instrument carries two sets of inter-calibrated L u (λ, θ view , φ), L s (λ, θ view , φ), and E d (λ) channels (380-950 nm, re-sampled to 1 nm step size), which faced north-northeast (NNE) and north-northwest (NNW), respectively.Viewing angles θ view , θ view of down-and upwelling radiance sensors were fixed to 40 • from zenith and nadir (Fig. 1), respectively.On 11 October 2015 ten measurements from each channel were recorded every ten minutes in the period 10:30 am to 4:00 pm (local time).The resulting 680 measurements of L s (λ), L u (λ), and E d (λ) (2 channels × 10 measurements per cycle × 34 cycles) was not subjected to additional quality control procedures to allow a realistic assessment of sun glint and surface-reflected sky radiance under changing environmental conditions.

Accompanying measurements
For every Baltic Sea station, light absorption of coloured dissolved organic matter (CDOM) a CDOM (λ) was determined from water samples filtered through a 0.2 µm membrane filter and analysed with a dual beam spectroradiometer (190-800 nm) against a reference of ultrapure water.The absorbance measurements were converted to units of absorption.Measurements from cuvettes of different lengths (1, 4, 5, or 10 cm) were combined to optimise signal quality in both the ultraviolet and visible to near infra-red range.The exponential model a CDOM (λ) = a 440 nm CDOM • exp(−S • (λ − 440 nm)) + k was optimized in a least-squares minimization for the wavelength range 350-600 nm to derive CDOM absorption slope S and offset k [e.g.23].The offset compensates for residual noise or scattering [24].Wind speed v wind (in m s −1 ) was retrieved at each Baltic Sea station from the ship weather station and averaged to 5.7 ± 2.9 m s −1 (0.4 -12.5 m s −1 ).Solar azimuth θ sun was calculated based on GPS location and time using the open-source (LGPL) Python library PyEphem.For the Paterswoldsemeer data set no accompanying measurements were taken.

Models
The 3C model corrects measurements of L u (λ) E d (λ) with respect to sun glint and surface-reflected sky light, resulting in remote-sensing reflectance R r s (λ).The following sections elaborate on the underlying analytical model (section 2.2.1), model inversion (section 2.2.2), and model evaluation (section 2.2.3).
The GC90 [18] parametrisation of downwelling direct and diffuse irradiance was adapted to spectrally resolve the offset parameter ∆(λ) in Eq. ( 6), as follows: with ρ dd and ρ ds as the reflectance factors for direct and diffuse downwelling irradiance, respectively.Based on the components of GC90, irradiance ratios were calculated according to Gege and Groetsch [32]: with Rayleigh transmission coefficient T r (λ), and aerosol transmission coefficient T as (λ).The wavelength dependency of τ a (λ) is approximated by the Ångström law with turbidity coefficient β at λ a = 550 nm, and Ångström exponent α.Further details on the calculation of atmospheric path length M, atmospheric path length corrected for air pressure M , aerosol forward scattering probability F a , and aerosol single scattering albedo ω a are given in Gege [19].The satisfactory performance of this model was demonstrated in Gege and Groetsch [32], where it was able to accurately reproduce several hundreds of E d (λ)-normalized sky radiance measurements collected under various cloudless atmospheric conditions.Influence of clouds are not covered by ∆(λ) but taken into account with sky radiance measurements in eq. ( 6).Examples for derived Rayleigh and aerosol transmission coefficients, sky reflectance, and sun glint are also shown in Gege and Groetsch [32].All parameter abbreviations are summarized in Table 1.Fig. 3. Spectrally resolved weighting factors that were used in 3C model optimization.

Model inversion
Least-squares minimization using the constrained limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS-B) [33] of the residual-sum-of-squares (RSS) between observed λ) and modelled was carried out with the bounding parameters listed in Table 1 and weighting vector W(λ).The 675-750 nm spectral region was weighted at 10 % because chlorophyll-a fluorescence is not accounted for in the bio-optical model.Similarly, spectral weighting (10 %) was assigned to the range 760-775 nm to lessen the influence of the deep oxygen A-band rotational absorption lines near 765 nm on the fit result.In contrast, the spectral region up to 500 nm was weighted with a factor of 5 to add extra weight to accurate modelling of Rayleigh-scattered radiance.Weights were determined by experimentation and are shown in Fig. 3.The RSS was calculated as follows: ∆(λ, α, β, ρ dd , ρ ds , θ sun ) was subsequently calculated from the derived model parameters α, β, ρ dd , ρ ds , leading to the definition of R r s (λ) from Eq. ( 4) extended with spectrally dependent ∆(λ): The same optimisation procedure was carried out using the L10 model and parameters as indicated in Table 1 to derive the scalar offset factor ∆, which was then used instead of ∆(λ) in Eq. ( 21) to calculate R r s (λ).The minimisation procedure and its modelling components (GC90, 3C, inversion scheme) were implemented in Python and are publicly available under open license xyz (free for non-commercial use).The version of this package used to produce the results of this paper is available on Zenodo (http://doi.org/10.5281/zenodo.293851).Future versions will be available on gitlab (https://gitlab.com/pgroetsch/rrs_model_3C).GC90 and 3C are also implemented in the public domain software "Water Color Simulator" (WASI) [17] (http://www.ioccg.org/data/software.html).

Model evaluation
For each Baltic Sea station individual R r s (λ) spectra were derived with the spectral optimization procedure of section 2.2.2 and parameters as indicated in Table 1.A pre-fit was carried out on the mean-averaged L u (λ) E d (λ) spectrum of each station to determine optimal parameter starting values for the individual fits.Subsurface reflectance R − (λ) (section 2.1.2) was translated to above-surface reference reflectance R with water-to-air radiance-divergence factor ζ = 0.518, water-to-air reflectance for isotropic irradiance γ = 0.48 [31], the ratio of upwelling irradiance to upwelling radiance Q, and a scalar offset factor δ. Q equals π for an isotropic distribution, with typical values from 4 to 5 sr [31].
The scalar offset δ accounts for shortcomings in model inversion (e.g.due to ambiguities between bio-optical and surface reflection model) or data acquisition (e.g due to instrument self-shading, white caps and foam, or ship tilt and roll).Despite not being of primary interest for the present study, these differences in scaling (Q) and offset (δ) usually need to be accounted for to allow unbiased comparison of R − (λ) and R r s (λ).Analytical models for the Q-factor [e.g.25,14] do not account for wave-induced variability.Therefore, both factors were treated as bounded parameters when minimizing the root-mean-square-error (RMSE) between R r s (λ) from above and below surface observations, with geometric parameter Q bounded between 3 and 8 sr (starting value: 5 sr) and offset δ bounded between -0.01 and 0.01 sr −1 (starting value: 0).Accounting for δ ensures that error metrics (RMSE, and normalized RMSE, nRMSE) are indicative of residual spectral differences only, which are attributed to suboptimal correction for sun glint and reflected sky light in R r s (λ).nRMSE was calculated by normalizing RMSE with the mean R r e f r s (λ) value.Daytime measurements at Paterswoldsemeer provided two R r s (λ) channel pairs, facing NNW and NNE.During each observation, ten individual measurements were recorded over typically less than 10 seconds, depending on integration time.Observations from both channels were processed to R r s (λ) with 3C and L10 according to the inversion method in section 2.2.2.No accompanying a CDOM (λ) measurements were available and thus CDOM absorption slope S was set to a value of 0.012 nm −1 in the optimization.This value was estimated from model optimization with S as a free parameter.Bio-optical parameter starting values were adjusted to levels typical for relatively turbid, eutrophic lakes: C chl = 30 mg m −3 , C S P M = 5 g m −3 , a 440 nm CDOM = 5 m −1 [e.g 34].

Results
Quality control of above-surface Baltic Sea observations (section 2.1.3)removed 383 out of 3149 observations and 1 out of 45 stations.For each of the remaining stations R r s (λ) spectra were derived with 3C and L10 (section 2.2.2) and matched to corresponding sub-surface reflectances (section 2.2.3).If reflectances at a station could not be matched within physically reasonable boundaries for the Q-factor (3-8 sr), or spectral optimization failed (RSS > 1e − 4 sr −1 ), the station was discarded (3C: 5 stations, L10: 8 stations).Four stations are presented in detail in Figs. 4 (measurements) and 5 (modelled results).These stations were manually selected to illustrate results for a range of typical cloud cover conditions indicated by L s (λ) E d (λ) -ratios of varying shape, intensity, and variability.
1. Clear sky.This case is characterized by low variability in E d (λ) and L s (λ), and a L s (λ) E d (λ)ratio dominated by blue sky light at relatively low intensity ( L s E d (750 nm) 0.1sr −1 ), indicating few or no clouds.Low wind speed (5.4 m s −1 ) and low variability in L u (λ) and r s (λ) shows that surface-reflected radiance predominantly affected the blue part of the spectrum.Normalized root-mean-square-errors (nRMSE) of 4.67 ± 1.08 % and 16.32 ± 3.56 % between R r s (λ) and R r e f r s (λ) were obtained using 3C and L10, respectively.2. Haze.This station is characteristic of hazy skies and rough water surfaces due to higher wind speed (9.5 m s −1 ).Compared to the clear sky case, E d (λ), L s (λ), and L s (λ) E d (λ) -ratios were more variable, indicating an inhomogeneous sky radiance distribution.This, and the presence of longer wave periods could have caused considerable variability in L u (λ) E d (λ)ratios.L s (λ) E d (λ) was clearly dominated by blue light, yet at higher absolute intensities ( L s E d (750 nm) 0.1sr −1 ) compared to the clear sky case.This suggests considerable aerosol-scattered light contribution.nRMSEs of 7.97 ± 1.10 % and 16.72 ± 2.38 % for 3C and L10, respectively, are higher than for the clear sky case.

Scattered clouds. Inhomogeneous cloud cover conditions are indicated by spectrally neutral
L s (λ) E d (λ) -ratios close to π −1 at concurrent high variability in all three channels.The ratio L u (λ) was elevated towards the blue compared to R r e f r s (λ).This suggests that Rayleigh-scattered light that originated from outside the projected area of L s (λ) was reflected into L u (λ) observations, which is less likely to be corrected for with a spectrally flat offset ∆ in L10 (nRMSE 14.09 ± 1.02 %).Modelled diffuse radiance allowed 3C to compensate for this suboptimal measurement condition, which resulted in a residual error comparable to the clear sky case (nRMSE 4.96 ± 0.91 %).
4. Fully overcast.This case is representative of fully and homogeneously overcast conditions.
Figure 6 summarizes the 3C processing results for a daytime radiometric measurement cycle in two fixed directions, NNE (CH 1) and NNW (CH 2).The main source of variability in modelled R r s (λ) was due to variable scaling: spectrally averaged remote-sensing reflectances R r s (λ) varied slightly throughout the measurement period (Fig. 6, panel B), with relative standard deviations (RSD) of 6.6 % (CH 1) and 5.3 % (CH 2).The constant offset between channels might have been caused by residual dirt layers on cosine correctors, calibration errors, or adjacency effects.To compare the spectral shape of R r s (λ) over the course of the measurement period,  λ) , and the ratio λ) .Remote-sensing reflectance R r e f r s (λ) was derived from sub-surface reflectance (assuming a Q-factor of 5), which contains no surface reflections, and plotted in the last row for visual comparison with (λ) , which contains surface reflections.For further information on parameters and units, please refer to Table 1.
Variability of ρ ds within measurement cycles of typically less than 10 seconds was accounted for in the bio-optical model inversion (8.53 % and 9.84 % mean RSD for CH 1 and CH 2, respectively), and could not have been achieved by a L r (λ) correction that is based only on wind speed that is considered constant throughout a measurement cycle.Throughout the day, ρ ds varied by RSD 29.4 % (CH 1) and 36.5 % (CH 2), with clear maxima for when the sun was located directly behind (180 • relative azimuth) the respective radiance sensor pair.This can be observed for both channels and is likely due to Rayleigh scattering being highest in the direction  Fully overcast case (L10, #58) ∆ : 0. 001 ± 0. 00038 sr −1 NRMSE (Rrs, R ref rs ): 5.91±1.18% Fig. 5. Processing results for the four cases depicted in Fig. 4, processed with 3C (left panels) and L10 (right panels).Modelled (λ) (blue) to derive remote-sensing reflectances R r s (λ) (thick grey lines, error bars indicate standard deviations).Reference remote-sensing reflectances R r e f r s (λ) were derived from sub-surface reflectances R − (λ) and are depicted as black lines (error bars indicate standard deviations).For further information on parameters and units, please refer to Table 1.The dashed blue vertical line indicates ±135 • from the sun in the azimuthal plane, which is considered optimal for remote-sensing reflectance observations from above water.For further information on parameters and units, please refer to Table 1.

Discussion
A spectrally resolved correction factor ∆(λ) aids the processing of in situ hyperspectral radiometry to accurate R r s (λ).The 3C model can account for conditions where sky radiance measurements in the specular direction to the sea-viewing sensor do not exist, or are not representative of water surface-reflected radiance.Average sun glint contributions were generally small, which was at least in case of the Baltic Sea observations likely due to the sun glint-minimizing viewing geometry maintained by the Rflex system.However, sun glint contributions were highly variable, which highlights the importance of this correction component.The correction for diffuse sky radiance ρ ds was also typically small compared to the surface reflectance factor ρ s , indicating a small correction to the well-established methodology behind Eq. ( 3).Nevertheless, spectrally resolving this correction term with 3C consistently improved bio-optical model fits and R r s (λ) retrievals, regardless of sky condition, and as the authors of the L10 model suggested it should.A reliable correction for surface reflections could not be achieved with L10 under variably clouded sky conditions.Such conditions are rather common, which certainly contributes to typically high quality control rejection rates for automated shipborne R r s (λ) observations [e.g.10,11].
Remote-sensing reflectances were derived with 3C at comparable spectral accuracy for clear, mixed, and overcast sky conditions, which is concluded to be a major improvement over a scalar offset correction.Remote-sensing reflectance is an apparent optical property -changes in light field distribution have an impact on observed R r s (λ).The most obvious dependence is on θ sun (see e.g.Eq. ( 11)), which is why normalized radiances and reflectances ρ w (λ, θ sun = 0, θ view = 0) are often favoured in satellite remote sensing applications [e.g 1].For in situ observations, wave-induced and thus potentially short-lived changes in the Q-factor may have a scaling effect on R r s (λ).To focus on spectral differences rather than changes in scaling over the course of a day, derived remote-sensing reflectances were equalized to a common mean intensity in the presented daytime measurement cycle.Equalized remote-sensing reflectances from both channels in the fixed-angle measurements at Paterswoldsemeer matched closely and showed nearly no temporal variability, which could not be achieved with a scalar offset correction.These are strong indications that 3C accounted reliably for variable water surface reflections.Correction for diffuse sky light was highest for observations that were recorded with the sun exactly behind the respective radiance sensors, which is where Rayleigh scattering is most efficient.From this we conclude that a spectrally resolved offset correction for diffuse sky light, such as validated in this study, is crucial to account for the strong directional dependency of Rayleigh-scattered radiance.Sun glint contributions were also small but highly variable in this time series.Elevated values of ρ dd were observed in both channels when relative azimuth angles to the sun and sun elevation were lowest.Based on ray-tracing simulations, a relative viewing azimuth angle to the sun of approximately ∆θ = 135 • was recommended by Mobley [2,6] to minimize the chance of sun glint contamination.The present results empirically support that ∆θ = 135 • offers a good compromise between low intensity of Rayleigh-scattered diffuse sky light and a small chance of sun glint contamination.However, if this viewing geometry is not realisable, e.g. in airborne or satellite remote sensing, due to self-shading, or for fixed position instruments, 3C can compensate for sun glint and increased Rayleigh-scattered sky light reflections.
Reflectance derived from depth-profiles was used as a reference.Intense fluctuations of downwelling irradiance in water are induced by wave focusing effects [35,36,37] and low data quality of upwelling irradiance can be caused by instrument noise due to low light levels at depth.In order to account for these effects, multiple casts were averaged, relatively more measurements were taken near the surface, and depth ranges were manually adjusted.Both effects can nevertheless not be ruled out to have caused weak spectral distortion.Casts were recorded several meters away from the ship, on the sun-exposed side to avoid shading.Such close proximity to the ship can introduce adjacency effects, which would explain variable offsets but also affect the best fit of the Q-factor.Free-falling profilers that can operate at a greater distance to the ship avoid such contamination, but these are not particularly suitable for the highly absorbing waters of the Baltic Sea where the time taken by the profiler to collect a hyperspectral measurement would equate to a depth similar to the euphotic zone itself.Residual influence of foam and white-caps may be present in above-surface observations even after quality control, and is not accounted for in 3C.Variable ship pitch and roll can also alter the viewing geometry such that sky radiance is geometrically not representative of surface-reflected radiance.It is likely that these fundamentally different sampling methodologies have an impact on offsets δ between reflectances sampled above and below surface, such that they had to be accounted for in our statistical analysis, but bear no further interpretation.
The GC90 atmospheric model inherent to 3C is defined for cloudless maritime atmospheres.Despite the limitation of GC90 to clear sky conditions, 3C can to an extent also correct observations recorded under partially or fully overcast skies.This apparent contradiction can be resolved based on the Baltic Sea case 3: sky radiance was not spectrally representative of water surface-reflected light, probably because cloud occupied the complete sky radiance sensor field-of-view in an otherwise dominantly blue sky.Wind speeds close to 5 ms −1 roughened the water surface such that light from a large fraction of the sky (including blue parts) was reflected into the L u (λ) sensor field-of-view [see e.g. 2, Fig. 2].Since clouds appear spectrally neutral (white), at least in the visible spectral domain, the measurement of L s (λ) E d (λ) contained no spectral information that could be used to correct for Rayleigh-or aerosol-scattered radiance.The spectrally resolved offset parameter ∆(λ) in 3C can compensate for such observational shortcomings in sky radiance measurements with modelled radiances as observable under cloudless maritime atmospheres.However, compositions of aerosols may exist that cannot be parametrised accurately with the Ångström approximation of scattering in GC90, e.g.due to urban or industrial air pollutants, bush and forest fires, and sand storms.Under such conditions air mass type might require adjustment and 3C model optimization is more likely to produce a poor fit, which could not be verified with the data set analysed here.Average photon path lengths in the atmosphere can differ significantly for L r (λ) and E d (λ), especially when clouds are present.This effect introduces atmospheric absorption features to L u (λ) E d (λ) -ratios, which cannot be corrected for with 3C.However, such absorption features surface also in L s (λ) E d (λ) -ratios (e.g.Fig. 4, case 3 and 4) and can be used to correct for the effect in R r s (λ) [10] if L s (λ) happens to be representative of L r (λ), e.g. for a flat or very calm water surface.Variable atmospheric path length, cloud transmission, and a wider range of aerosol optical properties may be added to 3C to account for these effects, and further research in this direction is recommended.
Model optimization should not be expected to result in modelled L u (λ) E d (λ) -ratios that perfectly match observations, because this would require a finely tuned bio-optical and aerosol model and perfect instrument calibration.The aim is to realistically separate water leaving signal from residual signal that is attributed to L r (λ) [12].This separation is dependent on a suitable model parametrisation, e.g.realistic choice of specific inherent optical properties (sIOPs) and parameter ranges.Absorption by CDOM was recorded in the Baltic Sea data set and used to derive CDOM slope S, which is required in the bio-optical model inversion.When such measurements are lacking, S can be added to the inversion as a free fit parameter, on the condition that reasonable starting values and fit ranges can be provided, or be set to a realistic average value.The latter option even slightly improved the average nRMSE for Baltic Sea observations, and was chosen for processing the daytime measurement sequence.The bio-optical model by Albert and Mobley [25] is specified to simulate reflectances for a wide range of water types, and performed equally well for dark Baltic Sea waters and turbid inland lake reflectance in this study.In case of oceanic observations, a much simpler model could be applied to limit the number of fit parameters, and thus potential ambiguities in fit results.The choice of a suitable bio-optical model and its parametrisation for the water body in question remains with the user, as does the interpretation of fit residuals.Nevertheless, there is ample opportunity to determine the optical water type and associated starting values for the bio-optical model from the radiometric signals without first modelling R r s (λ) to perfection.Further validation studies should be carried out to define suitable protocols for a wider range of water types and sky conditions.

Fig. 1 .
Fig.1.Measurement geometry for water leaving radiance L w , upwelling radiance L u , and sky radiance L s (adapted from Simis and Olsson[10]).A minus sign (-) indicates sub-surface.Parameter dependencies were omitted for brevity and for further information on parameters and units, please refer to Table1.
with the Fresnel reflectance factor ρ f (θ view , n w ), which is a function of viewing zenith angle θ view and refractive index of water n w [e.g.1].The bio-optical model by Albert and Mobley[25] relates modelled remote-sensing reflectance R m r s (λ, a, b b , θ sun , θ view ) to absorption a(λ) and backscattering b b (λ) properties of water (w), chlorophyll-a (chl), CDOM, and suspended particulate matter (SPM): with water absorption a w (λ) and backscattering b b,w (λ), chlorophyll-a concentration C chl and concentration-specific absorption a * chl (λ), and SPM concentration C S P M and dry-weight specific absorption b * b,S P M (λ), CDOM absorption at 440 nm a 440 nm CDOM and CDOM absorption slope S. Proportionality factors f (λ, ω b , θ sun ) and f r s (λ, ω b , θ view , θ sun ) were taken from Albert and Mobley [25, Eqs. ( λ) using the analytical model proposed by Lee et al.[30, Eq. (24)]:

1 ] 1 ] 1 Fig. 4 .
Fig.4.Each column depicts observations that belong to one of four cases representing clear, hazy, scattered cloud, and fully overcast sky conditions.Corresponding station locations are indicated in Fig.2.In the first row, individual downwelling irradiance E d (λ) spectra are plotted together with their mean spectrum.In the same fashion, the following rows depict downwelling radiance L s (λ), upwelling radiance L u (λ), the ratio

Fig. 6 .
Fig. 6.A daytime measurement cycle of radiometric observations that was recorded 11 October 2015 from 10:30 to 16:00 (local time, sun zenith angle lower than 70 • ) in Paterswoldsemeer (The Netherlands) by an Ecowatch instrument that features two remotesensing reflectance R r s (λ) channels, facing NNE (CH 1) and NNW (CH 2).Panel A depicts equalized R r s (λ) that were derived with 3C for each channel over the course of the day, with standard deviations plotted as error bars.Panel B shows spectrally averaged L u (λ) E d (λ) -ratios (upper curves) and spectrally averaged R r s (λ) (lower curves) as a function of azimuth angle difference between sun and sensor.Panels C to F depict 3C model parameters as a function of azimuth angle difference: diffuse and direct reflectance coefficients ρ ds and ρ dd , respectively, Ångström exponent α, and turbidity coefficient β.Error bars in panel B to F indicate standard deviations over ten measurements that were recorded per cycle.The dashed blue vertical line indicates ±135 • from the sun in the azimuthal plane, which is considered optimal for remote-sensing reflectance observations from above water.For further information on parameters and units, please refer to Table1.

Table 2 .
Average distribution statistics (mean and standard deviation (std)) of 3C model optimization results (upper section), and regression statistics between derived and reference remote-sensing reflectances (lower section), for Baltic Sea stations recorded under the following sky conditions: clear skies ( L s E d (750 nm) < 0.1, 17 stations), mixed sky conditions (0.1 >= L s E d (750 nm) < 0.3, 11 stations), and overcast skies ( L s E d (750 nm) >= 0.3, 9 stations).For further information on parameters and units, please refer to Table1.

Table 3 .
Average distribution statistics (mean and standard deviation (std)) of L10 model optimization results (upper section), and regression statistics between derived and reference remote-sensing reflectances (lower section), for Baltic Sea stations recorded under the following sky conditions: clear skies ( L s E d (750 nm) < 0.1, 15 stations), mixed sky conditions (0.1 >= L s E d (750 nm) < 0.3, 9 stations), and overcast skies ( L s E d (750 nm) >= 0.3, 9 stations).For further information on parameters and units, please refer to Table1. .ρ dd resulted in very small absolute values (mean 1.62e − 12 for CH 1 and 4.50e − 12 for CH 2) that were highest for the smallest relative azimuth angles to the sun and corresponding low sun elevation.Extremely hight variability in ρ dd (164 % and 145 % mean RSD for CH 1 and CH 2, respectively) underlines the importance of correcting for sun glint even if the average contribution is small.Ångström exponent α from both channels resulted in similar values between approximately 12:00 and 14:30 hours (local time).High variability towards the afternoon in CH 2 indicates developing cloud formations or haze in that direction of the sky.This is supported by an increase in estimated turbidity coefficient towards the afternoon as observed by CH 2, but not by CH 1.