1. Introduction

The radiative transfer equation (RTE) describes the propagation of radiation in turbid media phenomenologically [1, 2, 3, 4, 5, 6]. Recently, the rederivation of the RTE from the electromagnetic wave theory establishes a new interpretation of the radiative transfer theory [7]. The RTE is a differential-integral equation which generally has no analytical solution for realistic medium configurations and boundary conditions. Therefore, efficient numerical solutions are crucial for the applications of the radiative transfer theory. A non-complete list of these methods includes the discrete-ordinate (DISORT) method[1, 8, 9, 10, 11, 12, 13, 14, 15], the invariant embedding method [1, 16, 17, 18, 19, 20], the adding and doubling method [3, 21, 22], the matrix operator method [23, 24, 25, 26, 27, 28, 29], the spherical harmonics method [30, 31], the multicomponent method [32, 33], the spherical harmonics discrete ordinate method [34, 35], the FN method [36, 37], the Successive Order of Scattering (SOS) method [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48], and the Monte Carlo method [49, 50, 51, 52, 53, 54, 55, 56].

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 multi-components 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 ², 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.

2. 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 I; π = 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 ₁ and i ₂ 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 e 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(τ ₀, μ < 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.

3. 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 ₁ in Eqs. (7). In other words, Eqs. (5) and (6) only contains the diffuse light in the turbid medium. In Eqs.(7), E ₀ = (E ₀,0,0,0)^T is the irradiance vector for the solar radiation; The matrices r and t are the Fresnel reflection and transmission matrices, respectively. The readers are refered to [57] for the exact expressions (see Eqs. (12–18) of [57]). μ ₀ = cos(θ ₀) and φ ₀ are the cosine of the solar zenith angle and the solar azimuthal angle, respectively. μ ₀ < 0 because θ ₀ > 90°. μ′₀ = cos(θ′₀). θ′₀ is the refracted angle for the solar incident angle θ ₀, i.e. sin(π - θ ₀) = n_w sin(π - θ′₀) and n_w is the refractive index of the water relative to air. The factor of |μ ₀/μ′₀| takes the different projection areas of the solar beam at the two sides of the ocean surface into account. The radiance vectors at the boundaries are:

where r ^*(μ,φ,μ′,φ′) is the bidirectional reflectance matrix for the ocean bottom. For a Lambertian surface, r ^* _1,1 = ρ is a constant, and r ^* _i,j = 0 for i ≠ 1 or j ≠ 1; ρ is the single scattering albedo of the reflective surface. μ′ = cos(θ′) and sin(θ) = n_w sin(θ′) in Eqs. (8c) and (8e). The factors of n_w ² in Eqs.(8c) and (8e) are due to the n-square law of the radiance transmission through a dielectric interface [4].

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 et al. have recently used the so-called generalized Legendre functions as the expansion basis, which are real by definitions [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(ϕ - ϕ ₀)] and sin[m(ϕ - ϕ ₀)] 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

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) 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 ² _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.

4. 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,diff (τ,μ) 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 the atmosphere:

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.(14) 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. Thus 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 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.

5. 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, τ ₁ = τ ₀ + δτ,…,τ 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, k = K_a+1,K_a+2, …,K_a+K_o+1 with τ _Ka+1 = τ _o, τ _Ka+2 = τ _Ka+1+δτ, … , τ _Ka+Ko+1 =τ ^* _o.

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 ₁ (τ _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 ₁ (τ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 τ _Ka+Ko+1 to τ _Ka+1 for the ocean and from τ _Ka to τ ₀ 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 _^mn (τ_k,μ_p > 0) according to Eq. (11b) for the ocean. For the interface, only the first terms in Eqs. (13c) and (13e) are added to L ^m_n (τ _k,μ_p) for the purpose of 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 ² _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:

6. Validation and results

The validation of the SOS code is done with the Monte Carlo code by Zhai et al. [65]. This Monte Carlo code uses biased techniques and estimation schemes to increase the efficiency 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 ^-½ _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:

 

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 θ ₀ = 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.

Download Full Size | PDF

The ocean phase function is a Henyey-Greenstein (HG) scattering phase function [81]:

where μ = cos(θ) and g is the asymmetry parameter defined by:

 

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

Download Full Size | PDF

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.

Table 1. Average CPU time for the SG and DG schemes.

View Table | View all tables in this article

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°. 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 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.

 

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.

Download Full Size | PDF

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 smaller than the ocean Legender expansion order of L = 36. The following sensitivity studies of scattering N use the DG scheme with P = 36 and M = 8. The results with the scattering orders of N = 40 are regarded as the true solution.

 

Fig. 4. The effects of the orders of the Gaussian quadrature and Fourier series orders

Download Full Size | PDF

Table 2. Scattering orders needed for convergence.

View Table | View all tables in this article

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.

7. 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.