Monte Carlo Simulations of Radiative Transfer in Cloudy Atmosphere over Sea Surfaces

A rigorous and modular Monte Carlo radiative transfer model MCRT has been developed to compute radiance for realistic cloudy atmosphere over sea surfaces in the visible and near-IR spectra. The simulation methods of the sample’s trajectories are consistently described in the whole simulation spatial domain, which consists of a number of voxel cells for atmosphere and square cells for the sea surface. An inverse transformation method is applied in an atmospheric scattering event simulation and a rigorous accept/reject method is formulated for simulating a reflection event at Cox-Munk (CM) anisotropically and Nakajima-Tanaka (NT) isotropically roughened sea surface. Both methods have been implemented within the model and were tested as two individual modules achieving high accuracy. The model as a whole system was also validated by other codes. The mean difference is about 0.084% and 0.528% in comparison to the Spherical Harmonics code (SHARM) and Monte Carlo Atmospheric Radiative Transfer Simulator (MCARaTS) respectively for 1D atmosphere and the NT model. The comparisons of the Spherical Harmonic Discrete Ordinate Method (SHDOM) show that agreements are obtained in the sun glint regions for 1D atmosphere and the CM model and the mean differences are below 0.478% for a 3D cloud field with lambertian. In general, the major advantage of MCRT is that it could simulate more accurate reflection direction at sea surfaces in a shorter time and hence more accurate radiance in complex cloudy atmosphere over sea surfaces with high numerical efficiency especially on a small PC.


