Retrieving the Availability of Light in the Ocean Utilising Spectral Signatures of Vibrational Raman Scattering in Hyper-spectral Satellite Measurements

The availability of light in the ocean is an important parameter for the determination of phytoplankton pho-tosynthesis processes and primary production from satellite data. It is also a useful parameter for other applications, e.g. the determination of heat fluxes. In this study, a method was developed utilising the vibrational Raman scattering (VRS) effect of water molecules to determine the number of photons available in the ocean water, which is expressed by the depth integrated scalar irradiance E 0. Radiative transfer simulations with the SCIATRAN fully coupled ocean–atmosphere radiative transfer model (RTM) show clearly the relationship of E 0 with the strength of the VRS signal measured at the top of the atmosphere (TOA). Taking advantage of VRS structures in hyper-spectral satellite measurements, a retrieval technique to derive E 0 in the wavelength region from 390 to 444.5 nm was developed. This approach uses the weighting function differential optical absorption spectroscopy (WF-DOAS) technique, applied to TOA radiances, measured by the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY). Based on the approach of Vountas et al. (2007), where the DOAS method was used to fit modelled spectra of VRS, the method was improved by using the weighting function of VRS (VRS-WF) in the DOAS fit. This was combined with a look-up table (LUT) technique, where the E 0 value was obtained for each VRS satellite fit directly. The VRS-WF and the LUT were derived from calculations with the SCIATRAN RTM (Rozanov et al., 2014). RTM simulations for different chlorophyll a concentrations and illumination conditions clearly show that low fit factors of VRS retrieval results correspond to low amounts of light in the water column and vice versa. Exemplarily, 1 month of SCIAMACHY data were processed and a global map of the depth integrated scalar ir-radiance E 0 was retrieved. Spectral structures of VRS were clearly identified in the radiance measurements of SCIA-MACHY. The fitting approach led to consistent results and the WF-DOAS algorithm results of VRS correlated clearly with the chlorophyll concentration in case-I water. Comparisons of the diffuse attenuation coefficient, extracted by VRS fit results, with the established GlobColour K d (490) product show consistent results.


