A geometry-dependent surface Lambertian-equivalent reflectivity product for UV–Vis retrievals – Part 2: Evaluation over open ocean

Satellite-based cloud, aerosol, and trace-gas retrievals from ultraviolet (UV) and visible (Vis) wavelengths depend on the accurate representation of surface reflectivity. Current UV and Vis retrieval algorithms typically use surface reflectivity climatologies that do not account for variation in satellite viewing geometry or surface roughness. The concept of geometry-dependent surface Lambertian-equivalent reflectivity (GLER) is implemented for water surfaces to account for surface anisotropy using a Case 1 water optical model and the Cox–Munk slope distribution for ocean surface roughness. GLER is compared with Lambertian-Equivalent reflectivity (LER) derived from the Ozone Monitoring Instrument (OMI) for clear scenes at 354, 388, 440, and 466 nm. We show that GLER compares well with the measured LER data over the open ocean and captures the directionality effects not accounted for in climatological LER databases. Small biases are seen when GLER and the OMI-derived LER are compared. GLER is biased low by up to 0.01–0.02 at Vis wavelengths and biased high by around 0.01 in the UV, particularly at 354 nm. Our evaluation shows that GLER is an improvement upon climatological LER databases as it compares well with OMI measurements and captures the directionality effects of surface reflectance.

1 Introduction

Satellite retrievals of clouds, aerosols, and trace gases rely on the accurate representation of surface reflectivity. Many modern satellite ultraviolet (UV) and visible (Vis) trace-gas algorithms use the mixed Lambertian-equivalent reflectivity (LER) model which assumes the measured top-of-atmosphere (TOA) radiance is a combination of the clear-sky and cloudy-sky radiances weighted by an effective cloud fraction (ECF) (Koelemeijer et al., 2001; Seftor et al., 1994; Stammes et al., 2008). While many databases of LER currently exist, they are typically climatological LER datasets based on the minimum or the mode of the LER distribution. Therefore, they do not account for the directionality effects due to the variation of satellite viewing and solar illumination geometries. Further, they do not capture day-to-day change in the roughness of the water surface due to changing wind speed which can impact surface reflectance. Additionally, climatological LER databases may be affected by residual aerosol and cloud contamination. These databases include Kleipool et al. (2008) using data from the Ozone Monitoring Instrument (OMI) for wavelengths 328–499 nm, Koelemeijer et al. (2003) from the Global Ozone Monitoring Experiment (GOME) for wavelengths 335–772 nm, and Tilstra et al. (2017) from GOME-2 for wavelengths between 335 and 772 nm.

The bidirectional reflectance distribution function (BRDF) describes the variations of the surface reflectivity for all illumination and viewing geometries. Given a specific solar and viewing angle from a satellite and the BRDF of the surface, a quantity known as the bidirectional reflectance factor (BRF) can be derived for that surface. BRF is defined as the ratio of the radiant flux reflected by a surface to the radiant flux reflected into the identical beam geometry by an ideal diffuse Lambertian surface, irradiated under the same conditions as the sample surface (Schaepman-Strub et al., 2006). A peak in the BRDF arises due to sunglint phenomena when the Sun and satellite viewing angles are similar and oriented in a forward scattering geometry. Sunglint is strongest for smooth water surfaces that permit nearly perfect Fresnel reflection of direct light from the ocean surface (Kay et al., 2009; Cox and Munk, 1954; Thomas and Stammes, 1999). Additionally, there is an increase in radiance near the outer edges of the satellite swath viewing geometry, when TOA radiances increase due to increased effects of Rayleigh scattering. The increase in diffuse sky reflection from the ocean surface becomes more significant at these longer path lengths, and the relative contributions from the water-leaving radiance to the intensity reaching the satellite is reduced (Vasilkov et al., 2017).

Vasilkov et al. (2017) introduced a concept known as geometry-dependent LER (GLER), where for a specific viewing geometry, TOA radiances are simulated over a non-Lambertian surface using the BRDF. GLER can be easily implemented into trace-gas algorithms by simply replacing the currently used LER climatologies. For Vis wavelengths, the change in surface reflectivity associated with the implementation of GLER was found to directly decrease the OMI NO₂ air mass factor (AMF) over land and oceans by as much as 15 % and change the OMI NO₂ AMF indirectly by an additional −22 % to 13 % through changes to retrieved cloud properties (Vasilkov et al., 2017). In the UV, Ahn et al. (2014) noted that unrealistic surface albedo is one of the main causes for uncertainty in aerosol retrievals, while Torres et al. (1998) reported that surface albedo errors can lead to as much as a 5 % error in retrieved aerosol optical depth (AOD) for weakly absorbing aerosols.

Qin et al. (2019) evaluated the OMI 466 nm GLER product over land generated using BRDF data from the Moderate Resolution Imaging Spectroradiometer (MODIS). They compared GLER with OMI-derived LER globally and found GLER was biased low by 0.01 to 0.02 relative to OMI. The difference was attributed to several factors, including small calibration differences between MODIS and OMI and possible residual cloud and/or aerosol contamination in the OMI data that was not completely filtered out.

Here we evaluate the GLER product generated for the OMI instrument for ocean scenes. OMI is a hyperspectral UV–Vis (270–500 nm) imager on board the NASA Aura satellite, which was launched in July 2004. The high spectral resolution (0.42–0.63 nm) of the UV (270–370 nm) and Vis (350–500 nm) channels enables retrievals of many important atmospheric constituents including O₃, NO₂, SO₂, and aerosols. OMI has a spatial resolution of 13 km × 24 km at nadir with a field of view (FOV) of 0.8^∘ in the flight direction and 115^∘ across the swath. Prior to 2008, OMI provided global coverage daily with a repeat cycle of 16 d. The OMI row anomaly affects data coverage starting in mid-2007 (Schenkeveld et al., 2017), but a substantial amount of high-quality global data remain available thereafter.

In this work we focus on evaluation of GLER at UV and Vis wavelengths over oceans. In the UV, 354 and 388 nm are chosen for evaluation as both are important in the OMI retrieval of aerosol properties (Torres et al., 2007), and additionally OMI Raman cloud retrievals are performed at 354 nm (Vasilkov et al., 2008). For Vis wavelengths, 440 and 466 nm are chosen as they are important for O₂–O₂ cloud retrievals (Vasilkov et al., 2018; Veefkind et al., 2016) as well as NO₂ retrievals (Krotkov et al., 2017; Lamsal et al., 2014).

Whereas the land product described in Qin et al. (2019) used a model of BRDF with input from MODIS, GLER for ocean scenes is produced solely by modeling of water-leaving radiance and surface reflection. These surface-leaving radiance contributions are geometry-dependent, and the anisotropic nature of light backscattered by the ocean has been studied in many papers (see, e.g., Gordon , 1989; Morel and Gentili, 1991, 1993, 1996; Park and Ruddick, 2005; Lee et al., 2013). For GLER, water-leaving radiance is simulated using a Case 1 water model (Morel, 1988; Morel and Maritorena, 2001) that depends only on chlorophyll concentration. This model is described in detail in Sect. 2.2. Reflection from the ocean surface is modeled using the Cox–Munk slope distribution (Cox and Munk, 1954) and is further described in Sect. 2.3.

The algorithms and approaches described in this paper are relevant to NASA's future Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission. PACE is currently planned to launch in 2022–2023. Global PACE observations will provide data to monitor oceanic and atmospheric variables important for Earth's ecosystem, carbon cycle, and climate studies. The PACE Ocean Color Instrument (OCI) is designed as a wide-swath imaging spectrometer with a 1 km ground nadir resolution, a 5 nm spectral resolution between 345 and 890 nm, and several shortwave infrared bands (Werdell et al., 2019). As compared with the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), the Moderate resolution Imaging Spectrometer, and the Visible Infrared Imaging Radiometer Suite (VIIRS), the OCI will additionally measure TOA radiances in the UV to help identify phytoplankton composition and harmful algal blooms. PACE's spectral coverage from the UVA to green wavelength region and in the red to near-infrared will enable unparalleled evaluation of ocean ecosystem properties in optically complex waters and in regions of increasing eutrophication (Cetinic et al., 2018).

2 Data and methods

2.1 VLIDORT radiative transfer model

