Assessment of Remote Sensing Reflectance Glint Correction Methods from Fixed Automated Above-Water Hyperspectral Radiometric Measurement in Highly Turbid Coastal Waters

1. Introduction

Remote Sensing Reflectance (R_rs(λ), sr^-1) is a crucial quantity in ocean color remote sensing tasks, which is necessary to develop water quality algorithms and evaluate the efficiency of the atmospheric correction algorithms [1,2,3]. In the field, R_rs(λ) is calculated as the ratio of water-leaving radiance (L_W(λ), mW m^-2 nm^-1 sr^-1) to downwelling irradiance (E_d, mW m^-2 nm^-1) just above the air-sea interface [4]. L_W(λ) contains information of inherent optical properties (IOPs) such as absorption (a(λ), m^-1) and backscattering (b_b(λ), m^-1) which are controlled by types and concentrations of Water Constituents Concentrations (WCCs) in the water column [5]. In-situ measurements of L_W(λ) are challenging due to specular reflections caused by surface-reflected light causes on the water surface [6]. The skylight-blocked approach (SBA) method has been utilized to measure L_W(λ) in the field. The SBA can measure L_W(λ) with high precision across different water types, but it is susceptible to uncertainties caused by self-shading, floating structure, and marine bio-foaling [6,7,8]. In contrast, L_W(λ) can be easily determined from above-water measurements using low-cost instruments and over long periods of time [9,10]. The above-water measurements include downwelling irradiance, upward radiance, and reflected sky radiance, which are measured simultaneously. In this way, the resulting L_W(λ) is influenced by sensor viewing angle (θ_V), sensor viewing azimuth (φ_V), solar zenith angle (θ_S), solar azimuth (φ_s), wind-induced surface roughness, and sky conditions (e.g., cloudiness) [4,11,12]. Both above and in-water measurements have their advantages and disadvantages for accurate estimations of R_rs(λ) [13,14]. In the following, we focus on above-water measurements.

Automated above-water hyperspectral radiometry in a fixed station is a cost-effective solution that provides long-timeseries of measurements. This type of measurement records the solar planar downwelling irradiance (E_d(λ), mW m^-2 nm^-1), and sky (L_S(λ), mW m^-2 nm^-1 sr^-1) and surface (L_T(λ), mW m^-2 nm^-1 sr^-1) radiances in short time steps (i.e., minutes or seconds) simultaneously. However, atmospheric conditions have a significant impact on E_d(λ) and L_S(λ), and water surface roughness affects the LT(λ). Both factors introduce considerable uncertainty to R_rs(λ). To overcome this problem, researchers have developed various methods, which can be estimated through bio-optical modelling [15,16] or through a similarity spectrum approach [17]. In this way, R_rs(λ) is calculated from simultaneous measurements of E_S, L_S, and L_T as follow [4,18]:

$$R_{rs}\left(\lambda ,{\theta }_v,{\Delta \phi }_v\right)={\displaystyle \frac{L_T(\lambda ,{\theta }_v,\Delta \phi )}{E_d(\lambda )}}-\rho \left({\theta }_v,\Delta \phi \right){\displaystyle \frac{L_S\left(\lambda ,{\pi -\theta }_v,\Delta \phi \right)}{E_d\left(\lambda \right)}}-\Delta L$$

(1)

where the geometry factor of (θ_v, Δφ) indicates the directions seen by radiometer. ρ(θ_v,Δφ) (dimension less, dl) is the effective surface reflectance that rescales the sky radiance to the reflected sea surface radiance. It is assumed that it is independent of the wavelength [4,19,20,21] or dependent on the wavelength [2,22,23]. ΔL (sr^-1) is the reflectance offset due to the total glint contribution and/or variabilities of E_d(λ) and L_S(λ) that are not consistent with L_T(λ). For flat surface, ρ(θ_v, φV) equals the Fresnel reflectance of the sea surface. Various approaches have been suggested for estimating ρ, including: fixed ρ values [20,24,25], sky-conditions dependent ρ values [26], lookup tables of ρ(θ_v,Δφ) for various directional factors and wind speed [4,27], and wavelength-dependent of ρ(λ,θ_v,Δφ) for various directional factors, wind speeds, and aerosol loads [22]. However, a sensor geometry with θ_v = 40° and φ = 90°−135° is considered optimal to minimize glint perturbations [4], so in practice these geometries are usually fixed. Similarly, researches have considered different approaches to estimate the ΔL estimation, including a fixed value equal to the R_rs(λ) at Near Infra-Red (NIR) [28,29,30,31], values of the similarity spectrum at NIR [17], the difference between bio-optical models and simulated R_rs(λ) [15], and estimating the direct surface reflected radiance to simulate ΔL [19,21]. The last approach substitutes a wavelength-dependent correction factor of ΔL(λ) in Equation (1). Since automated above-water radiometry is subject to unavoidable varying environmental conditions, application of above-water radiometric measurements introduces additional challenges. Consequently, there is a need for a thorough review and assessment of methods because the estimates ρ and ΔL may not be accurate enough for long timeseries of measurements in all environmental conditions.

The analytical methods for deriving IOPs from R_rs(λ) comprise inverse modeling that cannot be exactly reduced to an analytical equation or expression and relies on empirical and/or semi-empirical relationships [32]. However, The Radiative Transfer (RT) simulations and the quasi-single scattering approximation (QSSA) forward models have been developed to estimate R_rs(λ) as a function of spectral IOPs [5]. Several forward QSSA analytical expressions have been developed that estimate the R_rs(λ) as a function of a(λ), m^-1) and b_b(λ) [5,33,34]. Similarly, RT models simulate R_rs(λ) as a function of user-input IOPs and environmental parameters [4,33,35]. Therefore, RT and QSSA forward models explicitly relate IOPs to R_rs(λ) free of sun and sky glint. The resulting R_rs(λ) can be used as a reference to evaluate the above-water R_rs(λ) from Equation (1). In this regard, the lack of IOPs is a critical factor that significantly impacts the outcomes of this study. The validated regional bio-optical models can be utilized to overcome the lack of IOPs data [5].

