Radiative transfer code SHARM-3D for radiancesimulations over a non-Lambertian nonhomogeneoussurface: intercomparison study

Results of an extensive validation study of the new radiative transfer code SHARM-3D are described. The code is designed for modeling of unpolarized monochromatic radiative transfer in the visible and near-IR spectra in the laterally uniform atmosphere over an arbitrarily inhomogeneous anisotropic surface. The surface boundary condition is periodic. The algorithm is based on an exact solution derived with the Green’s function method. Several parameterizations were introduced into the algorithm to achieve superior performance. As a result, SHARM-3D is 2–3 orders of magnitude faster than the rigorous code SHDOM. It can model radiances over large surface scenes for a number of incidenceview geometries simultaneously. Extensive comparisons against SHDOM indicate that SHARM-3D has an average accuracy of better than 1%, which along with the high speed of calculations makes it a unique tool for remote-sensing applications in land surface and related atmospheric radiation studies. © 2002 Optical Society of America OCIS codes: 290.4210, 010.1300, 010.1310, 010.1320.


Introduction
In two recent papers 1,2 Lyapustin and Knyazikhin presented a new solution to the three-dimensional ͑3-D͒ problem of radiative transfer over a nonhomogeneous and anisotropic surface. We derived both an exact solution in operator form by using the Green's function method and an accurate parameterized solution designed for fast calculations with arbitrary surface conditions. In this paper, the new code SHARM-3D is validated against code SHDOM 3 (spherical harmonic discrete ordinate method), which is widely recognized in atmospheric radiative transfer communities. The SHDOM is a rigorous numerical code that handles the full 3-D problem with a lateral variability of both surface and atmospheric optical properties. This remarkable code can be used to generate the benchmark radiance and flux data. Created for the purpose of realistic simula-tions in clouds, SHDOM is not optimized for clearskies conditions when the atmosphere can be considered horizontally homogeneous. For these conditions, which are important in land remote sensing, an efficient 3-D solution 2 is implemented in code SHARM-3D. To establish the benchmark reference for the required one-dimensional (1-D) Green's function calculations, the discrete ordinate 4 (DISORT) and successive orders of scattering 5 (6S) codes were also used.

Intercomparison Study: One-Dimensional Case
For a better understanding of 3-D simulations and to establish the reference point, the study was performed in four steps of increasing level of complexity: ͑1͒ standard 1-D solution for the path radiance ͑zero surface reflectance͒, ͑2͒ 1-D solution with a non-Lambertian surface, ͑3͒ and ͑4͒ 3-D solution for several test models of surface for Lambertian and non-Lambertian reflection, respectively.
The surface bidirectional reflectance distribution function (BRDF) was modeled with the Rahman-Pinty-Verstraete 6 (RPV) function, which is recognized as one of the most versatile and accurate models and which is implemented in code SHDOM.

A. Models of the Atmosphere
The new code was tested for a wide range of atmospheric conditions in terms of opacity and scattering properties. The atmosphere was modeled differently in the visible and in the near-IR ͑NIR͒ parts of spectrum. In the NIR range the atmosphere was represented by a uniform aerosol layer with a singlescattering albedo a ϭ 0.95 in clear ͑ a ϭ 0.2͒ and hazy ͑ a ϭ 0.8͒ conditions. The height of the layer was H ϭ 2 km for clear conditions and 5 km for hazy conditions; these two models of purely aerosol atmosphere are referred to here as NIR Clear and NIR Hazy. I used the scattering function corresponding to the continental aerosol model of Elterman 7 at 0.75-m wavelength. This phase function ͑Fig. 1, left͒ has a strong forward-scattering peak ͓͑0°͒ ϭ 155.95͔ and an asymmetry of scattering of g ϭ 1 ͞3 ϭ 0.59 ͓ min ͑115°͒ ϭ 0.248͔, and it suits well the purpose of testing the accuracy of the radiative transfer codes.
In the visible spectrum I used a stratified threelayer atmosphere bounded at H ϭ 0, 2, 10, 30 km. The aerosol distribution within layers was ⌬ a ϭ 0.2, 0.04, 0.01. The aerosol-scattering function, represented by the Elterman model, and the singlescattering albedo ͑ a ϭ 0.95͒ remain constant with altitude, though the corresponding atmospheric parameters change because of the changing weight of molecular scattering. The latter was modeled by use of the 1962 U.S. Standard vertical profile with an integral optical thickness 0.2207 at wavelength 0.45 m. The atmospheric model described is referred to further as a Vis. Clear model. In this study I also used a Vis. Hazy model, which has a tripled aerosol content ͑⌬ a ϭ 0.6, 0.12, 0.03͒.
Testing for a purely molecular atmosphere was performed with a model called Molec. Scatt. with ϭ 1, ϭ 0.3444, and a Rayleigh phase function.