For radiative transfer calculations we use the Vector Linearized Discrete Ordinate Radiative Transfer (VLIDORT) model. VLIDORT is a multiple scattering radiative transfer model that can simulate Stokes vectors at any level in the atmosphere and for any scattering geometry with a Lambertian or non-Lambertian underlying surface (Spurr, 2006; Spurr and Christi, 2019). VLIDORT can simulate attenuation of solar and line-of-sight path radiances in a spherical atmosphere. In this study, we correct for the effects of atmospheric sphericity for both incoming solar and outgoing viewing directions based on a regular pseudospherical geometry calculation. This is important for large solar and viewing zenith angles. We also include polarization using the vector mode, because to neglect it can lead to considerable errors for modeling backscattered radiances in the UV–Vis wavelength range.

2.2 Water-leaving radiance implementation

VLIDORT has a supplement (VSLEAVE) for the generation of surface-leaving radiances for use as inputs to the main radiative transfer calculation in the atmosphere. This supplement can be used for either simulations of solar-induced fluorescence or water-leaving radiances from Case 1 waters. Our Case 1 water model accounts for the bidirectional effects following Morel and Gentili (1996). In this paper we do not account for vibrational Raman scattering in ocean water (Morel et al., 2002; Vasilkov et al., 2002). The common Case 1 water model developed for the Vis (Morel, 1988) was extended to the UV using a parameterization of the particulate matter absorption coefficient from Vasilkov et al. (2002, 2005). The model requires as input several quantities that affect absorption and scattering properties of seawater and its constituents. Extinction coefficients for water absorption are taken from Lee et al. (2015), chlorophyll absorption coefficients from Vasilkov et al. (2005) below 400 nm and Lee et al. (2005) at longer wavelengths, and colored dissolved organic matter (CDOM) absorption from Morel and Maritorena (2001). The water scattering coefficients we use are from Morel et al. (2007), and those for chlorophyll scattering are from Morel and Maritorena (2001). A detailed description of these parameterizations is provided in Appendix A2.

The computation of emerging water-leaving radiance L_w depends not only on the optical properties of marine constituents and radiative processes in the near-surface ocean, but also on the total atmospheric direct and diffuse downwelling transmittance T_atm of atmospheric light through the air–water interface. This complicates separation of the water-leaving calculation and the calculation of atmospheric radiance propagation. Additionally, T_atm will in general depend on the surface-leaving contribution and hence on marine constituents. VLIDORT and its supplement VSLEAVE are therefore coupled. This coupling can be treated formally with a coupled ocean–atmosphere radiative transfer model such as that described in Spurr et al. (2007). Here, however, we have developed a simple coupling scheme for VLIDORT that ensures the value of L_w used as a input at the ocean surface will correspond to the correct value of the downwelling flux reaching the surface interface. The first applications of this new water-leaving model were presented in Vasilkov et al. (2017) and Sayer et al. (2017). The coupled model approach is described further in Appendix B.

2.3 Cox–Munk BRDF implementation

A supplement (VBRDF) is implemented in VLIDORT to account for the reflection of the water surface using the Cox–Munk slope distribution (Cox and Munk, 1954). We use the full form of the Cox–Munk distribution in which the facet slope variance depends on both wind speed and wind direction. Polarization at the ocean surface is accounted for using a full Fresnel reflection matrix as suggested by Mishchenko and Travis (1997). Additionally we account for contributions from oceanic foam that can be significant for high wind speeds using work by Frouin et al. (1996).

2.4 Ancillary data for water model

As mentioned in the introduction, modeling of the water-leaving radiance requires information on the chlorophyll concentration, and the modeling of Cox–Munk surface roughness depends on wind speed and wind direction. These inputs are not available directly from the OMI satellite and so other sources of ancillary inputs are required.

The wind speed measurements for GLER come from a pair of satellite microwave imagers. The first is the Advanced Microwave Scanning Radiometer – Earth Observing System (AMSR-E) instrument on board the NASA Aqua satellite with a spatial resolution of 0.25^∘ (Wentz and Meissner, 2004). The AMSR-E instrument, however, ceased operations in October 2011 due to an issue with the spinning mechanism (Wentz and Meissner, 2007). After October 2011, data are taken from the Special Microwave Imager/Sounder (SSMIS) with a spatial resolution of 0.25^∘. SSMIS is on board the Air Force Defense Meteorological Satellite Program (DMSP) satellite F16 (Wentz et al., 2012). While the SSMIS instrument has been operating since October 2003, the AMSR-E instrument was chosen for the first half of the OMI mission due to the small difference in Equator crossing times of 7–15 min between the Aqua and Aura satellites, whereas the F16 satellite crossing times range from 6 h behind Aura in 2005 to 2 h behind Aura currently. In future work, we plan to replace the SSMIS F16 wind speed with the AMSR-2 wind speed data. AMSR-2 is on board the Global Change Observation Mission Water Satellite 1 (GCOM-W1), which has an Equator crossing time more similar to OMI (Imaoka et al., 2012). Gaps in the wind speed data due to extreme glint or rain are filled by the Global Modeling Assimilation Office (GMAO) Goddard Earth Observing System Model Forward Processing for Instrument Teams (GEOS-5 FP-IT) near-real-time assimilation with a spatial resolution of 0.625^∘ longitude by 0.5^∘ latitude and a temporal resolution of 3 h (Lucchesi, 2013). Wind direction data are also from the GEOS-5 FP-IT model.

Monthly chlorophyll data from the MODIS instrument, which is on board the NASA Aqua satellite, were used in modeling of the water-leaving radiance (Hu et al., 2012). These chlorophyll data have a spatial resolution of 4 km. The MODIS daily chlorophyll data are not used due to large gaps caused by clouds and aerosols. Some gaps still exist in the monthly data and are filled by other datasets from the MODIS team. These include monthly climatological and yearly chlorophyll datasets which are also at a 4 km resolution. The benefit of using the monthly chlorophyll data instead of the climatological data comes from the ability to capture interannual trends due to phenomena such as the El Niño–Southern Oscillation (ENSO).

2.5 Calculation of LER and GLER

Using the equation from Dave (1978), LER (R) can be calculated from TOA radiance (I_comp) by inverting the following:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{I}_{\mathrm{comp}}& (\mathit{\lambda },\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\varphi },P_{\mathrm{s}},{\mathrm{BRF}}_s)=I_{\mathrm{0}}(\mathit{\lambda },\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},\mathit{\varphi },P_{\mathrm{s}})\\ & +{\displaystyle \frac{RT(\mathit{\lambda },\mathit{\theta },{\mathit{\theta }}_{\mathrm{0}},P_{\mathrm{s}})}{\mathrm{1}-R{S}_b(\mathit{\lambda },P_{\mathrm{s}})}},\end{array}\end{array}$$

where λ is wavelength; θ is the viewing zenith angle (VZA); θ₀ is the solar zenith angle (SZA); ϕ is the relative azimuth angle (RAA); P_s is the surface pressure; I₀ is the path scattering radiance by the atmosphere, calculated as the TOA radiance for a black surface; T is total (direct + diffuse) solar irradiance reaching the surface and reflected back to the satellite, multiplied by the transmittance; S_b is the diffuse flux reflectivity of the atmosphere; and R is LER.

In order to calculate GLER, we use VLIDORT to simulate I_comp for a clear sky over a non-Lambertian surface with the water-leaving radiance model described in Sect. 2.2 along with the Cox–Munk slope distribution for surface roughness. A lookup table (LUT) approach is used to calculate TOA radiance operationally, as running VLIDORT can be computationally expensive. Details on the LUT approach can be found in Appendix C. Given TOA radiance from VLIDORT, we can then calculate GLER using Eq. (1).

2.6 OMI data and selection criteria