This research aims to evaluate the performance of various sun, sky, and surface glint correction models from a fixed automated above-water radiometric measurements in highly turbid coastal waters. We will assess the strengths and limitations of the proposed models in various environmental conditions. Finally, we provide two showcases to further confirm the significance of the total glint correction models for above-water radiometry.

2. Field Data

2.1. Above-Water Hyperspectral Radiometric Data

We utilized a time-series of above-water hyperspectral measurements spanning from 2006-2023 at the Royal Netherlands Institute for Sea Research (NIOZ) jetty station (NJS) located near the Marsdiep tidal inlet of the Dutch Wadden Sea (Figure 1a). NJS has been equipped with a set of 6 TRIOS-RAMSES radiometers, comprising 1 Ramses-ACC for measuring solar downwelling irradiance (E_d, mW m^-2 nm^-1), 1 Ramses ACC-UV for measuring solar downwelling irradiance in the ultraviolet region (E_S-UV, mW m^-2 nm^-1), 2 pairs of Ramses-ARC for measuring sky radiance (L_S, mW m^-2 nm^-1 sr^-1), and total upwelling radiance from water surface (L_T, mW m^-2 nm^-1 sr^-1). The spectral data was recorded every 10 minutes at 1 nm intervals, ranging from 350 nm to 950 nm. The ARC sensors were set up with viewing angle θ_V = 35° and direction φ_V = 135° (south-east) and φ_V = 225° (south-west) (Figure 1b–d). A set of quality control tests were applied to the hyperspectral data to flag the invalid or suspicious spectra [36]. The quality-controlled flagged spectra were excluded from further analysis. In this study, we used the quality-controlled E_d(λ), L_S(λ), and L_T(λ) covering 6:00 AM to 6:00 PM (UTC).

2.2. Meteorological Data

Hourly meteorological observations from the De Kooy station (Figure 1), located approximately 8 km southwest of NJS, covering the period from 2006-2023, were obtained from The Royal Netherlands Meteorological Institute (KNMI) database (https://dataplatform.knmi.nl). Measurements included cloud cover (CC, coverage of the upper sky in eighths, with ‵9‵ indicating the upper sky invisible), atmospheric pressure (hPa), relative humidity (%), and wind speed (m s^-1).

2.3. IOPs Data

Reference [25] measured IOPs at 156 stations across the entire Wadden Sea between 2006 and 2007. Of these, 22 stations collected at NJS during March, June, and July 2006, were used in this study (hereafter referred to HA09). The HA09’s IOPs dataset included laboratory spectrophotometry of total absorption (a(λ), m^-1), Colored Dissolved Organic Matter (CDOM) absorption (a_CDOM(λ), m^-1), particle absorption (ap(λ), m^-1), Non-Algal Particle (NAP) absorption (a_NAP(λ), m^-1), and phytoplankton absorption (a_Phy(λ), m^-1). They have measured particle absorption with a spectrophotometer (Ocean Optics) using the filter pad method. In addition, beam attenuation (c(λ), m^-1) and total absorption were measured in-field using an AC9 (Wet Lab, 25 cm cuvette, linearity range 0.001-30 m^-1), hereafter c_AC9(λ) and a_AC9(λ). The AC9 operated at nine standard wavelengths (412, 440, 488, 510, 555, 630, 650, 676, and 750 nm) ([25], their table 2). Both in-situ and ex-situ measurements were conducted based on standard protocols [37], corrected for blank measurements on MilliQ water, and corrected for temperature-salinity and scattering [38], respectively. Only the AC9 data at depth ≤ 1.0 m were utilized, along with spectrophotometer data. The particle scattering was calculated as b_p(λ) = c_AC9(λ) - a_AC9(λ) for in-situ AC9 data. The light backscattering coefficients of particles were determined as follow [35,39]:

$$b_{bp}\left(\lambda \right)=b_p(\lambda )\times B$$

(2)

where B is the backscattering probability (or backscattering ratio), assumed to be spectrally constant [40,41]. The values of B ranged between 0.024 to 0.042 for surface waters of the coastal area of the English channel [42]. The value of B has been assumed to be 0.026 in the Wadden Sea (Please see supplementary, Section S1).

The HA09 dataset also included supplementary in-situ data of wind speed (m s^-1), water temperature (C°), Ph, salinity (PSU), Secchi Disk Depth (SDD, m^-1), and wave height (m).

2.4. Concentrations of Chla and SPM

A total of 1015 surface water samples were collected at NJS during 2006-2020 for laboratory Chla and Suspended Particle Matter (SPM, g m^-3) measurements. Of these, 648 samples included measurements of both Chla and Total Suspended Matter (TSM, g m^-3). The Chla concentrations of the NJS dataset included two groups of High-Performance Liquid Chromatography (HPLC) and Fluorometry measurements. The Chla concentrations of HA09 and NJS-HPLC samples were measured using standard HPLC measurements on samples filtered through GF/F filters then frozen and transported to laboratories for analysis. The Chla concentrations of NJS-Fluorometry samples were measured using a laboratory Fluorometer on filtered samples using standard methods by 43]. TSM concentrations were measured using the standard filtration method described by [37] using pre-ashed, rinsed, and weighted 47 mm GF/F filters. Readers are referred to [25,24] for more details on these measurements.

3. Methods

Since the temporal coverage of measured IOPs was limited to HA09 datasets, the bio-optical models were used to simulate the IOPs using the average of HA09 dataset and measured WCCs at NJS. Subsequently, the R_rs(λ) from above-water radiometric measurements were evaluated against glint-free R_rs(λ) simulated using IOP-to-R_rs(λ) forward models.

3.1. Bio-Optical Models