B. Path Radiance
The theory behind the numerical algorithm of the code SHARM-3D is founded on accurate knowledge of 1-D atmospheric path radiance and Green's function.
Therefore it is important to evaluate the accuracy of 1-D code SHARM 8 with respect to other recognized and validated radiative transfer codes. To this end, the path radiance was compared among the codes SHARM, DISORT, SHDOM, and 6S; the last-named code was included because of its particular popularity within the land remote sensing community.
The benchmark level of calculations was established with high-order solutions of codes SHARM ͑P 128 ͒, DISORT ͑N streams ϭ 128͒, and SHDOM ͑N ϭ 128, N ϭ 128͒. For SHDOM, 10 -14 and 24 vertical layers were used to reach convergence in the vertically homogeneous and nonhomogeneous cases, respectively.
The 6S code does not permit the angular resolution of integration ͑order of solution͒ to be changed easily, so this code was used as is. To ensure consistency with other codes, I turned off the depolarization factor in code 6S and used the node wavelengths in calculations to avoid errors of spectral interpolation within 6S.
The results of top-of-the atmosphere ͑TOA͒ path radiance intercomparison are presented in Fig. 2. The data are plotted in the 21 Gaussian quadrature points of the view zenith angle ͑VZA͒ for 2 azimuthal angles, 0°͑solid curves͒ and 180°͑dashed curves͒, and 5 solar zenith angles ͑SZA͒, 0°, 30°, 45°, 60°, and 75°. Figures 2͑a.1͒-2͑a.3͒ show excellent agreement of the codes SHARM and DISORT, to Ϯ0.05% and better, except for most grazing-view angles. An exception occurs in directions of the aureole ͑␥ ϭ 0°͒ and of the backscattering ͑␥ ϭ 180°͒; the latter are represented by spikes at ϭ Ϫ0.5 at SZA ϭ 60°. This discrepancy comes from the single-scattering term D 1 , which is calculated differently in SHARM and in DISORT: Whereas SHARM uses an exact tabulated phase function, subject only to errors of interpolation between the nodes, DISORT uses a reconstructed phase function from the Legendre expansion of order L͓ L ͑␥͔͒, which has insufficient accuracy in these particular directions even for high orders of L. This well-known problem of orthogonal expansion of the discontinuous function ͑see, for example, Refs. 9 and 10͒ is illustrated in the middle and right-hand parts of Fig. 1, which show the error of reconstructed phase function L ͑␥͒ for orders L ϭ 32, 64, 128, 256, 900. The maximal error occurs on the boundaries of the interval 0°-180°. Even at the L ϭ 900 expansion terms used in DISORT calculations, the phase function error in the backscattering direction is still as high as 0.656%.
The downwelling path radiance at the bottom of atmosphere, which is not shown here, agrees with the same high accuracy for codes SHARM and DISORT, except for the aureole direction.
A comparison of SHARM and SHDOM shows a slight underestimation of the path radiance by SHDOM in NIR Clear and NIR Hazy conditions ͓Figs. 2͑c.1͒ and 2͑c.2͔͒ within 0.1-0.5% at Ͻ Ϫ0.2 ͑VZA Ͻ 78°͒. In the vertically nonhomogeneous conditions ͓model Vis. Clear, Fig. 2͑c.3͔͒, models have an excellent nonbiased agreement. An error spike in the backscattering direction in SHDOM data has the same origin as in DISORT, as both codes use the Nakajima-Tanaka 11 truncation approximation with L correction of the solution in the single scattering.
In contrast to the codes discussed, the 6S code has an accuracy that is notably lower, 0.6% in the case of pure molecular scattering ͓Fig. 2͑b.3͔͒ and 1-2% in calculations with a selected aerosol phase function ͓Figs. 2͑b.1͒ and 2͑b.2͔͒. It should be mentioned that all the above codes use the plane-parallel geometry that is valid only for zenith angles of up to 78 -80°for the description of the radiative transfer in the Earth's atmosphere.