Introduction
Sunlight is the source of radiation which propagates through the ocean and drives the main biological and physical processes in the water.As pointed out by many studies (e.g.Morel, 1978Morel, , 1991)), the knowledge of the availability of light, which is determined by the scalar irradiance at depth E 0 (z), is required to quantify processes of photosynthesis, primary production and heat transfer.
The determination of the availability of solar light and accordingly radiant energy in the ocean from satellite remote sensing data is still a challenging task.The widely used approach to capture the amount of light in the water column from satellite data is to determine optically relevant parameters such as the diffuse attenuation coefficient K d and the light field below the water surface, and to calculate the scalar irradiance E 0 (z) at depth by using radiative transfer models T. Dinter et al.: Vibrational Raman scattering in oceanic water (e.g.Liu et al., 2006).Disregarding assumptions about the parameter of the ambient light, small errors in the determination of K d can result in high errors of E 0 (z).For instance, 30 % inaccuracy in K d can lead to a factor of 2 error in E 0 (z) (Lee et al., 2005b).
In this study, we used a new approach to retrieve the availability of light in the ocean by exploiting the spectral signatures of vibrational Raman scattering (VRS), detected in hyper-spectral satellite data.VRS is an inelastic scattering effect by the water molecules themselves.Energy is transported from the photon to the molecule during the scattering process and the water molecule is vibrationally excited.The emitted photon has another energy, i.e. a different wavelength.Rayleigh (and Mie) scattering is an elastic scattering process, where no energy is transferred between the scattering molecule or particle and the photon.In contrast, inelastic scattering occurs if the scattering molecule changes its state of excitation during the scattering process.Some of the photon's energy is then passed from the photon to the molecule (Stokes lines) or vice versa (anti-Stokes).VRS in liquid water is an inelastic scattering process, which provides a mean wave-number shift of ν = 3357 cm −1 (Walrafen, 1967).This leads to a wavelength shift of 40-120 nm in the wavelength region of 350-600 nm ( λ = − ν • λ 2 ).It involves two fundamental O-H stretch vibration modes of the water molecule that are further modified by hydrogen bonding and a librational fine structure.These interactions induce a broad band of emissions around the excitation wave-number shift, so that due to water, VRS radiation is re-emitted over a band of 30-50 nm depending on the wavelength (Haltrin and Kattawar, 1993) (see Fig. 1).
The VRS signal results from the scattering events of solar electromagnetic radiation with water molecules in the ocean and can contribute significantly to the marine upwelling radiance field in the UV and visible wavelength region (e.g.Westberry et al., 2013).The number of scattering events determines the strength of the VRS signal, detected at the top of the atmosphere.Changes in absorption, scattering, ocean surface or ambient light conditions directly affect the number of scattering events on water molecules, which produce a change in the VRS signal accordingly.Thus the VRS signal corresponds directly to the amount of light propagating through the ocean.Therefore, the strength of the VRS spectral signature is related to the light availability.Following the terminology of Morel (1991), we define the light availability as the scalar irradiance integrated over the water column, denoted in the following as the depth integrated scalar irradiance E 0 .
As a transpectral process, VRS contributes to the fillingin of solar Fraunhofer lines (Vasilkov et al., 2002;Vountas et al., 2003;Joiner et al., 2004) significantly.These spectral features of VRS were clearly identified in hyper-spectral satellite measurements, as shown for data of the satellite sensors GOME and SCIAMACHY, and were correlated with the chlorophyll a concentration (chl a) for at least case-I seawa- as measured by the SCIAMACHY satellite sensor.Incoming radiation in the VRS excitation wavelength region (390-444.5nm) leads to a filling-in in the VRS emission wavelength region (450-524 nm) as calculated with SCIATRAN (brown) according to Eq. ( 9).For three specific excitation wavelengths, their appropriate redistributions in the emission range are shown (in magenta, green, and cyan).
The jump of the solar spectrum between 390 and 400 nm leads to an equivalent jump in the VRS spectrum at 445 to 465 nm.
ter (Vasilkov et al., 2002;Vountas et al., 2007).In Bracher et al. (2009), the VRS signal was exploited as a proxy for the light path length in the water in order to derive chl a of different phytoplankton groups from slant column concentrations of their specific absorption spectra.
In this study we directly connect the strength of the VRS signal in SCIAMACHY measurements to the depth integrated scalar irradiance E 0 .Theoretical background information on definitions and the weighting function approach, used in the differential optical absorption spectroscopy (WF-DOAS) retrieval technique, are given in Sect. 2. To determine and to verify the relationship between the VRS signal and the in-water scalar irradiance, in Sect. 3 extensive radiative transfer calculations, carried out by the SCIATRAN coupled ocean-atmosphere radiative transfer model (RTM), are described.In order to account for this redistribution of photons due to the VRS effect, a spectroscopic model of VRS coupled with an adequate description of the interaction between light, seawater and atmosphere is required.Therefore the SCIATRAN coupled ocean-atmosphere RTM is used in this study (Rozanov et al., 2002(Rozanov et al., , 2014;;Blum et al., 2012).
On the basis of these simulations, a relationship between E 0 and the VRS signal was derived, parameterised and stored in a look-up table (LUT), appropriate for use in the retrieval algorithm.To derive VRS and thus inverted E 0 from SCIA-MACHY data, the application of the WF-DOAS method is shown in Sect. 4. Comparisons of the SCIAMACHY results to an established satellite product, the diffuse attenuation coefficient K d (490) (Mueller, 2000;Werdell and Bailey, 2005), are shown in Sect. 5.The advantages and limitations of the retrieval method and results as compared to other similar satellite data products are discussed in the last section.Here an outlook for further applications is also given.

The relationship between the intensity of the radiation field and the amount of light in the ocean
To determine a formulated relationship between the measured intensity of the radiation field at the top of the atmosphere (TOA) and the amount of light in the ocean, we start with defining the amount of the radiation energy within the ocean.For this, we use one of the main characteristics of the radiation field, i.e. the density of radiation energy which is given by e.g.Sobolev (1972): where c is the speed of light, z is the depth changing from z = 0 at the surface to z = H at the bottom of the ocean, the variable := {µ, ϕ} represents a pair of angle variables, where µ ∈ [−1, 1] is the cosine of the polar angle, and ϕ ∈ [0, 2π ] is the azimuthal angle, I λ (z, ) is the intensity of the radiation field, and the integration is performed over the unit sphere.Taking into account that the scalar irradiance is often defined as follows (e.g.Mobley, 1994), we obtain the following relationship between the density of radiation energy u λ (z) and the scalar irradiance E 0 (λ, z): To describe the amount of radiation energy in the ocean, we integrate u λ (z) over the entire ocean depth and over spectral range [λ 1 , λ 2 ]: The introduced value E 0 characterises the abundance of light in the vertical column [0, H ] and in the spectral range λ = λ 2 − λ 1 , and is also called the depth integrated scalar irradiance.It follows that E 0 with the unit W m −1 is directly related to the amount of radiation energy u in the vertical column by the speed of light.
In this study the intensity of the radiation field in the ocean, I λ (z, ), and consequently E 0 , is calculated under given conditions employing a coupled ocean-atmosphere radiative transfer model including VRS processes.To account for the fact that the intensity of the radiation field depends on numerous atmospheric and oceanic parameters, we assume that perturbations of E 0 caused by changes in these parameters can be described in a linear approximation as where E 0 is the perturbed depth integrated scalar irradiance, composed of E 0 and the perturbation Here, p i and q i are relevant atmospheric and oceanic parameters.p i and q i are their variations, and N a and N o are the number of considered atmospheric and oceanic parameters.
It follows that the computation of E 0 according to Eq. ( 5) requires the global information on numerous atmospheric and oceanic parameters.To obtain the required information, hyper-spectral satellite measurements of the backscattered earth shine radiation were used.Considering that the intensity at the TOA depends on the same parameters as E 0 , the logarithm of the TOA intensity can be represented as follows: ln where ln I λ is the logarithm of perturbed intensity at the wavelength λ.This formulation is restricted to the linear term of the Taylor series expansion of the perturbed logarithmic intensity.
Taking into account that the measurement of the intensity at the TOA is performed in a spectral range where the contribution of trans-spectral processes due to VRS are not negligible, the inelastically scattered radiation can be presented in the following form: where I − λ is the intensity calculated excluding the VRS process, and the additive component V λ , which will be called the VRS reference spectrum, is introduced as where describes the variation of the VRS reference spectrum caused by the variation of atmospheric and oceanic parameters.By introducing the so-called weighting functions W p i (λ) and W q i (λ) for elastically scattered radiation as Eq. ( 10) can be rewritten as follows: ln This equation constitutes the required relationship between the intensity at the TOA and the variation of the atmospheric and oceanic parameters.If all parameters can be obtained by solving Eq. ( 13), the variation of E 0 can be calculated according to Eq. ( 6).

The weighting function DOAS technique
In this section we introduce the weighting function DOAS (WF-DOAS) technique, which is used to obtain the solution of the formulated equation (Eq.13) derived in the previous section.The WF-DOAS algorithm is, analogous to the standard DOAS algorithm, a linear least squares algorithm, which yields typical species-dependent information only from differential absorption features.Weighting functions are the derivatives of the radiation field with respect to ocean, atmosphere, or surface parameters and are used in the retrieval of atmospheric trace gases.The WF-DOAS algorithm was originally applied to the retrieval of vertical amounts of strongly absorbing trace gases from GOME and SCIAMACHY data (Rozanov et al., 1998;Buchwitz et al., 2000;Coldewey-Egbers et al., 2005).Here we extend this method for the retrieval of oceanic parameters by optimising the wavelength region and adding appropriate weighting functions (see Sect. 4).All broadband contributions that affect the radiance are compensated for by using a low-order polynomial in the fit routine.As an result, the retrieval is relatively insensitive to aerosols, optically thin clouds, surface reflectivity, and other broadband absorption features.
The polynomial subtraction also reduces the sensitivity to any broadband residual radiometric calibration errors.Since we can not retrieve all parameters needed to calculate E 0 according to Eq. ( 6), we had to formulate an adequate approximated solution.Thus we assume that variations of the VRS reference spectrum, V λ , given by Eq. ( 11), are caused only by the variation of oceanic parameters.Assuming that the main driver for the oceanic inherent optical properties (IOPs) in case-I waters, typically encountered in the open ocean, is the phytoplankton (Morel and Prieur, 1977), which is characterised by its chl a, C, we can rewrite Eq. ( 11) as follows: where C is the variation of chl a, and comprises the contribution of all parameters other than chl a.
Introducing the effective parameter q V , which describes the variation of the VRS reference spectrum as the resulting expression, describing the logarithm of the perturbed intensity at the TOA given by Eq. ( 13), can be rewritten as follows: ln Here, the weighting function for VRS, is the derivative of the VRS reference spectrum with respect to chl a.The effective parameter q V comprises the variation of all relevant parameters.In this generalised form the introduced parameter q V depends on the wavelength λ.However, the wavelength dependence of the fit parameter (e.g.slant column in atmospheric evaluations) is a typical case for the standard DOAS retrieval algorithm (Burrows et al., 2011).Using this approximation, we are able to derive in accordance with Rozanov et al. (1998) and Buchwitz et al. (2000) the general formulation of the WF-DOAS approach with respect to the spectral impact of the target species VRS in TOA radiance measurements.In this case I λ is a modelled reference intensity for a fixed set of oceanic and atmospheric parameters and I λ in Eq. ( 17) can be replaced by a measured intensity, I meas λ , plus a low-order polynomial: The parameters q V , q i , and p i and the polynomial coefficients a k are obtained by adjusting the right-hand side of Eq. ( 19), which represents the model spectral intensity, to the measured spectral intensity (left-hand side).The unknown parameters ( q V , q i , and p i ) are determined by a (weighted) linear least squares minimisation procedure using a Levenberg-Marquardt fitting routine.We note that although the target parameter is q V , the atmospheric p i and ocean water q i parameters should be retrieved simultaneously because wrong estimations of these parameters, when resulting in differential spectral structures, may disturb the retrieval of q V or lead to high fit residuals.Such behaviour indicates a mismatch between the adjusted linear model and the measurement.
Assuming further that the variation of E 0 , given by Eq. ( 6), can be considered as a function of the effective parameter q V , the perturbed E 0 , given by Eq. ( 5), can be represented as follows: We note that the function F is defined in this way so that F (0) = E 0 and represents the reference point with no perturbation.This non-linear relationship is obtained by performing retrievals of numerous modelled spectra under known conditions and is established in the form of a look-up table (see details in Sect.3).The main differences between the standard DOAS (Platt and Stutz, 2008) for ocean applications used in previous publications (Vountas et al., 2003(Vountas et al., , 2007;;Bracher et al., 2009;Sadeghi et al., 2012) and the WF-DOAS technique are 1. the introduction of a modelled reference intensity spectrum for I λ , which corresponds to a reference point of a fixed oceanic, atmospheric and surface state, instead of using an extraterrestrial solar irradiance spectrum, and 2. the use of wavelength-dependent weighting functions (∂ ln I − λ /∂p i , ∂ ln I − λ /∂q i ) instead of pseudo and specific absorption reference spectra.
In line with standard DOAS, the logarithm of the intensity rather than the intensity itself is modelled (note that ∂ ln I /∂x = I −1 • ∂I /∂x).The WF-DOAS approach requires a radiative transfer model for the accurate simulation of the TOA radiance and its derivatives, considering oceanic, atmospheric and surface parameters.

Radiative transfer simulations
In this study the SCIATRAN version 3.2 was used for radiative transfer simulations (Rozanov et al., 2014).In contrast to earlier studies, where the downwelling irradiance at the sea surface and the light field in the ocean have mostly been calculated with separate RTMs (e.g. Lee et al., 2005a), this SCI-ATRAN version provides combined ocean-atmosphere calculations in one package.Thus, feedback (coupling) effects between ocean and atmosphere were included in the calculations.Besides coupling effects, SCIATRAN allows one to model inelastic scattering processes such as rotation Raman scattering (RRS) in the atmosphere and vibrational Raman scattering in water.The details of implementation and verification of RRS in the SCIATRAN software package are given by Vountas et al. (1998) and Rozanov and Vountas (2014).The verification of VRS was performed by comparing the VRS reference spectra with other model data and with VRS spectra obtained from hyper-spectral shipborne measurements of the solar radiation backscattered from the ocean (Kattawar and Xu, 1992;Peters, 2013).
The input data to perform radiative transfer calculations are an extraterrestrial solar spectrum, oceanic inherent optical properties (IOPs), and atmospheric and oceanatmosphere interface (ocean surface) optical properties.In particular, simulations were performed for the following scenario and will hereafter be referred to as the reference scenario.
-An extraterrestrial solar spectrum measured by the SCIAMACHY instrument was used (see Fig. 1, solid blue line) (Skupin et al., 2005).
-A cloud-and aerosol-free Rayleigh atmosphere including ozone absorption was assumed.
-The vertical profiles of temperature, pressure, and ozone concentration were set according to the mid-latitude standard atmosphere model (Sinnhuber et al., 2009).
-The absorption cross section of ozone was used according to Bogumil et al. (2003).
-The approach of Buiteveld et al. (1994) Kopelevich (1983).The concentrations of small and large particles were parameterised on the basis of Haltrin (1999), where a oneparameter model of seawater optical properties on the basis of chl a was presented.
-The total absorption coefficient was calculated based on the case-I water model of Morel and Maritorena (2001) as a combination of clear water, CDOM, and phytoplankton absorption, whereas the spectral phytoplankton absorption is calculated with the Prieur and Sathyendranath (1981) model.The clear water absorption coefficient is a merged spectrum of Sogandares and Fry (1997) for 340 to 380 nm and of Pope and Fry (1997) for 380 to 725 nm.
-The ocean-atmosphere interface transmission and reflection properties were approximated on the basis of Cox and Munk (1954) and followed by the Nakajima and Tanaka (1983) approach, where a typical wind speed of 4.1 m s −1 was assigned.
-A 500 m deep homogeneous mixed ocean layer with a black bottom albedo was assumed.
To study the influence of underwater light field conditions on satellite measurements, combined simulations of TOA radiances and radiation fluxes into the ocean were performed.This framework of in-water flux simulations in combination with the modelling of radiation for satellite geometry at the TOA delivered insight to infer several findings.Since the VRS effect is a transpectral process, different (excitation and emission) wavelength regions have to be taken into account.The amount of radiation in the excitation wavelength region of VRS, λ ex (in our study from 390 to 444.5 nm), leads to a filling-in in the emission wavelength region from 450 to 524 nm.Fig. 1 illustrates this relationship.Photons are transported from the excitation into the emission wavelength region and are redistributed as a sum of four Gaussian functions (Haltrin and Kattawar, 1993).The relative amount of light due to VRS emission is represented in Fig. 1 by the VRS reference spectrum according to Eq. ( 9).
Therefore subsurface light field simulations have to be combined with TOA radiance simulations.In particular, a direct relationship between E 0 in the excitation region and the VRS signal in the emission wavelength region detected at the TOA has to be established.This will be described in the next subsections.

Modelling of the depth integrated scalar irradiance and the averaged diffuse attenuation coefficient
The subsurface radiation fluxes were calculated based on the reference scenario, described above, for 23 different chl a concentrations (0-30 mg m −3 ) and for 73 solar zenith angles (SZA) in the range of 0 to 89 • .The downwelling, E d (z), upwelling, E u (z), and scalar irradiance, E 0 (z), were modelled for 19 different water depths from 0.001 up to 500 m accordingly.
The simulated scalar irradiance E 0 (z, λ ex ) in Fig. 2a, as in all other figures in this section, is plotted for a fixed SZA of 40 • and shows an exponential behaviour with depth for different chl a as expected.The scalar irradiance, integrated over the depth and over a spectral range λ ex (see Eq. 4), where λ 1 and λ 2 were set to 390.0 and 444.5 nm, respectively, results in E 0 ( λ ex ), which is illustrated in Fig. 2b.E 0 (dependence on λ ex is omitted in the following) shows an inverse S-shaped relation to the log-scaled chl a.Low chl a leads to high E 0 and vice versa.By using the SCIATRAN coupled ocean-atmosphere RTM, the following relationship between E 0 and chl a was obtained: To complete the study and for the purpose of comparisons to multispectral data products (see Sect. 5), the attenuation depth z 90 and the diffuse attenuation coefficient K d were obtained from the subsurface downwelling flux simulations The attenuation depth z 90 is defined as the depth at which the downwelling flux is 1/e times smaller than the subsurface downwelling flux E d (0 − ) (Gordon and McCluney, 1975).The depth z 90 depends on the wavelength λ and is widely used as the penetration depth of light into ocean water.The spectral behaviour of z 90 for all wavelengths between 340 and 600 nm is shown in Fig. 3a.This figure shows that z 90 has a strong dependence on chl a and reaches almost 100 m at a wavelength of 430 nm in the absence of phytoplankton absorption.Figure 3b shows the dependence of z 90 ( λ ex ) on chl a, where the downwelling fluxes were integrated over the VRS excitation wavelength region λ ex of 390-444.5 nm.
The averaged remote sensing diffuse attenuation coefficient at wavelength λ can be determined following Lee et al. (2005b) as Then a simple relationship between K d and z 90 can be obtained: The spectra of K d (λ) were calculated using Eq. ( 23) and are shown in Fig. 4a for different chl a.The dependence of K d ( λ ex ) in the VRS excitation wavelength region λ ex on chl a is drawn in Fig. 4b and can be expressed in the following form: This function shows a non-linear exponential relationship between the diffuse attenuation K d ( λ ex ) and the chlorophyll a concentration C. SCIATRAN scalar irradiance E 0 (∆λ ex ) at depth

Weighting functions for VRS and oceanic parameters
In order to account for the variation of all relevant oceanic parameters in the WF-DOAS retrieval algorithm, one needs to calculate the derivatives described in Eqs. ( 12) and ( 18).Assuming case-I water conditions, the one-parameter model of seawater optical properties on the basis of chl a described in Sect. 3 was applied, which also reduced the number of derivatives to calculate the weighting functions to 1.
To extract the changes in the TOA radiance due to the modification of the VRS signal, the calculation of the derivative of the VRS reference spectrum with respect to chl a, according to Eq. ( 18), was required.Taking into account that this derivative can not be derived analytically, the following finite difference approximation was used: Here, intensities I λ and I − λ are calculated including and excluding VRS processes, respectively, with chl a equal to C and C + C.
Figure 5 shows the derivative of the VRS weighting function, calculated according to Eq. ( 25), with a change in chl a from 0.1 to 0.11 mg m −3 ( C = 0.01 mg m −3 ) for a SZA of 40 • .This weighting function was used in the DOAS retrieval algorithm, described in the following sections.Additionally, another VRS-WF is shown in Figure 5, which was calculated for a chl a of 0.5 mg m −3 with the same change a) b)  of C = 0.01 mg m −3 .By scaling, both spectra show very similar differential spectral features.To account for the variation of all other (except VRS) relevant oceanic parameters in the WF-DOAS retrieval, the calculation of the derivatives of the intensity I − λ , according to Eq. ( 12), with respect to these parameters was required.To simplify the independent variations of these parameters, it was convenient to combine these in one weighting function.From a mathematical point of view, the sum of variations of all oceanic parameters given by the third term on the right-hand side of Eq. ( 19) can be rewritten as follows: where the weighting function W Oc (λ) is given by This weighting function comprises contributions of all other oceanic parameters (except VRS) changing with chl a, as described in the beginning of Sect. 3. The weighting function W Oc (λ) was calculated using a finite difference approximation as follows: • ln where the intensities I − λ (C) and I − λ (C+ C) were calculated as described before with chl a equal to C = 0.1 mg m −3 and C + C = 0.11 mg m −3 , respectively.The spectral shape of the W Oc (λ) is plotted in Figure 5 and results, as expected, in a combination of water and phytoplankton absorption.