InTRoduCTIon
Because it is quite time-consuming to solve radiative transfer in realistic cloudy atmosphere, homogenous planeparallel clouds are still widely applied in radiative transfer models (Mayer 2009).Nevertheless, there is no reason to assume cloudy atmosphere as a series of horizontally homogeneous layers.In recent years, a variety of three-dimensional (3D) radiative effects have been reported, such as the albedo reduction due to cloud inhomogeneity (Macke et al. 1999), roughening and smoothening of radiative fields (Iwabuchi and Hayasaka 2002;Várnai and Marshak 2003) and enhanced clear sky reflectance near clouds (Kobayashi et al. 2000;Nikolaeva et al. 2005;Wen et al. 2007;Várnai and Marshak 2009).3D radiative transfer model is the only tool which can examine these effects and can be solved by Monte Carlo method, which estimates interest quantities from a random process.Many papers (O'Hirok and Gautier 1998; Barker et al. 2003;Iwabuchi 2006;Mayer 2009;Pincus and Evans 2009;Cornet et al. 2010;Buras and Mayer 2011;Wang et al. 2011) have discussed the Monte Carlo radiative transfer in 3D cloudy atmosphere in detail.However, the surface was usually assumed to be lambertian and simulations with realistic surfaces were rarely reported in detail.In some sense, simulation results may be distorted for realistic surfaces such as sea surfaces.The CM model (Cox and Munk 1954) has been evaluated as the best one of sea surface models and the NT model (Nakajima and Tanaka 1983) is very close to the CM model (Zhang and Wang 2010).In this study, a forward Monte Carlo radiative transfer model MCRT is presented to solve for cloudy atmosphere with both sea surface models.The advantage of the model is that 3D cloudy atmosphere could couple with sea surfaces consistently.It could realistically simulate sun glitter and clouding simultaneously, which are helpful in interpreting 3D radiative effects and correcting biases due to the homogenous plane-parallel assumption.The paper is organized as follows: section 2 describes basic simulation procedures; sections 3 and 4 further discuss simulation methods of interaction with atmospheric particles and sea surface; validations and conclusions are presented in sections 5 and 6.

bASIC SIMulATIon pRoCeduReS
The MCRT has been designed to solve for 3D radiative transfer equation in cloudy atmosphere over sea surfaces.The model simulates trajectories for a large number of samples, illustrated in Fig. 1.These samples, which are emitted from a light source (e.g., the sun), walk randomly with a series of collision events in the atmosphere over sea surfaces until they are either completely absorbed or escape at the top of atmosphere (TOA).The model's spatial domain is assumed to be a cellular structure with a cyclic boundary condition in horizontal space, with N X × N Y × N Z cubic cells for the atmosphere and N X × N Y square cells for the sea surfaces.The optical properties are constants within each cell and include the following: the extinction coefficient σ, the single scattering albedo s and the phase function P(u 0 , u 1 ) for each atmosphere component; the parameters of bidi-rectional reflectance distribution function (BRDF) R(u 0 , u 1 ) and planar albedo α(u 0 ) (Marshak and Davis 2005) for sea surface.Here, u 0 is the original propagation direction and u 1 is the new propagation direction after interaction with the atmospheric component or sea surface.
The flowchart of the MCRT is shown in Fig. 2. We initialize the sample's status at the source.For the sun, each sample with a weight w = 1 is randomly distributed at the TOA with a propagation direction , , sin cos sin sin cos u 0 0 0 0 0 i z i z i = ^h (i.e., the sun illumination direction), where , 0 0 i z ^h are its zenith and azimuth angle.Simulation procedures primarily include two parts: propagation and interaction.During propagation, the sample is tracking voxel by voxel to integrate the path optical thickness ( ) p i x for each voxel until it approaches a free optical thickness f x (see Fig. 1). ( ) x for the voxel with index i is given by where l ( ) p i is the path length the sample has traveled in the voxel; f x is selected randomly from free path distribution (Marshak and Davis 2005) and given by where R is a uniform random number in (0, 1).Hence the interaction position could be determined in this way in the atmosphere.In addition, the sample may interact with a sea surface when it arrived at the sea surface but the integrated path optical thickness was less than f x ; or it may exit the TOA, whereby we should simulate a new sample.The interaction could be simulated as follows: the kind of interaction is determined from a random number R according to the fraction of the total scattering coefficient for each atmospheric component (O'Hirok and Gautier 1998) or the fraction of the total albedo for each sea surface model; the sample's weight, w, is multiplied by the total single scattering albedo s in the atmosphere or total planar albedo α over the sea surfaces; a new propagation direction u is sampled at interaction point according to P(u 0 , u 1 ) or R(u 0 , u 1 ); finally, the local estimate method (Marchuk et al. 1980;Antyufeev 2000;Marshak and Davis 2005) is used to estimate radiance contribution.The detailed simulation procedures are presented in Fig. 2. We could repeat these procedures until it exits the TOA.In order to accelerate the convergence rate of simulation, several variance reduction methods have been used in the MCRT, such as multiple-scaling methods for anisotropic phase function and optically thin regions (Wang et al. 2011).In the following sections, we discuss the methods of direction selection after interaction with atmospheric component and sea surface.

InTeRACTIon In The ATMoSpheRe
The inverse transformation method (Robert and Casella 2004) has been applied to simulate a new propagation direction from the normalized phase function P i ^h for various atmospheric components (i.e., gas molecules, aerosol particles, water droplets, and ice crystals etc).We first integrate the phase function P i ^h to obtain the cumulative probability distribution function (PDF) # and then invert it to obtain scattering polar angle p R where R is a uniform random number in (0, 1).However it only handles a few phase functions analytically, such as polarized Rayleigh phase function PRayleigh i ^h (Vermote et al. 2006) and Henyey-Greenstein phase function PHG i ^h.They are defined as follows where the depolarization factor d = 0.0279, and the asymmetry factor g. The cosines of scattering polar angle are given by where u b a b + , the parameters are summarized as follows: a @.However, the inverse cumulative PDF p R 1 -^h could not be solved analytically for the complex phase functions, such as those computed for water droplets and ice crystals.An alternative is to make a lookup table (LUT) to tabulate i for a set of values in (0, 1) in advance and then interpolate i for any given R.Moreover, the scattering azimuth angle is chosen randomly from (0, 2π).So the new propagation direction can be calculated using scattering polar and azimuth angle according to space analytic geometry (Marchuk et al. 1980).