C. Anisotropic Surface
Of the many different types of surface studied, I show here results for the grass 12 in the NIR spectral region ͑Fig. 3͒, which is rather representative of the other vegetative land cover types in the NIR. The NIR albedo of grass is close to 0.5, so nearly half of the incident solar energy is reflected back into the atmosphere in each instance of reflection. The same orders of solution as in calculations of path radiance were used for SHARM and SHDOM. A lower-order N streams ϭ 80 was used for DISORT because the current beta version still estimates poorly the high-order Fourier components of the surface bidirectional reflectance. The next version of DISORT ͑to be re- leased͒ will fix this problem and significantly reduce computational time. 13 As in the previous case, the codes SHARM and DISORT show excellent agreement, to 0.01-0.02% at ͉͉ Ն 0.2, except for SZA ϭ 45, and 75°for the NIR Hazy atmospheric model. We can speculate that the DISORT solution did not yet achieve convergence at N streams ϭ 80 at these solar angles. This point of view is supported by the absence of similar error peaks in SHARM-SHDOM diagrams and by excellent agreement ͑to Ϯ0.02%͒ of SHARM radiances with independent calculations by the rigorous 1-D Green's function method ͓Ref. 1, Eqs. ͑18͒-͑21͔͒, which are not shown here.
Code SHDOM agrees consistently with codes SHARM and DISORT, to an accuracy of approximately 0.2-1%, slightly underestimating the TOA radiance.
The 6S code has a typical pattern of error in viewing angle at considerably lower accuracy. Its average error for different surface types and atmospheric conditions is in the range of 4 -6%, although in cases of high surface anisotropy it can be as high as 9 -10%. Table 1 allows one to compare the speed of solution by the codes SHARM, DISORT, and 6S in calculations of path radiance and TOA radiance with the anisotropic surface boundary for the vertically homogeneous atmosphere. Handling the geometry is the largest factor in the computer time difference. For example, code SHARM solves the radiative transfer problem for all the solar-view angles simultaneously, DISORT yields a separate solution for each solar zenith angle, and 6S can deal only with one set of incidence-view angles at a time.
The numbers for the SHDOM code are not presented, as it has to use as many as 15 layers to achieve convergence even in the vertically homogeneous case. However, I found that the per-layer efficiency of SHDOM is one of the highest among the codes considered.

A. Variable Lambertian Surface
Two surface models were used in the 3-D studies: ͑1͒ a model of a circle, and ͑2͒ a bounded cascade model. 14 The first model represents a square or rectangular area of albedo q 1 , with a contrasting cir-  cle of radius R and albedo q 2 centered in the middle of area. This model is convenient for testing the general correctness of a solution and for studying the error pattern. The second model allows the scenes to be rendered with realistic spatial correlation.
I studied scenes of 64 by 64 pixel size with 10 vertical layers at the order of SHDOM solution N ϭ 24, N ϭ 48. At this angular resolution, SHDOM underestimates radiance over a Lambertian surface by ϳ0.35%. For this parameter configuration the 0.5 Gbyte RAM of the PC that is used just meets the SHDOM memory requirements without swamping the hard drive. Figure 4 shows the relative difference between the SHDOM and the SHARM-3D solutions for the circle model with q 1 ϭ 0.40 and q 2 ϭ 0.02 in NIR Clear atmospheric conditions. The high contrast of the scene ensures that the errors of the linearized SHARM-3D solution are maximized. 2 The radius of the circle is 1.5 km. At a resolution of 0.1 km, the total size of the scene was 6.3 km ϫ 6.3 km. Results are shown for SZA ϭ 60°and ϭ Ϫ0.7071 ͑VZA ϭ 45°͒. The relative azimuth was ϭ 0°, which corresponds to the position of an observer at the east when the Sun is at the west. One can see that the mean difference lies between Ϫ0.5% and Ϫ1%, of which Ϫ0.35% is a systematic contribution of SHDOM. The largest discrepancy occurs in 1-2 pixels on the rim of the circle, in the area of the contrast transition. These outliers are explained by the difference in handling the surface grid. While SHARM-3D calculates radiance in the grid nodes, SHDOM bilinearly interpolates the surface albedo between the nodes and calculates radiance in the center of a pixel bounded by these nodes. When SHARM-3D calculations are averaged over four grid points to simulate the radiance in the center of the pixel, the maximal difference decreases to 2.5%.
Simulations show that the maximal error is observed on the x transsect in the middle of the image. Figure 5 shows comparisons along two x transsects, at y ϭ 3.2 km ͑middle of image͒ and at y ϭ 1.8 km ͑edge of the circle͒, that are shown by the horizontal lines in Fig. 4. The SHARM-3D data are shown for two view angles, 45°͑squares͒, and 0°͑circles͒, at SZA ϭ 60°and ϭ 0°for NIR Clear ͑top͒ and NIR Hazy ͑bottom͒ atmospheric models. The respective SHDOM data are shown by the nearest dashed curves. Also, the horizontal dashed lines give the respective 1-D solutions over the bright background and over the dark circle for each view angle. The comparison shows very good general agreement between the two codes, within 4.2% for Hazy and within 0.6% for Clear conditions. In the nadir direction the agreement is consistent within several tenths of a percent.
The difference between two codes drops almost linearly as contrast decreases, so in the natural environment the errors will typically be lower. To evaluate the accuracy of a code in more-realistic conditions I used the bounded cascade model with a resolution of 1 km, a mean albedo of 0.2, and transfer factors of the model of p1 ϭ p2 ϭ 0.2. The albedo image, rendered with characteristics q min ϭ 0.066, q max ϭ 0.51, and q ϭ 0.218, is shown in the left part of Fig. 6. The middle and the right-hand figures  show the relative difference SHDOM Ϫ SHARM-3D for clear and hazy conditions, respectively, for ϭ Ϫ0.7071. At other view angles used, the error is smaller. For example, for clear conditions the error ranges from Ϫ1.54% to 0.28% for ϭ Ϫ0.4472 and is within Ϯ0.9% for ϭ Ϫ1. One can see that in these tests the relative difference does not exceed Ϯ2.3% of the total TOA radiance. As in the previous case, the maximal discrepancies are found on the boundary of the high-contrast areas. After the SHARM-3D calculations are normalized to the center of the SHDOM pixels, the total range of differences drops below 1-1.5%.
The effectiveness of the new code can be judged based on its performance. The total SHARM-3D solution consists of three steps. Obtaining the highorder solution for the 1-D Green's function and path radiance takes only several seconds ͑see Table 1͒. Computing the atmospheric optical transfer func-tion 15 ͑3-D Green's function with isotropic source͒ is the longest part of the solution, which may take as long as 10 -20 min. With atmospheric functions precalculated, it takes only 0.5 s for the 64 ϫ 64 image and ϳ36 s for the 1024 ϫ 1024 image to find radiance for any given surface scenario. For comparison, SHDOM's execution times for an image 64 ϫ 64 pixels are 1512 and 2010 s for the NIR Clear and the NIR Hazy cases, respectively.

