1. Introduction

Vibrational Raman scattering (VRS) in water is an inelastic scattering process occurring when photons interact with water molecules. It leads to filling-in of Fraunhofer lines [1] and a general spectral redistribution of the radiation in the ocean and in water-leaving radiance. VRS needs to be taken into account to accurately determine the underwater light field and to exploit this information in oceanic remote sensing applications, especially if the spectral resolution of measurements is high. Not considering oceanic Raman scattering in radiative transfer modeling affects inversion models on hyperspectral in situ data for determining inherent optical properties (IOPs) [2] and atmospheric trace gas retrievals on hyperspectral top-of-atmosphere (TOA) radiances [3]. With upcoming ocean colour satellite missions planned to record radiances at higher spectral resolution (e.g., PACE, https://pace.gsfc.nasa.gov/), VRS is going to play an increasingly important role for satellite-based ocean color retrievals.

The VRS signal in water-leaving radiances carries information about the number of inelastic scattering events in the water and consequently the number of photons in the ocean. Since the VRS contribution from the atmosphere can be considered negligible [4], the VRS signal in TOA radiances can be exploited to estimate light levels in the ocean [5]. Global data sets with good spatio-temporal resolution that characterize the underwater light field are needed for estimating primary productivity from satellite derived chlorophyll-a concentrations [6] and heat budget [7]. At present, planar photosynthetically available radiation (PAR) above surface, ultraviolet (UV) irradiance above surface, and diffuse attenuation coefficient at one wavelength are operationally provided from earth observations [8]. An independent global data set could be provided by exploiting the VRS signal in satellite-measured TOA radiances to spectrally derive the light availability and diffuse attenuation coefficient in the ocean.

Several hyperspectral satellite sensors designed for the retrieval of atmospheric trace gases have been used not only for quantifying the VRS contribution to TOA radiances and assessing its retrieval-influence but also for inferring information about environmental parameters. Using VRS, information on cloud pressure [9] and ocean parameters linking VRS to chlorophyll-a concentration and light penetration depth [9–11] have been derived.

These earlier studies investigated VRS in the UV whereas more recent publications analyzed VRS in the visible wavelength region to directly investigate the light used for photosynthesis. Dinter et al. [5] derived the light availability in the ocean from VRS as the depth and wavelength integrated scalar irradiance in the wavelength range from 390 to 444.5 nm. Through coupled ocean-atmosphere radiative transfer modeling, they related the VRS signal in TOA radiances to light availability in the ocean which was built into a look-up-table (LUT). Using Differential Optical Absorption Spectroscopy (DOAS), VRS was first retrieved from TOA radiances measured by the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) and converted to the amount of light in the ocean using the LUT. As evaluation of their light availability product, the diffuse attenuation coefficient is additionally determined using radiative transfer modeling. For one month of data, it was compared to the diffuse attenuation coefficient at a wavelength of 490 nm K_d(490) from the merged ocean color product GlobColour (http://www.globcolour.info/).

The connection between phytoplankton biomass and the VRS signal is more straightforward in the visible than in the UV, since this is the spectral region where chlorophyll-a absorbs. VRS retrieved in the visible with DOAS from TOA radiances measured by SCIAMACHY has been used as a proxy for the light penetration depth to calculate chlorophyll-a concentrations of different phytoplankton functional types (PFTs) [12].

SCIAMACHY measurements are only available until April 2012. Newer satellite sensors with improved temporal coverage and spatial resolution can be explored. The Ozone Monitoring Instrument (OMI) [13] provides data with significantly improved spatial resolution and daily global coverage from 2004 to date. The Global Ozone Monitoring Experiment 2 (GOME-2) [14] provides global coverage within one and a half days and long data time series being an operational mission with the first sensor launched in 2006 on Metop-A, the second sensor in 2012 on Metop-B, and the third sensor in 2018 on Metop-C. So far, the VRS signature has not been explicitly investigated in TOA radiances measured by OMI or GOME-2, especially not in the visible wavelength region. For both sensors, spectral structures caused by VRS in TOA radiances were pointed out to be interfering with atmospheric signals in various trace gas DOAS retrievals although often not explicitly accounted for [15–18]. Merely, vibrational and rotational Raman scattering have been simultaneously included in a model in the UV to retrieve cloud pressure for OMI [19]. Like SCIAMACHY, OMI and GOME-2 can potentially be used for phytoplankton functional type retrievals. DOAS retrievals on hyperspectral satellite data generally require an estimation of the light penetration depth to calculate PFT chlorophyll-a concentrations. Adequate retrieval of VRS from OMI and GOME-2 are therefore required. The approach developed by Dinter et al. [5] to derive the amount of light in the ocean from VRS can also be transferred to other sensors such as OMI and GOME-2. Their study was a feasibility study showing one month of results. Combining VRS from different satellite sensors could give a time series of global light availability since 2002.

This study therefore has two objectives. First it is investigated, if VRS can be found in TOA radiances recorded by OMI and GOME-2 and be adequately retrieved using DOAS. Second it is investigated, if and how a VRS derived product, in this case the diffuse attenuation coefficient, can be merged from different hyperspectral satellite sensors. The study focuses on the diffuse attenuation coefficient as VRS derived product instead of the light availability. For this quantity, better possibilities exist to compare with established ocean color products from multispectral sensors. It is also easier to intercompare diffuse attenuation coefficients from different satellite sensors with differing overpass times and viewing geometries, since the diffuse attenuation coefficient is regarded as a quasi-inherent optical property of ocean water [20].

Diffuse attenuation coefficients are calculated from VRS retrieved from TOA radiances measured by SCIAMACHY, OMI, and GOME-2 closely following the approach presented by Dinter et al. [5]. The quality of the resulting diffuse attenuation coefficient product is evaluated by comparison with the diffuse attenuation coefficient at 490 nm provided by the Ocean Colour Climate Change Initiative (OC-CCI) data set for nearly one year of data. The same time series data are explored in different Longhurst biogeochemical provinces in order to illustrate their broader comparability and possibility for a merged product.

2. Instrumentation and methods

2.1. Satellite sensors OMI, SCIAMACHY, and GOME-2

VRS signals are retrieved from three different hyperspectral satellite sensors, the Ozone Monitoring Instrument (OMI), the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY), and the Global Ozone Monitoring Experiment 2 (GOME-2) on Metop-A. All three sensors are on satellites in sun-synchronous orbits and have the mission goal to monitor atmospheric trace gases. GOME-2 and SCIAMACHY have a more similar instrument design compared to OMI, since they are both whisk-broom instruments whereas OMI is a push-broom sensor.

OMI [13] onboard the satellite Aura operated by NASA was launched on 15 July 2004 and is still in operation. It has a local overpass time at the equator at 1:45 pm in ascending node. OMI is designed as a nadir-viewing imaging spectrometer with 60 across track ground pixels that are measured simultaneously. Ground pixel size is 13 km by 24 km at nadir becoming significantly larger at the swath edges. With a swath of 2600 km, OMI achieves daily global coverage. It records spectra in the UV/Vis from 264–504 nm. Level 1b data (collection 003) from the visible channel (350 to 500 nm) were used for the VRS retrievals. The spectral resolution of this VIS channel is 0.63 nm. OMI measures the extraterrestrial irradiance spectrum and the upwelling earthshine spectrum. Daily measured OMI irradiance spectra are of poor quality due to low signal. For the OMI retrieval, a solar irradiance spectrum averaged over all irradiance measurements recorded in the year 2005 was used as reference [21].

