Monte Carlo simulation of halos in the crystal clouds

In this paper, we try to answer the question: how the multiple scattering, the sun elevation, shape and orientation of ice crystals in the cirrus clouds affect a halo pattern. To study the radiation transfer in optically anisotropic clouds, we have developed the software based on Monte Carlo method and ray tracing. In addition to halos, this software enables one to simulate “anti-halos”, which above the cloud layer can be seen by observers. We present the visualization of halos and anti-halos generated by the cirrus clouds for different shapes and orientations of ice crystals.


Introduction
The ice crystal clouds (cirrus, cirrostratus and cirrocumulus) regularly cover about 20-30% of the Earth and play an important part in the radiation balance of the Earth's atmosphere [1, p.328], [2], [3]. There are many scientific problems related to the atmospheric optics and radiative transfer in the ice clouds. Monte Carlo method provides an efficient and flexible technique to simulate the radiation transfer process taking into account the multiple scattering of light [1], [4]. The simulation of the radiative transfer in clouds is usually based on the assumption of the optical isotropy of a scattering medium, but for the ice clouds this assumption often seems to be far from reality. Ice crystals in the cirrus clouds often exhibit various shapes and size, as is shown from different in situ measurements and high-resolution digital images [1], [5], [6], [7]. In order to study the light scattering and absorption characteristics by ice crystal particles from the prospects of microphysics, a hexagonal model is conventionally employed. Its rationality and effectiveness are arguably justified based on the chemical fundamental growth theory of ice crystals [1, p.17], [8] and comparisons of mathematical simulation results and observations from the physical remote sensing experiments [1], [5], [9]. A special attention must be given to the orientation of ice crystals, especially, when the research is related to the atmospheric halos. The existence of oriented ice particles has been unequivocally confirmed, for example, by the LIDAR measurements [2], [9] and observations of numerous atmospheric halos [10], [11]. The ice crystals usually fall with their largest area perpendicular to the direction of movement [11, p.6], [12], [13]. In this study, we use the following terms [11]: thin (short) and thick (long) prisms are called plates and columns, respectively. The plate crystals tend to fall with their basal horizontal faces (plate orientations). Column crystals tend to fall with their horizontal axes (column orientations), but it should be noted that ice crystals are usually randomly (random orientations) oriented in space due to atmospheric instability.
The halo optical phenomenon occurs as a result of the sunlight or moonlight scattering on the ice cloud crystals. There are many halo forms [10], [11]: annular 22° and 46° halos, subheliac arcs, tangent arcs, parhelic circles, subsuns, etc. The variety of shapes and types of halos is associated with anisotropic scattering of light by crystals of various shapes, sizes and orientations. Halo phenomena can already be mathematically explained and can be simulated on a computer, see, e.g. [10], [11], [14], [15], [16], but most simulations do not take into account the multiple scattering of light and thickness of clouds. According to the cloud classification by the International Satellite Cloud Climatology Project (ISCCP), the optical thickness of the high ice clouds can reach 23 (for cirrostratus) and higher values [17]. In our paper, we study the influence of the light multiple scattering in optically thick ice clouds on the halo shape.
The objective of our paper is to understand how the cloud optical thickness, multiple scattering, sun elevation, shape and orientation of ice crystals affect the radiation transfer processes and shapes of halos and anti-halos. In Section 2, we briefly describe a mathematical model of the radiative transfer in an optically anisotropic medium and the ray tracing procedure in ice crystals. Results of Monte Carlo simulation we present in Section 3. Here we consider several optical models of the cirrus clouds and analyze the dependence of integral radiation properties, shapes of halo and anti-halo patterns on different parameters.