B. Variable Anisotropic Surface
In this test the Landsat-7 subimage of Oklahoma in NIR band 4 ͑left-hand image in Fig. 7͒ was used to generate a realistic surface reflectance. The area of 64 ϫ 64 pixels outlined by the white square was used to generate the spatial distribution of parameter of the RPV model from the TOA reflectance. At the time of acquisition ͑4 April 2000͒ the scene's primary composite elements were the plowed fields ͑dark ar- Fig. 6. Inhomogeneous surface albedo rendered with the bounded cascade model ͑left͒, and a relative difference ͑SHDOM Ϫ SHARM͒͞ SHDOM ϫ 100% for NIR Clear and NIR Hazy conditions ͑right͒ for VZA ϭ 45°, SZA ϭ 60°, and ϭ 0°. eas͒ and spring grasses ͑bright areas͒. Based on this, the shape of the BRDF ͑the other two parameters of the RPV model͒ was prescribed to every pixel as follows: the plowed field BRDF was assigned to the dark regions ͑image reflectance, Յ0.2͒, the shape of the grasses was assigned to the bright regions ͑image reflectance, Ͼ0.31͒, and some intermediate shape was assumed for the intermediate pixels.
The nadir TOA radiance for the NIR Clear atmospheric model calculated by SHARM-3D at a resolution of 30 m is shown at the middle of Fig. 7. The image that was generated almost reproduces the originally measured radiance. The top image at the right of Fig. 7 shows the relative difference between SHARM-3D and SHDOM calculations. One can see that the accuracy of the new solution is better than 1.2% for more than 90% of the pixels, with a maximal difference with SHDOM of Ϫ2.2%. If we take into account the systematic underestimation of radiance by SHDOM ͑0.3-0.5%͒ and the difference in gridding assumptions discussed in Subsection 3.A, the maximal difference reduces to 1.5%. For comparison, I also calculated radiance by using the independentpixel approximation ͑IPA͒ with an accuracy of 0.05%. This approximation gives the conventional 1-D solution for each pixel, and, with rare exceptions, it represents the common basis for land surface and atmospheric analysis of the remote sensing data. The relative difference SHDOM Ϫ IPA, shown at the bottom right in Fig. 7, has a range of errors of Ϯ17%.
The purpose of the next test was to study the accuracy of the SHARM-3D solution in the highcontrast conditions that provide more rigorous validation than the relatively benign conditions. The bounded cascade model was used to render the spatial distribution of parameter of the RPV model. The other two RPV parameters were generated randomly for each pixel in the range 0.4 Ͻ k Յ 0.85, Ϫ0.33 Ͻ g Ͻ 0. These boundaries of variation were inferred from the analysis of RPV fit data for a broad collection of field BRDF measurements and model simulations. 16 Care was exercised to keep the surface albedo below 1 by additionally restricting parameter k at high values of .
The albedo for generated scene is shown in Fig. 8 Fig . 8. Image of the surface albedo at SZA ϭ 30°in NIR Clear conditions. The spatial distribution of the surface reflectance was generated by the cascade model ͑parameter ͒, and the BRDF shape was generated randomly within the limits measured over land or predicted by the BRDF models. The scale of the axes is in kilometers. for SZA ϭ 30°at a resolution of 1 km. The albedo varies from ϳ0.1 to 0.75 at a mean value of 0.3. This rather high surface contrast can be encountered, for example, in conditions of a partial snow cover. The range of SHARM Ϫ SHDOM differences is shown in the two leftmost sets of error histograms ͑Fig. 9͒ for two zenith angles ͑0°and 63.4°͒ at SZA ϭ 30°and ϭ 0°. These results represent NIR Clear conditions and were not corrected for the gridding difference. After correction, the total range of SHARM Ϫ SHDOM differences did not exceed 2-3%. In Hazy conditions SHARM-3D calculations provide even better agreement with SHDOM, to Ϯ2-3%, with 85-90% of points within 1% accuracy.
The middle and rightmost parts of Fig. 9 show error histograms for the independent pixel approximations for Clear and Hazy conditions, respectively. One can see that in the high-contrast conditions the new method has a dramatically higher accuracy than the conventional IPA method, even at a resolution of 1 km, which traditionally is believed to belong to the 1-D radiative transfer regime.
The computer time for 64 ϫ 64 scenes was ϳ7 h and 27 m for SHDOM compared with 3-5 min for SHARM-3D.