The ESA instrument SCIAMACHY [22] onboard ENVISAT collected data from 2002 to April 2012. Local overpass time at the equator was 10:00 am in descending node. The instrument had a spatial resolution of 30 km by 60 km for nadir ground pixels. With a swath width of 960 km and alternating nadir-limb mode, global coverage was reached within six days. Level 1b data (version 7.04.) from cluster 15 (424 to 525 nm) of channel 3 with a spectral resolution of 0.47 nm were used for the VRS retrievals. SCIAMACHY daily recorded solar spectral irradiances.

The first GOME-2 [14] is hosted by the Metop-A satellite operated by EUMETSAT. It was launched in 2006 providing data from January 2007 to date. The spatial resolution was 40 km by 80 km for nadir ground pixels until 15 July 2013 when it was changed to 40 km by 40 km at reduced coverage. The sensor has a swath width of 1920 km and reaches global coverage in 1.5 days. Local overpass time at the equator is 9:30 am in descending node. Level 1b data (version 5.3) from visible channel 3 (397 to 604 nm) were used for the retrieval which have a spectral resolution of 0.51 nm. Solar spectral irradiances are recorded daily.

2.1.1. Ground pixel selection

Only ground pixels with solar zenith angle (SZA) smaller than 70° were used. Reflectance thresholds and operationally provided cloud products were used to filter cloud, glint, and ice affected ground pixels. An average reflectance R over the wavelength window used in the VRS fit was calculated as R = πI_TOA/(F₀ · cos(SZA)) using the spectrally averaged TOA intensity I_TOA and spectrally averaged solar irradiance F₀ as input. Reflectance thresholds were set to around 0.3. OMI and GOME-2 data were additionally cloud filtered using the O2-O2 [23] and the FRESCO [14] cloud products, respectively. Cloud fractions were set to 0.05. Backscans were excluded in SCIAMACHY and GOME-2 data. OMI suffers from a row anomaly which varies with time (http://projects.knmi.nl/omi/research/product/rowanomaly-background.php). It is assumed to be caused by a blocking object in front of the nadir port leading to changes in radiation level and wavelength shift. It first appeared in June 2007 affecting CCD rows 53 and 54 (0-based). From May 2008 on, the number of affected rows increases. Operationally provided flagging was used to filter affected ground pixels [24].

2.1.2. Observation period

The year 2007 is chosen as time period for our investigations. This year appears to be optimal for an intercomparison of the three hyperspectral sensors since all three sensors were measuring and OMI measurements are not severely affected by the row anomaly yet. GOME-2 measurements do not exist for the full month of January. For simplicity, the entire analysis for all three sensors was performed for 01 February 2007 to 31 December 2007.

2.2. Vibrational Raman scattering

Vibrational Raman scattering occurs in liquid water when vibrational modes of the water molecules are excited by inelastic scattering with photons. The wavelength of the photon is shifted in this process. In this work, we considered only the shift to a longer wavelength (Stokes line) around which radiation is emitted in a broad band. Based on experimental work by Walrafen [25], Haltrin and Kattawar [26] formulated a function to describe the broad band VRS emission by a combination of four Gaussian functions. In terms of wavelength, the shift and broad band width is wavelength dependent. In the spectral range from 400 nm to 700 nm, which phytoplankton use for photosynthesis, the light that has been Raman scattered in water originates from the wavelength region 352 to 567 nm and the emission band width ranges between 31 and 41 nm. This study evaluated the VRS emission in a spectral window from 450 to 497 nm which corresponds to an excitation wavelength range from 390 to 426 nm.

VRS leads to filling-in of Fraunhofer lines. It can be treated as a pseudo-absorber. Following Vountas et al. [3], a VRS pseudo-absorption, also called VRS compensation spectrum σ_VRS, can be defined as

2.3. Differential Optical Absorption Spectroscopy

Differential Optical Absorption Spectroscopy (DOAS) is a technique commonly used for the retrieval of atmospheric trace gases by distinguishing their high frequency absorption features [27]. The DOAS method has been extended for investigating oceanic parameters. The amount of VRS [10], light availability [5], chlorophyll-a concentration of different phytoplankton functional types [11, 28], and chlorophyll-a fluorescence [29, 30] have been successfully retrieved from SCIAMACHY and partly (chlorophyll-a fluorescence) from GOME-2 measurements.

Starting from Beer-Lambert’s law, the optical depth τ is separated into broadband and high frequency components by introducing a low-order polynomial.

In this study, the diffuse attenuation coefficient in the ocean was derived from the VRS fit factor by comparison with theoretical VRS fit factors from modeled data. Details on the comparison and radiative transfer modeling are elaborated below. Settings for the VRS DOAS fit were chosen based on the latest SCIAMACHY VRS DOAS fit results as in Losa et al. [12]. Cross sections for following components were included in the DOAS VRS fit

Although no weighting function DOAS (WF-DOAS) was used in this study, the ocean weighting function from Dinter et al. [5] was kept as cross section in the retrieval. Fit factors for wf_oc were not used for deriving any oceanic quantities, so their magnitude was not of importance in this study. The advantage of using wf_oc is that it combines the effect of absorbing oceanic constituents and liquid water in a single spectrum, avoiding the fit of multiple correlating absorption cross sections.

Pseudo-absorption cross sections for RRS and VRS as well as the ocean weighting function were calculated (see section 2.5) for each OMI CCD row separately and fitted for each row specifically. The row-specific solar spectrum was used as background spectrum for each CCD row. Daily recorded solar spectra were used as background spectrum for GOME-2 and SCIAMACHY (Azimuth Scan Mirror).

As in Losa et al. [12], instrumental artifacts were characterized by analyzing residual structures over the Sahara region [37]. A DOAS fit was made over a small region in the Sahara (2°E to 6°E, 23.5°N to 34.5°N) where no VRS signal is expected. In this pre-fit, only atmospheric components were considered and the resulting residuals were analyzed using a principal component analysis (PCA). The first two eigenvectors from this PCA were calculated on a monthly basis and included as cross sections in the DOAS fit. This method was only used for SCIAMACHY and GOME-2. No correction of instrumental effects through eigenvectors was used for OMI. As a test, the first two eigenvectors from a PCA over the Sahara region were calculated as for SCIAMACHY and GOME-2 for each CCD row separately and included in the VRS fit, however, no significant improvement of the fit was obtained. The global VRS fit factors calculated with and without eigenvectors had a correlation of r=0.99. To reduce the amount of retrievals and computation time, we decided not to include eigenvectors in the VRS fit for OMI.

The wavelength window for the fit was set to 450 nm to 497 nm which is slightly larger than in Losa et al. [12] and significantly shorter than in Dinter et al. [5]. The upper boundary is limited by OMI’s visible channel. The fit window was set the same for all three sensors, since this simplifies comparison between the sensors and merging of the data sets.

The only significant difference in the DOAS settings used here in comparison to Losa et al. [12] is that a second order polynomial was used instead of a third order polynomial. Due to the modifications in the RTM settings for calculating the pseudo-absorption cross sections (see section 2.5), a second order polynomial quickly proved to be more optimal.

2.4. Connection between VRS, light availability, and diffuse attenuation

This study exploits the VRS signal to derive the light availability in the ocean. An increased number of photons in the ocean leads to an increased number of inelastic scattering events. Since the VRS contribution from the atmosphere is negligible [4], the VRS signal in the TOA radiances is closely linked to the light availability in the ocean. This section summarizes the approach presented in Dinter et al. [5] on how to derive light availability and diffuse attenuation coefficient from VRS. They formulated a relationship between the change in VRS signal at TOA and the amount of radiation energy in the ocean. For a given wavelength region, the amount of radiation energy in the ocean ū is the depth and wavelength integrated scalar irradiance Ē₀ divided by the speed of light c given as

As pointed out by Dinter et al. [5], it is difficult to evaluate the quality of the resulting Ē₀ data set, since no comparable ocean color satellite product for the depth integrated scalar irradiance in the investigated wavelength range from 390 to 426 nm exists. However, the diffuse attenuation coefficient at 490 nm K_d(490) is often retrieved from ocean color satellite data. It is closely linked to the light availability in the ocean and was used, similar to Dinter et al. [5], for comparison in this study.

The diffuse attenuation coefficient can be calculated as a mean value over distant depths z₁ and z₂ from the change in downwelling irradiance E_d(z) [39]

2.5. Radiative transfer model simulations

Modeled TOA radiances were used to calculate VRS pseudo-absorption cross sections and to perform comparative DOAS retrievals whereas underwater fluxes were used to calculate the light availability and diffuse attenuation coefficient. The model should accurately describe radiative transfer processes, especially inelastic processes, in the atmosphere and in the ocean at high spectral resolution matching the spectral resolution of the satellites of about half a nanometer. We used the ocean-atmosphere coupled RTM SCIATRAN [41–43] version 4.0.8 to model various case I scenarios. The theoretical description of VRS is based on the above described formulation of VRS by Haltrin and Kattawar [26]. Correct implementation of the VRS in SCIATRAN has been evaluated by comparison with other radiative transfer models and experimental data from satellite, ship-based, and underwater instruments [43].

TOA radiances were modeled for 23 different case I scenarios with chlorophyll-a concentrations ranging between 0 and 30 mg/m³. A standard case I model was used based on Morel and Maritorema’s [44] parameterization for chlorophyll and CDOM absorption. Model settings follow Dinter et al. [5] with some changes:

Instrument-specific extraterrestrial solar spectra were used for each sensor for the TOA radiance calculations, since spectral alignment is very important for the DOAS retrieval. For SCIAMACHY, a solar irradiance measured by SCIAMACHY was used [52]. For OMI, an average over solar irradiance spectra recorded by OMI in the year of 2005 was used as mentioned above [21]. For GOME-2, SCIATRAN calculations were made at high spectral resolution of 0.01 nm using a Kurucz solar spectrum [53]. Resulting radiances were convolved to instrument spectral resolution using the GOME-2 slit function. Underwater fluxes were calculated instrument-unspecific at a spectral resolution of 0.5 nm also using the Kurucz solar spectrum since they should be insensitive to the exact spectral resolution.

Geometry settings were chosen as follows:

DOAS fit settings for the retrieval of theoretical VRS fit factors from the modeled TOA radiances were the same as for the retrieval on satellite radiances (see section 2.3) except for atmospheric cross sections. Only ozone was included as trace gas in the modeled atmosphere, so only the ozone absorption cross section was used when fitting modeled TOA radiances.

2.6. Retrieval sensitivity

The LUT used to convert VRS fit factors to K_d(λ_ext) was kept simple with SZA as only additional input parameter. In the following, the sensitivity of our retrieval with respect to model assumptions on AOD and CDOM absorption was investigated. Four new sets of model simulations were calculated. In two of these sets, AOD was increased or reduced to 0.2 or 0.05 at 550 nm, respectively. The other two sets were calculated with increased or reduced CDOM absorption. The coefficient for CDOM absorption in the Morel and Maritorema parameterization [44] was changed from 0.2 to 0.4 or 0.1, respectively. VRS retrieval as described in section 2.3 was performed for each of these four synthetic TOA radiance sets. Retrieved VRS fit factors were then converted to K_d(λ_ext) using the LUT. For each set, resulting K_d(λ_ext) for each scenario were compared to the expected $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ in this scenario calculated from flux simulations with the four modified settings, i.e. higher or lower AOD and higher or lower CDOM absorption. Figure 1 shows the relative deviation of derived K_d(λ_ext) from expected $K_d^{\text{exp}}({\lambda }_{\text{ext}})$, $(K_d({\lambda }_{\text{ext}})-K_d^{\text{exp}}({\lambda }_{\text{ext}}))/K_d^{\text{exp}}({\lambda }_{\text{ext}})$, as function of $K_d^{\text{exp}}({\lambda }_{\text{ext}})$. The analysis is shown in a $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ range from 0 to 0.15 m⁻¹, since our retrieval is only applicable in this range due to too low VRS signals for higher attenuating waters as discussed in section 4. For most SZAs, an increase in AOD leads to an overestimation of $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ whereas a reduction leads to an underestimation. The retrieval is especially sensitive to AOD when $K_d^{\text{exp}}({\lambda }_{\text{ext}})<0.01\hspace{0.17em}\hspace{0.17em}{\text{m}}^{-1}$. For these extremely clear waters, deviations are significantly larger than in the range 0.01 m⁻¹ to 0.15 m⁻¹. Here, deviations from $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ are nearly constant and less than −5% in case of reduced AOD and less than +20% in case of increased AOD. Higher or lower CDOM absorption in the RTM simulations leads to underestimated or overestimated $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ for most scenarios, respectively. Deviations increase with increasing $K_d^{\text{exp}}({\lambda }_{\text{ext}})$. In the investigated $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ range, derived K_d(λ_ext) deviate from $K_d^{\text{exp}}({\lambda }_{\text{ext}})$ by less than +15% in case of a reduced CDOM absorption and between −10% and +5% in case of an increased CDOM absorption.

 

Fig. 1 Retrieval sensitivity analysis: Relative deviation of derived K_d(λ_ext) from expected $K_d^{\text{exp}}({\lambda }_{\text{ext}})$, $(K_d({\lambda }_{\text{ext}})-K_d^{\text{exp}}({\lambda }_{\text{ext}}))/K_d^{\text{exp}}({\lambda }_{\text{ext}})$, in modeled open ocean scenarios with reduced (solid lines) and increased (dashed lines) (a) AOD and (b) CDOM absorption. Compared to the average global open ocean scenario used for the retrieval and described in section 2.5, AOD and CDOM absorption were reduced and increased by a factor of two. Deviations are shown for a range of SZAs.

Download Full Size | PDF

The retrieval sensitivity to AOD and CDOM absorption is rather low, but not negligible. Regarding AOD, the method could be improved by adding another dimension to the LUT and taking AOD from a data base as third input parameter. Wolanin et al. [29] calculated the average polynomial, resulting from the DOAS fit, to evaluate the brightness of each ground pixel. They used this estimate of all broad band effects as input parameter in a LUT. This approach could be implemented for the method presented here to reduce the sensitivity to broad band effects such as AOD and CDOM absorption. We decided to keep the LUT simple, since it turned out that the noise introduced by mismatches between model scenarios and observed scenes is a secondary effect. As discussed in section 4, spectral features caused by instrumentation have a larger effect on the results.

Other model settings additionally influence the retrieval. Further sensitivity analyses on the choice of phytoplankton spectra and vertical profiles of chlorophyll-a concentration can be found in Dinter et al. [5]. The approaches are not one-to-one comparable, but the sensitivity analyses presented there give an impression what to expect in our particular approach.

2.7. OC-CCI as reference data set

Diffuse attenuation coefficients were taken from the release version 3.1 of the OC-CCI [54] of the European Space Agency ESA (http://www.esa-oceancolour-cci.org/) for evaluating the quality of the diffuse attenuation coefficients derived from VRS fit factors. OC-CCI provides the diffuse attenuation coefficient at 490 nm merged from three different multispectral sensors MODIS, MERIS, and SeaWiFS.

Since OC-CCI diffuse attenuation coefficients are given at a wavelength of 490 nm $K_d^{\text{OC}}(490)$ diffuse attenuation coefficients from the OC-CCI data set were converted to the wavelength range of our retrieval output. SCIATRAN flux calculations were performed for wavelengths up to 500 nm to investigate K_d’s wavelength dependence. The modeled diffuse attenuation coefficient in the excitation wavelength region $K_d^{\text{mod}}({\lambda }_{\text{ext}})$ was related to the modeled diffuse attenuation coefficient at 490 nm $K_d^{\text{mod}}(490)$. $K_d^{\text{mod}}(490)$ was calculated by integrating fluxes between 485 nm and 495 nm. A 10 nm integration width was chosen here to best resemble the band width of multispectral ocean color sensors used to calculate $K_d^{\text{OC}}(490)$. A linear relationship between $K_d^{\text{mod}}(490)$ and $K_d^{\text{mod}}({\lambda }_{\text{ext}})$ was found,

3. Results

Here, we present the results of the DOAS fits and diffuse attenuation coefficient derivation. A separate discussion follows in section 4.

3.1. Global vibrational Raman scattering DOAS fit

Fit residuals and global distributions of VRS fit factors were analyzed for May 2007 as a first quality check of the DOAS fit. Examples of the differential optical depth (second term in Eq. (2)) for VRS and the ocean weighting function accounting for absorption by phytoplankton and liquid water are shown in Fig. 2. Results are shown for ground pixels from the South Pacific Gyre where a strong VRS signal is expected. The residual has about the same magnitude as the reference indicating that the signals from VRS and other ocean parameters are found moderately well by the DOAS fit in the TOA radiances. The mean absolute residual of all cloud-free pixels in a small region in the South Pacific Gyre (22–30°S, 115–125°E) during the month of May 2007 amounts to 0.0009 for SCIAMACHY, 0.0006 for GOME-2, and 0.0013 for OMI.

 

Fig. 2 Resulting differential optical depth for (a) VRS and (b) other oceanic parameters from the spectral fit on OMI (upper panel), SCIAMACHY (center panel), and GOME-2 (lower panel) TOA radiances. The solid lines are the references which are the differential cross sections multiplied by the retrieved fit factor for the corresponding ground pixel. The dashed line is the retrieval result showing the reference plus the overall fit residual. Results are from ground pixels in the South Pacific Gyre obtained on 19 May 2007 (OMI: 24.29°S, 235.03°E; SCIAMACHY: 24.25°S, 235.03°E; GOME-2:24.68°S, 235.34°E).

Download Full Size | PDF

Figure 3 shows monthly gridded global VRS fit factors for May 2007 from (a) SCIAMACHY, (b) GOME-2, and (c) OMI. For comparison, Fig. 3(d) shows the diffuse attenuation coefficient $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ from the OC-CCI data set for the same month. Due to the definition of the VRS pseudo absorption cross section, VRS fit factors are usually negative. We defined them as positive here for easier understanding by multiplying with −1. Now, high VRS fit factors correspond to high VRS signals. VRS fit factors show similar global patterns between the different sensors. These patterns are in agreement with regions generally known to have high and low amount of light in the water. Clearly, there is an anti-correlation between high/low VRS fit factors corresponding to regions with high/low diffuse light attenuation for all three instruments.

 

Fig. 3 Global VRS fit factors retrieved using DOAS on TOA radiances from (a) SCIAMACHY, (b) GOME-2, (c) OMI for May 2007. (d) $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ from OC-CCI. (e) Global normalized frequency of VRS fit factors excluding land and coast pixels for May 2007 for OMI (blue, horizontally striped), SCIAMACHY (green, dotted), and GOME-2 (red, diagonally striped). For OMI, only ground pixels of the two nadir rows in the middle of the CCD are included (row 29 and 30 starting to count at 0). f) Average VRS fit factor for OMI’s two nadir CCD rows (see (e)) as retrieved from synthetic TOA radiances simulated with SCIATRAN for different diffuse attenuation coefficients K_d(λ_ext). Only results for every second simulated SZA are shown for clarity.

Download Full Size | PDF

Figure 3(e) shows the normalized frequency of global VRS fit factors for the month of May 2007 excluding land and coast pixels. Only ground pixels from OMI’s nadir CCD rows 29 and 30 were considered roughly enabling comparison with distributions from GOME-2 and SCIAMACHY. VRS fit factor ranges are significantly different between the three sensors which is expressed in a difference in mean value of 0.81 for OMI, 0.44 for SCIAMACHY, and 0.49 for GOME-2 and standard deviation of 0.35 for OMI, 0.31 for SCIAMACHY, and 0.22 for GOME-2. For all three instruments, negative VRS fit factors are present in the global distribution. Among the three sensors, SCIAMACHY has the highest amount of negative values. About 4% of all ground pixels have negative VRS fit factors. Figure 3(f) shows the VRS fit factor $S_{\text{VRS}}^{\text{mod}}$ as function of K_d(λ_ext) as calculated from SCIATRAN simulations for different SZAs. Since SCIATRAN simulations were calculated instrument-specific, i.e., for each sensor individually, this relationship varies among the sensors. Displayed is the averaged result for OMI’s two center CCD rows 29 and 30 (0-based) for every second SZA that was simulated. The relationship $K_d({\lambda }_{\text{ext}})\leftrightarrow (S_{\text{VRS}}^{\text{mod}},\text{SZA})$ looks very similar to the one for OMI’s nadir rows if SCIAMACHY and GOME-2 settings were used. VRS fit factors retrieved on instrument-specific simulated TOA radiances for SCIAMACHY and GOME-2 deviate by less than 3% from those in Fig. 3(f) for the lowest simulated diffuse attenuation coefficients.

VRS fit factors are not directly comparable between the sensors, mostly, because they depend on SZA which differs among the sensors due to different local overpass times, and, minorly, because of different instrument design. A one-to-one comparison of the VRS fit factors between the different hyperspectral sensors might not be meaningful. But, comparing the theoretically expected VRS fit factor range as presented in Fig. 3(f) with the global VRS fit factor frequencies in Fig. 3(e), differences in fit factor frequencies are expected, but appear to be large. As seen in Fig. 3(f), the model predicts only positive VRS fit factors over the entire simulated K_d(λ_ext) range, but negative VRS fit factors are seen in the satellite data. This disagreement between retrieval results on model and satellite data indicates a correlation between cross sections as discussed in section 4.1.

3.2. Agreement between model and satellite data

A thorough comparison of the VRS fit performance between the three hyperspectral sensors and quality estimation of each single sensor can only be made when VRS fit factors are converted to a physical quantity, here, the diffuse attenuation coefficient in the excitation wavelength region. A correct conversion from VRS fit factor to diffuse attenuation coefficient requires that the VRS fit on modeled data accurately represents the VRS fit on real satellite data. The robustness of the conversion can be estimated from comparing the hyperspectral satellite derived diffuse attenuation coefficient to an independent data set.

Starting from the retrieved VRS fit factors, the diffuse attenuation coefficient for the excitation wavelength range K_d(λ_ext) was derived using a LUT built up by SCIATRAN calculations (details in section 2.5) relating $S_{\text{VRS}}^{\text{mod}}$, SZA, and K_d(λ_ext) based on the relationship shown exemplarily in Fig. 3(f). Figure 4 shows retrieved K_d(λ_ext) from (a) SCIAMACHY, (b) GOME-2, and (c) OMI in comparison to the wavelength-converted $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ using $K_d^{\text{OC}}(490)$ from OC-CCI as input (details in section 2.7). The scatter plot contains data from 01 Feb 2007 to 31 Dec 2007. The comparison was calculated on a daily basis using data gridded to 0.5°. K_d(λ_ext) derived from SCIAMACHY and GOME-2 are much larger than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ from the OC-CCI data set. Linear total least squares regression shown by the dash-dotted black line yielded a slope of 3.02 for SCIAMACHY and 2.05 for GOME-2 which indicates that their K_d(λ_ext) are roughly 3 and 2 times, respectively, larger than OC-CCI values. OMI shows a different behavior. K_d(λ_ext) derived from OMI is smaller than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ over the entire data range expressed in a slope of 0.30. Regression statistics are summarized in Table 1.

 

Fig. 4 Scatter plots show K_d(λ_ext) from VRS obtained with the DOAS fit settings as in Losa et al. [12] from (a) SCIAMACHY, (b) GOME-2, and (c) OMI in comparison to $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ which was wavelength converted from OC-CCI $K_d^{\text{OC}}(490)$. The solid red line is the 1:1 line. The dash-dotted black line is a linear total least squares regression through the global data set. Daily gridded data with 0.5°grid cell size were used as input. The comparison contains all data from 01 Feb 2007 .

Download Full Size | PDF

Table 1. Slope, intercept, and Pearson correlation of the total linear least-squares regression between K_d(λ_ext) and $K_d^{\mathbf{OC}}({\lambda }_{\mathbf{ext}})$ comparing 11 months of diffuse attenuation coefficients derived with the original DOAS settings on the full data set (left column: orig., full), with the modified DOAS settings on the full data set (middle column: mod., full), and with modified (SCIAMACHY and GOME-2) or original (OMI) DOAS settings on the reduced data set (K_d(λ_ext) and $K_d^{\mathbf{OC}}({\lambda }_{\mathbf{ext}})<0.15\hspace{0.17em}\hspace{0.17em}{\mathbf{m}}^{-1}$) (right column: mod.,red. or orig.,red.).

View Table

In summary, the slopes found in comparison with the OC-CCI data differ significantly from one and indicate a disagreement between the VRS fit factors found by fitting global satellite data and those found by fitting model data. High slopes for SCIAMACHY and GOME-2 express a lower VRS fit factor range in satellite data than in model data. Conversely, underestimation of K_d(λ_ext) by OMI indicates a higher VRS fit factor range than the model predicts.

3.3. Modification of the VRS fit

Modifications of the VRS fit were tested to improve the agreement between VRS fit factor ranges from fits on model and satellite data. For SCIAMACHY, including zeta and eta functions was tested which were taken from calibration key data [60] for correcting polarization structures. For GOME-2 and SCIAMACHY, the number of eigenvectors accounting for instrumental artifacts were reduced from two to one eigenvector in the DOAS fit. Since no eigenvectors were used in the DOAS fit on OMI data, no further modifications of the DOAS fit were tested for OMI in this study. We merely tested if reducing the number of CCD rows to the inner 20 CCD rows closest to nadir improves the results. Looking at the 11 month comparison between K_d(λ_ext) and $K_d^{\text{OC}}({\lambda }_{\text{ext}})$, a slightly different slope was found, but no change in correlation of the data set. Since no significant improvement was achieved, diffuse attenuation coefficients from all CCD rows were kept in the data set.

For GOME-2, reducing the number of cross sections from two to one eigenvector significantly improved the derived diffuse attenuation coefficient product. Changing this fit setting, shifted the globally retrieved VRS fit factor distribution towards higher values. For example for May, the global distribution has now a higher mean value of 0.56 while the standard deviation stayed the same.

For SCIAMACHY, best results were found when two eigenvectors were included for correcting instrumental effects and additionally the zeta function was fitted. VRS fit factors generally increased and the number of negative VRS fit factors decreased. Similar behavior to GOME-2 was found for the example month May. The global distribution has a similar standard deviation of 0.30, but a higher mean value of 0.51 compared to the old DOAS settings without zeta function.

Figure 5 shows K_d(λ_ext) derived from VRS fit factors resulting from the modified DOAS fit on (a) SCIAMACHY and (b) GOME-2 in comparison to the wavelength-converted OC-CCI diffuse attenuation coefficient product $K_d^{\text{OC}}({\lambda }_{\text{ext}})$. Regression results can be found in Table 1 (middle column). Considering all global ground pixels, linear total least squares regression yielded slopes of 1.73 for SCIAMACHY and 1.24 for GOME-2. For large values, the spread in the scatter plots becomes large for both instruments. Especially for SCIAMACHY, two pronounced wings are observed in the scatter plot, where values from OC-CCI are very large when values derived from SCIAMACHY are small and vice versa. We decided to restrict the retrieval of K_d(λ_ext) from hyperspectral satellite data to waters with diffuse attenuation coefficients below 0.15 m⁻¹ given by $K_d^{\text{OC}}({\lambda }_{\text{ext}})$. Additionally, values of K_d(λ_ext) above 0.15 m⁻¹ were discarded. This filtering left only data outside the grey area in Fig. 5. Linear total least squares regression is shown as dashed magenta line. Reducing the data set increased the correlation between K_d(λ_ext) and $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ from 0.27 to 0.66 for SCIAMACHY and from 0.25 to 0.66 for GOME-2. The slopes changed for SCIAMACHY from 1.73 to 1.34 and for GOME-2 from 1.24 to 1.25. The same filtering approach was applied to results for OMI displayed in Fig. 5(c). Correlation increased here from 0.35 to 0.65 and the slope increased from 0.30 to 0.77. Regression results for the reduced data sets can also be found in Table 1 (right column).

 

Fig. 5 Scatter plots showing K_d(λ_ext) from VRS obtained with the modified DOAS fits from (a) SCIAMACHY, (b) GOME-2, and (c) OMI in comparison to $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ wavelength converted from OC-CCI $K_d^{\text{OC}}(490)$. The solid red line is the 1:1 line. The dash-dotted black line is a linear total least squares regression through the global data set and the dashed magenta line is a linear total least squares regression through the global data set where values above 0.15 m⁻¹ were filtered out (grey area). Daily gridded data with 0.5°grid cell size were used as input. The comparison contains all data from 01 Feb 2007 until 31 Dec 2007.

Download Full Size | PDF

All results presented below were obtained using modified VRS fit settings and data filtering as described in this section. The hyperspectral satellite K_d(λ_ext) was filtered in waters where $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ predicts values above 0.15 m⁻¹ by gridding both data sets to the same grid of 0.5° gridsize.

3.4. Derived global diffuse attenuation coefficients

Global K_d(λ_ext) maps derived from VRS fit factors obtained with modified DOAS settings are shown in Fig. 6 for (a) SCIAMACHY, (c) GOME-2, and (e) OMI after K_d(λ_ext) was filtered as describe in section 3.3. Systematic differences between the three sensors can already be deduced from these global maps. Diffuse attenuation coefficients derived from OMI VRS fit factors are globally lower than those derived from SCIAMACHY and GOME-2 VRS fit factors. Higher values in the northern hemisphere are the most prominent difference between SCIAMACHY and GOME-2 derived K_d(λ_ext).

 

Fig. 6 Global distribution of the diffuse attenuation coefficient at the excitation wavelength region K_d(λ_ext) derived from (a) SCIAMACHY, (c) GOME-2, and (e) OMI VRS fit factors gridded for the time period 01 Feb 2007 to 31 Dec 2007. Data was gridded to 0.5° grid cell size. Globally mapped distribution of the difference between the OC-CCI and the hyperspectral satellite data sets $K_d({\lambda }_{\text{ext}})-K_d^{\text{OC}}({\lambda }_{\text{ext}})$ for (b) SCIAMACHY, (d) GOME-2, and (f) OMI. In all six images, K_d(λ_ext) was filtered below 0.15 m⁻¹ and restricted to waters where $K_d^{\text{OC}}({\lambda }_{\text{ext}})<0.15\hspace{0.17em}\hspace{0.17em}{\text{m}}^{-1}$.

Download Full Size | PDF

To identify regional differences between the hyperspectral satellite and the OC-CCI products, the difference $K_d({\lambda }_{\text{ext}})-K_d^{\text{OC}}({\lambda }_{\text{ext}})$ is shown as global map in Fig. 6 for (b) SCIAMACHY, (d) GOME-2, and (f) OMI. As the pronounced difference between GOME-2 and SCIAMACHY in the northern hemisphere suggested, a difference between northern and southern hemisphere is observed for SCIAMACHY. In the northern hemisphere, K_d(λ_ext) is smaller than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ for most ground pixels, whereas in the southern hemisphere K_d(λ_ext) is larger than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ for most ground pixels. The observed spread in the scatter plot in Fig. 5(a) as two wings, where $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ is large when K_d(λ_ext) is small and vice versa, is mainly caused by a north-south trend in K_d(λ_ext) originating from a north-south trend in SCIAMACHY VRS fit factors. For GOME-2, no significant north-south trend is observed. In low latitudes, $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ is larger than K_d(λ_ext) in the Mauretanian, Benguela, and Malvinas upwelling regions where diffuse attenuation coefficients are the highest. For open ocean situations such as the Southern Ocean, GOME-2 derived attenuation coefficients tend to be smaller than the OC-CCI product. However, especially in the high latitudes in the north, spatial patterns are not so clear. As seen from the global K_d(λ_ext) map, averaged over one year, OMI derived K_d(λ_ext) are generally smaller than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$. Largest differences between the two data sets can be found in northern high latitudes.

3.5. Diffuse attenuation coefficient time series

Figure 7(a) shows mean diffuse attenuation coefficients in the excitation wavelength range for SCIAMACHY, GOME-2, OMI, and OC-CCI in five different Longhurst provinces [61] over the course of 2007, starting 01 Feb 2007. Mean K_d(λ_ext) were calculated as 7-day composites. Relative differences $(K_d({\lambda }_{\text{ext}})-K_d^{\text{OC}}({\lambda }_{\text{ext}}))/K_d^{\text{OC}}({\lambda }_{\text{ext}})$ between hyperspectral satellites and OC-CCI in the five provinces are shown in Fig. 7(b). In general, largest differences between K_d(λ_ext) and $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ are observed in the NADR and ARCT region where a significant amount of data points deviates by up to ±60%. Positive and negative deviations are observed for GOME-2 in all five provinces showing a seasonal behavior. E.g., large positive deviations up to +120% are observed in the winter months February, March, and December. During the other months, GOME-2 deviates from $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ by less than ±30% in the ISSG, ARAB, and BENG provinces and less than ±40% in the NADR and ARCT provinces. OMI and SCIAMACHY diffuse attenuation coefficients are mostly smaller than $K_d^{\text{OC}}({\lambda }_{\text{ext}})$. This trend becomes more striking with increasing latitude. Although SCIAMACHY and OMI generally have a significant bias, they tend to be more stable with time than GOME-2. E.g., for the ARCT province, OMI time series has a large bias, but is relatively stable in time in comparison to $K_d^{\text{OC}}({\lambda }_{\text{ext}})$. Deviations range between −40% to −65% for most data points while they range between −10% to −60% for SCIAMACHY and −40% to +40% for GOME-2 (excluding winter months).

 

Fig. 7 (a) Mean K_d(λ_ext) from SCIAMACHY, GOME-2, OMI, and OC-CCI and (b) relative deviation of SCIAMACHY, GOME-2, and OMI $K_d^{\text{sat}}({\lambda }_{\text{ext}})$ from OC-CCI $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ in five different Longhurst provinces [61] over the course of the year 2007. 7-day composites were calculated starting 01 February 2007. Naming convention follows the original naming by Longhurst et al. [61]. (i) ISSG - Indian S. Subtropical Gyre Province, (ii) ARAB - NW Arabian Upwelling Province, (iii) NADR - N. Atlantic Drift Province, (iv) BENG - Benguela Current Coastal Province, (v) ARCT - Atlantic Arctic Province.

Download Full Size | PDF

4. Discussion

4.1. Correlating cross sections in DOAS

The spectral signature of VRS was identified in all three satellite sensors OMI, SCIAMACHY, and GOME-2. Global VRS maps can be retrieved from TOA radiances measured by the three sensors and resulting global patterns agree with regions of high and low diffuse light attenuation. However, a disagreement in VRS fit factor range between model and satellite data was found. For nadir viewing geometry, VRS fit factors retrieved from satellite data were lower than predicted by the model in case of SCIAMACHY and GOME-2, whereas they were slightly higher than predicted by the model in case of OMI. This discrepancy indicates a correlation between cross sections in the DOAS fit.

The difficulty of separating VRS from liquid water absorption spectral features in DOAS is known. For circumventing this problem in phytoplankton DOAS retrievals, Bracher et al. [11] performed a PCA over clear water regions and included the first two eigenvectors in the DOAS fit instead of including a separate VRS cross section. Peters et al. [18] tried to experimentally determine a VRS cross section to improve DOAS NO₂ fits, but were not able to separate liquid water absorption from VRS. They obtained the best NO₂ results over water regions when using a combined spectrum. Since we are particularly interested in the VRS signal, we need it as separate cross section in the DOAS fit. In principle, a combined spectrum of VRS and liquid water absorption could also be used as a proxy for the diffuse attenuation coefficient in the ocean. However, the interpretation of such an approach is more difficult, since the VRS excitation is in a different wavelength region than the DOAS fit window where liquid water absorbs. It is not straightforward for which wavelength range the light availability should be calculated.

To a certain extent, the comparison with model results accounts for the correlation between the spectra, since the same DOAS fit is performed on modeled TOA radiances. A theoretical correlation between cross sections also impacts retrieved fit factors from modeled data. The VRS fit factors retrieved from satellite data are directly related to those retrieved from modeled data to derive the diffuse attenuation coefficient in our approach. Correlating spectra affect the model retrieval and the satellite retrieval in the same way and the direct comparison of the two should compensate for this effect. However, in real world applications with satellite spectra that are noisy and show spectral features caused by the instrumentation, model and satellite results are not one-to-one comparable.

Differences in VRS fit factors between the three sensors suggest that not only the correlation between cross sections, but primarily the instrument itself influences the fit results. Major differences are that SCIAMACHY and GOME-2 are whisk-broom sensors whereas OMI is a push-broom sensor, that SCIAMACHY and GOME-2 provide higher quality daily sun spectra than OMI for which a fixed sun spectrum has to be used, and that SCIAMACHY and GOME-2 have no polarization scramblers. All these differences in instrument design plus differences in light path in the instrument interior cause different spectral signatures in the recorded TOA radiances and sun spectra, e.g., caused by straylight in the instrument [62] or structures on the diffuser used for recording sun spectra [63]. The eigenvectors that were used in the DOAS fit in this study to account for instrumental spectral artifacts that are not removed during the calibration differ between GOME-2 and SCIAMACHY. The eigenvectors could also correlate with the VRS pseudo-absorption cross section and change the theoretically expected VRS fit factor. In general, SCIAMACHY and GOME-2 have a more similar instrument design in comparison to OMI. This may explain, why the retrieved VRS fit factor ranges are more similar for these instruments and why including eigenvectors in the DOAS fit for OMI hardly affects the result.

We want to note that the discrepancy found here between VRS results from simulations and measurements determined with the DOAS settings as in Losa et al. [12] does not mean that their results on SCIAMACHY measurements suffer from the same issues. For converting their fit factors into PFT chlorophyll-a concentrations, they use a ratio of fit factors for two quantities, i.e. PFT and VRS, while here only one fit factor is used in combination with model results. Ratios of fit factors are more robust providing that instrumental artifacts affect the two quantities similarly.

4.2. Implications of DOAS setting modifications

Modifications of the DOAS fit for GOME-2 and SCIAMACHY were tested and significantly improved resulting diffuse attenuation coefficients. Evaluating the results by comparison with wavelength-converted diffuse attenuation coefficients from OC-CCI data sets yielded slopes close to one in the scatter plots (Fig. 5). Through modifying the DOAS fit for GOME-2 and SCIAMACHY, we were able to bring together fit performance on modeled radiances with that on satellite-recorded radiances.

Including zeta and eta functions for polarization corrections in DOAS fits has been suggested by McLinden [64] and shown to improve DOAS fit results on SCIAMACHY radiances for various atmospheric trace gas retrievals [65]. In the VRS DOAS fit studied here, only including the zeta function worked best. For other retrieval targets, other combinations of zeta and eta function were found beneficial. E.g., for BrO vertical profile retrievals from SCIAMACHY limb measurements, including the eta function improved the fit, but including the zeta function destabilized the fit [66].

For OMI, DOAS fit settings were not modified. No significant improvement was achieved through CCD row filtering for the 20 middle CCD rows which are closest to nadir and have more optimal viewing geometries for ocean observations than the outer CCD rows. It can be concluded that for OMI the viewing angle dependence in the model accurately describes the situation found in measurements.

4.3. Restricting diffuse attenuation coefficient data sets

Large spreads in the scatter plots were observed for high diffuse attenuation coefficients, especially for OMI and SCIAMACHY, leading to wings in the scatter plot. Such a butterfly-shape has been already observed by Dinter et al. [5] when comparing one month of derived diffuse attenuation coefficients from SCIAMACHY to K_d(490) from GlobColour. They suggested the following reasons for the wings. The horizontal wing is possibly caused by cloud-contaminated SCIAMACHY ground pixels leading to an overestimation of K_d(λ_ext) and the vertical wing is caused by deviations of vertical chlorophyll-a profiles and specific chlorophyll-a absorption from those set in the model. They restricted their statistical analysis of the comparison to diffuse attenuation coefficients below 0.06 m⁻¹ and improved the correlation from 0.42 to 0.69. One-to-one comparison for SCIAMACHY results found here with those published in Dinter et al. [5] is not possible since DOAS fit settings, observed time period, and multispectral reference product used for the comparison differ. However, it can be stated roughly that, when comparing correlations found here and in their study, a similar degree of agreement was found, when K_d(λ_ext) from SCIAMACHY is compared to diffuse attenuation coefficients from multispectral sensor data sets.

We similarly restricted our retrieval for all three hyperspectral sensors to $K_d^{\text{OC}}({\lambda }_{\text{ext}})<0.15\hspace{0.17em}\hspace{0.17em}{\text{m}}^{-1}$ to improve data quality justifying and adding general suggestions here for reasons behind the observed spread in the scatter plots. Waters with high diffuse attenuation coefficients have a low VRS signal. The sensitivity of the VRS DOAS fit in these waters is therefore low and the fit is expected to be less accurate there [10]. OC-CCI diffuse attenuation coefficients by the Lee method [39] have better accuracies for waters with high diffuse light attenuation [67,68]. The VRS method is expected to be more sensitive to waters with low diffuse light attenuation.

Additionally, derived diffuse attenuation coefficients from the hyperspectral sensors were cut-off at 0.15 m⁻¹ even if OC-CCI predicts values below 0.15 m⁻¹. This filtering was necessary to improve the data quality sufficiently and was easy to apply. From a physical point of view it is more difficult to understand though: Considering OC-CCI values as correct, the VRS fit gives wrong results for filtered out ground pixels, although the VRS fit sensitivity should have been high. Partially this can be caused by noise in the recorded TOA radiances. Cloud-contamination of ground pixels can significantly reduce the VRS signal as already pointed out by Dinter et al. [5]. Other possible causes are a mismatch between average ocean scenario used in the model simulations and the true state of the ocean during the measurement, e.g. extremely high wind speeds or aerosol loads. By restricting derived diffuse attenuation coefficients to waters where OC-CCI predicts values below 0.15 m⁻¹ and additionally cutting off all remaining values above 0.15 m⁻¹, reduces the data sets by less than four percent. It is unfortunate that these scenes have to be removed from the data set, since biological activity is high and significantly contributes to the overall primary production in open oceans, however, valuable information on open ocean primary productivity may still be drawn from the remaining data points.

4.4. Spatial and temporal biases in diffuse attenuation coefficients

Despite the DOAS fit modifications, issues in the VRS DOAS fit remained. Possible origins of these issues and impacts on the quality of the derived diffuse attenuation coefficient became apparent when looking at spatial and temporal patterns of the derived diffuse attenuation coefficient.

For SCIAMACHY, a clear north-south trend was found in the data which is most likely caused by spectral features arising from polarization sensitivity of the instrument. Scatter geometries differ from North to South causing differences in spectral features due to polarization [69]. If these spectral features, included as eigenvectors in the fit, correlate with the VRS pseudo-absorption cross section, the north-south trend in these spectral features is passed on to the VRS fit factors. The zeta function was included to account for polarization. Including the zeta function in the fit reduced but did not fully remove the north-south trend.

For OMI, agreement with OC-CCI becomes worse with increasing latitude. This finding suggests that the model predicts a different SZA dependence of the VRS fit factor than found in the DOAS retrieval results on measured TOA radiances. Detailed investigations on this behavior have not been performed in this study.

For GOME-2, spread does not arise from spatial trends, but mainly from temporal variations in the fit factors. Seasonal variations as observed here have been identified in trace gas columns retrieved with DOAS from GOME and attributed to the diffuser plate causing spectral artifacts in the irradiances that correlate with absorber cross sections [63]. This effect of seasonal variation has also been observed for the weak signal of glyoxal in DOAS columns from GOME-2 measurements [16]. Using a constant solar spectrum instead of daily solar spectra could improve the situation for 2007, but is not a solution for the whole GOME-2 mission period due to significant instrumental drifts [17,70]. Another option as background spectrum (e.g. [70], [71]) is to use daily earthshine spectra measured by GOME-2 over a region where no VRS signal is expected, e.g., over the Sahara. We have quickly tested this option for the month of December 2007. Resulting diffuse attenuation coefficients were closer to OC-CCI compared to results obtained with daily solar spectra but still too high with a slope in the monthly scatter plot of 0.57 (filtered below 0.15 m⁻¹). Other correlating spectral features seem to be introduced that disturb the fit results. The seasonal variation could also be empirically corrected, e.g., through adjusting retrieved fit factors by forcing the fit factors in a certain region to a known value [72]. The seasonal variation dominates the time series in the five Longhurst provinces. To build up a consistent time series over the entire GOME-2 mission period and analyze trends in diffuse attenuation coefficients from this series, a correction for this variation is required.

Although we clearly saw that biases and seasonal variations in the three hyperspectral sensors dominate the time series in the Longhurst provinces, also a difference between OC-CCI data set and hyperspectral data sets arises due to a difference in data coverage in the five Longhurst provinces. Coverage was calculated for all five provinces by dividing the number of unique gridded ground pixels in one Longhurst province within one 7-day composite by the total amount of gridded ground pixels fitting into the Longhurst province. Largest differences in coverage between OC-CCI and the hyperspectral sensors were found for NADR and ARCT which are strongly impacted by cloud occurrence (e.g. [73]). In combination with a spatial inhomogeneity in biomass in these provinces due to the occurrences of phytoplankton blooms, this may explain, why the diffuse attenuation coefficients from hyperspectral sensors compare less well with the $K_d^{\text{OC}}({\lambda }_{\text{ext}})$ data set in these two provinces than in the other three provinces.

Due to the large ground pixel sizes of hyper-spectral satellite sensors and the associated inhomogeneity within a ground pixel, in combination with the rather complex retrieval method presented in this study, we expect diffuse attenuation coefficients derived from hyperspectral sensors to be associated with larger uncertainties and regard OC-CCI K_d(490) as a solid data set to evaluate our retrieval performance. However, also the OC-CCI diffuse attenuation coefficients are associated with uncertainties. $K_d^{\text{OC}}(490)$ was determined using the semianalytical method by Lee et al. [39]. This method is a two step process. First, inherent optical properties are derived with a Quasi-Analytical Approach from remote sensing reflectances [74]. Extensive radiative transfer simulations have been made to relate these inherent optical properties, i.e., absorption and backscattering, in a second step, to an apparent optical property, the diffuse attenuation coefficient. Model runs were made with a whole range of different inherent optical properties found in different water classes. Other conditions were set to standard values describing the average global scenario, though, limiting the accuracy of the approach. Similar settings as chosen for the simulations in this study were used, e.g., a mean wind speed of 5 m/s, vertically homogeneous waters, and clear sky conditions. The OC-CCI data set has been validated in comparison with in situ matchups. On regional scales, the performance quality of the Lee method can vary [68,75]. The DOAS-based method presented here and the method by Lee et al. [39] differ fundamentally in determining diffuse attenuation coefficient, but both approaches are rather analytical with a similar degree of assumptions that have to be made. One advantage of the DOAS-based method is that it effectively removes the atmospheric signal from the TOA radiances. No atmospheric correction is needed as additional step before the VRS retrieval.

As discussed above, the time series analysis showed that biases and seasonal variations are significant in the five Longhurst provinces. As the diffuse attenuation products are now, the agreement between the three hyperspectral sensors and overall data quality are not robust enough for merging the three data sets. Several studies have pointed out how intermission differences affect trend analysis on merged ocean color time series [76–78]. Based on the findings in this study, the intersensor consistency of the diffuse attenuation coefficient product should now be increased, e.g., through empirical corrections of the seasonal variations in GOME-2 as mentioned above, statistical bias correction [79], better calibration of the level 1 data [80], or a combination of several correction schemes.

5. Conclusion

We evaluated the possibility of exploiting the VRS signal in radiances measured by hyperspectral satellite sensors for assessing the diffuse light attenuation over long time scales building up on the work by Dinter et al. [5]. VRS signals were retrieved and diffuse light attenuation coefficients were derived for nearly one year of earth observations from three different hyperspectral satellite sensors, SCIAMACHY, GOME-2, and OMI.

VRS signals were detected in radiances from all three hyperspectral sensors. Starting with the latest published DOAS VRS retrieval settings, VRS fit factors were found to be sensitive to spectral structures in the measured radiances caused by instrumentation. For example, SCIAMACHY VRS fit factors were highly influenced by polarization features. Modifying the DOAS retrieval for SCIAMACHY and GOME-2 by changing how instrumental effects are accounted for in the retrieval, resulted in significantly improved diffuse attenuation coefficients. For SCIAMACHY, the zeta function was included to better correct for spectral features caused by polarization sensitivity of the sensor. For GOME-2, the amount of cross sections accounting for instrumental effects determined through a PCA on fit residuals over the Sahara were reduced from two to one eigenvector.

Through these DOAS fit setting modifications, VRS fit factor ranges from DOAS fits on satellite measurements were brought in agreement with VRS fit factor ranges predicted by fits on synthetic radiances calculated by radiative transfer modeling. In comparison to the established diffuse attenuation coefficient product from OC-CCI, robust correlations and slopes close to unity were found for all three sensors when only diffuse attenuation coefficients below 0.15 m⁻¹ were used. Highest correlation was found for GOME-2 data.

For further evaluating the quality of the derived diffuse attenuation coefficients from GOME-2, SCIAMACHY, and OMI, spatial and temporal deviations from the OC-CCI product were investigated by assessing global maps and time series of 7-day composites in five different Longhurst provinces. SCIAMACHY and OMI both showed latitudinal trends in derived diffuse attenuation coefficients. A north-south trend was found for SCIAMACHY indicating that spectral structures caused by polarization sensitivity could not be completely corrected despite including the zeta function in the DOAS fit. For OMI, biases with respect to the OC-CCI product become larger with increasing latitude indicating a disagreement between the SZA dependence of the VRS fit factor found in measurements in comparison to that predicted by radiative transfer modeling. This effect was not further investigated here. GOME-2 shows spatially a more homogeneous agreement with OC-CCI, but exhibits periodic changes in VRS fit factors over time which are probably caused by spectral structures in recorded solar measurements changing over the course of one year due to altering measurement geometries.

This first intercomparison of the performance of DOAS ocean retrievals on radiances from three different hyperspectral satellites revealed how challenging it is to accurately calibrate hyperspectral sensors to reliably and consistently retrieve small ocean signals. Merging the data sets is difficult since significant biases between the sensors were found. This study provides a good basis for further improving overall data quality and intersensor consistency in order to create a long-term diffuse attenuation coefficient time series in the blue spectral range. Based on the results of this study, it is unlikely that diffuse attenuation coefficients derived from measurements of the particular hyperspectral sensors presented here could compete in quality with the diffuse attenuation coefficient from multispectral sensors. However, since the method presented here is independent and offers assessing diffuse attenuation in a different spectral range, it is worth further developing this approach, especially with respect to upcoming satellite missions. Better calibration of hyperspectral satellite sensors with respect to ocean applications could improve this situation. Building up on this study, empirical corrections could be introduced to reduce biases and stabilize the time series of hyperspectrally derived diffuse attenuation coefficients. A more complex LUT could be build to reduce noise as shown by the presented sensitivity analysis. A thorough uncertainty estimation is desired. Comparison with in situ measurements becomes possible for the most recent or upcoming sensors that provide measurements at much higher spatial resolution. For example since October 2017, TROPOMI (TROPOspheric Monitoring Instrument) [81], similar in mission goal and design as OMI, provides global hyperspectral radiances at a spatial resolution of 7 km by 3.5 km. Validating 7 km by 3.5 km large pixels with in situ point measurements would still be difficult, but becomes more feasable for in situ measurements from towed systems yielding transects of in situ data. With increased spatial resolution, DOAS fit quality is also expected to improve due to increased homogeneity in satellite scenes.