brief Introduction of Sea bRdf
Sea surface can be regarded as a collection of capillary wave facets, each randomly tilted with respect to the local horizon.Figure 3 shows the coordinate system of light reflection at the sea surface with a tilted facet.The direction, U n , pointing normal to the tilted facet, makes the zenith angle n i and azimuth angle n z ; the direction, U 0 , pointing toward the incident light, makes 0 i and 0 z while the direction, U, pointing toward the reflected light, makes i and z.The wind azimuth angle is w z .The BRDF of sun glitter is given by For CM model without the Gram Charlier series, p(z x , z y ) is expressed as where c v and u v are the crosswind and upwind root mean square slope parameters respectively.Aerial photographs of sun glint were used to evaluate these statistical parameters as a function of wind velocity W (Cox and Munk 1954).The results are given by .
. W 0 003 0 00192 For the NT model, p(z x , z y ) was formulated independently on wind direction (Nakajima and Tanaka 1983) and given by The mean square slope parameter is given by A shadowing factor G(μ 0 , μ) (Nakajima and Tanaka 1983), which is interpreted as a probability that an individual facet is visible, is included in the BRDF of sun glitter for the NT model.

formulation of the Simulation Method
It is difficult to directly sample the reflection direction from the BRDF [see Eq. ( 7)] as we should usually develop the cumulative joint probability table for BRDF in advance (Cornet et al. 2010), which consumes a large amount of memory especially when considering inhomogeneous surface.In this study, an accept/reject method (Robert and Casella 2004) is formulated to select the reflection direction for both the CM and NT models according to the corresponding bidirectional reflectance PDF (BRPDF), The basic idea of the method is to simulate direction from instrumental PDF , , , pins 0 0 n z n z ^h rather than BRPDF, and then to accept or reject the direction randomly according to the accept/reject (AR) function, , , ,  , , ,   , , ,  p p , , , pins 0 0 n z n z ^h should be chosen to match , , , pBR 0 0 n z n ^ zh as closely as possible to have a high acceptance probability.We could find that p(z x , z y ) in the facet slope coordinate is a good choice.In other words, , , , pins 0 0 n z n z ^h in the spherical coordinate is the product of p(z x , z y ) and the Jacobian determinant J of transformation from the slope coordinate (z x , z y ) to the spherical coordinate , n z ^h, and is given by , , , , ; , , where the factor, J, is calculated as By substituting Eqs. ( 7), ( 16), ( 18), and (19) into Eq.( 17), we obtain , , , , cos The numerical efficiency could be evaluated by the maximum value , p max AR 0 0 n z ^h of , , , pAR 0 0 n z n z ^h for the given incident direction , 0 0 n z ^h, since it represents the average sampling times from , , , pins 0 0 n z n z ^h required by sampling one reflection direction from the BRPDF (i.e., the average iterations for one reflection direction simulation, see detailed examples in section 5.3) (Robert and Casella 2004).
On the basis of the above discussion, we could present simulation procedures in detail.First, the CM model is considered.The Box-Muller algorithm is used to sample z x and z y from Eq. ( 11) in the slope coordinate where both R 1 and R 2 are uniform random numbers in (0, 1).By substituting Eq. ( 10 Then rearranging Eqs. ( 23) and ( 24) according to Eqs. ( 8) and ( 9), we obtain U n , U 0 and U are coplanar according to reflection law, thus U is give by The same procedures could be also applied to the NT model if the AR function is slightly modified as , , , cos However, we could make several simplifications for the NT model due to the independence of wind direction, such as u nx and u ny could be reduced as follows In addition, it should be noted that the factor 1 0 n a 6 , 0 0 n z ^h @ could be discounted in Eqs. ( 20) and ( 29), since the factor is constant for a given incident direction.At same time, we could calculate the LUT of the albedo for a set of incident direction in advance and then interpolate it for any given incident direction to update sample's weight w.It is worth noting that LUT for NT is a 1D array while the LUT for the CM is a 2D array as it is not only dependent on the incident zenith angle but also on wind azimuth angle.Although the LUT should be made for both models, it could also save substantial memory because of LUT independence of the reflection direction.Whitecaps as a lambertian are also implemented into the MCRT model.Elsewhere, simulation for the lambertian has been discussed in much detail (Marshak and Davis 2005).

VAlIdATIonS
The validation of MCRT was divided into two parts.First, the selection of a new propagation direction after interaction with atmospheric particles and sea surface was tested as individual modules.Second, the model as a complete system was validated with several 1D and 3D cases.

Selection of Scattering polar Angle in Atmosphere
A sharply peaked scattering phase function was calculated for water droplets with gamma size distribution using the Mie theory, shown in the inset of Fig. 4. We used 10 8 samples to select scattering polar angles i from the phase function and then estimated portion of light scattered at each i, which corresponds to the PDF for the phase function.It, in theory, is the product of the phase function and the sine of the associated i. Figure 4 clearly shows that the selected i perfectly follows the proper PDF.It demonstrates that the LUT is an effective tool for scattering polar angle selections.However, LUTs should be crafted with high accuracy and fine angular resolution.Various tests show that LUTs made with 0.01° resolution by cubic spline interpolation are reasonable for water droplets with an effective radius less than 20 μm at visible wavelength and it is highly accurate to select i from LUTs by linear interpolation.

Selection of Reflection direction at Sea Surface
A reflection direction simulation at sea surface for both the NT and CM models has been tested as another individual module.As an example shown in Fig. 5, we estimated the PDF for both models [i.e., the BRPDF defined in Eq. ( 16)] using 10 8 samples when the incident direction was , 0 0 i z = ^h (20°, 0°) (see the parameters in Fig. 5).Both the simulated and theoretical PDF for both models are almost the same over the upper atmosphere above the sea surface.
The errors are always less than 10 -3 .The results in sun glint regions are superior to the results far from sun glint regions.We could improve the PDF accuracy by introducing more samples.By comparing (a) and (b) in Fig. 5, effects of wind direction on the pattern shape of sun glint have been properly achieved for the CM model.Moreover, we have successfully simulated a large number of PDF for both models with various wind parameters in the visible and near infrared spectral range.

1d and 3d Atmosphere Cases
Radiance simulation has been performed for both 1D and 3D cases in the steps of increasing levels of complexity to validate MCRT model as a complete system.Radiance was normalized as cos I F0 0 r i ^h , where I is radiance, F 0 is incident solar irradiance.
For 1D cases, which were used primarily for checking sea surface models, we only varied the wind velocity while other parameters remained constants.The maritime aerosols were chosen with a visibility of vis = 23 km; the standard US atmosphere was used to consider Rayleigh scattering.The models were tested among three codes: SHARM (Lyapustin 2005), MCARaTS (Iwabuchi 2006) and SHDOM (Evans 1998).Because of only the NT model was implemented in SHARM and MCARaTS and only the CM model was implemented in SHDOM, we used SHARM and MCARaTS to check the NT model and SHDOM to check CM model.
Figure 6 presents the results of TOA radiance intercomparison for the NT model without whitecaps and shows excellent agreement of the MCRT and SHARM for entire viewing zenith angles along the principal plane.The mean   difference is 0.084%; the mean maximum difference is 0.291%.Even for large viewing zenith angle v i , the difference is below 0.398% for all cases.This minor discrepancy arises from the slow rate of convergence at a large v i as few samples are scattered in horizontal direction for the MCRT.However the comparison of the MCRT and MCARaTS shows slight biases by the MCARaTS.The maximum difference is 2.21% and the mean difference is about 0.528%.The reason is largely due to negligence of the Jacobian determinant J [see Eq. ( 19)] during derivation of the AR function (Iwabuchi and Kobayashi 2008).As an example, Fig. 7 shows that the simulated PDF (similar to the PDF in Fig. 5) without the factor J largely deviates from the theoretical standpoint for the NT model with wind velocity 12 m s -1 (i.e., the simulated PDF incorrectly increase at large reflection zenith angle and decrease at small reflection zenith angle in sun glint regions).It clearly explains the positive biases at large v i and negative biases at small v i for the NT model with wind velocity 12 m s -1 in Fig. 6.Even so, the radiance biases for MCARaTS are very small in Fig. 6 as the radiance contribution for the interactions occurring after the first interaction with sea surfaces (i.e., the contribution biased by incorrectly simulating the reflection direction) is too small for the clear 1D atmosphere.However, for cloudy atmosphere, such as cumulus, this contribution could become important due to multiple-interactions and thus radiance could be heavily biased.More importantly, since the denominator on the right side of Eq. ( 19) could be very small especially for the large 0 i , the factor can significantly decrease the value of p max AR (i.e., the average iterations for one reflection direction simulation, shown in Fig. 8) Fig. 7.A simulated PDF with and without Jacobian determinant J and theoretical PDF along the plane with a reflection azimuth angle of z = 180° at sea surface with the incident zenith angle 0 i = (20°, 50°, 80°) and azimuth angle 0 z = 0° for the NT model with wind velocity 12 m s -1 .
Fig. 8.The average iterations for one reflection direction simulation as functions of the incident zenith angle 0 i for the NT model with a wind velocity at 2, 8 and 12 m s -1 with and without Jacobian determinant J. according to Eqs. ( 17) and ( 18), thus bringing about a higher numerical efficiency.Conversely, due to the lack of J, the iterations are so large especially when W = 12 m s -1 and 0 i > 80° that the iteration limitation should be made for MCARaTS to accelerate simulation (Iwabuchi and Kobayashi 2008).Hence, the major advantage of the MCRT over MCARaTS is that it could rapidly simulate more accurate radiances for the NT model.
Figure 9 presents the results of TOA radiance intercomparison for CM model with whitecaps.We have removed Gram-Charlier terms and water-leaving BRDF in SHDOM to maintain the consistency with MCRT.The results show the significant overestimation of radiance by SHDOM for large viewing zenith angle, although they indicate good agreement in sun glint regions, which partially validates the CM model in the MCRT.The discrimination comes from the lack of angular resolution though N μ = 64 discrete zenith angles used for SHDOM.We could improve the accuracy for SHDOM by increasing N μ , but the amount of memory required for SHDOM may be a large constraint for a small PC computer (Pincus and Evans 2009).
For the 3D case, used primarily for checking MCRT in realistic cloudy atmosphere, we took a cloud field from large-eddy simulations, a marine cumulus (see section 4 of Pincus and Evans 2009), having a square domain 4.8 km × 4.8 km with a horizontal resolution of 50 m, vertical depth of 3.0 km with a vertical resolution of 20 m, filled with some scattered clouds.Aerosols were included in the field and surface was assumed to be lambertian with an albedo α = 0.05.Wavelength was 2.13 μm.The SHDOM as a suitable tool to solve cloudy atmosphere problem was used to validate the MCRT.The simulation was performed with the sun in the direction of , 0 0 i z = ^h (40°, 180°).Reflected radiance from the field was computed at seven viewing angles symmetric about the nadir, at v i = (0°, 15°, 30°, 45°) and v z = (0°, 180°).As an example, the pixel radiance field with 4 × 10 8 samples in a direction of , ) is illustrated in Fig. 10.MCRT provides excellent agreement with SHDOM, always less than 0.922% (see Table 1), for domain radiance calculation.The slight difference is chiefly due to the different truncation approximation for the phase function between SHDOM with the delta-M scaling method (Evans 1998) and MCRT with one geometric truncation approximation (Wang et al. 2011).On the whole, SHDOM is much more efficient than MCRT for pixel radiance calculations, although SHDOM consumes huge amounts of memory during a simulation in cloudy atmosphere (Pincus and Evans 2009) especially with sea surface.Thus, the MCRT, in other words, is a better choice to solve for cloudy atmosphere with sea surfaces on a small PC computer.
Up to now, we have vastly presented MCRT's abilities in simulating interaction with sea surfaces and cloudy atmosphere.Finally, the simulation for the cloud field and the CM model with whitecaps was performed at the highest level of complexity.Figure 11   60°), which were viewed at altitudes of 5 km.The radiation interaction between the atmosphere and sea surface was simulated at highly vivid level.Due to no direct illumination from sun at large 0 i , cloud sides pointing to the viewer become darker with an increasing 0 i .Moreover, sun glint and cloud shadow upon the glint are also clearly shown in Fig. 11.Under this combined configuration of atmosphere and sea surface, more 3D effects could be obtained with a larger 0 i .Thus we can use the MCRT to analyze various 3D effects for complex atmosphere and sea surface conditions.