Non-linear relationship between VRS fit factor and chl a concentration
In Sect.3.1 the non-linear relationship between the depth integrated scalar irradiance E 0 and chl a as given by Eq. ( 21) was established.Taking into account that the WF-DOAS retrieval algorithm provides the VRS fit factor q V from TOA radiance simulations, the relationship between q V and chl a was derived.Combining then both functions, we obtain the required resulting relationship between E 0 and q V (see Eq. 20).For this purpose, the TOA radiance simulations were performed for the same scenarios as the subsurface radiation calculations in Sect.3.1.Assuming that the TOA radiance is affected by variations of chl a Eq. ( 19  ln where C i for i = 1, 2, . .., 23 are perturbations of chl a from a priori C 0 = 0.1 mg m −3 .In our retrieval, no weighting function for atmospheric trace gases but the ozone cross section σ Oz and the slant column density SC Oz were used.This is a valid assumption for weak absorbing trace gases (Rozanov and Rozanov, 2010).Solving this equation for all C i by using the WF-DOAS technique described in Sect.2.2, the following non-linear relationship was derived: Figure 6a and b shows the fit results of the VRS weighting function against the modelled TOA radiance in the wave-length region of 450-524 nm.This fit window was chosen because 1. the appropriate excitation wavelength region of VRS (390-444.5nm) (Fig. 3b) shows for lowest chl a the highest penetration depth of radiation into the ocean and therefore the largest variation of E 0 ; 2. at 455 nm (see Fig. 5), there is a large step in the VRS reference spectrum and weighting function, which leads to distinct differential structures.These structures are very different compared to the atmospheric Ring effect, which is also an inelastic scattering effect but with a much narrower distance between the excitation and emission wavelength regions (Grainger and Ring, 1962); 3. it allows validation with comparable data products of different satellite sensors (e.g.K d (490)) (see Sect. 5); and 4. it is situated within the range of SCIAMACHY cluster 15 (424-525 nm) of channel 3, which was used for our retrievals.
The scaling of the VRS weighting function is expressed by q v and is plotted in Fig. 7 against a log-scaled chl a.It shows the non-linear relationship from Eq. ( 30) and behaves in an asymptotic S-shaped characteristic.The fit factor of 0 is retrieved under the reference conditions of chl a with C 0 = 0.1 mg m −3 .Therefore, the function f q , given by Eq. ( 30), satisfies the following condition:   2012), respectively) and profiles (profile-1 for a stratified and profile-2 for a mixed water profile according to Uitz et al. (2006)).(b) VRS fit factor q V vs. chl a for the same scenarios as in (a).