The measured LER data used to evaluate ocean GLER in this study were retrieved from OMI Collection 3 Level 1b Vis channel radiance data by inverting Eq. (1) where I_obs is used in place of I_comp. The OMI radiances are normalized using the OMI day-1 solar irradiance spectrum adjusted for variation in Earth–Sun distance when radiance measurements were collected. The GLER product is designed to characterize the magnitude and the angular variability of the Earth's surface reflectance in a Rayleigh atmosphere, and therefore several aspects of instrument calibration must be considered. The absolute radiometric calibration error will introduce bias and inconsistency across the measurement swath for LER at any single wavelength. Dobber et al. (2008) estimated that the uncertainty in radiometric calibration of OMI Collection 3 Sun-normalized radiances is under 2 % and that the relative viewing angle dependence is also less than 2 %. Schenkeveld et al. (2017) evaluated long-term changes in the absolute radiometric response of the OMI instrument and estimated degradation of approximately 1 %–1.5 % per decade in the wavelength region used in this study. Since we compare results at 354, 388, 440, and 466 nm in this study, the spectral dependence of the OMI calibration is also an important consideration. Little work has been published on this topic, although the study by Jaross and Warner (2008) compared the OMI Sun-normalized radiances to radiative transfer model simulations over Antarctic ice and showed that the spectral dependence of the OMI calibration is within the absolute radiometric uncertainty.

Beginning in mid-2007, OMI experienced an anomaly known as the “row anomaly” that has affected the L1b radiance data. There have been several impacts from the row anomaly including decreased radiances due to possible blockage, increased signal due to sunlight being reflected into the instrument, a wavelength shift due to a change in the slit function, and earthshine radiances from outside the FOV that are reflected into the nadir port. The row anomaly is further explained in Schenkeveld et al. (2017). For this reason, after 2007 we focus only on rows 1–21, which are not affected by the row anomaly.

Absorption by O₂–O₂ and O₃ were accounted for at 440 and 466 nm but neglected for 354 and 388 nm. Since GLERs were simulated for a Rayleigh-only atmosphere, pixels with absorbing aerosols are removed using the OMAERUV absorbing aerosol index AI ($\left|\mathrm{AI}\right|>\mathrm{0.5}$ are removed) (Torres et al., 2007). We compare the evaluation for two independent cloud screening methods to determine which will better represent the GLER evaluation. The MODIS geometrical cloud fraction (GCF) is retrieved from the 15 µm CO₂ absorption region (Menzel et al., 2007) and is colocated to the OMI FOV in the OMMYDCLD product (Joiner, 2014). The OMI Raman cloud product contains cloud pressure based on rotational Raman scattering and ECF calculated at 354 nm using the Cox–Munk distribution to model ocean surface reflectivity (Vasilkov et al., 2008).

We compare cases with and without sunglint separately because the reflection of light in each case is quite different. For comparisons excluding sunglint scenes, we screen out data with a sunglint angle of less than 20^∘ in which sunglint can occur. For the comparisons with sunglint, while the sunglint angle of 20^∘ is again used to choose the sunglint region, additional screening based on the spectral dependence of the measured LER is performed to remove clouds within the sunglint region. The reason for this is that cloud fraction retrievals are affected by sunglint. The difference in LER occurs because of a spectrally dependent error in the underestimation of the Rayleigh scattering of diffuse light when one assumes a Lambertian ocean surface, when the reflectance is in fact specular. We select sunglint scenes when the difference between the measured LER at 354 and 388 nm is less than −0.05. We note that some weaker sunglint has an LER difference that is not below this threshold, but here we focus on stronger glint that has no cloud contamination.

In addition to the OMI-derived LER, we compare with the Kleipool LER climatology (Kleipool et al., 2008), since a number of current operational algorithms use these LER data as input. There are two LER datasets available from the Kleipool data, one representing the monthly minimum LER and another determined through interpretation of LER histograms. Both are shown in our evaluation as each is used in some algorithms.

3 Results and discussion

3.1 Global comparison of GLER and OMI-derived LER

First we compare GLER with the OMI-derived LER globally for January and July 2006 at four wavelengths in Figs. 1 and 2. To determine a cloud screening method for the evaluation of sunglint-free data, in Fig. 1 we compare GLER with the OMI-derived LER using the cloud screening methods introduced in Sect. 2.6. We note there is a spectral dependence in the difference between GLER and the OMI-derived LER. At 354 nm, GLER is biased high compared to the OMI-derived LER, whereas no bias exists at 388 nm. For longer wavelengths of 440 and 446 nm, GLER is biased low compared to the OMI-derived LER. As shown in Table 2, the GLER and the OMI-derived LER compare best in January at 388 nm where R² is 0.76 and the bias is 0.002 with the Raman-based ECF (using the MODIS GCF R² is 0.60 and the bias is 0.007).

Figure 1 Scatterplots of OMI-derived LER vs. GLER for two months in 2006 (January on left and July on right) with possible sunglint removed at four wavelengths (354, 388, 440, and 466 nm). Data are for deep ocean (based on OMI Level 1b ground pixel quality flags) and have been screened for aerosols (OMAERUV |AI| < 0.5). Clouds are screened through two different methods which are described in Sect. 2.6. Only data with a solar zenith angle less than 70^∘ are included as cloud shadowing at high view angles can decrease the OMI-derived LER.

Figure 2 Scatterplots of OMI-derived LER vs. GLER for two months in 2006 (January on left and July on right) for pixels with sunglint at four wavelengths (354, 388, 440, and 466 nm). Data are for deep ocean with no screening for clouds or aerosols.

Overall the comparison is better when using the Raman-based ECF cloud screen than when using the MODIS GCF. This is expected given that there is a 15 min window between the Aqua and Aura overpass times in 2006 (becomes 7 min in 2009) leading to some change in cloud cover. It is also worth noting, since OMI has a wider swath than MODIS, cloud retrievals are not available from MODIS for pixels on the edge of OMI swath (these pixels are not shown in Fig. 1). For these reasons the Raman-based ECF will be used for cloud screening in the rest of the paper.

In Fig. 2, the comparisons of GLER and OMI-derived LER are presented for data with sunglint. These data were not screened with cloud or aerosol retrievals, as strong glint can lead to artificial classification of aerosol or clouds in the retrieval algorithms. As shown in Table 3, the bias between GLER and the OMI-derived LER is smaller for the data with sunglint than the sunglint-free data. For sunglint pixels, the bias is smaller at UV wavelengths where the water-leaving radiance contributes the least. For the brighter glint data, there is much more uncertainty and GLER is biased high compared with the OMI-derived LER. If the measured wind speed is too low, the model will overestimate the LER of the glint, whereas when the measured wind speed is too high, the model will underestimate the LER of glint. This sensitivity to the wind speed will be further evaluated in Sect. 3.5. The small bias is possibly caused by aerosols which scatter or absorb the direct light causing a small dimming effect in the OMI glint data. This issue will be evaluated in Sect. 3.3.

Table 1 Statistical analysis of GLER vs. OMI LER for nonsunglint scenes. January (July) 2006 deep ocean only with Raman-based ECF = 0.0 (corresponds to the second and fourth column of Fig. 1).

Table 2 Statistical analysis of GLER vs. OMI LER for nonsunglint scenes. January (July) 2006 deep ocean only with MODIS GCF = 0.0 (corresponds to the first and third column of Fig. 1).

Table 3 Statistical analysis of GLER vs. OMI LER for sunglint-only scenes. January (July) 2006 deep ocean only.

3.2 Angular behavior of GLER

Figures 3 and 4 show a comparison of the cross-track dependence of GLER and OMI-derived LER along with the Kleipool LER climatology for a few solar zenith angle ranges screened for clouds using the Raman-based ECF and screened for absorbing aerosols with the OMAERUV AI. GLER follows a similar cross-track pattern at various solar zenith angles to the OMI-derived LER that varies with wavelength. However, there is a bias between GLER and the OMI-derived LER which varies by wavelength. For Vis wavelengths, the OMI-derived LERs are biased high by as much as 0.01–0.02 compared to the GLER, whereas for UV wavelengths the bias is nearly zero at 388 nm and at 354 nm GLER is biased around 0.01 higher than the OMI measurements. In the UV channels, especially at 354 nm, the bias varies both with cross-track and solar zenith angle.

Figure 3 January 2006 LER as a function of VZA for select SZA ranges at 354 and 388 nm over the Pacific Ocean (only deep ocean pixels considered). The blue error bars represent the 2 % calibration uncertainty of OMI (Dobber et al., 2008). Negative VZAs represent the west side of the OMI swath (forward scattering), whereas positive VZAs represent the east side of the OMI swath (backward scattering). Data are screened for clouds with the Raman-based ECF and aerosols are removed with the OMAERUV AI.

Figure 4 Same as Fig. 3 but for 440 and 466 nm.

