Impact of vertical stratification of inherent optical properties on radiative transfer in a plane-parallel turbid medium

The atmosphere is often divided into several homogeneous layers in simulations of radiative transfer in plane-parallel media. This artificial stratification introduces discontinuities in the vertical distribution of the inherent optical properties at boundaries between layers, which result in discontinuous radiances and irradiances at layer interfaces, which lead to errors in the radiative transfer simulations. To investigate the effect of the vertical discontinuity of the atmosphere on radiative transfer simulations, a simple two layer model with only aerosols and molecules and no gas absorption is used. The results show that errors larger than 10% for radiances and several percent for irradiances could be introduced if the atmosphere is not layered properly. ©2009 Optical Society of America OCIS codes: (010.1310) Atmospheric scattering; (010.5620) Radiative transfer References and links 1. K. Stamnes, S. C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete ordinate method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27(12), 2502–2509 (1988). 2. A. Berk, L. S. Bernstein, and D. C. Robertson, “MODTRAN: A moderate resolution model for LOWTRAN 7,” GL-TR-89–0122, Phillips Laboratory, ADA214337 (1989). 3. K. F. Evans, and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transf. 46(5), 413–423 (1991). 4. E. F. Vermote, D. Tanre, J. L. Deuze, M. Herman, and J.-J. Morcette, “Second simulation of the satellite signal in the solar spectrum, 6S: an overview,” IEEE Trans. Geosci. Rem. Sens. 35(3), 675–686 (1997). 5. Q. Min, and M. Duan, “A successive order of scattering model for solving vector radiative transfer in the atmosphere,” J. Quant. Spectrosc. Radiat. Transf. 87(3-4), 243–259 (2004). 6. M. Duan, and Q. Min, “A polarized radiative transfer model based on successive order of scattering method,” Adv. Atmos. Sci. Doi: 10.1007/s00376-009-9049-8 (in print). 7. G. E. Thomas, and K. Stamnes, Radiative transfer in the atmosphere and ocean. (Cambridge, 1999), P160. 8. J. Xie, and X. Xia, “Long-term trend in aerosol optical depth from 1980 to 2001 in north China,” Particuology 6(2), 106–111 (2008). 9. P. M. Teillet, “Rayleigh optical depth comparisons from various sources,” Appl. Opt. 29(13), 1897–1900 (1990). 10. L. Harrison, and J. Michalsky, “Objective algorithms for the retrieval of optical depths from ground-based measurements,” Appl. Opt. 33(22), 5126–5132 (1994). 11. L. Harrison, J. Michalsky, and J. Berndt, “Automated multifilter rotating shadow-band radiometer: an instrument for optical depth and radiation measurements,” Appl. Opt. 33(22), 5118–5125 (1994). 12. Q. Min, and L. Harrison, “Cloud Properties Derived From Surface MFRSR Measurements and Comparison With GOES Results at the ARM SGP Site,” Geophys. Res. Lett. 23(13), 1641 (1996). 13. Q. Min, E. Joseph, and M. Duan, “Retrievals of thin cloud optical depth from a multifilter rotating shadowband radiometer,” J. Geophys. Res. 109(D2), D02201 (2004), doi:10.1029/2003JD003964. 14. M. Alexandrov, A. Marshak, B. Cairns, A. A. Lacis, and B. E. Carlson, “Automated cloud screening algorithm for MFRSR data,” Geophys. Res. Lett. 31(4), L04118 (2004), doi:10.1029/2003GL019105. 15. M. Alexandrov, A. Lacis, B. Carlson, and B. Cairns, “Remote sensing of atmospheric aerosols and trace gases by means of multifilter rotating shadowband radiometer. Part I: Retrieval algorithm,” J. Atmos. Sci. 59(3), 524–543 (2002). 16. H. R. Gordon, and D. J. Castaño, “Coastal zone color scanner atmospheric effects and its application algorithm: multiple scattering effects,” Appl. Opt. 26(11), 2111–2122 (1987). #121457 $15.00 USD Received 14 Dec 2009; revised 19 Feb 2010; accepted 1 Mar 2010; published 4 Mar 2010 (C) 2010 OSA 15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 5629 17. M. Wang, and H. R. Gordon, “Retrieval of the columnar aerosol phase function and single-scatting albedo from sky radiance over the ocean: simulations,” Appl. Opt. 32(24), 4598–4609 (1993). 18. F. Kuik, J. F. De Haan, and J. W. Hovenier, “Benchmark results for single scattering by spheroids,” J. Quant. Spectrosc. Radiat. Transf. 47(6), 477–489 (1992). 19. M. Duan, Q. Min, and J. Li, “A fast radiative transfer model for simulating high-resolution absorption bands,” J. Geophys. Res. 110(D15), D15201 (2005), doi:10.1029/2004JD005590.