A mathematical model and photon trajectory simulation in the media with optical anisotropy with respect to the zenith angle
We assume that the scattering medium of the cloud layer is optically anisotropic, and optical properties depend not on the azimuth direction but only on the zenith angle of the photon direction. In this case, the radiation transfer equation may be written down in the form Here we use the following notation: ( , ) is the radiance in the direction = ( , , ) at the point , 0 ( , ) is the source radiance, ′ = ( ′ , ′ , ′ ), ′ is the cosine of the zenith angle for the direction ′ , ′ and are the azimuth angles (with respect to the horizontal axis OX, for instance) for the directions ′ and , respectively. The cosine of the zenith angle ′ =< ′ , > is defined by the scalar product of ′ and the unit vertical vector . The extinction coefficient and the single scattering albedo depend only on the cosine of the zenith angle. The scattering phase function ( ′ , ′ , , ) is the probability density to scatter in the direction with the cosine of the zenith angle and the azimuth angle if the direction before scattering has cosine of the zenith angle ′ and the azimuth angle ′. The phase function satisfies the normalization condition: for all ′ , ′ . Moreover, it is natural to assume that the phase function depends only on the difference − ′ and it is symmetric with respect to the azimuth change: ( ′ , ′ , , ) = ( ′ , ′ , , − ) . It means that the phase function ( ′ , ′ , , ) actually depends only on the absolute value ∆ = | − ′ |.
In this case, the extinction coefficient and the single scattering albedo do not depend on a direction, while the phase function depends on the cosine μ =< ′ , > of the angle between the directions before and after scattering, The relations between the phase functions ( ) and ( ′ , ′ , , ) are the following: The Jacobian of the transformation ( , ) ↔ ( , ) is equal to unit [18]. There are some obvious specific features of the radiation simulation in the described optically anisotropic medium, because the single scattering albedo and distribution of the photon free path length depend on the cosine of the zenith angle of the photon direction. The procedure of the scattering simulation is specific as well. Let us denote the directions of the photon before and after scattering by Here we assume that | ′ | ≠ 1 (if | ′ | = 1 , then (a, b) is an isotropic random vector). For the scattering simulation, we need to know a family of distributions ( ΄,·), ( ΄, ,·) for a set of the parameters ΄[0,1], [−1,1]. If 0 , 0 are probability densities of the distributions and , then we have the following representation of the phase function in terms of the marginal and conditional densities: In the numerical experiments carried out, we assumed that the scattering medium in a cloud is spatially homogeneous and the single scattering albedo is equal to 1 (absorption in the visible wavelength range may be neglected). The distributions and were computed by the ray tracing technique and approximated by piecewise constant probability densities.

The ray tracing in ice crystals
When the ice crystal size in the cirrus clouds (from a few to thousands of micrometers [1]) is much larger than the incident wavelength (a typical human eye will respond to wavelengths from about 0.380 to about 0.750 micrometers), the geometric optics method, or the ray tracing method, is valid for the computation of the light scattering in crystals. In addition, the geometric optics method provides a  [1].
In the process of the ray tracing, the sunlight is treated as a bunch of individual parallel straight lines. Mathematically, the equation of a straight line may be determined by a point and a direction vector. If a photon hits a crystal, it may be reflected away from the crystal or refracted into the crystal. When a photon is refracted into the crystal, it is traced until it escapes from the crystal. The direction of the reflected and refracted ray and its intensity can be calculated by means of the Snell law and the Fresnel formulas. After finishing the ray tracing process, one can count the number of photons or add up the energy (or weights), which could be attached to each photon. More specifically, the ray tracing algorithm for a hexagonal column (with the height ℎ and the side length of a basal face) can be described as follows: where , , are the rotation angles. All the vertex points are rotated around the axes of the global coordinate system. To simulate the sun elevation, initial points and direction are rotated about the axis Y.
Step 3: Determine whether a random incident ray, defined by and 0 , intersects with the crystal. In this study we use the approach proposed in [19], which is suitable for all three-dimensional convex crystals no matter how complex a crystal may be.
Step 4: If there is no intersection, then go to step 1. Otherwise, the coordinates of the intersection point are calculated according to the equations of straight lines (ray) and planes (faces of a crystal).
Step 5: The new direction of the photon 1 = ( 1 , 1 , 1 ) and reflection coefficient  Step 6: Go to step 1, when a photon escapes from the particle. The loop stops until the number of sampled photons reach a pre-given value.