The Kleipool LER climatology compares well with the OMI-derived LER near nadir but does not capture any of the BRDF effects seen in both GLER and OMI-derived LER since it is a climatology that averages all viewing geometries. We note that there is a slight increase in Kleipool LER at higher viewing zenith angles, but this is simply due to the sampling used for this analysis. Compared to the OMI-derived LER, the Kleipool data have a cross-track-dependent difference of as much as 0.04 outside of the sunglint and 0.1 in the sunglint for UV wavelengths. The difference between Kleipool and the OMI-derived LER (0.01) is smaller for Vis wavelengths outside of the sunglint. At the higher solar zenith angles, the Kleipool monthly data appear to be adversely affected either by ice or residual cloud that was not fully removed from the Kleipool LER climatology. This is evident by the large bias compared with OMI-derived LER as well as the large variability coming from spatial sampling.

3.3 Simulating GLER with aerosols

As noted in earlier sections, aerosols can have an impact on measured LER. For this reason, a simulation was performed to calculate GLER including the effect of aerosols. Aerosol-related input parameters to VLIDORT (layer AOD; single scattering albedo, SSA; and scattering matrix) are from MERRA-2 aerosol reanalysis data which are produced using the GEOS-5 atmospheric model and data assimilation system. MERRA-2 assimilates radiance data from a variety of satellite sensors which are then used to train a neural network to produce AOD which is calibrated to the Aerosol Robotic Network (AERONET) direct-Sun point measurements (Randles et al., 2017; Gelaro et al., 2017). The species-specific aerosol scattering matrices are characterized by six sets of matrix elements generated from the Mie theory.

Figure 5 Map showing the results from GLER aerosol simulation on 10 April 2006 for orbit 9229. Panel (a) is GEOS-5 470 nm AOD used in the simulation and panel (b) is change in GLER when the aerosol contribution is added to the simulations at 466 nm.

Figure 6 Comparison of GLER versus OMI-derived LER for orbits on 10 April 2006 over the Pacific Ocean. Panel (a) is GLER compared to the OMI-derived LER at 466 nm, whereas panel (b) is GLER with contribution from aerosols vs. OMI-derived LER for 466 nm. All data are filtered to remove clouds with the Raman-based ECF.

The simulation was performed for 10 April 2006 over the central Pacific Ocean in order to determine the aerosol effect for a general background oceanic aerosol case. Figure 5 shows the MERRA-2 AOD and the LER change due to added aerosols for orbit 9229 where AOD ranges from around 0.05 in the South Pacific Gyre to larger than 0.4 in the northern Pacific. The impact of the aerosols on GLER ranges from a decrease in LER of 0.01 in sunglint to an increase in LER of greater than 0.04 outside of the glint, especially at higher viewing angles. Figures 6 and 7 show the comparison between GLER and the OMI-derived LER as well as a comparison between GLER with a contribution from aerosols and the OMI-derived LER.

In Fig. 6 the addition of aerosols at 466 nm increases the LER by about 0.01 over darker surfaces, whereas the brighter regions which have some sunglint contribution show a reduction in LER of around 0.01. The combination of these changes improves the regression slope from 1.16 before considering aerosols at 466 nm to 1.0 after aerosols are introduced. After accounting for aerosols, OMI-derived LER is still around 0.01 higher than GLER at 466 nm for darker surfaces. The brighter surfaces, however, which have some sunglint contribution, have little to no bias after accounting for aerosols. Figure 8 shows that aerosols increase GLER generally by 0.01–0.02, with the largest increase at large forward scattering viewing zenith angles.

Figure 7 Same as Fig. 6 but for 354 nm.

Figure 8 LER as a function of VZA at 354 and 466 nm for orbits on 10 April 2006 over the Pacific Ocean. As in Fig. 3 negative VZAs represent the west side of the OMI swath whereas positive VZAs represent the east side of the OMI swath. Since only simulated data are shown, no cloud or aerosol screening is performed.

At 354 nm the aerosol impact is not as significant as that seen for 466 nm with little to no change in the bias for darker surfaces and a reduction of around 0.01 for the brighter surfaces. In Fig. 8 there is a cross-track dependence in the aerosol contribution to GLER at 354 nm. For geometries with forward scattering (negative VZAs), the aerosol contribution can effectively increase the derived LER by 0.01 or more, whereas for backward scattering geometries (positive VZA) there is a small decrease in LER. This view angle dependence of the aerosol impact would remove the view-angle-dependent bias between GLER and OMI-derived LER seen in Fig. 3. Therefore, applying the aerosol impact to Fig. 3 would result in GLER being approximately 0.01 higher than the OMI-derived LER at 354 nm for all view angles.

In this case study, we note that an AOD of 0.1–0.15 increased the LER by as much as 0.01–0.02 at 466 nm, with the largest increase being in the forward scattering direction. At 354 nm, however, similar AOD values slightly decrease LER in the backward scatter direction but can increase LER by as much as 0.01 in the forward scattering direction. While this analysis was for only a specific case study, we note that the aerosol contribution likely accounts for some of the difference between GLER and the OMI-derived LER.

3.4 Interannual variability of LER

Surface LER over the ocean can change day to day with the chlorophyll concentration affecting the water-leaving radiance contribution and changes in the wind speed altering the roughness of the water surface. There is also a seasonal variation in LER due to the changing viewing geometry of satellite measurements as the SZA changes through the year. In Figs. 9 and 10, GLER and OMI-derived LER are shown for the duration of the OMI mission to evaluate the ability of GLER to capture these variations. In this figure, data were selected for a region in the South Pacific Gyre (180–120^∘ W, 30^∘ S–Equator) as this region is generally free of strong aerosols and has relatively low cloud fractions. OMI row 10 is evaluated to avoid sunglint as well as to avoid the OMI row anomaly which impacts many of the OMI rows starting around 2007 (Levelt et al., 2018; Schenkeveld et al., 2017). While this relatively small bias between GLER and the OMI-derived LER is evident, GLER does follow the same general trend as the OMI measured data.

Figure 9 Monthly mean of GLER and OMI-derived LER in the equatorial Pacific Ocean (180–120^∘ W, 30^∘ S–0) for four wavelengths (354, 388, 440, and 466 nm) as a function of time for the OMI mission from row 10. Corresponding chlorophyll and wind speed measurements along with sunglint angle used for GLER are shown in the bottom three panels. Data are filtered to remove clouds with the Raman-based ECF.

Figure 10 Trend in the difference between the monthly mean GLER and the OMI-derived LER for the same location as in Fig. 9 at UV and Vis wavelengths for row 10. Data are filtered to remove clouds with the Raman-based ECF.

Figure 10 shows there is about a 0.01 seasonal variation in the GLER and OMI-derived LER difference for all wavelengths which could possibly be due to seasonal changes in aerosols or seasonality in cirrus clouds that are too thin to be retrieved by the Raman-based ECF. We note that in Fig. 10 there is a small downward trend in the difference between GLER and OMI-derived LER of at most 0.005 in LER. At 354 nm a change of 0.005 LER corresponds to approximately 1 % TOA radiance, which is close to the 1 %–1.5 % TOA radiance degradation noted by Schenkeveld et al. (2017).

3.5 Sensitivity to chlorophyll and wind speed

Table 4 Sensitivity to wind speed.

Table 5 Sensitivity to chlorophyll.

Table 6 Sensitivity to wind speed and chlorophyll.

In order to determine the uncertainty in GLER calculations, a sensitivity test was performed based on the inputs of chlorophyll and wind speed. The wind speed measurements were perturbed ±1 m s⁻¹ as Wentz and Meissner (2000) note that is the uncertainty in their wind speed algorithm. The MODIS Ocean Color Team notes that the chlorophyll uncertainty varies regionally but can possibly be as high as 35 % (Moore et al., 2009). We perturb chlorophyll by this 35 % in order to gather an absolute bound of GLER to the chlorophyll input. We place an absolute lower bound on the wind speed in our calculations of 0.4 m s⁻¹ and a lower bound of 0.01 of milligrams per cubic meter (mg m⁻³) on the chlorophyll datasets as measurements of these input below these lower bounds are unrealistic and could lead to large errors in calculation of GLER. In Tables 4, 5, and 6 the magnitudes of the mean and max differences are reported in units of LER.