Introduction
The radiative transfer (RT) equation is an integro-differential equation.Even in onedimensional plane-parallel cases, it is usually solved by numerical approximation [1][2][3][4][5][6].In most numerical RT models, the atmosphere is often divided into several layers.Each layer is assumed to be homogenous but the inherent optical properties (IOPs) are allowed to vary from layer to layer in order to resolve the vertical variation in the IOPs.This plane-parallel configuration artificially introduces vertical discontinuities of the atmospheric IOPs at layer boundaries that lead to an unphysical discontinuity in the radiation field.In this study, we will investigate the effects of the vertical discontinuities on simulations of radiances and irradiances.

Derivation of the discontinuity and its underlying physics
The radiative transfer equation for a plane-parallel scattering and absorbing medium can be expressed as [1]: Here µ is the cosine of zenith angle (Fig. 1), positive for downward and negative for upward directions, φ is azimuth angle, τ is the optical depth, I is intensity (radiance) of the radiation field, and J is source function given by: 0 / 0 0 0 The first and second terms on the right side of Eq. ( 2) are source functions due to single scattering and multiple scattering, respectively.ω is the single scattering albedo, P is the scattering phase function due to single scattering, πF 0 is the solar irradiance at the top of the atmosphere (TOA), and µ 0 and φ 0 are the cosine of the solar zenith angle and solar azimuth angle, respectively.The optical depth, τ, is given through integration of extinction coefficient β: Formal solutions of Eq. ( 1) for radiances in the downward and upward directions at optical depth τ can be written as: where b is the total optical thickness of the atmosphere, the symbols ↑ and ↓ denote the upward and downward radiance, respectively.In the Eqs.( 4) and ( 5), the values of µ are always positive, that isµ = |µ|>0.Physically, the radiance I ↑ and I ↓ for the horizontal direction (µ = 0) at τ should have exactly the same value if the atmospheric optical parameters vary continuously with optical depth.Therefore, the following equality must be satisfied:  = (6) However, since the whole atmosphere is often divided into many layers and each layer is assumed to be homogeneous, the IOPs, such as single scattering albedo ω and phase function P, differ between neighboring layers.Thus, the IOPs are discontinuities across layer boundaries.In this situation, Eq. ( 6) may not be satisfied.To provide direct insight into the consequence of this artificial discontinuity, a two-layer model (Fig. 1) is used to simulate the radiances and irradiances.The single scattering albedo and scattering phase function are set to be ω x , and P x for layer x, where x = 1, 2 (for details, see section 3).

Single Scattering Approximation
For the single scattering case, the source function J in Eqs. ( 4) and ( 5) are replaced by: Because both layers are homogeneous, insertion of Eq. ( 7) into (4) and ( 5), leads to the following expressions for the single scattering radiance: For µ = 0, the single scattering radiance in the horizontal direction becomes: ( , 0, ) ( , 0, ) ) Therefore, the downward radiance in Eq. ( 8) and the upward radiance in Eq. ( 9) give different results in the horizontal direction at the interface between layers 1 and 2. The difference depends on the difference between the products ω 1 P 1 and ω 2 P 2 .