Look-up table for the relationship between E 0 and VRS fit factor
By combining the model results of the subsurface scalar irradiance calculations given by Eq. ( 21) with the model results of TOA radiance simulations given by Eq. ( 30), the non-linear relationship between the depth integrated scalar irradiance and the VRS fit factor is obtained as follows: Comparing the derived expression with Eq. ( 20), the function F ( q V ) introduced above can be obtained as a combination of the functions f c and f −1 q .The relationship between E 0 and q V given by Eq. ( 32) was established in the form of a look-up table (LUT) and is shown by the solid magenta line in Fig. 9.
Previous research shows that the relationship between chl a and the absorption of phytoplankton is complex (e.g.Morel and Bricaud, 1981;Sathyendranath et al., 1987;Babin et al., 1993;Bricaud et al., 1995).Packaging effects and different pigment compositions lead to different specific (chl a normalised) absorption and may differ up to 1 order of magnitude.In addition the assumption of a vertically homogeneous mixed water body is not a realistic scenario.To investigate the impact of different phytoplankton compositions (i.e.different pigment compositions and specific absorptions) and different vertical chl a distributions on the relationship between E 0 and q V in the reference LUT (see Fig. 9) and Eq. ( 32), especially with respect to the different wavelength regions of excitation and emission of VRS, additional scenarios with the absorption characteristic of three different phytoplankton types and two realistic phytoplankton profiles were investigated.The calculations of the functions f c (C) Resulting relationship between the VRS fit factor q v and E 0 , derived according to Eq. ( 32) for different phytoplankton types and profiles as in Fig. 8.The solid magenta line is a fitted third-order polynomial to the reference scenario (magenta points) with a SZA of 40 • and is used as a LUT for the satellite data retrieval.
and f q (C) were repeated in the same way as for the reference scenario.
In particular, the specific absorption spectra of diatoms and cyanobacteria were taken from Bracher et al. (2009), and of coccolithophores (Emiliania huxleyi) from Sadeghi et al. (2012).These spectra were not normalised to their absorption at 440 nm as the phytoplankton spectrum of Prieur and Sathyendranath (1981) to produce simulations with strong differences in absolute values and spectral shape.Thus, for example, the maximum of the diatom spectrum at 440 nm has a value of 0.015 as compared to 1.0 in the reference scenario.Additionally, typical phytoplankton profiles for mixed and stratified waters according to Uitz et al. (2006) were included in this sensitivity study.The mixed vertical profile has a constant concentration up to 30 m depth and decreases exponentially underneath, whereas the stratified profile is characterised by a smooth ascent to a chl a maximum at 60 m and also by an exponential decrease underneath (see Uitz et al. (2006), Fig. 5 profile M 2 and Fig. 4 profile S 4, respectively).
The obtained functions f c (C) and f q (C) for the five additional scenarios, introduced above, are presented in Fig. 8a  and b.To facilitate the comparison, the functions for the reference scenario are given in the same figures.It can be seen that the relationship between E 0 and chl a as well as between q V and chl a strongly depends on the specific absorption coefficient and stratification of the ocean layer.Thus, it follows from Fig. 8a that using the specific diatom absorption spectrum at chl a of 0.1 mg m −3 leads to 2 times larger E 0 as compared to the reference scenario.
The strong dependence of the VRS fit factor on the different absorption coefficients and the chl a profiles is shown in Fig. 8b.The combination of the functions f c (C) and f −1 q ( q V ) as given by Eq. ( 32) results in a significant weaker dependence of the relationship E 0 ↔ q V on the specific absorption coefficient and the vertical chl a distribution (see Fig. 9).The concentration of particles and the particle size distribution were changed according to the one parameter model (chl a) of Haltrin (1999), which provides appropriate parameters for the Kopelevich (1983) scattering model.By considering different specific phytoplankton absorption spectra, the relationship between scattering and absorption also changed extremely.This shows that the fit of the VRS weighting function reacts similarly to E 0 with respect to changes in the seawater optical conditions and is able to compensate for unknown variations of IOPs.Nevertheless, small deviations can be explained: changes in specific absorptions and chl a profiles lead to spectral deformations in the VRS weighting function and its differential form, and can not be compensated for completely within the DOAS algorithm.A third-order polynomial is fitted to the reference scenario (magenta line and points in Fig. 9) and is used as LUT for the satellite data retrieval.It follows from the LUT in Fig. 9 that for all considered additional scenarios the deviations of a specific absorption coefficient and chl a profile from the reference scenario may lead to small errors up to ∼ 10 %.