As seen in Table 4 the average wind speed sensitivity is quite small at around 0.0015, but the maximum sensitivity can become as high as 0.3 due to the high sensitivity of sunglint to wind speed. This is because for extreme sunglint cases a small roughness in the ocean surface will lead to increased scattering of light which will significantly diminish the strength of the glint. It is worth noting that while such roughness decreases the strength of the glint, it will increase the overall size of the region affected by sunglint due to the scattering of light at the ocean surface. The wind speed sensitivity decreases with decreasing wavelength because the fraction of the direct solar light which is responsible for sunglint decreases for shorter wavelengths where the contribution of the diffuse light increases due to Rayleigh scattering.

In Table 5 the chlorophyll sensitivity is shown to be much smaller than the wind speed sensitivity. In contrast with the wind speed sensitivity, the chlorophyll sensitivity is largest at UV wavelengths because the CDOM absorption, which rises with increasing chlorophyll according to the Case 1 water model, exponentially increases for shorter wavelengths.

To determine the combined effect, we additionally calculated GLER perturbing both the wind speed and chlorophyll for the four possible combinations. Table 6 shows the mean difference from the combined sensitivity analysis in GLER is similar to that obtained by only perturbing the chlorophyll. The maximum difference from the combined sensitivity test, however, is similar to that of the wind speed perturbation. This is because while the wind speed has a significant impact on sunglint, only a small fraction of OMI pixels are impacted by glint.

3.6 Additional sources of uncertainty

We note that, in addition to the sensitivities from the input data, uncertainties also result from the modeling of the GLER as well as the OMI data used for the evaluation. One possible source of uncertainty is the water-leaving radiance model being used in the calculations. Here we implement a Case 1 water model that is assumed to be representative of the global oceans with dependence only on chlorophyll. Szeto et al. (2011), however, showed that the world's oceans are optically variable and that optical parameters such as chlorophyll absorption vary even for Case 1 open ocean. Work by Lee et al. (2006) similarly showed that large deviations exist from the presumed Case 1 water model due to uncertainty in the optical properties used to parameterize the models. They determined that the SeaWiFS remote sensing reflectance retrievals at 555 nm for Case 1 water have a deviation of ±50 % from that of a Case 1 bio-optical property model.

Our simulations do not include any vibrational Raman scattering effects which can increase the water-leaving radiance. Westberry et al. (2013) show that Raman scattering can impact the water-leaving radiance in Vis wavelengths as much as 4 %–7 % for low chlorophyll concentrations.

In Sect. 3.3 it is shown that aerosols can increase the LER derived from OMI. Uncertainty in AOD used for this analysis could have an appreciable impact on the evaluation results. Randles et al. (2017) compare the MERRA-2 550 nm AOD with the Maritime Aerosol Network (MAN) and find the MERRA-2 AOD to be biased low by 0.009 with a large spread in the comparison which they note could be due to the uncertainty of the MAN AOD of ±0.02. Even though the AOD uncertainties are small, we have shown that even a 0.05–0.1 AOD increase can increase 466 nm LER by 0.005–0.01.

With regard to the OMI measurements, uncertainty could arise from the cloud screening of the OMI measurements as retrieval of cloud properties can become difficult for thin cirrus clouds, which are especially prevalent over the western Pacific (Nazaryan et al., 2008). It has been shown that such contamination can actually increase MODIS AOD retrievals by 0.015–0.025 (Kaufman et al., 2005). As previously mentioned, there is also up to a 2 % TOA radiance absolute calibration uncertainty with the OMI measurements (Dobber et al., 2008) which can lead to an uncertainty of around 0.01–0.02 in LER in the UV and up to 0.005 LER in the Vis. Several OMI algorithm teams apply adjustments to remove the residual error in viewing angle dependence in the OMI measurements by looking at Earth radiances over land. When these corrections were applied in our evaluation of GLER, we found the difference between GLER and the OMI measurements in the UV and Vis was reduced.

Given the uncertainties listed above, we believe differences between GLER and OMI measurements are likely caused by some combination of these factors. It is possible that these factors vary with wavelength. For example, we have shown that aerosols have a greater impact at 440 and 466 nm, whereas other factors such as chlorophyll uncertainty are more important at 354 and 388 nm where the water-leaving radiance contributes more than the direct reflectance.

4 Conclusion

Previous work Qin et al. (2019) introduced the GLER product for land surfaces based on BRDF input from MODIS. In this paper we evaluate a surface LER product called GLER which accounts for the ocean surface BRDF effects at UV and Vis wavelengths. Surface roughness is modeled using the Cox–Munk slope that depends on wind speed measurements from the AMSR-E and SSMIS instruments. A contribution of water-leaving radiance is also included which is based on chlorophyll input from MODIS.

We evaluated the GLER product over water by comparing with OMI-derived LER at UV and Vis wavelengths. The BRDF effect in the OMI-derived LER is captured well with the GLER data at Vis wavelengths. At UV wavelengths, however, the difference between OMI-derived LER and GLER varies some with viewing angle. This UV effect could possibly be explained by the scattering of aerosols. There is, however, a bias between GLER and the OMI-derived LER which is no more than 0.01 after accounting for the effect of aerosols.

We note the GLER data capture the seasonality and interannual variability seen in the OMI measurements that may be caused by variations in the viewing angles as well as changes in the chlorophyll and wind speed data due to meteorological phenomena. The bias between GLER and OMI-derived LER could be caused by a combination of factors including small calibration errors in the OMI-derived LER, deviations in the OMI measurements from the Case 1 water model, and residual thin clouds that are difficult to screen.

There are several possible applications for the GLER product. It can be used to replace climatological LERs currently used by many cloud, trace-gas, and aerosol algorithms. Additionally, GLER can be used as a tool to evaluate satellite calibration to detect possible drift or striping in instrument retrievals. In future missions such as PACE, GLER can be adapted to perform retrievals of water-leaving radiance at UV and Vis wavelengths.

In future work we plan to implement a Case 2 water model for turbid and coastal waters as well as replace our LUT approach with a neural network approach to reduce the computational time to produce the GLER product.

Appendix A: Description of the VSLEAVE water-leaving radiance model

In this Appendix, we give details of the water-leaving radiance scheme included in the VSLEAVE supplement to VLIDORT version 2.8. Section A1 of this appendix deals with the basic water-leaving formulation, while Sect. A2 deals with the ocean optics model. In particular, the material in Sect. A2 is based on the work of Sayer et al. (2010), which has a comprehensive review of the semiempirical marine optics formula, and a companion paper (Sayer et al., 2017), with the latter containing important updates to the optics model. The treatment is for Case 1 waters.

A1 Water-leaving radiance model

Here we summarize the computation of water-leaving radiances using the VSLEAVE supplement to VLIDORT. A full description of the VLIDORT VSLEAVE supplement used here may be found in Spurr and Christi (2019). Water-leaving output from VSLEAVE consists of three terms which are Sun-normalized radiances. The first is a direct term $L_{\mathrm{w},\mathrm{direct}}(\mathit{\mu },{\mathit{\mu }}_o,\mathit{\varphi })$ which is the water-leaving radiance for solar illumination angle θ₀ and cosine μ₀=cos (θ_o) going into the viewing direction with zenith angle θ and μ=cos θ and the relative azimuth angle ϕ between the solar and viewing directions.