Multiple Scattering Radiance
For the second and higher scattering orders, the source function J in Eqs. ( 4) and ( 5) is replaced by: where the subscript "n" stands for the n th order of scattering.Inserting Eq. ( 13) into Eq.( 4) and letting µ tend to 0, we find that the first term of Eq. ( 4) tends to zero, and the second term tends to ( , 0, ) n J τ φ ↓ [7].Therefore, the horizontal radiance for the n th order scattering resulting from the expression for the downward radiance (see Fig. 1) can be written as: and, similarly, the horizontal radiance resulting from the expression for the upward radiance becomes: Therefore, if 1 1 2 2

P P
ω ω ≠ , we have: ( , 0, ) ( , 0, ) ) Thus, again different horizontal radiances are obtained at the interface between the two layers if the IOPs are different in the two layers.Equation ( 16) also holds for the total radiance including all orders of scattering if we replace J n by the total source function J in Eqs. ( 14) and (15).Comparing Eq. ( 16) with ( 12), we easily see that Eq. ( 12) is a special case of Eq. ( 16) for the single scattering case for which n = 1.

Case study
To illustrate the effect of the vertical discontinuity of the atmospheric IOPs on radiative transfer simulations in a plane-parallel atmosphere, we use a simple two-layer model with only aerosols and molecules.The coefficients for aerosol extinction and Rayleigh scattering are assumed to decrease exponentially with height.The scale height of atmospheric density is about 8 km for the US standard atmosphere.For rural aerosols, the scale height is about 1.7 km and a value of 2 km is used in our case study.Thus, we adopt the following expressions for aerosol extinction and molecular scattering (for simplicity we ignore molecular absorption in this study): ,0 exp( / 8) where the subscripts a and m, respectively, stand for aerosol and molecule, and "0" means the coefficients at height z = 0km.The column integrated optical depths due to aerosol extinction and molecular (Rayleigh) scattering are τ a = 2β a,0 , τ m = 8β m,0 .
In the case study, we assume a total aerosol optical depth of 0.4, which is the mean value over China [8].The molecular scattering optical depth for the US standard atmosphere is given by [9]: where λ is the wavelength in micrometers (µm).We use a total Rayleigh scattering optical depth of 0.316, which corresponds to λ = 0.415µm, which is one of the wavelengths used in the Multi-Filter Rotating Shadowband Radiometer (MFRSR) designed for retrieval of aerosol and cloud optical depths [10][11][12][13][14][15].It is also close to the 0.413µm, which is adopted in instruments used in ocean color satellite remote sensing [16,17].The aerosol scattering phase function P a (Fig. 2) for randomly-oriented prolate spheroids with aspect ratio a/b = 4.0, size parameter x = 10.079368, and a refractive index m = 1.55-0.01i is used [18], and P a and the aerosol single scattering albedo ω α are assumed to be constant in each layer.As already mentioned, we further assume that there is no gas absorption.If the whole atmosphere is separated into several layers, and each layer is assumed to be homogeneous, then the optical depth ∆τ x , single scattering albedo ω x , and phase function P x for a layer x located between z 1 and z 2 can be written as: , , x a a x a m x m x x where τ a,x and τ m,x , respectively, are the contributions of aerosol and Rayleigh scattering.
Unless otherwise stated, we use the DISORT [1] radiative transfer code, to investigate the effects of the vertical IOP discontinuities on simulated radiances and irradiances.The whole atmosphere is separated into two layers at z = 2 km, as shown in Fig. 1.Thus, the aerosol optical depth of layers 1 and 2 are 0.147 and 0.253 respectively, and the Rayleigh scattering optical depth of layers 1 and 2 are 0.246 and 0.07, respectively.The solar zenith angle is set to be 53.13degrees (µ 0 = 0.6).All the results shown in this paper are normalized by multiplying by 1/F 0 .