Solar zenith angle-dependent look-up table
In order to extend the LUT for different illumination conditions, additional in-water and TOA radiation simulations for seven SZAs between 20 and 80 • were performed.The same procedure, described in the last section, to build up the LUT of E 0 and the VRS fit was repeated and fitted third-order polynomials were applied to each function.By this technique a linear three-dimensional interpolation scheme was established.The resulting three-dimensional LUT, shown in Fig. 10, was used to calculate an SZA-dependent relationship between E 0 and VRS fit for the satellite data retrieval described in the next section.

SCIAMACHY satellite data retrieval
The SCIAMACHY instrument (Scanning Imaging Absorption Spectrometer for Atmospheric CHartographY) (Bovensmann et al., 1999)  in the UV-Vis region.Fully calibrated Level-1C data from channel 3, cluster 15 (424-525 nm) were used in this study.This spectrometer's cluster provided a spatial resolution of 30 km × 60 km per pixel with an integration time of 0.25 s.The swath width for the nadir scan measurement was 960 km.With a continuous changing of nadir-limb viewing direction and a resultant interrupted nadir along-track orbit, the sensor reached a global coverage time within roughly 6 days.The main objective of the SCIAMACHY mission was the determination of the abundances of atmospheric trace gases.Nevertheless, large parts of the solar radiation penetrate the ocean surface and are influenced by absorption and scattering effects of seawater and its constituents, which are included in the backscattered radiation detected by the instrument.The global fits were performed up to a SZA limit of 80 • due to signal-to-noise issues.Furthermore, a TOA reflectance threshold (R thresh = I TOA • π/I sol / cos(SZA)) of 0.29 was applied to remove cloud, aerosol and glint contaminated pixels.Here, I TOA and I Sol are the TOA radiance and the extraterrestrial solar irradiance measurement spectrally averaged over the whole wavelength range of SCIAMACHY cluster 15 (424-525 nm).