The other two water-leaving radiance outputs may be written $L_{\mathrm{w},m}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})$ and $L_{\mathrm{w},m}({\mathit{\mu }}_i,{\mathit{\mu }}_{\mathrm{0}})$, where ${\mathit{\mu }}_i(i=\mathrm{1},\mathrm{\dots }{N}_{\mathrm{d}})$ are the discrete-ordinate polar cosines, and m is the Fourier component index, $m=\mathrm{0},\mathrm{1},\mathrm{\dots }\mathrm{2}{N}_{\mathrm{d}}-\mathrm{1}$. These are diffuse-term contributions: $L_{\mathrm{w},m}({\mathit{\mu }}_i,{\mathit{\mu }}_{\mathrm{0}})$ is required for the inclusion of surface leaving in the diffuse-scattering boundary condition at the surface, while the term $L_{\mathrm{w},m}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})$ is required for postprocessing of the discrete-ordinate solution. Fourier terms arise from cosine-azimuth expansions of the full functions: $L_{\mathrm{w},\mathrm{direct}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi })=L_{\mathrm{w},\mathrm{0}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})+\mathrm{2}{\sum }_{m=\mathrm{1}}^{\mathrm{\infty }}{L}_{\mathrm{w},m}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})\cos \left(m\right(\mathit{\varphi }\left)\right)$. In the discrete-ordinate approximation with N_d streams, we can only use 2N_d−1 components in this sum. In the postprocessing, it is more accurate to use the complete term $L_{\mathrm{w},\mathrm{direct}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi })$ itself in place of the (less-accurate) Fourier-series truncation, and this exact-term correction is the default in VSLEAVE. In this case, Fourier terms $L_{\mathrm{w},m}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)$ are not needed. Note that we will always need the Fourier components $L_{\mathrm{w},m}\left({\mathit{\mu }}_i,{\mathit{\mu }}_{\mathrm{0}}\right)$ for the diffuse-field calculation. However, when there is no azimuth dependence, only $L_{\mathrm{w},\mathrm{0}}\left({\mathit{\mu }}_i,{\mathit{\mu }}_{\mathrm{0}}\right)$ for m=0 survives. In this study we consider an anisotropic distribution of the water-leaving contribution, but the model can also generate it as an isotropic term L_w,iso(μ₀) which depends only on the incoming solar direction (no azimuth dependence, all outgoing directions equal), in which case $L_{\mathrm{w},m}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})=\mathrm{0}(m\ge \mathrm{1})$ and $L_{\mathrm{w},\mathrm{0}}(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}})=L_{\mathrm{w},\mathrm{iso}}\left({\mathit{\mu }}_{\mathrm{0}}\right)$ for all outgoing polar directions μ , and also $L_{\mathrm{w},\mathrm{direct}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi }\right)=L_{\mathrm{w},\mathrm{iso}}\left({\mathit{\mu }}_{\mathrm{0}}\right)$. The current default for VSLEAVE is for an unpolarized azimuth-independent formalism. Thus only the intensity component of the water-leaving Stokes vector is nonzero, and there is no azimuthal dependence.

Water-leaving radiance may be written as

$$\begin{array}{}\text{(A1)}& L_{\mathrm{w}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi }\right)=L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi };\mathrm{Chl},V\right)T_{\mathrm{atm}}\left({\mathit{\mu }}_{\mathrm{0}}\right),\end{array}$$

for any given combination of angles, where the transmittance T_atm(μ₀) depends only on the solar angle, and $L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi },\mathrm{Chl},V\right)$ is computed from the marine optical properties using a semiempirical model which depends explicitly on the chlorophyll concentration and wind speed V. The ocean-optics model for the determination of $L_{\mathrm{w}}^{\ast }$ is described below.

A2 Ocean-optics model

Though the longest wavelength in the GLER product is presently 466 nm, we describe our model of ocean optics over a wider spectral range below as the model has applications in aerosol and ocean color studies where the longer visible and near-infrared wavelengths are typically used. The water absorption α_W(λ) coefficients have been linearly interpolated from a table of values at every 5 nm from 200 to 900 nm constructed from a number of literature sources. These are Quickenden and Irvin (1980), interpolated with Lee et al. (2015) between 325 and 345 nm; 350–550 nm from Lee et al. (2015); Pope and Fry (1997) for 555–725 nm; and Hale and Querry (1973), Table 1, for 725–900 nm (the latter with 25 nm increments linearly interpolated to 5 nm values). Table entries provided as the imaginary index of refraction k_W are converted using ${\mathit{\alpha }}_{\mathrm{W}}\left(\mathit{\lambda }\right)=\frac{\mathrm{4}\mathit{\pi }{k}_{\mathrm{W}}}{\mathit{\lambda }}$, where wavelengths are in meters (m).

The chlorophyll absorption α_ph(λ) comes from two sources. The first source (in the range 300–400 nm) relies on linear interpolation of two sets of coefficients {a₁(λ),b₁(λ)} given at 10 nm intervals in this range (Vasilkov et al., 2005). The absorption is given by

$$\begin{array}{}\text{(A2)}& {\mathit{\alpha }}_{\mathrm{Ph}}\left(\mathit{\lambda },\mathrm{Chl}\right)=\mathrm{Chl}\cdot a_{\mathrm{1}}\left(\mathit{\lambda }\right)\cdot {\mathrm{Chl}}^{-b_{\mathrm{1}}\left(\mathit{\lambda }\right)},\end{array}$$

where Chl is the chlorophyll concentration and λ is wavelength. The value at 300 nm is used for all λ<300 nm. The second source (over the range 400–720 nm) is based on linear interpolation of two sets of coefficients a₂(λ),b₂(λ) at 10 nm intervals (Lee et al., 2005). The absorption formula in this regime is given by

$$\begin{array}{}\text{(A3)}& {\mathit{\alpha }}_{\mathrm{Ph}}\left(\mathit{\lambda },\mathrm{Chl}\right)=\left[a_{\mathrm{2}}\left(\mathit{\lambda }\right)+b_{\mathrm{2}}\left(\mathit{\lambda }\right)\ln \left(a_{\mathrm{440}}\right)\right]a_{\mathrm{440}},\end{array}$$

where $a_{\mathrm{440}}=\mathrm{0.06}\cdot {\mathrm{Chl}}^{\mathrm{0.65}}$ (Morel and Maritorena, 2001). The value at 720 nm is used for all λ>720 nm. The CDOM absorption is given by Morel and Maritorena (2001):

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{CDOM}}\left(\mathit{\lambda },\mathrm{Chl}\right)=& \mathrm{0.2}\cdot \left({\mathit{\alpha }}_{\mathrm{w}}\left(\mathrm{440}\mathrm{nm}\right)+\mathrm{0.06}\cdot {\mathrm{Chl}}^{\mathrm{0.65}}\right)\\ & \exp \left[-\mathrm{0.014}\left(\mathit{\lambda }-\mathrm{440}\right)\right],\end{array}\end{array}$$

where α_w(440 nm)=0.00635 and λ is in nanometers (nm). The complete absorption is then

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)=& {\mathit{\alpha }}_{\mathrm{W}}\left(\mathit{\lambda }\right)+{\mathit{\alpha }}_{\mathrm{Ph}}\left(\mathit{\lambda },\mathrm{Chl}\right)\\ & +{\mathit{\alpha }}_{\mathrm{CDOM}}\left(\mathit{\lambda },\mathrm{Chl}\right).\end{array}\end{array}$$

We use the following formula for the backscattering coefficient, assuming it is half of the scattering coefficient for pure water Rayleigh scattering, as described in Morel et al. (2007):

$$\begin{array}{}\text{(A6)}& b_{\mathrm{W}}\left(\mathit{\lambda }\right)=\mathrm{0.0028}{\left({\displaystyle \frac{\mathrm{420}}{\mathit{\lambda }}}\right)}^{\mathrm{4.3}},\end{array}$$

with λ in nanometers. For the particulate matter backscattering coefficient, we use the following from Morel and Maritorena (2001):

$$\begin{array}{}\text{(A7)}& {\displaystyle }{b}_{\mathrm{Ph}}\left(\mathit{\lambda },\mathrm{Chl}\right)=b_{\mathrm{Pb}}\left(\mathrm{Chl}\right){\mathit{\beta }}_{\mathrm{bbp}}\left(\mathrm{Chl},\mathit{\lambda }\right),\text{(A8)}& {\displaystyle }\begin{array}{rl}{b}_{\mathrm{Pb}}\left(\mathrm{Chl}\right)=& \mathrm{0.416}{C}^{\mathrm{0.766}};{\mathit{\beta }}_{\mathrm{bbp}}\left(\mathrm{Chl},\mathit{\lambda }\right)=\mathrm{0.002}\\ & +\mathrm{0.01}\left[\mathrm{0.5}-\mathrm{0.25}{\log }_{\mathrm{10}}\mathrm{Chl}\right]{\left({\displaystyle \frac{\mathit{\lambda }}{\mathrm{550}}}\right)}^V,\end{array}\end{array}$$

where the exponent V=0 for Chl>2, and $V=\mathrm{0.5}\left[{\log }_{\mathrm{10}}\mathrm{Chl}-\mathrm{0.3}\right]$ for Chl≤2. The complete backscattering is then

