A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method

: A vector radiative transfer model has been developed for coupled atmosphere and ocean systems based on the Successive Order of Scattering (SOS) Method. The emphasis of this study is to make the model easy-to-use and computationally efﬁcient. This model provides the full Stokes vector at arbitrary locations which can be conveniently speciﬁed by users. The model is capable of tracking and labeling different sources of the photons that are measured, e.g. water leaving radiances and reﬂected sky lights. This model also has the capability to separate ﬂorescence from multi-scattered sunlight. The δ - ﬁt technique has been adopted to reduce computational time associated with the strongly forward-peaked scattering phase matrices. The exponential - linear approximation has been used to reduce the number of discretized vertical layers while maintaining the accuracy. This model is developed to serve the remote sensing community in harvesting physical parameters from multi-platform, multi-sensor measurements that target different components of the atmosphere-oceanic system.

The solution for the RTE in a coupled Atmosphere -Ocean System (AOS) is particularly of interest in the ocean color remote sensing, the general circulation model, and many other studies.There are a few radiative transfer codes currently available for an AOS [26,27,32,40,57,58,59,60,61,62,63,64,65].Based on the matrix operator method, [26] and [27] solve the scalar RTE in an AOS with wind induced rough ocean surface.The multicomponents method developed by Zege et al. [32,33] solves the RTE in an AOS including polarization.The model by Chami [40] is based on the SOS method which also includes polarization.A flat ocean surface is assumed and the circularly polarized component is ignored in this model.[60,63] solve the scalar RTE using the DISORT method.Bulgarelli et al. has developed finite element method for solving the scalar RTE in an AOS with a smooth ocean surface [61].Another important work is the Hydrolight code based on the invariant embedding method by Mobley [4,59].Polarization is not included in the Hydrolight.The model by Chowdhary et al. solves the vector radiative transfer equation in an AOS based on the adding method [64].The rest of the list are mainly based on the Monte Carlo method [54,57,58,62,65].The Monte Carlo technique is slow for remote sensing purposes because it requires an accumulation of statistics.In our classification, the models in [20,45,46,47] do not solve the RTE in an AOS because they either treat the ocean surface as bidirectional reflectance surface or use a semi-analytical model for the ocean [66], which is not a real simulation of light propagation in ocean.
Most discrete ordinate methods (excluding the Monte Carlo method) apply the Fourier transformation to the azimuth dependence of the Stokes vector in RTE in order to reduce the dimension of differential equations, as the computation time for eigen solution is proportional to the third power of the dimension.For each Fourier component, Gaussian quadrature scheme is used to discretize the cosine of the zenith angle μ = cos θ for numerical integration of the scattering term in RTE.At an ocean surface, the refracted light in the ocean from the atmosphere is confined within the Fresnel cone, which only cover part of the downward solid angle.To match the boundary condition at the ocean surface, two sets of quadrature points are usually selected for ocean [32,60,63].One corresponds to the refracted directions of the quadrature points in the atmosphere.The other covers the region outside the Fresnel cone.While the two sets of Gaussian quadratures in water are required in most models that treat the RTE as eigenvalue problem, it is not entirely clear whether that is a neccessity for the successive order of scattering approach.While less accurate, relaxing to one set of Gaussian quadrature may help reducing computational time and memory requirement.
For optically thin media (such as clear sky with moderate aerosol loading) and for highly absorbing media (such as ocean), the SOS method has advantages when compared with other methods.First, each order of scattering has clear physics interpretation.This results in efficient programing.For example, because of the fact that the source term associated with the refracted light into the water can be formulated analytically, two sets of quadrature points for the ocean are not mandatory.The same quadrature points can be used in the ocean just as in the atmosphere.In other words, the quadrature points in the ocean do not have to match Snel's law [67] with the quadrature points in the atmosphere, as long as the angle change associated with refraction is properly represented in the scattering source term (scattered light at arbitrary quadrature angle).The use of one set of Gaussian quadrature points could decrease the memory requirements for the ocean by a rough factor of 1/4 because a major part of the memory in the computer code is allocated to store the phase matrices, whose number of elements are proportional to P 2 , where P is the number of quadrature points.The CPU time reduction, however, depends on the specific program implementation.
This paper documents the physics principle and detailed mathematical derivation of a Fortran 90 computing package for solving the vector RTE based on the SOS method.The current code is developed for a flat ocean surface.However, the technique can be extended to rough wave surface model [68].Both the single and double ocean quadrature versions of the code are available to the public science community by contacting the corresponding author.
Section 2 outlines the vector RTE, while section 3 describes the SOS formulas for solving the vector RTE.Section 4 presents the treatment of the air-sea interface boundary condition.Section 5 shows the numerical procedure used in the code.Section 6 provides the validation and simulation results.The final section, section 7, presents the conclusions.