Retrieval of the depth integrated scalar irradiance E 0
In order to derive E 0 , we used the WF-DOAS technique described in Sect.2.2.We denote in the following the retrieved depth integrated scalar irradiance E 0 from Eq. ( 20) as E 0 , i.e. omitting the superscript.The generalised form of the WF-DOAS equation given by Eq. ( 19) was solved in a least square minimisation problem applied to the wavelength range from 450 to 524 nm as where τ dif (λ) = ln(I meas λ /I λ ) (see also Eq. 19), W V and W Oc are the VRS and oceanic weighting functions introduced in Sect.3.2, q V and q Oc are the appropriate fit factors, R is a rotational Raman scattering reference spectrum, which compensates for the atmospheric Ring effect (Grainger and Ring, 1962), and σ i (λ) and SC i are the absorption cross sections and slant column densities of considered atmospheric trace gases, which were water vapour (H 2 O), ozone (O 3 ), glyoxal (CHOCHO), oxozone (O 4 ) and nitrogen dioxide (NO 2 ).We note that, in contrast to Eq. ( 19), we used in Eq. ( 33) the cross sections instead of weighting functions.This approximation is often used in the case of weak gaseous absorption (Rozanov and Rozanov, 2010).As discussed in Sect.2.2, the fitted third-order polynomial compensates for all additional broadband absorption and scattering features.
The parameters q V , q Oc , q R , and SC i and the polynomial coefficients a k are the unknown parameters and were obtained by computing the least square minimisation solution of Eq. ( 33).All DOAS algorithms require an accurate spectral calibration.In order to correct even for small spectral misalignments between the reference and the measured spectra, a non-linear shift-and-squeeze spectral correction algorithm is usually applied for standard DOAS retrievals, and is also applied for the WF-DOAS algorithm used in this study.Global VRS fit factors q V were obtained by applying the WF-DOAS technique described above to SCIAMACHY nadir measurements of October 2008.An example of a VRS WF-DOAS fit from a SCIAMACHY measurement over the South Pacific on 23 October 2008 is shown in Fig. 11.At this site (126.71• W, 34.85 • S) very low chl a and high E 0 values are expected (see Morel et al., 2007).In Fig. 11 the scaled VRS weighting function compares well to the overall fit residual, which equals the measurement subtracted by all components other than VRS.The results clearly show that the spectral structures of the VRS signal were extracted from the TOA radiance measurement.
Applying the reference LUT from Fig. 9 without SZA correction, a global distribution of E unc 0 is derived, which is shown in Fig. 12. Regions with high E unc 0 values, coloured in blue/violet, correspond to well-known large areas of oligotrophic waters.These regions also correspond to a high scaling of the VRS weighting function.The map shows similar variability and patterns to global satellite chl a from imaging spectrometers like SeaWiFS, MODIS and MERIS.This is coherent insofar as the main optical driver in open ocean case-I waters is phytoplankton, which here dominates the absorption of light.
Applying the SZA correction introduced in Sect.3.4, different dependencies are distinguishable (see Fig. 13): generally, there is low light availability at the high latitudes and high light availability at the low latitudes.However, in tropical regions (latitude < ±30 • ), patterns of phytoplankton abundance are observable, whereas for higher latitudes the SZA dependency dominates E 0 .