$$\begin{array}{}\text{(A9)}& b_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)=b_{\mathrm{W}}\left(\mathit{\lambda }\right)+b_{\mathrm{ph}}\left(\mathit{\lambda },\mathrm{Chl}\right).\end{array}$$

In the original formulation of water-leaving radiance in VLIDORT, the following formula was used to obtain the basic ocean-surface reflectance (Morel and Gentili, 1992):

$$\begin{array}{}\text{(A10)}& {\displaystyle }\begin{array}{rl}R\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)& =f\left({\mathit{\mu }}_{\mathrm{0}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)\\ & \equiv f\left({\mathit{\mu }}_{\mathrm{0}}\right){\displaystyle \frac{b_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}{a_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}},\end{array}\text{(A11)}& {\displaystyle }\begin{array}{rl}f\left(\mathit{\lambda },\mathrm{Chl},{\mathit{\theta }}_{\mathrm{0}}\right)& =d_{\mathrm{0}}-d_{\mathrm{1}}\mathit{\eta }-d_{\mathrm{2}}{\mathit{\eta }}^{\mathrm{2}}+\left(d_{\mathrm{3}}\mathit{\eta }-d_{\mathrm{4}}\right){\mathit{\mu }}_{\mathrm{0}};\\ \mathit{\eta }& ={\displaystyle \frac{b_{\mathrm{W}}\left(\mathit{\lambda }\right)}{b_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}}.\end{array}\end{array}$$

Here, $f(\mathit{\lambda },\mathrm{Chl},{\mathit{\theta }}_{\mathrm{0}})$ is given with five constants {$d_{\mathrm{0}},d_{\mathrm{1}},d_{\mathrm{2}},d_{\mathrm{3}},d_{\mathrm{4}}$} = {0.6279, 0.0227, 0.0513, 0.2465, 0.3119}, and μ₀=cos (θ₀) is the cosine of the solar zenith angle. In order to assign the water-leaving radiance, the complete reflectance term is given by

$$\begin{array}{}\text{(A12)}& R^{\prime }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)={\displaystyle \frac{R\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)}{\mathrm{1}-\mathit{\omega }R\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)}}.\end{array}$$

Here, albedo ω=0.485, using the value in Austin (1974). The isotropic water-leaving radiance is then obtained after passage through the air–ocean interface:

$$\begin{array}{}\text{(A13)}& S_{\mathrm{iso}}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)\approx {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathit{\pi }}}{T}_{\mathrm{Surf}}\left({\mathit{\theta }}_{\mathrm{0}}\right){\displaystyle \frac{R^{\prime }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\mu }}_{\mathrm{0}}\right)}{\mathrm{|}{n}_{\mathrm{w}}{\mathrm{|}}^{\mathrm{2}}}}.\end{array}$$

Here, n_w is the relative refractive index of water to air. For the flat surface case, the air–water boundary transmittance T_Surf(θ₀) is often set to 1.0. In practice we use Fresnel optics to compute this quantity; values are typically 0.96 or more, depending on the value of θ₀. In the rough-surface case, T_Surf(θ₀) may be computed using glitter calculations based on Gaussian probability wave-facet distributions characterized by wind speed and direction.

The above formulation does not account for the atmospheric transmitted flux T_atm(θ₀) at the ocean surface, a quantity which is propagated through the interface. In the previous formulation, the ratio $\frac{T_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{0}}\right)}{Q}$ was made implicit in the factor $\frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathit{\pi }}$ appearing in Eq. (A13). Also, we replace the $f\left(\mathit{\lambda },\mathrm{Chl},{\mathit{\theta }}_{\mathrm{0}}\right)$ calculation with the direction-dependent ratio $\mathit{\rho }\equiv f/Q$ from Morel and Gentili (1996) and Morel et al. (2002). The water-leaving radiance is then

$$\begin{array}{}\text{(A14)}& {\displaystyle }\begin{array}{rl}& S\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},\mathit{\mu },\mathit{\phi }\right)={\mathit{\mu }}_{\mathrm{0}}{T}_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{0}}\right)T_{\mathrm{Surf}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\\ & {\displaystyle \frac{R^{\ast }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},\mathit{\mu },\mathit{\phi }\right)}{\mathrm{|}{n}_{\mathrm{w}}{\mathrm{|}}^{\mathrm{2}}}},\end{array}\text{(A15)}& {\displaystyle }{R}^{\ast }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},\mathit{\mu },\mathit{\phi }\right)={\displaystyle \frac{\mathit{\rho }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},\mathit{\mu },\mathit{\phi }\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}{\mathrm{1}-\mathit{\omega }f\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}},\end{array}$$

where R_TOT(λ,Chl) is as defined in Eq. (A10), and ρ is the ratio f∕Q. We use a tabulated form of the ratio f∕Q in our calculations.

In order to obtain an isotropic surface-leaving radiance, we derive a quantity $\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)$ from the f∕Q tables by averaging over all outgoing zenith and relative azimuth angles, θ and ϕ, and then interpolating linearly with wavelength λ, followed by cubic spline interpolation and linear interpolation with the solar angle cosine μ₀ and with the logarithm of the chlorophyll concentration. Spline interpolation is necessary because we want smooth and continuous derivatives with respect to Chl when considering linearization, as discussed below. The quantity $\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)$ then defines the isotropic water-leaving contribution through

$$\begin{array}{}\text{(A16)}& {\displaystyle }{S}_{\mathrm{iso}}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)={\mathit{\mu }}_{\mathrm{0}}{T}_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{0}}\right)T_{\mathrm{Surf}}\left({\mathit{\theta }}_{\mathrm{0}}\right){\displaystyle \frac{\stackrel{\mathrm{‾}}{R^{\ast }}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)}{\mathrm{|}{n}_{\mathrm{w}}{\mathrm{|}}^{\mathrm{2}}}},\text{(A17)}& {\displaystyle }\stackrel{\mathrm{‾}}{R^{\ast }}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)={\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}{\mathrm{1}-\mathit{\omega }f\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}}.\end{array}$$

The azimuth dependence is very weak in the f∕Q tables, and we have omitted this dependence in the surface-leaving formulation. However, we can derive nonisotropic surface-leaving f∕Q values by interpolating table entries with the cosine of the outgoing angle μ. The resulting table extractions are then $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{v}}}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)$ and $\stackrel{\mathrm{̃}}{{\mathit{\rho }}_d}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_d\right)$ for each viewing angle μ_v and discrete-ordinate stream μ_d; these quantities are azimuth-averaged. We then have

$$\begin{array}{}\text{(A18)}& {\displaystyle }\begin{array}{rl}& S_{\mathrm{v}}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)={\mathit{\mu }}_{\mathrm{0}}{T}_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{0}}\right)T_{\mathrm{Surf}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\\ & {\displaystyle \frac{R_{\mathrm{v}}^{\ast }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)}{\mathrm{|}{n}_{\mathrm{w}}{\mathrm{|}}^{\mathrm{2}}}},\end{array}\text{(A19)}& {\displaystyle }{R}_{\mathrm{v}}^{\ast }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)={\displaystyle \frac{\stackrel{\mathrm{̃}}{{\mathit{\rho }}_{\mathrm{v}}}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}{\mathrm{1}-\mathit{\omega }f\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}}\right)R_{\mathrm{TOT}}\left(\mathit{\lambda },\mathrm{Chl}\right)}},\end{array}$$

and similarly for the discrete-ordinate directions.

In the rough-surface case, the above analysis for the ocean reflectance still holds, but now we need to generate glitter-dependent transmission terms through the water–air interface, for both (i) the incoming solar directions ${\overrightarrow{T}}_{\mathrm{aw}}\left({\mathit{\theta }}_{\mathrm{0}}\right)$ and (ii) the outgoing line-of-sight ${\stackrel{⟵}{T}}_{\mathrm{wa}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}})$ and discrete-ordinate ${\stackrel{⟵}{T}}_{\mathrm{wa}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_d\right)$ directions respectively. Thus, for instance, the rough-surface water-leaving term for a viewing angle μ_v is