The vector radiative transfer equation
A beam light of arbitrary polarization may be represented by the Stokes vector L = (I, Q,U,V ) T , where the four elements I, Q, U, and V are the Stokes parameters; the superscript T stands for the transpose of the vector.A formal definition of the Stokes parameters in a medium may be found in [7] (see Eq. 2.6.4).The mathematical definitions of the Stokes parameters depends on propagation direction, electric vector components, as well as the dielectric properties of the medium.The vector RTE which describes multiple light scattering in a plane-parallel medium can be written as : where τ is the optical depth defined as the product of the extinction coefficient β e and the path lengh l; μ = cos(θ ) and θ is the zenith angle; φ is the azimuthal angle.All the zenith angles in this paper are defined in terms of the upward direction.θ < 90 • means that the direction is pointing upward and θ > 90 • when pointing downward.Both μ and φ are pertained to the propagation direction of radiance L. S is the source matrix: where ω = β s /β e is the single scattering albedo and β s is the scattering coefficient.Mishchenko et al. have shown that Eq. ( 2) is only valid in a macroscopically isotropic and mirror symmetric medium [7].In this equation and hereafter the dot product is assumed between the matrix P and vector L. The same assumption applies if matrices are written next to each other.P is the phase matrix defined as: The rotation matrix R rotates the reference plane of the Stokes vector between the meridian plane and the scattering plane; i 1 and i 2 are the rotation angles [1,69]; F is the scattering matrix for the turbid medium; the scattering angle Θ is a function of variables μ, φ , μ , and φ .Equation (1) does not contain thermal emissions.The integral form of Eq.( 1) is needed to apply the SOS method.Standard manipulation leads to [1,5]: where τ l and τ u are the lower and upper limits of the medium under consideration, namely, τ l < τ < τ u .For the atmosphere, τ l = 0 and τ u = τ * a are the top and bottom of the atmosphere, respectively.For the ocean, τ l = τ o and τ u = τ * o are the top and bottom of the ocean, respectively.Please note that τ o = τ * a + ε is just below the ocean surface, where ε is an infinitesimal number.In other words, the ocean surface is assumed to be a degenerated layer with no geometric thickness.The only impact of the ocean surface to the system is through the refraction and reflection of the dielectric interface, which is governed by Snel's law and Fresnel formulas [57,65].L(0, μ < 0, φ ) is 0; and L(τ * a , μ > 0, φ ) includes the light reflected by the dielectric interface from the downward light in the atmosphere and the light transmitted through the interface from the ocean.L(τ o , μ < 0, φ ) includes the light reflected by the dielectric interface from the upward light in the ocean and the refracted light through the dielectric interface from the atmosphere; and L(τ * o , μ > 0, φ ) is the reflection of the ocean bottom.

The successive orders of scattering method
In the SOS method, the total radiance vector is expressed as a summation of contributions from photons scattered a number of times, from 0 to a maximum number of N [39,47,70,71], Substitution of Eq. ( 5) into Eqs.(4a) and (4b) leads to: where It should be noted that the direct solar light has been taken out of Eqs.(6).The source terms caused by the direct solar light have been integrated and written out explicitly as S 1 in Eqs.(7).
The next step is to expand the radiance vectors and the phase matrices into Fourier series of the azimuth.The expansion of the radiance vectors can be written as: with The expansion of the phase matrices can be written as: In the pioneer works [72,73,74], the analytical expressions of P m sin and P m cos have been found in terms of the generalized spherical functions (GSF) [75].Lenoble [47].Not trying to create our own convention, we find it is more convenient to use the Wigner d functions [76] as the expansion basis.The link between the generalized spherical functions and the Wigner d functions has already been suggested by Hovenier and van de Mee [74].
Previously, the Wigner d functions have been used as the basis to expand the scattering matrices by Mishchenko et al. [77].The Wigner d functions are also real and better documented than the generalized Legendre functions.More importantly, they can be calculated via stable up recurrence relations.Interested readers may refer to Appendix B of [77] for detailed properties of the Wigner d functions.The derivation of P m sin and P m cos in terms of the Wigner d functions is given in the Appendix A.
The standard procedure includes substituting Eqs.( 9) and ( 10) into Eqs.( 7) and ( 8), carrying out the integration over φ , and identifying the same terms containing cos[m(φ − φ 0 )] and sin[m(φ − φ 0 )] at both sides of the equations.The dependence of φ is then reduced to a set of equations with variables of τ and μ and order number m.The equations with different order m are completely decoupled and can be solved separately.For details, please refer to [47,73,74,78,79].After lengthy manipulations, the resultant formulas are: where and element is reduced to a constant, i.e., the single scattering albedo; other elements are zero.The dummy integration variable μ is used for later convenience.Equations ( 11), (12), and ( 13) are the equations to be used in the SOS method.From a mathematical point of view, the radiances with the factor of n 2 w in Eqs.(13) need to be treated with extreme care.At the ocean side, these terms are the contribution transmitted from the atmosphere, confined within the Fresnel cone.Thus extra integration quadrature points are generally needed to cover the region outside the Fresnel cone, which increases both the memory requirement and the CPU time for the computation.In the next section, a single quadrature scheme is presented based on the separation of directly transmitted sky light and the diffuse light in ocean.The double quadrature scheme is also implemented in the SOS code.Numerical comparison is made in section 6.