Models of the ice clouds
In the numerical experiments conducted, we simulated the radiative transfer process for a planeparallel source of visible light and homogeneous layers for several optical models of the ice crystal clouds. The single scattering albedo for the models of the ice clouds is equal to 1 for a visible wavelength. We consider several optical models of the ice clouds assuming that the ice crystals have the shape of hexagonal columns and plates. The shapes and orientations of the ice crystals that we use in our computations are presented in table 1. The ratio of the column length to diameter (of a circumscribed circle for the basal face structure) is equal to 2.5. For the hexagonal plates, this value is equal to 0.1. In addition, we consider the ice crystals with the trigonal geometry. The existence of this type of crystals has been confirmed by the ice chemistry theory, laboratory experiments and observations [8], [20]. But until now the trigonal crystals have rarely been used for optical models of the cirrus clouds. In Model 1, the axis of the hexagonal column is isotropically oriented in the horizontal plane. In Model 2, similar to Model 1, all crystals are isotropically oriented in the horizontal plane, and two lateral faces are parallel to the horizontal plane. It is called the Parry orientation. In Model 3, the axis of the hexagonal column is vertical. This orientation is rarely documented, even its existence in the atmosphere remains open [11]. Orientations of columns in trigonal Models 4, 5, 6 are the same as in hexagonal Models 1, 2 and 3. In Model 7 (Lowitz plates), ice crystals are hexagonal plates. They are uniformly rotating around one of its diameters, and diameter is rotating in the horizontal plane. Dependences of the extinction coefficient on the cosine of the zenith angle for Models 1-6 are shown in figure 2. The extinction coefficients here are divided by the corresponding extinction coefficients for the optically isotropic medium with the stochastic orientation of crystals.  Figure 2. Extinction coefficients as functions of the cosine of the zenith angle for Models 1 (left-red), 2 (left-green), 3(left-blue), 4(right-red), 5 (right-green) and 6 (right-blue). Diagrams present the ratio to the extinction coefficients of the corresponding models with isotropic orientation of crystals. Figure 3 represents Monte Carlo estimates of albedo of the cloud layers with optical thicknesses 1 and 5 for Models 1, 2, 3 with different orientations of ice columns and for the model with columns isotropically oriented in space. Numerical results show that the orientation of ice crystals may significantly affect the integral radiation characteristics of clouds. For thin ice clouds, the single scattering of light gives the main contribution to halos. For comparatively thick clouds, the contribution of light of higher scattering orders becomes greater, and it can introduce the new elements to halos (this effect was mentioned in [21], but certain numerical results presented in [21] were inaccurate). Figure 4 demonstrates how the multiple scattered radiation introduces additional elements to the halos generated by cloud Models 1 and 2 with optical thickness 5. The circles are the projections of the sky hemisphere onto the horizontal plane and represent the brightness of the corresponding regions of the sky. The zenith is at the center of the circle, while the points of the circle periphery correspond to the horizon. The distance from the center is proportional to the zenith angle. The colors of the images here are similar to colors of geographic maps: the minimum brightness level is shown in dark blue and a maximum is shown in brown. As one can see from the figures, additional rings and arcs appear in halos when the radiation of second and higher orders of scattering is taken into account.  With an increase in the optical thickness of the cloud layer, the halo image becomes more blurred and indistinct. Angular distributions of the radiation going through the ice clouds of different optical thickness are presented in figure 5 for Models 1, 2 and 7. The source zenith angle here is equal to 60 0 , and the values of the cloud optical thickness are 1, 5, 7 and 10.

Conclusion
We have developed software to study the halo formation and radiative transfer processes in the optically anisotropic ice clouds taking into account the multiple scattering of light. The software is based on the ray tracing technique and Monte Carlo method. Numerical experiments performed for the models of cloud layers with different optical thickness, shape and orientation of ice crystals, have revealed several peculiarities and eventual effects. Computations for the ice clouds with crystals of the same shape and different orientation have shown that orientation can essentially affect the albedo and transmittance of clouds. We have confirmed an observation made in [21] that the light of second and higher orders of scattering can introduce additional elements in a halo pattern. It means that for some types of clouds, the geometry of the halo formation depends on cloud thickness. For optically thick clouds, the angular distributions of transmitted radiation lose information about the source elevation. These distributions become more symmetric with respect to the nadir direction, but at the same time they depend on the shape and orientation of crystals. The analysis of the simulated anti-halo images shows that they are less dependent on the order of scattering and cloud thickness. It should be noted that an unusual behavior of albedo for the vertical column orientation in figure 3 means that albedo can be a non-monotonic function of the source zenith angle for optically anisotropic cloud layers.