$$\begin{array}{}\text{(A20)}& \begin{array}{rl}& S_{\mathrm{v},RS}\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)={\mathit{\mu }}_{\mathrm{0}}{T}_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{0}}\right){\overrightarrow{T}}_{\mathrm{aw}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\\ & {\displaystyle \frac{R_{\mathrm{v}}^{\ast }\left(\mathrm{Chl},\mathit{\lambda },{\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right)}{\mathrm{|}{n}_{\mathrm{w}}{\mathrm{|}}^{\mathrm{2}}}}{\stackrel{⟵}{T}}_{\mathrm{wa}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{v}}\right),\end{array}\end{array}$$

by analogy with Eq. (A18) and using Eq. (A19).

Appendix B: Coupling of VLIDORT and VSLEAVE

The simplest approximation to T_atm(μ₀) is the decoupled scenario where the transmittance has no dependence on ocean properties. In this case, we drop the T_atm(μ₀) term from the main VSLEAVE result in Eq. (A1) above and then reintroduce T_atm(μ₀) from an internal computation in the main VLIDORT model. The direct transmittance $T_{\mathrm{direct}}\left({\mathit{\mu }}_{\mathrm{0}}\right)=\exp [-{\mathit{\tau }}_{\mathrm{atm}}/{\mathit{\mu }}_{\mathrm{0}}]$, where τ_atm is the total atmospheric vertical optical depth; a closer value which includes a diffuse transmittance component is

$$\begin{array}{}\text{(B1)}& T_{\mathrm{atm}}\left({\mathit{\mu }}_{\mathrm{0}}\right)=\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{atm}}}{{\mathit{\mu }}_{\mathrm{0}}}}\right].\end{array}$$

This equation was adapted from a similar formula in Gordon and Wang (1994). Equation (B1) is easy to implement in VLIDORT. A more accurate expression may be obtained in certain cases by using a precalculated lookup table of T_atm(μ₀) values, computed offline with VLIDORT in a Rayleigh atmosphere over a 270–900 nm wavelength range, and for a number of θ₀. However, T_atm(μ₀) is still decoupled from the VSLEAVE water-leaving radiance output.

The coupling scheme works as follows. From Eq. (A1), we write

$$\begin{array}{}\text{(B2)}& L_{\mathrm{w}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)=L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)T^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right),\end{array}$$

where T^↓(μ₀) is the total (direct and diffuse) downwelling atmospheric transmittance at the ocean surface, and $L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)$ is the water-leaving radiance from VSLEAVE computed with unit transmittance. Here, μ₀ is the solar zenith cosine, and μ is any outgoing stream direction; we assume azimuth independence.

To find the coupling adjustment for T^↓(μ₀), we use an initial estimate $T_{\mathrm{0}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$ which could be the quantity in Eq. (B1) above; another value which we have tried is $T_{\mathrm{0}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)=\frac{\mathrm{3}}{\mathrm{2}}{T}_{\mathrm{Direct}}\left({\mathit{\mu }}_{\mathrm{0}}\right)$. With this starting value, we then have an adjusted water-leaving radiance $L_{\mathrm{0}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)=L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)T_{\mathrm{0}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$ which is then input to a Fourier-zero (azimuth independent) VLIDORT radiative transfer (RT) computation. From this RT computation we then derive an updated total downwelling transmittance $T_{\mathrm{1}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$, which in turn provides an updated water-leaving input $L_{\mathrm{1}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)=L_{\mathrm{w}}^{\ast }\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)T_{\mathrm{1}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$. We repeat the Fourier-zero VLIDORT radiative transfer calculation with this new input, yielding a new result $T_{\mathrm{2}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$ for the transmittance and a new water-leaving value $L_{\mathrm{w},\mathrm{2}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)$. This iteration is stopped when the relative difference in the value of T^↓(μ₀) between two iterations is less than some small convergence criterion. We have found that convergence is rapid: typically only three iterations are needed for convergence at the level of 10⁻⁶.

It is not necessary to carry out a full Fourier calculation for every step. The discrete-ordinate homogeneous solutions and particular integrals do not depend on the surface-leaving radiance, and they need to be established just once from the initial Fourier-zero computation. Also, the complete discrete-ordinate solution is determined through the linear-algebra boundary value problem (BVP) Ax=B, where matrix A is constructed entirely from the homogeneous solutions to the radiative transfer equation (RTE), x is the vector of unknown homogeneous-solution integration constants, and vector B is constructed from the layer particular integrals and also contains the surface boundary condition appropriate for water leaving. Once the matrix inverse A⁻¹ is found, the BVP solution is obtained through straightforward back substitution: $\mathit{x}={\mathbf{A}}^{-\mathrm{1}}\mathit{B}$. Thus, the first guess for water-leaving input L₀(μ,μ₀) will give rise to column vector B₀, with the corresponding solution ${\mathit{x}}_{\mathrm{0}}={\mathbf{A}}^{-\mathrm{1}}{\mathit{B}}_{\mathrm{0}}$. From the discrete-ordinate solution based on x₀, we then derive the next transmittance estimate $T_{\mathrm{1}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)$ and then form the next-guess water-leaving input $L_{\mathrm{w},\mathrm{1}}\left(\mathit{\mu },{\mathit{\mu }}_{\mathrm{0}}\right)$ and associated column vector B₁, from which we get the next solution ${\mathit{x}}_{\mathrm{1}}={\mathbf{A}}^{-\mathrm{1}}{\mathit{B}}_{\mathrm{1}}$, and so on. All column vectors B_p are similar and only the surface-leaving entries are different. Thus the coupling adjustment is tantamount to a series of back substitutions, and this represents very little extra computation load compared with the main radiative transfer equation, finding the inverse A⁻¹. A three-iteration calculation is approximately 2 % slower than a standard one.

Computation of the diffuse downwelling transmittance comes through the discrete-ordinate result:

$$\begin{array}{}\text{(B3)}& {\displaystyle }\begin{array}{rl}{T}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)=& T_{\mathrm{diffuse}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right)+T_{\mathrm{direct}}^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right);\\ & T_{\mathrm{diffuse}}^{\downarrow }={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{{\mathit{\mu }}_{\mathrm{0}}}}\sum _{\mathit{\alpha }=\mathrm{1}}^{n_d}{I}_{\mathit{\alpha }}^{\downarrow }{\mathit{\mu }}_{\mathit{\alpha }}{c}_{\mathit{\alpha }},\end{array}\text{(B4)}& {\displaystyle }{I}^{\downarrow }=\sum _{\mathit{\alpha }=\mathrm{1}}^{n_d}{L}_{\mathit{\alpha }}{Y}_{\mathit{\alpha }}^{-}{e}^{-k_{\mathit{\alpha }}\mathrm{\Delta }}+M_{\mathit{\alpha }}{Y}_{\mathit{\alpha }}^{+}+G^{\downarrow }\left({\mathit{\mu }}_{\mathrm{0}}\right).\end{array}$$

Here, ${\mathit{\mu }}_{\mathit{\alpha }},c_{\mathit{\alpha }},\mathit{\alpha }=\mathrm{1},\mathrm{\dots }{n}_d$ are the discrete-ordinate quadrature values; and I^↓ is the downwelling intensity field at the surface expressed in terms of homogeneous solutions $Y_{\mathit{\alpha }}^{\pm },k_{\mathit{\alpha }}$ in the lowest layer of the atmosphere, particular solutions G^↓(μ₀) in that layer, and integration constants L_α,M_α for that layer as determined from the BVP solution ${\mathit{x}}_{\mathrm{1}}={\mathbf{A}}^{-\mathrm{1}}\mathit{B}$. This flux computation does not require any postprocessing or any evaluations at other levels in the atmosphere.

Appendix C: Description of lookup tables

Processing of GLER using online radiative transfer calculations is not efficient for the OMI mission due to their computationally expensive nature. Instead, LUT interpolation is used for the OMI mission, which speeds up the calculations significantly. To calculate GLER, two separate LUTs were generated: one for the TOA radiances calculated to include geometry-dependent surface BRDF effects and the other to derive LER from these radiances using the quantities I_o, S_b, and T, as described in Sect. 2.5. I_o, S_b, and T depend only on viewing geometries and surface pressure, whereas the TOA radiance table additionally included dependencies on chlorophyll, wind speed, and wind direction. The LUT nodes for the TOA radiance table shown in Table C1 were chosen by analyzing TOA radiance as a function of each input parameter to keep the interpolation error below 0.5 %.

Table C1 LUT structure for TOA radiance calculations.