IOPs are expressed as the sum of absorption or scattering by water molecules and various dissolved and particulate water components [44]. CDOM, NAP, and phytoplankton are the representative of the natural water particles with distinct absorption and backscattering. The total absorption (a(λ), m^-1) and backscattering (b_b(λ), m^-1) coefficients of water can be expressed as follows [1,45]:

$$a\left(\lambda \right)=a_w\left(\lambda \right)+a_{Phy}\left(\lambda \right)+a_{CDOM}\left(\lambda \right)+a_{NAP}(\lambda )$$

(3)

$$b_b\left(\lambda \right)=b_{b,w}\left(\lambda \right)+b_{b,Phy}\left(\lambda \right)+b_{b,NAP}(\lambda )$$

(4)

where the subscripts w, Phy, CDOM, and NAP present contributions by pure water, phytoplankton, CDOM, and NAP constituents. Since CDOM backscattering is negligible, it is left out in Equation (4). The total non-water absorption (a_nw, m^-1) and particulate backscattering (b_bp, m^-1) are expressed as the sum of the non-water constituents: $a_{nw}\left(\lambda \right)=a_{Phy}\left(\lambda \right)+a_{CDOM}\left(\lambda \right)+a_{NAP}(\lambda )$, $b_{bp}\left(\lambda \right)=b_{b,Phy}\left(\lambda \right)+b_{b,NAP}(\lambda )$. Table 1 summarizes the algorithms for parameterizing absorption and backscattering of phytoplankton, CDOM, and NAP in coastal waters.

In this study, the absorption coefficients of water constituents were calculated (Eqs. 5-9) and the results were validated using the HA09 field measurements. The spectral slope of CDOM (S_CDOM, nm^-1) and NAP (S_NAP, nm^-1) were calculated by fitting the exponent of a_CDOM between 350-500 nm and 380-730 nm, respectively [54]. To calculate SNAP, the a_NAP(λ) in the 400-480 nm and 620-710 nm ranges were excluded to avoid pigment absorption residuals [54,55]. The absorption of pure water, a_w(λ), was taken from [56,57] in the ranges of 350-550 nm and 555-700 nm, respectively. According to Table 1, a_nw(λ) was calculated in two ways, one using the sum of Eqs. 6, 8, and 9 (a_nw,BR19x(λ) hereafter) and the other using the sum of Eqs. 7, 8, and 9 (referred to a_nw,LC99(λ) hereafter).

Reference [24,28] utilized the parametrization of b_b(λ) using Eqs. 10 and 12 to retrieve the WCC at the NJS. They showed that the measured SIOPs by [25] were valid for WCC retrieval. Therefore, we used the same parametrization for the backscattering model (b_b,BA18(λ) hereafter). In addition, we tested the parametrization of b_b(λ) using Equations (11) and (13) (b_b,GE05(λ) hereafter). The backscattering coefficient of pure water (b_b,w(λ)) was obtained from [39].

The Specific IOPs (SIOPs) (denoted by * in Table 1) were calculated as:

a*_Phy(λ) = aPhy(λ)/Chla (m² mg^-1), b*_SPM(λ) = b_bp(λ)/SPM (m² g^-1) , and a*_NAP(λ) = a_NAP(λ)/SPM (m² g^-1). Table 2 shows the statistical ranges of WCCs, absorption, and backscattering used for bio-optical model parametrization.

3.2. Generating the Glint-Free R_rs(λ)

The bio-optical models were used to obtain glint-free R_rs(λ). The bio-optical model could be represented as follows [58,59,60]:

$$R_{rs}\left(\lambda \right)={\displaystyle \frac{\alpha .r_{rs}\left(\lambda \right)}{1-\beta .r_{rs}\left(\lambda \right)}}$$

(14)

$$r_{rs}\left(\lambda \right)=g_w·u_w+\sum _{k=1}^4{g}_k{u}_p^k\left(\lambda \right)$$

(15)

$$u_w\left(\lambda \right)={\displaystyle \frac{b_{b,w}\left(\lambda \right)}{b_b\left(\lambda \right)+a\left(\lambda \right)}}$$

(16)

$$u_p\left(\lambda \right)={\displaystyle \frac{b_{b,p}\left(\lambda \right)}{b_b\left(\lambda \right)+a\left(\lambda \right)}}$$

(17)

where α and β are constant values retrieved from a lookup table based on θv, Δφ, θs, wind speed, and cloud cover presented in [35]. The gw and gk are the model parameters dependent on θ_S and sensor geometry, retrieved from a lookup table presented by [59]. This model has been used for high sediment load waters [59,61]. Equation (14) model was used to generate the glint-free R_rs(λ) spectra (R_rs,ref(λ) hereafter) using the average IOPs and measured WCCs.

3.3. Methods of ρ and ΔL Estimation from Above-Water Radiometry

In this study, the determination of ρ(θv,Δφ) and ΔL have been considered using methods presented in Table 2. Table 3 summarizes the methods of estimating the ρ(θv,Δφ) and ΔL factors. Regarding the ρ, the models were categorized into 4 groups, including: i) fixed values of ρ(θv,Δφ), ii) lookup tables for wind speed and sun-sensor geometries, iii) wavelength-dependent approaches and dependency of Aerosol Optical Thickness (AOT), wind speed, and sun-sensor geometries (HT23 and ZX17), and iv) using direct measurements of radiometric data. The Ku13 and JD20 models are not sensitive to the ρ(θv,Δφ) value, as they primarily focus on the shape of the water surface reflectance spectra. The ZX17 presents a lookup table of simulated ρ(λ) as a function of AOT, θs ≤ 60°, wind speed, and sun-sensor geometries in clear skies.

Since environmental conditions were not specifically considered in the models, we evaluated the performance of selected models in various environmental conditions, including sky cloudiness and AOT loads. The viewing geometry (θv,Δφ) and wavelength (λ) were omitted for brevity hereafter. In addition, the quality of R_rs(λ) was used to evaluate the results of models in different conditions using the Quality Assurance (QA) scores provided by [62]. The QA scores ranged from 0 (invalid) to 1 (perfect).