Horizontal radiance at layer boundary
The horizontal radiances at the interface (z = 2 km) between the two layers are shown in Fig. 3 as a function of the aerosol single scattering albedo ω a , and the errors of the horizontal radiances given by the 2-layer model are also illustrated."2km+" and "2km-" denote the horizontal radiances (µ = 0), which are calculated from the downward radiance in the upper layer ("2km+"), and from the upward radiance in the lower layer ("2km-") of the 2-layer model.The "true" value is calculated from the improved successive order of scattering method SOSVRT [5,6], the vertical variation of atmospheric optical properties is properly taken into account (detail in section 3.2).The "true" value is further validated against a DISORT calculation with 80 layers.With such a large number of layers, the discontinuity between layers can be ignored.The radiances for azimuth angles of 0, 90, and 180 degrees are illustrated in Fig. 3, which clearly shows that different horizontal radiances are obtained at the interface between the two layers from the downward radiance in the upper layer ("2km+") than from the upward radiance in the lower layer ("2km-").In the upper layer, there is more molecules than aerosols, and the molecules decrease more slowly, while in the lower layer, most of the scattering is due to aerosols.When the atmospheric layer is reconstructed as a vertically homogeneous layer, the vertical distribution of IOPs in the lower layer has a stronger variation than that in the upper layer, which explains why the radiances labeled "2km+" are more accurate than those labeled "2km-".
For radiances in directions other than the horizontal direction, the discontinuity is not easy to see because the radiances are forced to be the same by the layer interface condition: where the subscripts 1 and 2 denote layers 1 and layer 2 respectively.Why are different radiances in the same direction µ = 0 produced by the 2-layer model?The reason is that the atmosphere is broken into two different layers, but each layer is assumed to be homogenous.Thus, this 2-layer model is a very crude representation of the true vertical variation of the atmosphere.As a result, the atmospheric IOPs on one side of the interface between the 2 layers are different from the IOPs on the other side.This artificial discontinuity in atmospheric IOPs results in different horizontal radiances.Figure 3 shows that the differences increase with the increase in aerosol absorption (smaller single scattering albedo).Because we assume there is no molecular absorption in atmosphere, the smaller ω α results in bigger difference of the ω (or ωP) between the two layers and larger errors in the radiances.
Although the plane-parallel assumption is not applicable for direction near the horizon, the errors in the radiances at large zenith angles could introduce extra errors to radiances at small zenith angles, because the source functions due to multiple scattering are derived by integration over all polar angles.Such errors can be a problem for retrieval algorithms which use wavelengths in bands with strong gaseous absorption.For example, in high spectral resolution measurement of the oxygen A-band, the absorption optical depth varies from 0 to 10 or even 100 for ultra high spectral resolution [19], and the absorption coefficient of oxygen varies sharply with height, therefore, the single scattering albedo in a layer could be very small (close to 0), or close 1 due to the contribution of the oxygen absorption.
In real atmospheres, the atmospheric IOPs vary continuously with height.We ignore this continuity when we separate the atmosphere into a limited number of homogenous layers.Different concentrations of aerosols and molecules present in each layer, lead to different values of ω x and P x for the two layers, resulting in the difference in radiances.For example, if the lower layer contains absorbing aerosols (ω a = 0.5), whereas the upper layer contains molecules only, the single scattering albedo just below the interface of the two layers is 0.5, while it is 1.0 just above the interface.