The boundary condition at the ocean surface
The starting point is to evaluate the source function Eq.(12c) with the boundary conditions of Eqs.(13).In the following description, we assume that the source function in the ocean is to be found.Substituting Eq.(13e) in Eq. (11a) leads to: where the first term in Eq.( 14a) is the direct contribution from photons transmitted into the ocean from the atmosphere.This part of the radiance field is confined within the Fresnel cone.The total radiance field as shown in Eq. ( 14a) is not continuous because the first term jumps from a finite value at the edge of the Fresnel cone to zero just beyond the Fresnel cone.In the case of wind roughed ocean surface, this discontinuity is partly smoothed out but the rate of change of the radiance field near the critical angle will still be large.The term L m n,di f f (τ, μ) represents the diffuse radiance field.The factor r(π − θ ) also changes abruptly around the edge of the Fresnel cone because of the total internal reflection for incident angles larger than the critical angle.Both these two parts contribute to the discontinuity around the critical angle.The quadrature techniques, however, requires the function to be smooth in order to get an accurate integration.This property has prompted people to use two sets of quadrature points to do the integration in Eq. (12c).One set of μ covers the Fresnel cone and the other covers the rest of the region.
Further analysis finds that the integration of the first term of Eq. (14a) can be done by transforming the integral variables μ , which is below the ocean surface, to the variable μ , which is above the ocean surface.Specifically, sin(π − θ ) = n w sin(π − θ ) and cos(θ )dθ = n w cos(θ )dθ lead to the following equation: This relation is substituted into the source term caused by the direct sky light transmission from where μ c = cos(θ c ) and θ c = π − arcsin(1/n w ) is the critical angle which defines the edge of Fresnel cone; μ q and w q are the quadrature points and weights in the atmosphere, respectively.The summation in the last step is over half of the quadrature points, which covers the range of −1 ≤ μ ≤ 0. In the simulation, only P m (τ, μ p , μ q ) will be pre-calculated and stored in memory, where μ p,q are the quadrature points.In the single quadrature scheme, this means P m (τ, μ p , μ q ) will not be found because μ q is not one of the quadrature points in the ocean.The solution is simply to use linear interpolation between the two nearby quadrature points.Now it is assumed that the discontinuity caused by r(π −θ ) in the diffuse light is insignificant compared to the total radiance field.Thus the integration of the diffuse light in Eq. (14b) can be done using a single quadrature scheme: where μ p,q and w q are the quadrature points and the associated weights; P is the order of the quadrature scheme.The total source term is: Equations ( 16) and (17) only use one set of Gaussian quadrature points in the ocean, which leads to less memory usage and CPU time.By assuming insignificant contribution of r(π − θ ) in the diffuse light, the question remains as to how large the error will be.We have also implemented the double quadrature scheme in the ocean.Numerical comparison shows that the difference between the two schemes is negligible.
The source terms by the direct transmitted light from the atmosphere is taken care by Eq. ( 16).If the radiance field in the ocean needs to be exported at the end of the simulation, a detector array is defined to record the first term Eq.(14a) separately.The separation of the direct and diffuse contribution to the detector offers another benefit.If only the second term in Eq.(13c) is substituted in Eq.(11b), and the optical depth integration term is ignored, the resultant equation is the water leaving radiance: where the subscript wlr stands for the water leaving radiance.The water leaving radiance is important for the retrieval of the inherent optical properties (IOP) of ocean waters.ignores the multiple interaction between the diffuse light in the atmosphere and ocean.Equation (19), however, provides an exact solution for the water leaving radiance without causing any extra computational burden.