3.4. Identification of Environmental Factors

The environmental factors used in the previous sections were obtained as follows:

3.4.1. Sun Azimuth and Zenith Angle

The θ_S and φ_S were calculated for the date-time of observations using the method presented in the “Climate Data Toolbox” [63].

3.4.2. Aerosol Optical Thickness (AOT)

A spectrum of Ed(λ) is a sum of the direct and diffusive components of downwelling irradiance, which can be calculated based on the Gregg and Carder model (GC90 hereafter) [64]. The GC90 model has been developed primarily for cloudless maritime atmospheres. The RADTRAN solar irradiance model [65] has been developed based on the GC90 along with the use of spectral cloud transmission based on the delta-Eddington approximation of two-stream approach flowing Slingo [66]. An R-based implementation of this model is available at https://rdrr.io/github/sean-rohan-NOAA/trawllight/man/radtran.html. In this study, the direct and diffusive components of solar irradiance were calculated using the inverse RADTRAN model of E_d(λ) using the parametrization provided by [16] and corresponding meteorological data. Subsequently, AOT was calculated using the E_d(λ), direct and diffuse solar irradiance following the methods presented by [67].

The AOT(550) values were in the range of 0.01 to 1.2. In general, AOT(550) was found to be approximately 0.1 for clear atmospheres, 0.4 for moderately turbid atmospheres, and greater than 1.0 for very turbid atmospheres [22]. Here, we considered the AOT(550) ≤ 0.1 as ‘low AOT’, and AOT(550) ≥ 0.5 as ‘high AOT’ atmospheric conditions.

3.4.3. Sky Conditions (Clear, Scattered Clouds, or Overcast)

We have developed a method for identifying sky conditions as a function of cloudiness using hyperspectral observations of E_d(λ) and L_S(λ). This methods is based on the work of [68] that utilizes the spectral solar irradiance model developed by [64] for cloudless maritime atmospheres. Specific criteria were established to identify cloud presence based on time series data of E_d(λ) at 550 nm. This model differentiated cloudless skies, scattered clouds, and overcast conditions using RT simulations. The details of this method is available in [36]. Here, the hyperspectral data were collected over a broad range of environmental conditions, encompassing clear (11.8% of observations) and cloudy skies, low to high sun zenith angles (28-89°), and low to high wind speeds (1-14 m s^-1).

3.5. Statistical Metrics

The statistical parameters of standard deviations (std), normalized Root Mean Square Error (NRMSE), coefficient of determination (R2), Mean Absolute Percentage Error (MAPE), Unbiased Percentage Difference (UPD), normalized mean difference or bias (NBias), and mean ratio (MR) were used to quantify the consistency between models and observations.

$$NRMSE={\displaystyle \frac{\sqrt{\overline{{(x_{mod}-x_{obs})}^2}}}{\overline{x_{obs}}}}$$

(18)

$$MAPE={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{|x_{i,mod}-x_{i,obs}|}{x_{i,obs}}}\times 100$$

(19)

$$UPD=100\times {\displaystyle \frac{x_{mod}-x_{obs}}{x_{mod}+x_{obs}}}$$

(20)

$$NBias={\displaystyle \frac{\sum _{i=1}^n|x_{i,mod}-x_{i,obs}|}{\sum _{i=1}^n{x}_{i,obs}}}$$

(21)

$$MR={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{x_{i,mod}}{x_{i,obs}}}$$

(22)

where x_mod and x_obs represent the estimated and observed values, n is the number of observations. A statistical score was defined based on the values of NRMSE, MAPE, UPD, and R2 as follow [69]:

$$Total Score={(W}_{MAPE}\times {MAPE}_{nor})+{(W}_{RMSE}\times {RMSE}_{nor})+{(W}_{UPD}\times {UPD}_{nor})+{(W}_{R2}\times {R2}_{nor})$$

(23)