Comparison
In this section our WF-DOAS approach to derive E 0 is compared to another satellite product.Unfortunately, no appropriate and independent data sets of E 0 , to test the accuracy of our retrieval results, were available.However, global satellite data sets of the diffuse attenuation coefficient at 490 nm, K d (490), exist.For our comparison, the operational Glob-Colour K d (490) product was used (see Barrot, 2010).The GlobColour K d (490), an indicator of the turbidity of the water column, is computed directly from the CHL-1 chlorophyll product, merged from data of the three different satellite sensors SeaWiFS, MODIS-Aqua, and MERIS (Barrot, 2010).
As shown in Sect.3.1, we used the underwater downwelling flux simulations to calculate K d ( λ ex ) (see Fig. 4).In the same way as in Sect.3.4, the relationship between the VRS fit factor and K d ( λ ex ) was derived and expressed in a LUT.In particular, by combining Eqs. ( 24) and ( 30), the non-linear relationship between the VRS fit factor and K d ( λ ex ) was obtained as follows: and is shown in Fig. 14 for a SZA of 40 • .Additionally, identical to Sect.3.4.1 and Fig. 10, a threedimensional LUT was built up to consider the SZA dependency of the K d ( λ ex ) ↔ q V relationship.
For the comparison in Fig. 15, the K d ( λ ex ) values were calculated from the VRS retrieval for each SCIAMACHY pixel by applying the three-dimensional LUT for October 2008 and then averaging for each day on a 1 • × 1 • grid.The same spatial and temporal resolution was chosen for the GlobColour data and then matched to the SCIAMACHY daily grids.This comparison is shown in the density scatter plot of Fig. 15, which reveals an ambivalent behaviour: a good agreement between K d (490) and K d ( λ ex ) up to values of ∼ 0.06 m −1 is observable, whereas for higher values of K d , the scatter plot shows a butterfly shaped distribution.The comparison of K d in the range 0.02-0.06m −1 shows a small offset for small K d values, where the SCIAMACHY results are generally a bit lower than GlobColour.This is consistent with the wavelength dependence of K d (λ).Indeed, Fig. 4a shows that K d (λ) in the VRS excitation wavelength region is lower than K d (490) for low chl a and higher than K d (490) for high chl a.
Furthermore, the comparison shows that high GlobColour K d (490) values could not be reproduced by the SCIA-MACHY VRS retrieval results, and vice versa.We suppose that such a butterfly distribution can be explained firstly by the impact of some cloud-contaminated SCIAMACHY pixels.Indeed, if one SCIAMACHY pixel with a spatial resolution of 30 km×60 km is contaminated with clouds, this leads to an erroneous decreasing of the VRS signal and an overestimation of K d ( λ ex ) values in the daily grid cell.This can explain the horizontal wing of a butterfly distribution.
Secondly, the vertical wing of the butterfly distribution may be explained by variation in IOPs and in vertical chl a profiles.We like to point out that the calculation of K d (490) of the GlobColour product is performed according to the following expression (Barrot, 2010): where K w (490) = 0.0166 m −1 , χ (490) = 0.08349, e(490) = 0.63303, and chl is the chlorophyll a concentration expressed in mg m −3 .Here, chl a is the main parameter which is required to estimate the diffuse attenuation coefficient.Taking into account that the DOAS retrieval algorithm based on the VRS fit factor is not designed to retrieve chl a as an end product, its concentration can diverge significantly if e.g. the specific absorption of the phytoplankton differs highly from that which is considered in the reference scenario of the RTM.In such cases an underestimation of K d ( λ ex ) can explain the vertical wing of butterfly distribution.
The comparison of the two data sets, presented in Fig. 15, has a correlation coefficient of ∼ 0.42.However, removing the outliers by simply cutting the wings of the butterfly, nearly 10 000 (∼ 13 % of all data used in the comparison) pixels were excluded (see grey boxes in Fig. 15) and a correlation coefficient of 0.69 was derived.
We want to note that in this study we are not proposing a new method to retrieve a diffuse attenuation coefficient from space, but just use the GlobColour K d (490) product to verify our SCIAMACHY E 0 results.Nevertheless, the comparison shows, in regions of K d values lower than 0.06, where the variability of the IOPs is expected to be low, good agreement between the two data sets.This confirms our derived relationship between VRS signals and the availability of light in the ocean.