Numerical procedure
The atmosphere and ocean are discretized into K a and K o homogeneous layers, respectively.The optical thickness of each layer needs to be small enough to meet the required numerical accuracy.For the atmosphere, the optical depth at each layer is denoted as τ k with k = 0, 1, 2, ..., K a .τ 0 = 0, τ 1 = τ 0 + δ τ,...,τ k = τ k−1 + δ τ,...τ K a = τ * a .This means that the kth layer is between τ k−1 and τ k .Practically, δ τ does not have to be equal for each layer and can be adjusted for convenience.For the ocean, The zenith angle cosine μ is discretized into the Gaussian quadrature of even order P. μ p and w p are the quadrature points and weights with p = 1, 2, ..., P , respectively.P m (τ, μ p , μ q ) for each p and q are calculated and stored in memory before the SOS iteration.Even P m is a function of τ while many layers with different τ may share the same P m .Thus an index is created for each layer.The layers with the same scattering properties are labeled with the same index numbers.These index numbers will be used to refer to the correct P m when needed.
The source functions S m 1 (τ k , μ p ) for the atmosphere and ocean are initialized according to Eqs.(12a) and (12b), respectively.The SOS iteration begins with Eqs.(11) with known S m 1 (τ k , μ p ).At this point, the first terms in Eqs.( 11) are neglected temporarily.The integration for μ < 0 is done in the following way: The exponential -linear approximation has been implemented to evaluate the integration between τ k−1 and τ k [80].Equation (20) shows that L m n (τ k , μ p < 0) can be updated successively from k = 0 to K a for the atmosphere and from k = K a + 1 to K a + K o + 1 for the ocean.For μ p > 0, the procedure is: L m n (τ k , μ p > 0) will be updated successively from τ K a +K o +1 to τ K a +1 for the ocean and from τ K a to τ 0 for the atmosphere.Again the exponential -linear approximation has been implemented to evaluate the integration between τ k and τ k + 1 [80].After this step, all of the source terms should be reset to zero to prepare the source terms for the next order of scattering.
The next step is to consider the boundary conditions.The bottom reflection can be taken into account by adding Eq. (13f) to L m n (τ k , μ p > 0) according to Eq. (11b) for the ocean.For the interface, only the first terms in Eqs. ( 13c source integration in Eqs.(17).That is: where ⇐ is understood to evaluate the right hand side and give the resultant value to the left hand side; θ p = arccos(μ p ).After Eqs. ( 22), the single quadrature scheme uses Eqs.(17) to get the source function due to ocean diffuse light.The double quadrature scheme, however, integrates the source function in two separate regions, −μ c < μ < μ c , and |μ| > μ c .The integration of the second terms in Eqs.(13c) and (13e) involving n 2 w in both schemes is done using Eq.( 16) or its equivalent form for the transmission from the ocean to the atmosphere.The only difference is that the double quadrature scheme does not involve the interpolation of P m , because μ q in Eq.( 16) is one of the predefined quadrature points in the ocean.A separately allocated array L m n,det (τ k , μ p ) is used to store the radiance field for the final output: Equations ( 23) are done only at levels where the output is desired.Again, in the single quadrature scheme L m n (τ k , μ p ) is not directly available because μ p is not one of the quadrature points.The solution is also to interpolate the values at the two near points.This is also physically sound because L m n (τ k , μ p ) is generally smooth before interacting with the ocean interface.In addition, Eq. (13g) needs to be evaluated and stored to prepare the next order of scattering.After this step, L m n (τ k , μ p ) are set to zero.The steps from Eq.( 20) to here are repeated N times, where N is a preset value for the maximum order of scattering based on sensitivity studies.The N could also be determined based on convergence criterias during the run time [40].However, it is more attractive to use a preset value in terms of efficiency.After the iteration of the procedure for N times, the remainder is taken into account by a geometric series: Duan and Min have shown that this relation is satisfied with N ≈ 40 for most, if not all plane parallel systems [42].This requirement is smaller for turbid media with small asymmetry factors or single scattering albedos.This relation is assumed for all the four Stokes parameters unless L m n−1 (τ k , μ p ) = 0 for one of the elements Q, U, and V .Then the total radiance field is given by:

Validation and results
The validation of the SOS code is done with the Monte Carlo code by Zhai et   and accuracy.It can deal with both one-dimensional and three dimensional vector radiative transfer problems.Even it is well known that the Monte Carlo method converges slow, the error associated with the Monte Carlo method is easy to understand, i.e. ∝ N −1/2 photon , where N photon is the number of the photon bundles.As long as N photon is large enough, the code can produce reliable results.In the validation case, the atmosphere is chosen to be a conservative Rayleigh medium with optical depth of τ a = 0.1.The Rayleigh scattering matrix is: The ocean phase function is a Henyey-Greenstein (HG) scattering phase function [81]: where μ = cos(θ ) and g is the asymmetry parameter defined by:  The range of the asymmetry parameter is [−1, 1].g ≈ 1 sugguests the scattering function is strongly peaked in the forward direction.g ≈ 0 suggests the scattering function is nearly isotropic.In our calculation, we take g = 0.9185, which means the scattering function has a large forward peak.The HG phase function with g = 0.9185 has the same integrated backscattering coefficients b b as the Petzold phase function [82].The other scattering matrix elements for the ocean is determined by: where the superscript HG shows the matrix elements are with the HG scattering phase function while RAY shows that they are with the Rayleigh scattering phase function.In other words, the reduced scattering matrix elements of the ocean water are the same as those of the Rayleight scattering matrix [83].δ -fit technique for the phase function truncation has been extended to the full scattering matrix to fulfill the polarization calculation requirement in this code.The details are provided in Appendix B. The optical depth is 10.0 for the ocean.The single scattering albedo in ocean depends on the wavelength with a possible range of highly absorptive (ω ≈ 0.1) to nearly conservative (ω ≈ 0.9).The single scattering albedo has been chosen to be 0.5 as a representative value for the primary study in this paper.The bottom is a Lambertian surface with single scattering albedo of 0.5.The refractive index of the ocean water is 1.338.Figures 1 and 2 show the Stokes parameters I, Q/I, U/I, and V /I as functions of observation angles for the system at the top and bottom of the atmosphere, respectively.The observation angle is 0 • for the upward (from the bottom of the ocean to the top of the atmosphere (TOA) direction and 180 • for the downward (from the TOA to the bottom of the ocean) direction.The principle plane of the observation angles shown in the figures is obtained by rotating the solar incident plane by an azimuthal angle of 90 • clock-wise if looking in the zenith direction.The solid lines are calculated by the Monte Carlo method.The circle and cross symbols are computed by the SOS method with the single Gaussian (SG) quadrature and double Gaussian (DG) quadrature schemes, respectively.A solar source with an incident angle of 120 • is used.The principle plane for the observation angles is at an azimuthal angle of 90 • from the solar incident plane.For the Monte Carlo code, the total photon number used is 5000000.The Gaussian quadrature order of P a = 36 is used in the atmosphere.The ocean Gaussian quadrature order of P o = 72 is used for the DG case.The total number of scattering orders is N = 10.δ τ = 0.01 and 0.2 have been used for the optical depth integration of the ocean.δ τ = 0.01 has been used in the atmosphere if δ τ = 0.01 in the ocean.δ τ = 0.2, however, cannot be used for the optical depth integration of the atmosphere because it is larger than 0.1, the optical depth of the atmosphere.Thus the value δ τ = 0.05 is selected for the atmosphere in this case, which is the largest value one can select because there have to be at least two layers (K a = 3) in the atmosphere to solve for the three unknown coefficients in the exponential -linear approximation.The order of the Fourier expansion used is M = 8 in Eq. ( 9).The maximum order of the HG phase function expansion is L = 36 in Eq.(A-13).The δ -fit technique is generalized to the full scattering matrix for the polynomial expansion.The SOS code was run using an I-Mac with a 2.8 Ghz Intel processor.The fortran compiler is GFORTRAN with the −O3 compile option.Table 1 shows the average CPU times for both the SG and DG schemes for the cases shown in Figs. 1 and 2. Both the SG and DG quadrature schemes results agree well with the Monte Carlo results for all the four Stokes parameters at both the top and bottom of the atmosphere.Note that only the Stokes parameters within observation angles of 0 − 90 • are shown at the TOA because they are zeros for observation angles larger than 90 • .For both δ τ = 0.01 and δ τ = 0.2, the differences between the SG and DG are negligible, for the TOA and bottom of the atmosphere.For the TOA, the errors of the SG and DG scheme radiances with δ τ = 0.2 vary between 0.25% at 0 • to 1.0% at 82.6 • , if using the DG scheme with δ τ = 0.01 as a standard.The error at 87.5 • is -2.5%, which is slightly larger than smaller angles.The linear polarization components Q/I and U/I also have larger errors at larger observation angles.The differences of the circular polarization components V /I between the SOS code and Monte Carlo code are larger compared to those of I, Q/I, and U/I.The DG scheme has slightly better agreement with the Monte Carlo method.For the bottom of the atmosphere, the differences between the SG, DG, and the Monte Carlo results are smaller than those of the TOA.The errors of the SG and DG cases with δ τ = 0.2 are from -0.5% to 0.5%.Another interesting feature is that the ratios Q/I and U/I have smaller errors than the absolute radiance.Considering the large forward peak of the phase function, the results in Figs. 1 and 2 are quite impressive.
Figures 3 show the water leaving radiance for the same system as in Fig. 1 at the bottom of the atmosphere (just above the ocean surface).The radiance percentage differences between the SOS results and the Monte Carlo results are around 2%, which are larger than those of the total radiance plots.The reason is that the ocean has a large asymmetry factor g = 0.9185 and an optical depth of τ = 10, which is rather difficult for the Monte Carlo method to achieve convergence.The SOS method, however, virtually takes all order of scattering into account (see Eq. ( 25)).The better agreements between the two methods for the total radiance is due to the contribution of the atmospheric Rayleigh layer, which is thin and the phase function is symmetric around the scattering angle of 90  that the polarized components converge before the absolute radiance for both methods.The small value of the circular component makes it more sensitive to the optical depth integration step δ τ.Again, the difference between the SG and DG results are negligible.The water leaving radiance for δ τ = 0.2 has relative errors from −0.6% to −1.5%, if comparing to the results with δ τ = 0.01.The next case to be studied is designed to answer the following question.For a clear sky condition with an aerosol optical depth up to 0.3, how many orders of scattering is needed to achieve a radiance error smaller than 0.5% in an AOS?For this purpose, a layer of aerosol is added at the bottom of the atmosphere in the first case studied.The aerosol consists in a fine mode of spherical particles with log-normal size distribution [84].The modal radius is r m = 0.08μm and the variance is σ = 0.46.The refractive index is m = 1.45.The zero imaginary part of the refractive index leads to a single scattering albedo of 1.The incident wavelength is 0.443μm.The conservative Rayleigh layer of optical depth 0.1 is unchanged, located at the top of the aerosol layer.The aerosol layer at the atmosphere bottom does not mix with the Rayleigh layer.To make sure the standard converged results are obtained before the sensitivity study, an optical depth of τ sol = 0.3 is given to the aerosol layer.The ocean is also kept the same as the previous case.Figure 4 shows the effects of different Gaussian quadrature orders P and Fourier series orders M. The radiance and Q/I shown in the figure is at the TOA.A total scattering order of N = 40 is used for all the SOS calculation.The azimuthal angle for Fig. 4 is 180 • , because the Stokes parameter at this azimuthal angle contains minimum amount of sun glint, which is very important in the satellite remote sensing.The Gaussian order P shown in Fig. 4 refers to the Gaussian quadrature order used in the atmosphere.For the DG case, the corresponding Gaussian order used in the ocean is 2P.In all the cases, the differences between the SG and DG case are negligible.It is observed that the results with the Gaussian quadrature order of P = 18 and the Fourier series order of M = 8 are already converged to the results with P = 36 and M = 36.The Gaussian order P = 8 produces larger errors because P = 8 is much   Table 2 shows the order of scattering needed for convergence.A set of aerosol optical depth τ sol are used to test the sensitivity of scattering orders N. Other system configurations are the same as the previous cases.The radiance at the TOA are studied and the observation azimuthal angle is 180 • .Both the total radiance and water leaving radiance are tested.A tolerance level of ±0.5% is used in the calculation of Table 2.It is shown that scattering orders of N = 4 to N = 8 are needed to achieve errors smaller than ±0.5% for the total radiance.The scattering orders are from N = 8 to N = 10 are needed for the water leaving radiance calculation.The requirement for scattering orders N is higher for the water leaving radiance than the total radiance.The reason is that the total radiance consists of the contributions from both the atmosphere and the ocean.The atmosphere is a thin conservative medium in which a few order of scattering are enough for convergence.Thus the total radiance is easier to achieve the ±0.5% error tolerance level.However, the water leaving radiance comes only from the ocean.The ocean is a optically thick medium, which leads to more orders of scattering.The results with the scattering order of N = 10 have errors smaller than 0.1% for all the total radiance calculation.The maximum errors with N = 10 associated with the water leaving radiance are around 0.25%.
It should be emphasize that the SOS method will converge even for optically thick and near conservative media given a large enough order of scattering.Min and Duan have presented that the SOS method can converge for a turbid medium of optical depth up to 50, with a single scattering albedo of 0.99 [42].The limitations of the SOS code may be its computational speed, which can only be found by crosscheck with other methods.This will be an interesting work in the future.The last case we use 0.5% as a standard to study the how many orders of scattering are needed for converging a realistic atmosphere-ocean system.Applications of radiative transfer models for aerosol retrievals may require the polarized radiance measurement as accurate as 0.1% [85,86].As stated in the previous paragraph, the total polarized radiance converges within 0.1% at the TOA with the scattering order of N = 10 for the AOS studied.A higher single scattering albedo may lead to larger scattering orders.

Conclusion
This paper presents a vector radiative transfer model for an AOS.The model is based on the SOS Method.The SOS formulas have apparent physical significances and the errors are easy to be analyzed.The source terms in the RTE involves numerical integration with Gaussian quadratures.Previously published models normally use two sets of Gaussian quadrature points in ocean and one set in atmosphere to treat the interface boundary condition correctly.This paper has presented a Single Gaussian (SG) quadrature scheme in ocean in comparison with the Double Gaussian (DG) quadrature scheme.Although two sets of Gaussian quadratures is essential for discrete ordinate methods, the differences between the SG and DG schemes can be small for SOS in a typical atmosphere and ocean configuration.Both memory requirements and CPU time are smaller for the SG scheme, which can be helpful for remote sensing applications.The current SOS model (both SG and DG schemes) can differentiate water leaving radiance from light scattered by the atmosphere and by sea surface reflection of the sky light.The model results compare well against the Monte Carlo simulations.For a realistic case involving an ocean of optical depth τ = 10 and single scattering albedo ω = 0.5, the SOS converges with the order of scattering N = 10 and the Gaussian quadrature order P = 36.The δ -fit technique is generalized to the whole 4 × 4 scattering matrix to fulfill the polarization calculation requirement.This work is primarily of interest in aerosol and ocean color remote sensing, and other relevant scientific disciplines.Both the SG and DG scheme computing packages are available to the public science community by contacting the corresponding author.
where the order number s in [77] has been replaced with l; and the subscripts M are added to the coefficients in [77] for clarification.
In the visible range, scattering particles are often much larger than the incident wavelength.In turn, the scattering phase function a 1 usually has a strong forward peak, which leads to hundreds or thousands of expansion terms β l .In this case the truncation techniques are normally used to reduce the number of fitting coefficients significantly [87,88].In this paper the δ -fit technique is generalized to fit the scattering matrix [88], which is presented in Appendix B. If no truncation is needed, the method in [77] is used to calculate the expansion coefficients as in (A-12).Equations (10), which are also used in [47], are different from Eq. ( 6) of [73] by a factor of 2. Careful check finds that C m = 2P m cos and S m = 2P m sin , where C m and S m are the notations used in [73].These two factors of 2 cancel out and the whole formulas remain consistent.Now define P m = P m cos + P m sin D with D = diag{1, 1, −1.−1}.The expression of P m can be found by substituting Eqs.(A-6) and (A-3) into Eqs.(7a) and (7b) of [73], with the knowledge of which a factor of 2 has been canceled out: where L is the maximum order of l used for the expansion; d l m0 = d l m0 (μ), d l,+ m2 = d l,+ m2 (μ), and d l,− m2 = d l,− m2 (μ); similarly, d l m0 = d l m0 (μ ), d l,+ m2 = d l,+ m2 (μ ), and d l,− m2 = d l,− m2 (μ ).P m can be precalculated and stored in the computer memory before the SOS iteration, which will speed up the code significantly.The computation aspects and the recurrence relations of the Wigner d functions are documented in Appendix B of [77], which will not be repeated here.θ k is the scattering angle, w k is the weight for each scattering angle.By setting the weights for the forward-scattering angles to zero, the forward peak will be automatically truncated.The singular-value decomposition method is used to solve the equations ∂ ε/∂ β * l = 0.After a set of coefficients β * l are found, the renormalization factor is f = 1 − β * 0 .All the coefficients are divided by β * 0 to renormalize the fitted phase function.The δ -fit is proven to be accurate for radiance simulations [88].
This work deals with the polarized radiative transfer equation.Thus the δ -fit needs to be generalized to all the scattering matrix elements.To conserve the polarization properties of the radiance field, it is important to keep the reduced scattering matrix elements unchanged.Therefore the following procedure is proposed.
First, the expansion coefficients for the phase function are found by using the least squares fitting.The coefficients are substituted into Eq.(B-4) to produce the fitted phase function.The fitted scattering matrix elements are then obtained by: where F i j (θ k ) are the scattering matrix elements other than a 1 (θ k ) = F 11 (θ k ) in Eq. (A-8).For instance, F 12 (θ k ) = b 1 (θ k ), F 22 (θ k ) = a 2 (θ k ), ..., etc. Finally the expansion coefficients c l * i j for the elements F * i j (θ ) are found by solving: For mn = 00 and (i, j) = (4, 4), c l * i j correspond to the coefficients δ l in Eq.(A-9b).For mn = 20 and (i, j) = (1, 2), (4, 3), c l * i j correspond to γ l , and ε l in Eq.(A-10b), respectively.For mn = 22 and (i, j) = (2, 2), (3, 3) c l * i j correspond to ζ l and α l in Eq.(A-11b), respectively.Note that the weight is no longer needed because the truncation is already done by the least squares fitting of the phase function.To use the δ -fit coefficients, the optical depth and single scattering albedo for the layer have to be adjusted similarly as the Delta-M method:

)
In the previous equations, L m n = L m n,cos + L m n,sin and P m = P m cos + P m sin D with D = diag{1, 1, −1.−1}.The purpose of D is to flip signs for the upper right 2 × 2 sub-matrix of P m sin [47, 73, 74].r * m has been treated similarly as P m .For the Lambertian surface, its (1, 1) (C) 2009 OSA 16 February 2009 / Vol.17, No. 4 / OPTICS EXPRESS 2065

Fig. 1 .
Fig.1.The Stokes parameters at the top of the atmosphere.The solid lines are calculated by the Monte Carlo (MC) method.The circle and cross symbols are computed by the SOS method with the single Gaussian (SG) quadrature and double Gaussian quadrature (DG) schemes, respectively.The Gaussian quadrature order of P a = 36 is used in the atmosphere.The ocean Gaussian quadrature order of P o = 72 is used for the DG scheme.The Fourier series order used is M = 8.The solar incident angle is θ 0 = 120 • .The azimuthal angle of the observation is 90 • from the solar principle plane.The atmosphere is a conservative Rayleigh medium with an optical depth of 0.1.The ocean scattering function is a Henyey-Greenstein function with an asymmetry factor of 0.9185.The optical depth of the ocean is 10.0.The ocean bottom is a Lambertian surface.The single scattering albedos of both the ocean medium and bottom are 0.5.The total order of scattering is N = 10.

Fig. 2 .
Fig.2.The Stokes parameters at the bottom of the atmosphere (just above the ocean surface) for the same system as in Fig.1.

Fig. 3 .
Fig.3.The water leaving Stokes parameters at the bottom of the atmosphere (just above the ocean surface) for the same system as in Fig.1.

Fig. 4 .
Fig. 4. The effects of the orders of the Gaussian quadrature and Fourier series orders (C) 2009 OSA 16 February 2009 / Vol.17, No. 4 / OPTICS EXPRESS 2077 #105295 -$15.00USD Received 15 Dec 2008; revised 20 Jan 2009; accepted 26 Jan 2009; published 2 Feb 2009ε = ∑ k w k ( a * 1 (θ k ) a 1 (θ k ) − 1) 2 , The factor of |μ 0 /μ 0 | takes the different projection areas of the solar beam at the two sides of the ocean surface into account.The [40] an efficient and accurate model for water leaving radiance simulation is very useful.The paper by Chami et al. has tried to estimate the water leaving radiance by removing the sky light reflected from the ocean surface[40].That procedure is inefficient and less accurate in the sense that it (C) 2009 OSA 16 February 2009 / Vol.17, No. 4 / OPTICS EXPRESS 2067 #105295 -$15.00USD Received 15 Dec 2008; revised 20 Jan 2009; accepted 26 Jan 2009; published 2 Feb 2009

Table 1 .
Average CPU time for the SG and DG schemes.
•.It is obvious that the polarization components Q/I, U/I, and V /I agree well between the two methods.Another very interesting feature is

Table 2 .
Scattering orders needed for convergence.