where the subscript ‘nor’ indicates statistical normalization = [x – min(x)]/[max(x) – min(x), and W_x is the weight of each parameter, assigned an equal value of 0.25 for each parameter. The Total Score varies between 0 to 1, with higher values for higher statistical consistency between observed and modeled parameters.

4. Results

4.1. Parametrization and Validation of Bio-Optical Models

The bio-optical models are parametrized using the average values of measured a*_Chl(λ), a_CDOM(λ2), a*_NAP(λ), S_CDOM, S_NAP, b*_SPM(λ), and WCCs (Table 2). The normalized a_Phy(λ) shows low variability (CV = 0.06 at 410 nm and CV = 0.03 at 675 nm) indicating phytoplankton remain relatively consistent during the measurements. The normalized a_CDOM(λ) and a_NAP(λ) spectra for the entire dataset show a minimal difference across the spectral structure. As a result, the mean values of a*_Phy(λ), a_CDOM(λ2), a_NAP(λ), S_CDOM and S_NAP are used for parametrization of the absorption models in Eqs. 6, 8, and 9 (Table 3). Moreover, the mean values of b*_SPM(λ) = b_SPM(λ) / [SPM] m² g^-1 are used for the parametrization of the NAP backscattering coefficient (Equation (12)). Figure 2 shows the correlations between the measured and estimated absorption and backscattering coefficients. The values of MAPE, NRMSE, NBias, MR, and R² of a_Phy(λ) from Equation (6)\(7) against in-situ measurements are 15.0\24.2%, 0.15\0.26, 0.02\0.03, 0.88\0.75, and 0.91\0.84, respectively. The measured a_Phy(λ) spectra exhibit more consistency with the spectra calculated from Equation (6) than Equation (7). The exponential fits of a_CDOM(410) and a_NAP(410) show very good consistency with measurements (MAPE ≤ 3.6%, NRMSE ≤ 0.06, NBias ≤ 0.03, MR ≅ 1.0, and R² ≥ 0.95). The correlation between in-situ a_nw(440) and a_nw,BR19x(440) shows more consistency with measurements than a_nw,LC99(440) (MAPE = 4.3%, NRMSE = 0.11, NBias = 0.04, MR = 0.98, and R² = 0.86). Therefore, parametrization of Equation (3) using a_p,BR18x(λ) (Equation (6) + Equation (8) + Equation (9)) represents the optimum setup that correlates well with the in-situ measurements. Parametrization of b_bp(λ) using b_b,GE05(440) components show more consistency than those by b_b,BA18(440), with the values of MAPE = 18.8 & 23.8%, NRMSE = 0.12 & 0.16, NBias = 0.03 & 0.04, MR = 0.89 & 0.83, and R² = 0.92 & 0.87, respectively. All estimated IOPs fall within the boundaries of the 95% confidence level of the measurements. The findings indicate that the HA09 average SIOPs are valid for bio-optical model applications.

4.2. Evaluation of ρ and ΔL Estimation Methods

The R_rs,ref(λ) are calculated using Equation (14) utilizing the average of SIOPs (Table 2) and measured WCCs for the time periods when in-situ measurements are available. The R_rs(λ) spectra corresponding to the R_rs,ref(λ) are calculated using the quality-controlled E_d(λ), L_S(λ), and L_T(λ) spectra from south-east and south-west looking sensors base on the proposed methods outlined in Table 3. Figure 3 shows the average of modeled R_rs(λ) and R_rs,ref(λ) spectra in various environmental conditions. The corresponding error estimates of each case in terms of statistical Total Score (Equation (21)) between R_rs(λ) and R_rs,ref(λ) are shown in Figure 4. The full set of statistical parameters comparing R_rs(λ) and R_rs,ref(λ), and associated scatter plots are shown in the Supplementary Material Figures S1–S4. The environmental conditions are defined as follows.

Case I, high Δφ: Clear skies, θ_S ≤ 45°, wind speed ≤ 3 m s^-1, low AOT, and Δφ ≥ 60° (Figure 3a and Figure 4a). (N = 65). This case is characterized by low variability of R_rs(λ) spectra. The lowest deviations from R_rs,ref(λ) are achieved by 3C model (Total Score = 0.90, MAPE = 8.8%, UPD = -2.4, NRMSE = 16.1%) in the range of 400-720±5 nm, and the heist deviations by the Ru05 (Total Score = 0.76, MAPE = 32.9%, UPD = -9.7, nRMSE = 24.7%). Similar results are found in the range of λ=440±5 nm (blue), λ=560±5 nm (green), λ=680±5 nm (red), and λ=720±5 nm (NIR). The other methods show intermediate and relatively similar error estimates against Rrs,ref(λ), with an average of Total Score = 0.79, MAPE = 19-28%, UPD = -8-7, and nRMSE = 24-28%.

Case II, low Δφ: Clear skies, θ_S ≤ 45°, wind speed ≤ 3 m s^-1, low AOT, and Δφ ≤ 40° (Figure 3b and Figure 4b). (N = 83). This case is characterized by high deviations of the models that utilize lookup tables of ρ and BA18. The 3C model shows the best score with Total score = 0.91, MAPE = 8.7%, UPD = 2.9, and NRMSE = 15.8% in the range of 400-720±5 nm. Similar statistical results are observed for the blue, green, red, and NIR regions. The worst statistical scores are observed for MO99, MO15, BA18, HT23, and ZX17 methods with Total Score = 0.56, MAPE > 100%, UPD = -18, and NRMSE = 14-25%.

Case III, Overcast sky condition: Overcasts, 45° ≤ θ_S ≤ 65°, wind speed ≤ 3 m s^-1, low AOT, and 40° < Δφ < 60° (Figure 3c and Figure 4c). (N = 115). The MO99, MO15, BA18, and HT23 show a significant overestimations against R_rs,ref(λ) across all wavelengths, with Total Score = 0.54, MAPE = 91-97%, UPD = -14.6- -14.8, NRMSE = 20-24%. The lowest deviation belongs to the 3C model with Total Score = 0.93, MAPE = 7.4%, UPD = 1.5, and NRMSE = 11.8%. The visual comparison shows that the Ku13 and JD20 models underestimate significantly R_rs(λ) in the blue and NIR regions. The ZX17 is not available for this case.

Case IV, scattered cloud condition: Scattered clouds, 45° ≤ θ_S ≤ 65°, wind speed ≤ 3 m s^-1, low AOT, and 40° < Δφ < 60° (Figure 3d and Figure 4d). (N = 307). This case is characterized by high deviations and overestimations of MO99, MO15, BA18, and HT23 methods with R_rs(560) = 0.022-0.027 sr-1 and Total Score = 0.59, MAPE = 85-93%, UPD = -12.5 – 15.2, and NRMSE = 9.1-10.2%. The best simulation belongs to the 3C model with a Total Score of 0.92. The ZX17 is not available for this case.

Case V, strong wind speed: Clear skies, 45° ≤ θ_S ≤ 65°, wind speed ≥ 9 m s^-1, low AOT, and 40° ≤ Δφ ≤ 60° (Figure 3e and Figure 4e). (N = 69). Visual comparison of R_rs(λ) spectra indicates the relatively high variability and underestimations in the blue wavelength regions by JD20 and Ku13 methods, which utilize spectral shape for corrections. The 3C & Ru05 show the best performance with Total Score = 0.90 & 0.79, MAPE = 8.8% & 17.3%, UPD = 1.4 & -6.3, and NRMSE = 17.3% & 24.3%, respectively.

Case VI, high wind speed: Clear skies, 45° ≤ θ_S ≤ 60°, 4 ≤ wind speed ≤ 7 m s^-1, low AOT, and 40° ≤ Δφ ≤ 60° (Figure 3f and Figure 4f). (N = 211). This case is characterized by overestimation of the MO99, MO15, BA18, HT23, and ZX17 methods with an average of ΔR_rs(560) = 0.009 sr^-1 relative to R_rs,ref(λ), with Total Score = 0.58, MAPE = 92-98%, UPD = -16.8 - -16.1, and NRMSE = 9.8-11.2%. The 3C model shows the best performance with Total Score = 0.93, MAPE = 6.6%, UPD = 1.2, and NRMSE = 8.8%. The JD20 and Ru05 show the best performance next to the 3C model. The HT23 and ZX17 show similar result.

Case VII, high AOT: Clear skies, 45° ≤ θ_S ≤ 60°, wind speed ≤ 3 m s^-1, high AOT, and 40° ≤ Δφ ≤ 60° (Figure 3g and Figure 4g). (N = 32). This case is characterized by low variability of R_rs(λ) spectra. All models except Ru05 underestimate the R_rs(λ) spectra with an average of ΔR_rs(560) = 0.004 sr^-1 relative to the R_rs,ref(λ) spectra in the blue and green regions. The Ku13 and JD20 methods show the highest deviations in the blue and green regions (UPD = 11.5-18.8, and NRMSE = 37.6-60.1%). The HT23 and ZX17 show a good performance in the green-red regions (Total Score = 0.89, MAPE = 9.5%) and a relatively weaker performance in the blue region (Total Score = 0.75, MAPE = 29.3%). Overall, the 3C model show the best performance, with Total Score = 0.89, MAPE = 9.3%, UPD = 4.6, and NRMSE = 19.1%.

Case VIII, high sun-zenith angle: Clear skies, θ_S ≥ 80°, wind speed ≤ 3 m s^-1, low AOT, and 40° ≤ Δφ ≤ 60° (Figure 3h and Figure 4h). (N = 11). All methods show low variabilities. Apart from the 3C, the other methods underestimate the R_rs(λ) spectra in the blue-NIR regions. The MO99, MO15, Ru05, BA18, HT23, and JD20 show relatively similar performance with Total Score = 0.81-0.84, MAPE = 11.6-15.1%, UPD = -1.2- 4.6, and NRMSE = 30.1-32.4%. The Ku13 shows a very weak performance in the blue region with a Total Score of 0.36. The 3C model shows the best performance in this case with Total Score = 0.92, MAPE = 5.2%, UPD = 0.92, and NRMSE = 22.7%. The ZX17 is not available for this case.

4.3. Variability of ρ and ΔL

Figure 5a shows the variability of ρ values for the eight Cases. The ρ values are not highly sensitive to the environmental conditions and are in the range of 0.0271±0.0007 for MO99, 0.0309±0.0008 for MO15, 0.0265±0.0002 for Ru05, 0.03115±0.0006, and 0.0296±0.0008 for 3C. The 3C model exhibits the highest ρ = 0.0308 for Case-I and the lowest value of ρ = 0.0286 for Case-VI. Ku13’s fixed ρ of 0.02 significantly differs from other models. The fixed ρ = 0.0265 of BA18 indicates the minimum value of ρ in various conditions and is very close to the ρ values of Ru05. The Ku13 and JD20 models are not sensitive to the ρ value, as they primarily focus on the reflectance spectral shapes. The HT23 and ZX17 models demonstrate a wavelength-dependent in ρ values, decreasing from blue to NIR by an average value of 0.002 in all conditions. While the HT23 and ZX17 have been simulated based on similar environmental conditions (e.g., wind speed and AOT) the ρ(λ) values from these models are significantly different, with average of ρ(560) = 0.0296 and ρ(560) = 0.0252, respectively.

The models, which result in a wavelength-independent ΔL, exhibit low (≤0.005 sr^-1) and relatively uniform ΔL variations in all conditions (Figure 5b). In contrast, the models with wavelength-dependent ΔL (i.e., 3C and Ku13) exhibit higher and more variable ΔL in various conditions. The 3C model shows a wavelength-dependent ΔL with a minimum of 0.0080 sr^-1 in Case-VIII and a maximum of 0.0102 sr-1 in Case-IV. The values of 3C-derived ΔL at the blue to NIR regions are very close to the average values (std ≤ 0.001 sr^-1) in various conditions. The Ku13 model presents the wavelength-dependent ΔL with the highest values of 0.0091 sr^-1 for Case-VIII and 0.0141 sr^-1 for Case-II. The ΔL values of Ku13 decrease from blue to NIR regions by an average of 0.0021 sr^-1 in all conditions.

4.4. QA scores of Simulated Above-Water R_rs(λ)

Figure 6 displays the QA scores for the R_rs,ref(λ) and R_rs(λ) spectra. For the total R_rs,ref(λ) spectra, 55% of spectra achieve QA scores ≥ 0.9, and 88% QA scores ≥ 0.8. The minimum values of QA scores for R_rs,ref(λ) spectra range from 0.6 to 0.7, accounting for 6% of total spectra. The total spectra of the MO99, MO15, and BA18 models show relatively similar QA scores with an average of 25% of the spectra scoring greater than 0.9 and 24% scoring less than 0.6 (bad/unusable spectra). The total spectra of JD20 model show the lowest percentage of spectra with QA score greater than 0.9 (17%) and 23% of bad/unusable spectra (QA ≤ 0.6). Overall, significant portions of the total Rrs(λ) spectra using MO99 (25%), MO15 (27%), Ru05 (18%), BA18 (22%), HT23 (24%), Ku13 (22%), JD20 (32%), and ZX17 (21%) are suspicious or unusable (QA ≤ 0.5). The 3C model shows the best match of R_rs(λ) QA scores with the R_rs,ref(λ) spectra. The total spectra of the 3C model account for 77% of R_rs(λ) with QA scores greater than 0.8 and only 2.3% of spectra have QA scores less than 0.6.

For Case-I, 71% of R_rs,ref(λ), 74% of 3C, and 48-60% of MO99, MO15, Ru05, BA18, HT23, and ZX17 have QA scores ≥ 0.9. The Ku13 and JD20 models have 34% and 40% of spectra with scores ≥ 0.9. About 29% of JD20 spectra are unusable or suspicious (QA ≤ 0.6). The MO15 and BA18 have 9% and 6% of spectra with QA scores ≤ 0.6. For Case-II, III, V, and VI, the R_rs,ref(λ) have QA ≥ 0.7 which shows a relatively similar distribution of QA scores with the 3C model. In these conditions, less than 2% of 3C spectra are unusable/suspicious (QA ≤ 0.6), about 16%-28% of JD20 spectra have QA ≤ 0.6, and about 16%-22% of the spectrum from the other models are unusable/suspicious. For Case-IV, less than 3% of R_rs,ref(λ) spectra have a bad QA score. About 30% of JD20 spectra show a bad score, and about 14%-22% of the spectrum from the other models have QA ≤ 0.6. The 3C model shows the closest distribution of QA to the R_rs,ref(λ) with only 2.3% of spectra having QA ≤ 0.6. For Case-VII, low percentages of QA scores ≥ 0.9 are observed with 35% for R_rs,ref(λ) and 3C spectra, and less than 10% for the other models. In this case, more than 30% of all modeled spectra, except 3C, are unusable/suspicious. The lowest percentages of high-quality spectra (QA ≥ 0.9) are observed in Case-VIII, with 17% of R_rs,ref(λ), 33% of 3C, 17% of JD20, and zero for the other models. For Case-VIII, 50% of JD20 and Ku13, 34% of HT23, 17% of MO99, MO15, and BA18 are bad or suspicious spectra. Overall, the 3C spectra show the highest QA scores, with a distribution similar to that of R_rs,ref(λ), and less than 3% of bad/suspicious spectra under all conditions.

4.5. Showcases of R_rs(λ) Models

Two showcases are presented to illustrate the performance of R_rs(λ) models for Chla and TSM retrieval, respectively. We compare in-situ Chla and TSM with estimated values retrieved from glint-corrected R_rs(λ) based on methods in Table 3. An empirical algorithm developed by NIOZ [70] is used to estimate WCCs from the above-water radiometric measurements. High-quality R_rs(λ) spectra with QA ≥ 0.8 (Section 4.4) are selected for estimation of Chla and TSM. Figure 7 shows the correlations between in-situ and estimated Chla using different methods of R_rs(λ). The R_rs(λ), which utilize lookup tables of viewing direction and winds to estimate ρ (i.e., MO99, MO15, HT23, and Ru05), show relatively similar error statistical metrics, with an average Total Score = 0.63, MAPE = 44%, UPD = 21.5%, NRMSE = 0.36, and R² = 0.65. The ZX17 shows a better performance with Total Score = 0.74 and R2 = 0.84. The methods that utilize a fixed value of ρ (i.e., BA18, Ku13, and JD20) show more consistency with in-situ Chla, with an average Total Score = 0.67, MAPE = 53%, UPD = 18%, NRMSE = 0.37, and R² = 0.73. In contrast, the estimated Chla using the 3C model shows the best results, with Total Score = 0.81, MAPE = 27%, UPD = 8%, NRMSE = 0.22, and R² = 0.89. Figure 8 shows the second showcase of in-situ TSM vs. estimated TSM using different methods of R_rs(λ). The MO99, MO15, Ku13, BA18, and JD20 models show relatively similar results with an average Total Score = 0.69, MAPE = 21%, UPD = 7.8, NRMSE = 0.11, and R² = 0.73. The ZX17 shows a good performance with Total Score = 0.82 and R² = 0.83. The 3C model provides the better performance with MAPE = 15%, UPD = 3.25, NRMSE = 0.09, and R² = 0.79. The 3C model suggests that ~3% of Chla and TSM data falls outside the 95% confidence level, whereas other models indicate ~6-9% of data outside this level. Overall, there is no significant distinction between the models that utilize look-up tables and those that employ fixed ρ values for the estimation of Chla and TSM, with the exception of ZX17, which demonstrates better performance under clear sky condition.

5. Discussion

This study reviews various above-water R_rs(λ) glint-correction methods that use different approaches to estimate ΔL values. We consider the ΔL values of two groups of wavelength-independent and -dependent models, showing a significant difference between the groups (Figure 5b). Our results indicate that the wavelength-dependent ΔL values in the blue to NIR regions are not significantly variable. However, the uncertainties in R_rs(λ) spectra mostly depend on changes in the amplitudes of spectra rather than on the variability of ΔL values [20,71]. In addition, the results confirm the validity of wavelength-independent models for ΔL estimations (Figure 3 and Figure 4). For instance, the presented showcases (Figure 7 and Figure 8) show that the performance of the Ku13 model is relatively similar to the models that use wavelength-independent ΔL estimation (e.g., BA18). It can be inferred that the uncertainty in R_rs(λ) spectra primarily arises from the spectrally dependent shape of the surface-reflected skylight (ρ×L_S(λ)/E_d(λ), second term in Equation (1) and form accurate estimation of ΔL values, which may be either wavelength-dependent or wavelength -independent. The surface-reflected skylight spectra are highly variable under different sky and surface conditions, particularly at shorter wavelengths (e.g., 350-400 nm) [19]. Therefore, the performance of above-water R_rs(λ) models depends on the simulation of the surface-reflected skylight and accurate estimation of ΔL values.

The presence of aerosols significantly modifies the surface-reflected skylight spectra even for low values of aerosol optical thickness. The HT23 model, which estimates the ρ value as a function of directional factors, surface roughness, and AOT, should be used only under clear-sky conditions [22]. In this study, we estimated the AOT values using inverse modeling of solar downwelling irradiance where no direct measurements of aerosols are available (Section 3.5). It is important to acknowledge that uncertainty in the modeling of AOT affects the results of the HT23 model. However, the HT23 model shows a relatively good performance with a Total Score ≥ 0.8 in clear sky and low wind conditions (Case-I, Case-VII, and Case-VIII) and poor performance (Total Score ≤ 0.6) in roughened surface and cloudy conditions (Figure 4). A more detailed assessment of the HT23 model would involve direct retrieval of aerosols and distinctive measurements [72] or modeling [64,66] of direct and diffusive components of solar downwelling irradiance.

Although 3C is flexible and effective in a variety of environmental and measurement conditions, its performance is inherently tied to the bio-optical model used to simulate the water optical properties. Therefore, its dependency on bio-optical modeling introduces certain limitations and implications for R_rs retrieval. Since the 3C model requires predefined ranges of WCCs and IOPs for SIOP-based L_T/E_d modeling, if the assumed SIOPs do not accurately represent the actual optical properties of the water body, the retrieved R_rs may deviate from the true value. This limitation makes 3C less adaptable to waters with unknown or highly variable optical properties, such as coastal and estuarine waters where IOPs are spatially and temporally dynamic. The impact of this constraint can be significant when trying to retrieve low-magnitude R_rs values in highly turbid or oligotrophic waters (e.g., high solar zenith or low light conditions). This constraint indicates that 3C is not a universal solution for R_rs retrieval but rather an approach that works best in environments where the optical properties are well-characterized. Addressing bio-optical constraints through region-specific SIOPs datasets, improved uncertainty quantification, and validation efforts can enhance the robustness of 3C for different aquatic environments.

Moreover, inelastic scattering processes like Chlorophyll-a and CDOM fluorescence contribute significantly to uncertainties in R_rs(λ). Phytoplankton fluorescence at ~685 nm can contribute to a non-negligible fraction of R_rs(λ) that depends on Chlorophyll-a concentration, physiological state, and incident irradiance [73]. CDOM fluorescence at ~400–450 nm further complicates the interpretation of R_rs(λ) [74]. None of the R_rs(λ) glint correction methods (Table 3) account for inelastic scattering, treating it as part of the overall uncertainty. Some bio-optical models include a chlorophyll fluorescence term to correct R_rs(λ) in the 680–690 nm range (e.g., Fluorescence Line Height, FLH) [75]. CDOM fluorescence correction is more challenging, but hyperspectral techniques or fluorescence quantum yield estimates can help. The RT Simulations models, like Hydrolight, for inelastic scattering can be explicitly included in forward simulations to refine R_rs(λ) retrieval algorithms [76]. Given our expertise in this study, incorporating fluorescence into the proposed uncertainty analysis for the 3C might be a valuable next step.

Environmental conditions greatly impact the performance of above-water R_rs(λ) models and contribute to the uncertainty in the corresponding R_rs(λ) spectra. Our findings indicate that the sky cloudiness, solar position, and wind speed are the most important factors affecting fixed above-water radiometry. The above-water R_rs(λ) models are generally less sensitive under clear sky conditions compared to cloudy sky conditions. This is mainly due to smaller variations of ρ and ΔL values in clear sky conditions. Under scattered cloud conditions (Case-IV), models that simulate radiances (i.e., 3C) and account for wind speed to estimate ρ values, or models that utilize the spectral shape of R_rs(λ) to estimate ΔL values (i.e., JD20, Ku13), demonstrate good performance (Figure 4). In contrast, the other models that estimate the ρ and ΔL values regardless of cloudiness have poor performance. The MO99, MO15, and HT23 use a lookup table for ρ as a function of wind speed, assuming a Cox-Munk relationship between wind speed and surface roughness [77]. Nevertheless, these models show poor performance (Total Score < 0.6) during windy conditions (Case-V and Case-VI). The BA18, which uses a fixed ρ value, shows poor performance in cloudy and high wind speed conditions. The minimum discrepancies between models during clear sky conditions (Case-I) are observed at θs in the range of ~30° - 50° and Δφ > 40°. For Δφ ≤ 40°, the 3C, JD23, Ru05, and Ku13 perform relatively better performance, but it makes sense to omit these spectra. This can be attributed to the impact of the sun’s glint on the surface-reflected skylight, which makes the estimated ρ value from lookup tables invalid for estimating the ΔL values. For θs ≥ 80°, all models, with the exception of Ku13, demonstrate strong performance in the range of blue to red wavelengths. However, they exhibit limited effectiveness in the NIR region, with the notable exceptions of the 3C, Ku13, and JD20 models. This indicates that more than 60% of the modeled R_rs(λ) of any model are questionable or unusable, except for the 3C model, which shows the best performance overall.

This study indicates that the 3C-like methods account for marine environments with multiple scattering effects and varying sky and surface conditions, making it particularly useful in highly turbid waters. However, in estuarine environments where detrital organic matter and SPM significantly contribute to the optical properties, backscattering coefficients vary significantly [78]. This variability could influence the effectiveness of the 3C-like models if not properly parameterized. The 3C model assumes a particular range of backscattering conditions based on empirical data; therefore, modifications may be needed to optimize its performance in estuaries where NAP and resuspended sediments dominate. In addition, the contribution of CDOM is often significant in estuaries [79]. Since CDOM absorption impacts the water-leaving radiance, integrating inelastic scattering effects into the 3C framework could further improve accuracy. Therefore, a wavelength-dependent of ρ(λ) on estuarine water conditions might be necessary. Refinements to account for spectral variability in surface reflectance, Chla and CDOM fluorescence, and variable bbp contributions will enhance generalizability of 3C-like models across different coastal and estuarine environments.

6. Conclusions

This study reviews and evaluates the well-known bio-optical models and above-water R_rs(λ) glint-correction methods from a fixed automated above-water radiometry under various environmental conditions. The forward modeling of the IOP-to-R_rs transformation using the average of in-situ SIOPs provides a valid fiducial reference for evaluating the above-water R_rs(λ) models from radiometric measurements. Besides the performance of R_rs(λ) models, the quality of simulated R_rs(λ) should be considered for practical remote sensing applications. As a general result, the 3C model, which directly simulates the surface-reflected radiance rather than estimating the ρ and ΔL values, provides the best method for modeling R_rs(λ) from automated fixed above-water radiometric hyperspectral measurements during sub-optimal conditions.