Contribution of Raman scattering to polarized radiation field in Contribution of Raman scattering to polarized radiation field in ocean waters ocean waters

: We have implemented Raman scattering in a vector radiative transfer model for coupled atmosphere and ocean systems. A sensitivity study shows that the Raman scattering contribution is greatest in clear waters and at longer wavelengths. The Raman scattering contribution may surpass the elastic scattering contribution by several orders of magnitude at depth. The degree of linear polarization in water is smaller when Raman scattering is included. The orientation of the polarization ellipse shows similar patterns for both elastic and inelastic scattering contributions. As polarimeters and multipolarization-state lidars are planned for future Earth observing missions, our model can serve as a valuable tool for the simulation and interpretation of these planned observations.


Introduction
Scattering and absorption of light in ocean waters include both elastic and inelastic processes. The inelastic processes mainly include Raman scattering by pure liquid water and fluorescence by phytoplankton [1]. Raman scattering in bulk water transfers light energy from higher frequency to lower frequency regions and significant signals are observed with fixed frequency shift spectra centered around 3400 cm −1 , which are mainly attributed to the OH bond stretch mode [2]. It has been long recognized that the contribution of Raman scattering to the underwater radiation field is significant in the visible region for clear waters, both experimentally [3][4][5][6][7] and theoretically [8][9][10][11][12][13]. Studies have shown that Raman scattering has major impacts on the ocean color remote sensing algorithms [14,15].
The theoretical modeling of Raman scattering in ocean waters is mostly based on the Monte Carlo (MC) method [8][9][10][11][12][13]. Stavn and Weidemann studied the impacts of Raman scattering on the underwater irradiance [8]. Raman scattering is also known to be the cause of filling the solar Fraunhofer lines in ocean waters [9,10]. The angular resolved radiation field was reported in [9,[11][12][13]. Another known code for solving Raman scattering in ocean waters is Hydrolight, a commercial software based on the invariant embedding method that simulates the unpolarized radiation fields [16]. Semianalytical models have also been developed that approximate under water irradiance fields including Raman effects [17][18][19].
The modeling of Raman scattering in ocean waters mostly neglects the polarization of the radiation field due to Raman scattering. To date the only paper which studied the polarized radiation field of Raman scattering was by Kattawar and Xu based on the MC method [11]. The polarized radiation field reported in [11] was averaged over the viewing azimuth angle interval of 30 • in order to reduce the MC statistical noise. Though the MC method is a versatile tool for simulating radiative transfer in the ocean, it is also known that the convergence of the MC method is slow for cases with strongly forward-intensified phase functions [20], which prevail in ocean waters [1]. To improve the efficiency and accuracy of the vector radiative transfer simulation dealing with Raman scattering, a deterministic solution is desired for understanding the underwater polarized light field. In addition this kind of radiative transfer tools will be needed for the development of algorithms for retrieval of ocean inherent optical properties (IOP) using a polarimeter, which is a goal of next generation of ocean color missions such as NASA's Pre-Aerosol, Clouds, and ocean Ecosystem (PACE) mission [21].
In this paper we report an implementation of Raman scattering in a vector radiative transfer (VRT) code for coupled atmosphere and ocean systems (CAOS) based on the successive order of scattering (SOS) method [22,23]. The polarization of the radiation field due to both the elastic and inelastic (Raman) scattering is fully accounted for by employing Mueller matrices for both processes. All four Stokes parameters can be simulated for detectors located at flexible altitudes in the atmosphere or depths in the ocean with the dependence of the viewing zenith and azimuth angles fully resolved. This paper is organized in the following way: Sec. 2 describes the inelastic vector radiative transfer equation; Sec. 3 outlines the optical properties of ocean waters including Raman scattering parameters; Sec. 4 shows our computational scheme using the SOS method; Sec. 5 applies the SOS code to study the contribution of Raman scattering to the polarized radiation field in the CAOS systematically; and Sec. 6 summarizes the conclusions.

Inelastic vector radiative transfer equation
The inelastic scalar radiative transfer equation for a plane-parallel system without thermal emission [1,13] is generalized into a vector form: where L = (I, Q,U,V ) T is the Stokes vector which contains the four Stokes parameters; the superscript T means the transpose of the vector; z is the vertical depth location; µ = cos(θ ) and θ is the viewing zenith angle; φ is the viewing azimuthal angle. The zenith angle θ < 90 • means that the direction is pointing upward and θ > 90 • is for downward. Both µ and φ pertain to the propagation direction of radiance L. For clarity, Fig. 1 shows the definition of the solar zenith, viewing zenith, and viewing azimuth angles used in this work. The solar zenith angle θ s does not show up explicitly in Eq. (1), but it is implicitly included in the elastic and inelastic source matrix S and S R when L is decomposed into a sum of direct solar and diffuse scattered radiances. The symbol c is the extinction coefficient at wavelength λ , which is often called the attenuation coefficient in ocean optics.
The viewing geometry with the solar zenith, viewing zenith, and azimuth angles defined. The plane contains the x-axis and zenith is the principal plane, which also contains the sun.
The symbol S is the elastic source matrix: where b(z, λ ) is the elastic scattering coefficient; µ and φ denote the incident direction of a beam of light; P is the elastic phase matrix defined as: The two rotation matrix operations R(π − i 2 ) and R(−i 1 ) are to ensure the consistency of the reference plane of the Stokes vector before and after scattering; i 1 and i 2 are the rotation angles [24]; F is the scattering matrix for the turbid medium; the scattering angle Θ is a function of variables µ, φ , µ , and φ . The first element F 11 of F is the phase function. Note that the product b(z, λ )F 11 is referred to as the volume scattering function with the unit of m −1 Sr −1 in ocean optics and F 11 satisfies the following normalization condition: Note that in the atmospheric radiative transfer community this normalization is often set to 4π. If that convention were adopted, there would be an extra factor of 1/(4π) at the right hand side of Eq. (2). The last term S R in Eq. (1) is the inelastic Raman source matrix: where b R (z, λ , λ e ) is the spectral Raman scattering coefficient; λ e is the excitation wavelength. Note that the unit of b R (z, λ , λ e ) is m −1 nm −1 , while the elastic scattering and absorption coefficients are in the unit of m −1 . P R is the Raman phase matrix defined in a way similar to Eq. (3) with F replaced by F R . The {1, 1} element F 11,R of F R also satisfies the normalization condition similar to Eq. (4) The extinction coefficient c(z, λ ) can be expressed as: where a e (z, λ ) is the elastic absorption coefficient due to attenuation at wavelength λ , a R (z, λ ) = b R (z, λ , λ )dλ is the inelastic Raman absorption coefficient, where λ and λ are the excitation and emission wavelengths, respectively. The symbol a R is used to show the fact that the radiation energy is absorbed in the excitation wavelength λ . In practice, a e and a R are not separated in absorption measurements.
In this work we assume the inelastic scattering in the atmosphere is negligible, which means S R (z, µ, φ , λ ) is zero in the atmosphere. The only inelastic scattering process considered is Raman scattering in ocean.

Optical properties of ocean waters
There are two sets of optical properties involved in this problem: elastic and inelastic scattering properties. For the elastic scattering process, we need to determine the scattering and absorption properties for both the atmosphere and ocean. While the radiative transfer model presented here has the full capability of handling aerosols, clouds with vertical inhomogeneity, we have assumed the atmosphere is both aerosol and cloud free, i.e., only molecular scattering is considered, to simplify the discussion. Trace gas absorptions are also omitted without loss of generality. The Rayleigh scattering matrix with a depolarization of 0.0284 is used for molecular scattering [25,26]. The column molecular scattering optical depth is calculated with the method in [26] with molecular density determined from the 1976 US standard atmosphere [27]. The ocean surface is rough with its wave slope distribution correlated to wind speed [28]. The effect of wind direction on the ocean surface roughness is neglected which may have important implication in some applications [29]. This is however out of scope of this study.
The ocean can be vertically inhomogeneous in this model. It is, however, assumed to be vertically homogeneous in this study and the z dependence of the ocean inherent optical properties is dropped for simplicity. We have also adopted the scheme in which the inherent optical properties are parameterized in terms of chlorophyll a concentration [Chla] so that all the optical properties will have an additional dependence on [Chla]. The elastic optical properties of ocean waters are described in [30]. In open ocean waters, far from terrestrial influences, the inherent optical properties are often described by a three-component model that includes pure seawater, phytoplankton particles and other associated particles, and colored dissolved organic matter (CDOM). The total absorption coefficient at λ is assumed to follow: where the pure water absorption spectrum a w is from Pope and Fry [31]. Note that where a we (λ ) is the elastic part of water absorption. With Eq. (8) a e (z, λ ) can be written as: The CDOM absorption coefficient is [32]: The particle absorption takes the following form: where the coefficients A p and E p were provided by Bricaud [33]. The A p and E p table contains data only from the wavelength interval of 400 nm to 700 nm. In order to evaluate the inelastic source matrix in the visible, the IOPs need to be expanded into the UV. For this purpose we have used another set of particle absorption data provided by Lee [34] which spans from 350 nm to 750 nm. In order to keep the two sets of a p data consistent, we rescale Lee's data by a p (λ = 440, [Chla]) calculated from Eq. (11) to enforce that they are equal at λ =440 nm. In the end we only use Lee's scaled data in the UV region. For all other IOPs the parameterization equations summarized below (Eqs. 12 to 19) are used for both the visible and UV spectrum regions. It is assumed that the scattering coefficient of CDOM is zero. The total elastic scattering coefficient b is then the sum of two terms: The pure seawater scattering coefficient The salinity and temperature dependence of b w [35] is neglected for simplicity. The symbol b p is the scattering coefficient of phytoplankton and their covariant particles [30,36,37]: where The scattering phase function of the phytoplankton particle F 11,p (Θ) is chosen to be the Fournier-Forand function [38,39]. The algorithm in [16] is adopted to determine F 11,p (Θ) conveniently by the the backscattering fraction B bp ([Chla]): which is spectrally independent. The total elastic phase function is the pure water phase function: The elastic scattering matrix of the ocean water is F = F 11 ·F, whereF is the average reduced Mueller matrix for ocean waters from Voss and Fry [40]. The scattering range in [40] is only from 10 • to 155 • and we have used the parameterization by Kokhanovsky [41] to extend the scattering angle range to [0 • , 180 • ] The Raman absorption coefficient is [7]: The spectral Raman where f R (λ , λ ) is the Raman wavelength redistribution function. The function f R (λ , λ ) is a sum of four Gaussian functions with center frequency shifts all close to 3400 cm −1 [1,2]. The Raman scattering matrix is written as F R = F 11,R ·F R , where F 11,R is the Raman phase function [9,13]: where ρ = 0.17 is the average depolarization ratio for Raman scattering [42]. The reduced Mueller matrix for Raman scatteringF R is [11]: whereF other,R denotes any other element not listed otherwise.
To understand the polarization features of inelastic and elastic scattering, Fig. 2 shows the comparison of the elastic and inelastic reduced Mueller matrix elements used in this study. The elastic elements are from [40] and the inelastic elements are calculated using Eqs. 21. Among these the (1, 2) Mueller element is especially of interest in remote sensing because it directly transforms the unpolarized radiance into polarized components upon scattering. The resultant polarized field contains retrieval information on the underwater particle size [43]. Figure 2 shows that the peak value of elastic |F 12 | is higher than that of the inelastic (1, 2) element. This can be attributed to the depolarization ratio value of ρ=0.17 for Raman scattering, which is larger than that of the elastic scattering. Since the (1, 2) element transforms the unpolarized light into linearly polarized light, we would expect that the Raman scattering contribution should decrease the degree of linear polarization of the total radiance field. Also the (2, 2) element changes the Stokes parameter Q after a scattering event. The inelastic (2, 2) element is smaller than the elastic element, which also contributes to smaller Q values.

Computational scheme
To solve the inelastic radiative transfer equation, one needs to know the inelastic source term Eq. (5), which is an integration over excitation wavelengths weighted by f R (λ , λ e ). The function f R (λ , λ e ) satisfies the normalization condition [1]: For an emission wavelength λ at which we wish to evaluate the total radiation field, the range of the excitation wavelength λ e can be found with all λ e satisfying f R (λ , λ e ) = 0. The function f R (λ , λ e ) only depends on the frequency shift δ ν = ν e − ν, where ν = 10 7 /λ is the wavenumber with units of cm −1 and λ is in units of nm [1]. Furthermore, f R (δ ν) is a sum of four Gaussian function with center (peak) frequency shifts around 3400 cm −1 . We choose the frequency shift grid to consist of eight frequency points starting from 2960 cm −1 and ending at 3800.0 cm −1 with equal increments of 120 cm −1 between consecutive grid points. Numerical tests show that this option will lead to the normalization error of |1 − f R (λ , λ e )dλ e | smaller than 0.3%. For a desired emission wavelength λ , the excitation wavelengths are determined from the eight frequency shift grid points. The elastic vector radiative transfer is then solved at each of the excitation wavelengths using the SOS method for coupled atmosphere and ocean systems [22,23]. All the optical properties of ocean waters and atmosphere are determined by the spectral models of these properties as described above. Note that we have ignored the Raman scattering contribution to the excitation wavelengths, which has negligible effects at the final emission wavelength. After the polarized radiation fields at those excitation wavelengths are obtained, the inelastic source matrix S R is evaluated using the known excitation polarized radiation field and Raman scattering properties. Finally the SOS method is used again to solve Eq.

Results and discussion
In this section the inelastic vector radiative transfer solution described in Sec. 4 is used to simulate the polarized radiation field in ocean waters. The relative importance of Raman scattering The ratio of the inelastic scattering radiance to the total radiance as in Fig. 3 Fig. 3(a) shows the radiance in the nadir-viewing direction just below the ocean surface as a function of wavelength; and Fig. 3(b) shows the percentage ratio of the Raman scattering contribution (The radiance difference between the cases with and without Raman scattering) to the total radiance field. The portion of Raman scattering contribution to the total radiance shown in Fig. 3(b) increases as wavelength increases until λ ≈ 500 nm, beyond which the ratio is roughly 25% and mostly flat. The importance of Raman scattering decreases as [Chla] increases with a residue of a few percent for high [Chla] values. The radiance spectral behavior in Fig. 3(a) is generally consistent with the spectral shape of water leaving radiance in the literature, specifically, Fig. 1a of [15] and Fig. 3 of [13]. However, the ratio in Fig. 3(b) is higher than that in Fig. 1 in [15] but similar to those ratios reported in [13]. The results in Fig. 1 of [15] were based on simulations with Hydrolight 4.1 (while the rest of the results in [15] are based on Hydrolight 5 [45]). We attribute the discrepancy to different Raman scattering parameters in the two Hydrolight versions as well as different default chlorophyll-based bio-optical model used in the two Hydrolight versions. For the case of pure ocean water, the radiance percentage ratio largely agrees with [13] with an exception at 443 nm (see Fig. 3 of [13], showing a peak at 443 nm). Our results and the results based on Hydrolight does not have a peak between 400 nm and 450 nm (see Fig. 2b in [15]). The cause of the peak at 443 nm in [13] is unclear to us. Figures 4 show the radiance at 560 nm in ocean water as a function of the viewing zenith angle θ for three detector depths: 0 (just below the ocean surface), 20 meters, and 40 meters; and for four different [Chla] values. The viewing azimuth angle is zero, i.e., the radiance is plotted in the solar principal plane. The viewing zenith angles of θ = 0 • and θ = 180 • correspond to the nadir-viewing and zenith-viewing directions, respectively. For all four different [Chla] values, the angular features of radiance just below the ocean surface are similar, though the upwelling Raman scattering radiance is slightly higher for clearer waters (smaller [Chla] values). This is again consistent with and explained by Fig. 3(a)  rect solar radiances directly refracted by the wind-roughened ocean surface. Outside the direct solar region, the radiance is mostly isotropic for downwelling radiance, but slowly decreases as the viewing zenith angle approaches 90 • . At the viewing zenith angle of 90 • the radiance angular distribution shows a discontinuity due to the total internal reflection by the rough surface. The upwelling radiance variation is small in comparison with the downwelling radiance. As the detector depth goes deeper, these interesting surface phenomena gradually fade away and the direct solar peak is less pronounced due to attenuation. The radiance peak will migrate to θ = 0 • due to preferential attenuation of photons with longer pathlength. Note that the Raman contribution for the clear water case is larger than the elastic radiance by several orders of magnitude at depth. As [Chla] increases, the impact of Raman scattering decreases, especially for detectors at deeper locations, except for upwelling irradiance/radiance where the Raman signal could become dominant at wavelength of strong attenuation. Similar findings were published in [8] for irradiance. Figures 5 are the same as Fig. 4 except that they are for the wavelength of 660 nm. The near surface detector response is similar to those at 560 nm. As detectors go deeper, the radiance decreases faster than the case at 560 nm due to a larger water absorption coefficient. Correspondingly, the Raman scattering contribution is larger than the case of λ =560 nm, and the Raman scattering contribution is larger than the elastic scattering even for [Chla]=1 mg m −3 . For the detector at 40 meters, the total radiance field is completely isotropic.  scattering, this isotropy is not achieved for depth of 40 m, at which the radiance level is of order 10 −9 W m −2 µm −1 Sr −1 . Figure 6 shows the degree of linear polarization (DOLP) Q 2 +U 2 /I at 560 nm as a function of viewing zenith angle in the principal plane. For the detector just below the ocean surface, the DOLP shows variations due to the surface interaction with the radiation field. The upwelling DOLP is relatively smooth and the peak is at the viewing zenith angle around 60 • , which is close to the scattering angle Θ = 90 • relative to the directly transmitted solar radiance. This is consistent with Fig. 2 which shows the peak ofF occurs at Θ = 90 • . At the viewing zenith angle of 90 • , the slope of DOLP, as a function of viewing zenith angle, suddenly changes sign and shows a relatively flat region. This is due to the upwelling polarized radiance reflected by the rough ocean surface. If it were a flat surface, the angular distribution would be symmetric near the viewing zenith angle θ = 90 • due to total internal reflection. For this rough surface case, the slope of DOLP is flipped and the nearly constant DOLP is due to the depolarization of the rough surface. Note this happens in the Fresnel cone zone [46,47] determined by the critical angle, which is defined by 180 • − arcsin(1/n w ) ≈ 132 • and n w ≈ 1.34 is the refractive index of sea water. Outside the Fresnel cone, only the light totally reflected by the surface exists for the detector just below the ocean surface. Note that for a rough ocean surface the edges of the Fresnel cone is diffused so there is no clear distinction at the edge. Once the viewing angle moves within the Fresnel cone, the observed radiation field is a mix of sky light transmission through and ocean light reflected by the surface. Complicated patterns show up due to multiple   Fig. 6. The degree of linear polarization as a function of the viewing zenith angle at λ = 560nm in the principal plane for the system shown in Fig. 4 scattering. There are two peaks between the viewing zenith angle of 140 • and 180 • . The first peak between 140 • and 160 • is actually a continuation of the trend for the viewing zenith angle smaller 140 • , where the Stokes parameter Q decreases to zero and becomes negative. This sign information is lost in the DOLP plot. In addition, in this region the sign of Q for the direct solar transmission light and the multiply scattered light in the ocean is different. The sign of the total Q field depends on the relative magnitudes of these two components, which are in turn dependent on the viewing zenith angle. As the detector goes deeper, the surface influence on the radiance and DOLP gradually diminishes due to multiple scattering and absorption. Figures 6 show the DOLP for deeper ocean is also smoother. The Raman scattering contribution decreases the DOLP significantly for the low [Chla] cases. The influence is larger for greater ocean depth. The DOLP difference between the two cases (with and without Raman) is more than 0.2 for [Chla]=0.03 mg m −3 and an ocean depth of 20 m. The difference is observable even for [Chla]=0.3 mg m −3 . Raman scattering has slightly shifted the angular position of the DOLP peak, which was not observed in [11] due to their azimuth angle averaging schemes. Figures 6 and 7 show that Raman scattering has a greater influence on the DOLP at 660 nm than that at 560 nm. For [Chla]=1 mg m −3 , the difference between the cases with and without Raman scattering is pronounced at 660 nm, while the difference is much smaller at 560 nm for the same [Chla] value. The angular distribution of DOLP is anisotropic even for the detector depth of 40 meters. This means that the radiation field is not isotropic in terms of polarization, even though the radiance field indicates so. Moreover, the DOLP difference at the zenith is always small between cases with and without Raman  scattering. These features will have profound impacts for development of ocean water remote sensing algorithms using polarization [43,48]. Another important polarization parameter is the Orientation of the Polarization Ellipse (OPE) ϕ [49]. The OPE is normally defined as ϕ = arctan(U/Q)/2. The range of OPE is between 0 and 180 • [49], while the range of arctan(U/Q)/2 is between −90 • and 90 • . Ideally one needs to adjust ϕ values to be within [0 180 • ] according to the signs of Q and U. However, that adjustment could introduce spurious patterns of OPE due to the extremely small values of Q and U at depth, while the arctan is relatively stable and easy to understand. Thus we have just focused on ϕ = arctan(U/Q)/2 here. Figures 8 show the OPE of the upwelling radiation field at the detector depth of 40 m for the same case as in Fig. 7(a). The OPE is plotted in the polar coordinate system with the viewing zenith and azimuth angles indicated by the radial and angular coordinates, respectively. The viewing zenith angles are shown by the numbers along the horizontal radius; whereas the viewing azimuth angles are indicated by the numbers along the outmost semicircle. The center of the semicircles represents the zenith. The OPE is zero in the principal plane (the viewing azimuth angle of 0 • in Fig. 8). This is consistent with the fact that the system is symmetric about the principal plane, which leads to U = 0 at φ =0. The OPE patterns in Figs. 8(a) and 8(b) are largely similar though small differences are observed. At first glance this may be surprising, it is however understandable considering Fig. 2, in which the elastic and inelastic Mueller matrix elements have very similar behaviors as a function of scattering angle. Different from the DOLP, which shows the importance of the polarization relative to the radiance, the information of OPE is the relative magnitude of Q and U, which is largely determined by the reduced Mueller matrices in deep ocean. This similarity may help sea animal navigation using polarization [50].

Conclusion
This paper reports a solution of the inelastic vector radiative transfer equation using the SOS method. The inelastic process considered is Raman scattering. The elastic radiative transfer model is used to solve the polarized radiation fields at those excitation wavelengths, which are then integrated to find the Raman scattering source matrix at the emission wavelength. The elastic radiative transfer model is used again to solve the radiation field at the emission wavelength, with the known inelastic scattering source matrix. All four Stokes parameters representing the polarized radiation field are conserved in this procedure. The angular distribution of the radiation field is fully resolved in terms of both the viewing zenith and azimuth angles. The radiative transfer model is used to systematically study the contribution of Raman scattering to polarized radiation field in ocean waters. It is found that the impacts of Raman scattering is greater for clearer waters (small chlorophyll a concentration values) and longer wavelengths, which is consistent with the findings in the literature. At the surface, the radiation field is influenced by directly transmitted solar radiance, total internal reflection, rough surface, and other surface phenomena. As the detectors go deeper in the ocean, these surface phenomena gradually diminish and the angular distribution of the radiance field becomes smoother and more isotropic. In addition, Raman scattering causes the radiance field to be isotropic for a detector at 40 meters below the ocean surface at the wavelength of 660 nm. However, if polarization is considered, the isotropy is not achieved at the same level. Generally, Raman scattering decreases the DOLP in ocean waters due to a significant depolarization ratio for Raman scattering. The DOLP at the zenith is small, however, for cases both with and without Raman scattering. The OPE shows the relative magnitude of Q and U, which is largely determined by the reduced Mueller matrices. The OPE angular patterns are similar for elastic and inelastic scattering in deep ocean, due to the similar behavior of the elastic and inelastic reduced Mueller matrices. Analyses of polarized light fields in the ocean and from remote sensors have been expanding, in particu-lar in response to the success of polarized remote sensing of the atmosphere in constraining aerosol properties. It is expected that our model will be useful in understanding how hydrosols affect the above-and under-water radiometeric fields and their polarization properties, thereby similarly advancing hydrosol inversions from in-situ and remote measurements.