Conclusions
SHARM-3D is a new algorithm to model the unpolarized monochromatic outgoing radiance in the visible and near-IR spectral ranges over a realistic nonhomogeneous and anisotropic surface. As a byproduct, it calculates the instantaneous surface radiative budget, including surface albedo, and absorbed and reflected fluxes. This code is suitable for cloudfree or for homogeneous cloudy conditions when the atmosphere may be considered horizontally uniform.
The algorithm is a mixture of the exact Green's function solution and several parameterizations introduced to achieve superior efficiency of the code. The main approximations include ͑1͒ linearization of radiance in spatial variation of surface reflectance, ͑2͒ use of the maximum eigenvalue acceleration method along with the Lambertian approximation to model the multiple reflections of light from an anisotropic surface, and ͑3͒ modeling of the diffuse atmospheric transmission of the variation of the surface-reflected radiance in the Lambertian approximation. All the parameterizations have a strong physical basis and are supported by a theoretical or numerical analysis. As a result, SHARM-3D is 2-3 orders of magnitude faster than the rigorous code SHDOM. By comparison, SHDOM is considered to be a fast code relative to statistical Monte Carlo methods, which are usually used to model 3-D radiative transfer over the inhomogeneous surfaces.
In this paper I have presented only several examples from the intercomparison study, which was aimed at a comprehensive characterization of SHARM-3D in different view geometries and atmospheric and surface conditions. All the results obtained show an agreement between SHARM-3D and SHDOM to several tenths of a percent for the over-whelming majority of pixels of the analyzed images. For the pixels located on the boundaries of sharp contrast, the difference may be as high as several percent. The origin of the discrepancy between the two codes for these pixels is in the gridding assumption, as was discussed above. The normalization of the SHARM-3D radiances to the center of a pixel dramatically reduces these boundary effects, bringing the total range of differences to 2-3%.
One of the major advantages of the Green's function technique implemented in SHARM-3D is the separation of the atmospheric radiative transfer from the surface-reflective properties. As a result, atmospheric 1-D and 3-D Green's functions need to be calculated only once for a given optical state of the atmosphere. Then they can be used for calculations with various surface conditions. Conceptually, this method is ideally suited for the rigorous iterative atmospheric correction algorithms based on lookup tables. At present, this general concept is implemented operationally in the Multi-angle Imaging Spectro-Radiometer ͑MISR͒ atmospheric correction algorithm. 17 The code SHARM-3D will be released for public use in the near future.