Effects on calculation of radiances and irradiances
Generally, we are not interested in the radiance close to the horizontal direction in planeparallel atmospheres, but it is important to quantify the errors in radiances and irradiance incurred by the artificial discontinuity resulting from dividing the atmosphere into a small number of layers.Figure 4 illustrates the relative errors in radiances due to inadequate resolution of the vertical variation of atmospheric IOPs in the two-layer model.Only radiances for view zenith angles less than 75°, where the plane-parallel assumption is applicable, are plotted.This figure clearly shows that inadequate resolution of the vertical variation in the two-layer model can lead to significant errors in the radiances at small zenith angles, not only for radiance at the interface between layers, but for radiance of all levels including the TOA and the surface.The errors can be up to 10% or even larger, which becomes a serious concern in the development remote sensing algorithms.
Based on the two-layer model, the upward irradiance at the TOA and the downward irradiance at the surface are illustrated in Fig. 5 (left axis), relative errors are also plotted (relative to the right axis).For the irradiances at the TOA, errors of 20% are possible for aerosol with strong absorption and several percent for irradiances at other altitudes.To reduce the errors, more layers are needed in the radiative transfer simulation.To investigate how many layers are needed in the simulation, the whole atmosphere is divided into many layers based on Eq. ( 17) and each layer is assumed to be homogenous and has the same total optical depth.The optical parameters ω, ∆τ and P for each are given by Eqs.(20), (21).The "true" value is calculated from the improved version of SOSVRT by direct integration of the source function over τ with a small step size δτ; in the case study of this paper, we used δτ ≈0.002.For each integration step between τ i and τ i + δτ, the extinction coefficients β a and β m at three points, τ i , τ i + δτ/2 and τ i + δτ are computed through Eqs.(3), ( 17) and (18), and the single scattering albedos ω and phase functions P at the three points are given through Eqs. ( 21) and ( 22) by replacing optical depth ∆τ with extinction coefficient β.Allowing for this vertical variation within each layer in the SOSVRT is the main difference from a traditional n-layer model in which each layer is assumed to be homogeneous.Then the source functions can be calculated with Eq. ( 7) or ( 13) and integrated to compute the radiance by assuming it varies linear-exponentially with τ [5].Therefore, the vertical variation of atmospheric optical properties is properly taken into account.Figure 6 illustrates the maximum error in radiances and irradiances as a function of the number of layers.The results show that the stronger the aerosol absorption, the more layers are needed.For the case discussed in this paper, the number of layers needed to ensure 1% accuracy of radiances for zenith angles less than 75 degrees are 4, 5, 6, 7 for aerosol scattering albedo of 0.95, 0.8, 0.7 and 0.5, respectively, while 2, 4, 5 and 6 layers are needed to ensure 1% accuracy of irradiances.

Conclusion
In numerical simulations of radiative transfer, the atmosphere is often divided into many layers and each layer is assumed to be homogenous by neglecting the vertical variation of atmospheric inherent optical properties (IOPs) within each layer.This assumption introduces an artificial discontinuity of atmospheric IOPs and horizontal radiances at layer interfaces, and results in errors in radiances and irradiances also at small zenith angles and at all levels in the atmosphere.The bigger the difference in IOPs between layers, the larger the errors in radiances and irradiances.A simple two-layer radiative transfer model is used to investigate the impact of this discontinuity in atmospheric IOPs on radiances and irradiances.If the presence of strongly absorbing aerosols in atmosphere, the errors in radiances may be up to 10% for zenith angle less than 75°, and several percent in irradiances.Errors of this magnitude require serious consideration in the development of remote sensing algorithms and in climate modeling.For example, in radiative transfer simulations in absorption bands such as the oxygen A-band, which require high spectral resolution, and in atmospheres with strongly absorbing aerosols, such as particles including black carbon, the vertical variation in the IOPs must be properly taken into account.
The artificial separation of the atmosphere into several homogeneous layers with different IOPs, leads to errors in simulated radiances and irradiances.Although errors illustrated in this paper based on a two-layer model are quite large, they could be reduced or even eliminated by introducing a sufficient number of layers or by assuming a continuous variation of the IOPs.Therefore, radiative transfer codes such as DISORT, SOSVRT, PolRadtran, 6S, MODTRAN etc are expected to be accurate enough if the atmosphere is divided into a sufficient number of layers.

Fig. 3 .
Fig. 3. Horizontal radiances (left three panels) and relative errors (right three panels) at z = 2km computed from downward radiance in the upper layer (2km + ) and upward radiance in the lower layer (2km-) with the two layer model.The true value is derived from the improved SOSVRT using a continuous IOP profile.The aerosol optical depths for layer 1 and 2 are 0.147 and 0.253, respectively, and the Rayleigh scattering optical depths for layer 1 and 2 are 0.246 and 0.070, respectively.The IOPs of each layer is given by Eqs.(20)-(22), and the cosine of the solar zenith angle was µ0 = 0.6.

Fig. 4 .
Fig. 4. Relative differences in upward radiances at TOA and downward radiances at the surface produced by the 2-layer model for three different values of the aerosol single scattering albedo ωa.The optical properties of the two layers are the same as shown in Fig. 3.

Fig. 6 .
Fig. 6. , Maximum error of radiance and irradiances versus the number of layers used in the simulation. )