Discussion and outlook
The rate of photosynthesis of phytoplankton depends on several parameters: phytoplankton biomass, light availability, and the so-called P-I (photosynthesis-irradiance) parameters (e.g.Jassby and Platt, 1976;Morel, 1991).Thus, in order to develop realistic models of ocean primary productivity, and carbon and heat fluxes, global knowledge of the availability of underwater light is required.To date, most approaches to determine the amount of light in the ocean, incorporating satellite ocean-colour data, have been based on the determination of the diffuse attenuation coefficient K d (490) and the estimation of the photosynthetically active radiation at the ocean surface PAR(0 + ) using atmospheric corrected satellite reflectance measurements (e.g.Mueller, 2000;Werdell and Bailey, 2005;Lee et al., 2005b;Frouin and Pinker, 1995).
With the method developed in our study, we overcome the determination of K d and PAR(0 + ) by allowing the retrieval of E 0 for the wavelength region of 390-444.5 nm directly as a columnar value, by retrieving the VRS signal in hyper-spectral satellite data.As shown in our model study (see Sect. 3), the strength of the VRS signal correlates directly with the amount of light in the water.Changes in IOPs lead to adequate changes in the VRS signal, which is shown in a sensitivity study, where different phytoplankton absorption spectra and profiles were included (see Sect. 3.4).The retrieval is still very robust despite such variations in IOPs: its error is below 10 %.Additional testing of the retrieval with a different aerosol loading (here a maritime background aerosol with an optical thickness of 0.05 was considered) revealed low impact on the VRS fit.For clear water conditions this led to a deviation below 8 % for the retrieval of E 0 and to much lower deviation (< 1 %) for higher concentrations of water constituents.To parameterise the relationship of the VRS signal at the TOA with E 0 , extensive radiative transfer simulations were calculated.This relationship was then expressed in a three-dimensional LUT, including a SZA correction, and an adequate VRS weighting function was calculated.
The determination of the depth integrated scalar irradiance E 0 was derived from satellite SCIAMACHY measurements without former retrieval of certain IOPs, in particular the absorption and (back-)scattering coefficients, and the phase function of the water constituents.With the LUT and the VRS weighting function, E 0 at 390-444.5 nm was determined from hyper-spectral SCIAMACHY satellite data using an improved WF-DOAS method in a fitting wavelength region of 450-524 nm.We restricted the retrieval to this wavelength window, because the broader the fit window, the fewer spectral long-wave effects (compared to VRS structures) can be compensated for by an additional fitted polynomial, which can lead to deviations in the fit of the target species (W Oc and W VRS ).Also, spectral changes in absorption and scattering by water constituents in the ocean which are not considered by the W Oc weighting function and include long-wave structures are compensated for by the fitted third-order polynomial in the WF-DOAS retrieval.Uncompensated structures are left in the residual and lead to high fit errors and high chi-square (χ 2 ) values.However, the χ 2 values in the VRS fit retrieval were in the range of 10 −3 and show reliable fit quality.
In this study we focused on the wavelength region around 400-500 nm for the purpose of feasibility and comparison to other satellite products.The comparison for 1 month of data with an independent data set (GlobColour K d (490) product) and retrieved values of the diffuse attenuation coefficient from the VRS signal show consistent results.Nevertheless, this method is easily applicable to other wavelength regions to retrieve a wavelength-dependent E 0 .Former studies showed good VRS fit results from SCIAMACHY measurements in the UV wavelength region (Vountas et al., 2007;Bracher et al., 2009).As discussed above, the challenge is to find appropriate fit windows where 1. the spectral signature of VRS does not significantly correlate with other effects, such as filling-in of Fraunhofer lines by atmospheric rotational Raman scattering, 2. the fit window is not too wide, otherwise long wave effects (broadband structures) are not compensated for by the fitted polynomial, and 3. the instrument provides an appropriate spectral resolution (at least better than 0.5 nm) and a high signal-tonoise ratio (e.g.> 2000 for SCIAMACHY) to resolve the filling-in of Fraunhofer lines by weak VRS signals.Morel (1991) defines PAR as the photosynthetically available radiation where PAR(λ, z, t) = E 0 (λ, z, t) (Morel, 1991, p. 276, and Appendix 1).Generally, PAR is determined as an integrated value over the wavelength λ, depth z, and time t.Usually the integration ranges between the limits of the photosynthetic spectral domain (400-700 nm), the euphotic depth, and the time length of 1 day.RTM simulations show that the retrieval of E 0 at 390-444.5 nm considers from 10 % (high chl a) to 40 % (low chl a) of E 0 in the spectral domain up to 700 nm.Extending the retrieval of VRS in hyper-spectral satellite data to multiple wavelength regions has the potential to retrieve directly the availability of the photosynthetically active radiation PAR in the water column.But, the retrieval of E 0 for the UV-A and UV-B region is also possible.For SCIAMACHY data this task is limited by the changing of the integration time within channels 3 and 4, which changes the pixel size from 30 km × 60 km to 240 km × 60 km.Improvements in VRS fitting may be achieved by including the chl a product from multispectral satellite imaging: the appropriate VRS weighting functions determined for the correct chl a (from the chl a satellite product) should be calculated and used as a starting point for the WF-DOAS minimisation algorithm in order to fulfil the assumption of linearity around the reference point.In addition, changes from a homogeneous to a heterogeneous profile lead to consistent changes in the VRS signal and E 0 , which have been shown in the sensitivity study in Sect.3.4.But, for the purpose of refinements, more diverse scenarios have to be simulated by radiative transfer calculations, and different sets of weighting functions may be taken into account within the retrieval process.
Limitations to our VRS and E 0 results are the coarse spatial and temporal resolution of SCIAMACHY, which will be overcome by adapting this retrieval to similar hyper-spectral satellite missions, such as GOME-2 on the METOP series, OMI on AURA and the upcoming TROPOMI sensors on the Sentinel-5-Precursor, Sentinel-4 and Sentinel-5 missions (launches are scheduled for 2015, 2019 and 2020, respectively).All of these sensors have global coverage within 1 day and the latter have a pixel size of 7 km by 7 km.With this analytical and generic method, the establishment of unique long-term information on light availability in the ocean will be enabled.

Figure 1 .
Figure1.Incoming extraterrestrial solar irradiance spectrum (blue), as measured by the SCIAMACHY satellite sensor.Incoming radiation in the VRS excitation wavelength region (390-444.5nm) leads to a filling-in in the VRS emission wavelength region (450-524 nm) as calculated with SCIATRAN (brown) according to Eq. (9).For three specific excitation wavelengths, their appropriate redistributions in the emission range are shown (in magenta, green, and cyan).The jump of the solar spectrum between 390 and 400 nm leads to an equivalent jump in the VRS spectrum at 445 to 465 nm.

Figure 2 .
Figure 2. (a) Scalar irradiance E 0 (z, λ ex ) of the subsurface light field, integrated over the VRS excitation wavelength region λ ex from 390 to 444.5 nm, for different chl a and for a fixed SZA of 40 • .(b) Depth integrated scalar irradiance E 0 ( λ ex ) as a function of chl a. a) b)

Figure 6 .Figure 7 .
Figure 6.(a) Absolute spectral fit result (spectral optical density) of the VRS weighting function in the wavelength region of 450-524 nm for different chl a.(b) Differential spectral fit result (differential optical density) which is derived by subtracting a fitted third-order polynomial.

Figure 8 .
Figure 8.(a) E 0 vs. chl a for different phytoplankton types (using phytoplankton absorption spectra specific for diatoms (dia), cyanobacteria (cya) and coccolithophores (emi) taken from Bracher et al. (2009) and Sadeghi et al. (2012), respectively) and profiles (profile-1 for a stratified and profile-2 for a mixed water profile according toUitz et al. (2006)).(b) VRS fit factor q V vs. chl a for the same scenarios as in (a).

Figure 10 .
Figure10.Same relationship as for the reference scenario in Fig.9for different SZA in a three-dimensional LUT.

Figure 11 .
Figure 11.Example of the optical density (OD) of a VRS weighting function fit at an oligotrophic site over the South Pacific on 23 October 2008.The solid line is the scaled VRS weighting function from Fig. 5 and the dashed line is the measurement, where all components are subtracted besides VRS.

Figure 12 .
Figure 12.Global map of the depth integrated scalar irradiance E unc 0 at 390-444.5 nm from DOAS VRS weighting function fits in unit of photons s −1 m −1 without SZA correction.Red shows low and blue-magenta high E 0 values.Note the different colour table toFig.13.

Figure 13 .
Figure 13.Global map of the depth integrated scalar irradiance E 0 at 390-444.5 nm from DOAS VRS weighting function fits in unit of photons s −1 m −1 including the SZA correction from Fig. 10.Red shows low and blue-magenta high E 0 values.

Figure 14 .
Figure 14.Resulting relationship of the VRS fit factor with the diffuse attenuation coefficient for a SZA of 40 • , which is used as a LUT for the retrieval of K d ( λ ex ) from SCIAMACHY satellite data.

Figure 15 .
Figure 15.Density map of the scatter plot of K d ( λ ex ) retrieved from SCIAMACHY data and the K d (490) GlobColour data product, matched by daily grids with a spatial resolution of one degree.