SuMMARy And ConCluSIonS
We have developed a rigorous and modular Monte Carlo radiative transfer model (MCRT) for highly accurate radiance estimate in realistic cloudy atmosphere over the CM anisotropically and NT isotropically roughened sea surface with various wind parameters in the visible and near infrared spectral ranges.The paper briefly describes the model architecture (primarily the sample's propagation and interaction in the cellular spatial domain), and then put emphasis on simulation methods of the sample's interaction, i.e., the inverse transformation method to simulate scattering events in atmosphere and the accept/reject method formulated to simulate reflection events at the sea surface.The both methods as individual components of the MCRT have been validated with various tests.Moreover, the model as a complete system was also tested for various cases from simple to complex.For the 1D standard atmosphere with the NT model, the comparisons against SHARM indicate that the MCRT does not bias the result as the differences are only a few tenths of a percent (largely due to radiation noise); compared with the MCARaTS, it indicates that the accept/reject method, implemented in the MCRT, could simulate a more accurate reflection direction in a shorter time.The comparisons against SHDOM show good agreement in sun glint regions for the 1D standard atmosphere with the CM model and the domain radiance difference always less than 0.922% for one cumulus field with lambertian.The major advantage over SHDOM is that the MCRT is the better choice on a small PC because of the memory constraints for the SHDOM.Finally, in order to show the MCRT's abilities in various aspects, we have simulated a set of scenes, which depict pictures concerning the interaction between the atmosphere and sea surfaces realistically.In conclusion, the MCRT can be a useful tool to solve radiative transfer problems and examine various 3D radiative effects in complex cloudy atmosphere over sea surfaces.

Fig. 1 .
Fig. 1.Schematic of the sample's trajectories used in Monte Carlo simulations of radiative transfer in 3D atmosphere.
7)where p(z x , z y ) is the PDF of the facet slopes as functions of the slope components z x and z y , RF ~h is the Fresnel reflectance for unpolarized light as a function of incident angle ~ (

Fig. 3 .
Fig. 3. Coordinate system of light reflection at a sea surface.

Fig. 4 .
Fig. 4. A simulated (dotted line) and theoretical (solid line) PDF for the scattering phase function P i ^h as functions of a scattering polar angle i .The phase function was calculated for water droplets with an effective radius of r e = 10 μm and effective variance of υ e = 0.1 at wavelength λ = 0.67 μm, shown in the inset.

Fig. 5 .
Fig. 5. Simulated (left) and theoretical (middle) PDF over the upper atmosphere at sea surface with the incident direction , 0 0 i z = ^h (20°, 0°) for (a) NT with wind velocity 6 m s -1 and (b) CM model with a wind velocity 6 m s -1 and wind azimuth angle w z = 0° (right).The corresponding transects through the plane with various reflection azimuth angle z (see the legends), including principal and perpendicular plane.Dotted lines represent simulated values while solid lines represent theoretical values.
Fig. 6.Normalized TOA radiance and the corresponding relative difference I I I code MCRT MCRT h = -^h between the MCRT and the codes (SHARM and MCARaTS), which were computed for the standard US atmosphere with the NT model, as functions of viewing zenith angles v i along the principal plane.Atmospheric visibility vis = 23 km, wind velocity 2, 8 and 12 m s -1 at sea surface, the solar zenith angle 0 i = (0°, 40°, 70°), wavelength λ = 0.682 μm. 10 7 samples used for both MCRT and MCARaTS.The negative (positive) v i corresponds to forward (backward) viewing directions.
Fig. 10.Normalized radiance field reflected from top of marine cumulus computed by the MCRT model.The surface was assumed to be lambertian with an albedo α = 0.05.

Table 1 .
illustrates a set of local radiance fields with several solar zenith angles 0 Normalized domain radiance for both the MCRT and SHDOM and the corresponding relative difference i = (0°, 30°, 45°, v i