CrystalTrace: a Monte Carlo raytracing algorithm for radiative transfer in cirrus clouds with oriented ice crystals

We present CrystalTrace, a Monte Carlo raytracing algorithm designed for radiative transfer computations involving ice crystals with both random and preferred orientations. Integrated into the three-dimensional radiative transfer solver MYSTIC, which is part of libRadtran, CrystalTrace enables multiple-scattering simulations of absolute radiances for realistic atmospheric scenes, including clouds with different mixtures of ice crystals and water droplets, aerosols, molecules, and surface properties. It seamlessly extends MYSTIC's macroscopic raytracing down to the microscopic scale of individual ice crystals. By computing single-scattering properties online, CrystalTrace removes the need for precomputed look-up tables, which are also subject to limited angular resolution. The current version of CrystalTrace is implemented for the visible spectral range and supports individual hexagonal prisms with different sizes and aspect ratios. It currently does not account for diffraction and ice absorption. CrystalTrace, integrated with MYSTIC, fills a critical gap by providing a computationally efficient forward radiative transfer simulator for ice clouds containing oriented crystals, thereby enabling retrievals of ice crystal properties from ground-based, airborne, and satellite imaging observations. CrystalTrace is a Monte-Carlo raytracing algorithm for radiative transfer simulations of oriented ice crystals.

1 Introduction

Ice crystal orientation has significant effects on the global radiation budget (Noel and Sassen, 2005) and oriented crystals have been frequently observed by space-borne observations (Chepfer et al., 1999; Noel and Chepfer, 2004, 2010; Zhou et al., 2012; Saito et al., 2017). Ice crystal shape and orientation also determines the fall speed of crystals in ice clouds which strongly affects climate sensitivity in GCM (General Circulation Model) simulations (Sanderson et al., 2008; Mitchell et al., 2008). Better knowledge of ice crystal shape, surface roughness, and orientation are therefore essential for improving estimates of the radiative forcing of cirrus clouds as well as satellite retrievals of cirrus optical and microphysical properties (e.g. Yang et al., 2015; Liou and Yang, 2016). Remote sensing of cirrus cloud properties also depends crucially on an accurate representation of the ice crystal single scattering properties which are key to correctly model the enhanced backscattering from horizontally oriented ice crystal faces. In the solar spectrum these oriented crystal faces act as mirrors reflecting a significantly larger amount of radiation back to space compared to randomly oriented crystals.

Figure 1 HaloCam_JPG image from 1 April 2014 during the ML-CIRRUS campaign. Left: original HaloCam_JPG. Right: same image with an unsharp mask filter applied to better highlight the optical features.

From an airplane window or mountain top, cirrus clouds consisting of oriented crystals often form optical displays, such as a sub-sun, sub-sun parhelia, and a lower parhelic circle. In fact, frequently observed halo displays such as sundogs and upper tangent arcs can only be explained by the presence of oriented ice crystals (e.g. Wegener, 1925; Tricker, 1970; Tape, 1994; Tape and Moilanen, 2006). Continuous HaloCam observations during the ACCEPT campaign revealed that the majority (73 %) of observed halo displays were produced by oriented crystals (sometimes together with randomly oriented crystals) while the minority (27 %) were produced by randomly oriented crystals only (Forster et al., 2017). An example of halo displays produced by oriented ice crystals is shown in Fig. 1: An upper tangent arc and sundogs (22° parhelia) are visible, which are produced by oriented columns and plates, respectively. Above the tangent arc, another arc is visible which had been present for several minutes. This is most likely a rare Parry arc which forms by ice crystal columns with a specific horizontal orientation. This halo type is named after William Edward Parry who first described this type of halo display between 1819 and 1820 (e.g. Pernter and Exner, 1910; Tricker, 1970; Greenler, 1980) and is extremely rare to observe.

Radiative transfer simulations usually assume that scattering particles are randomly oriented within the atmosphere. This assumption greatly simplifies the parameterization of scattering because the phase function then depends only on the scattering angle. For liquid water clouds, random orientation is realistic. For ice crystals with preferred orientation, however, the single-scattering properties depend on the azimuth of the scattering plane as well as the zenith angle of the incident and outgoing ray (Liou and Yang, 2016). For a look-up-table-based approach, this implies storing the optical properties for each incident and scattered angle. Storing orientation-dependent phase matrices for all combinations of incident and scattered angles would increase storage by roughly a factor of 10⁷ (360 × 180 × 180) for 1° resolution compared to randomly oriented crystals. A practical alternative is to compute radiances directly using a geometric raytracing method to calculate the radiance field caused by the refraction and reflection of sunlight by ice crystals. This method was first developed by Wendling et al. (1979) for hexagonal ice crystal columns and plates and later advanced by using more complex particles and including absorption and polarization effects (Pattloch and Tränkle, 1984; Muinonen et al., 1989; Takano and Liou, 1990; Macke, 1994; Takano and Liou, 1995; Hess, 1996; Macke et al., 1996; Yang and Liou, 1998; Prigarin, 2009).

This study presents CrystalTrace, a computationally efficient raytracing model for ice crystals of different shapes with both random and preferred orientation. A key advantage of the geometric raytracing method is its compatibility with Monte Carlo methods which makes it straightforward to integrate into libRadtran's MYSTIC solver to allow simulating absolute radiances under realistic atmospheric conditions. This manuscript is organized as follows: Sect. 2 introduces CrystalTrace, Sect. 2.1 describes its implementation, and Sect. 2.2 its integration into MYSTIC. Section 2.3 presents CrystalTrace's validation, followed by a demonstration to simulate different halo displays in Sect. 2.4 including a real-world example.

2 CrystalTrace – the Monte Carlo Raytracing Algorithm

CrystalTrace is a Monte Carlo raytracing algorithm for simulating light scattering by non-spherical ice crystals in the geometric-optics limit. It supports mixtures of different ice crystal shapes with varying orientation angles, along different symmetry axes. The raytracing technique is computationally fast compared to rigorous electromagnetic methods like T-matrix or Discrete-Dipole Approximation (Draine and Flatau, 1994; Flatau and Draine, 2014) which is an important advantage for the use of radiative transfer simulations. We focus on the visible spectral range, where the wavelength is much smaller than typical ice crystal sizes resulting in size parameters of the order of 100. Raytracing is applicable when the particle size is much larger than the wavelength, i.e. for a size parameter χ≫1 (Liou and Yang, 2016). In this regime, incident light can be approximated as a bundle of localized rays undergoing external and internal reflections and refractions. Furthermore, it has to be assumed that the energy attenuated by the scatterer may be decomposed into equal extinction from diffraction and from Fresnel rays. The geometric optics approximation can be applied to crystals with particle sizes large compared to the wavelength of the incident radiation. For radiative transfer simulations of oriented ice crystals in the visible spectrum, this is well justified since the tendency of crystals to assume preferred orientations increases with size. According to Konoshonkin et al. (2017) the accuracy of the results will be at least 95 % for crystals with size parameters of 120 and larger.

2.1 Implementation of CrystalTrace

The implementation of CrystalTrace will be explained in the following which is building on the methodology described by Macke (1994). The current version of CrystalTrace does not account for the polarization state, diffraction, ice absorption and ice crystal surface roughness which could be incorporated in the future. The algorithm starts with:

Step 1 – Select a crystal shape. The default ice crystal shape implemented in CrystalTrace has a hexagonal (six-fold) symmetry. Hexagonal crystals are typical for the Earth's atmosphere and result from the non-uniform deposition coefficients across the crystal surface (e.g. Magono and Lee, 1966). The majority of halo displays observed on Earth can also be traced back to hexagonal prisms. CrystalTrace provides hexagonal prisms as shown in Fig. 2 with different aspect ratios ranging from columns (Fig. 2a) to plates (Fig. 2b), all of which are composed of two basal faces and six prism faces. The length of the c-axis is referred to by L, while the maximum diameter of the ice crystal top and base faces along the a-axis is denoted by 2a. A dedicated Python tool is provided with the CrystalTrace GitHub repository, which takes any combination of L and a (both in units µm) as input to compute the geometry of a corresponding hexagonal prism. Since diffraction is not implemented yet, note that differences in ice crystal size will have no effect on the result as long as the aspect ratio remains constant.

Figure 2 Definition of ice crystal symmetry axes and dimensions. The length of the c-axis is denoted by L. The maximum diameter of the ice crystal top and base is defined by 2a along the a-axis.

Figure 3 Ice crystal geometry and definition of Euler angles as used for CrystalTrace. The reference coordinate system is represented by the x, y, and z-axes. (a) Original geometry of the ice crystal rotated by the Euler angle α_Euler around its z-axis. The axes of the updated coordinate system are denoted by x^′, y^′, and z^′. The second rotation by β_Euler around the x^′-axis results in the crystal orientation shown in panel (b) with the updated coordinate system denoted by x^′′, y^′′, and z^′′. (c) Third rotation around the z^′′-axis by γ_Euler.

Step 2 – Determine Ice Crystal Orientation. The crystal's orientation is calculated by randomly sampling three Euler angles α_Euler (angle of crystal c-axis), β_Euler and γ_Euler. The Euler rotations around these three angles are illustrated in Fig. 3. Three configurations with different degrees of freedom (DOF) for the choice of Euler angles are usually distinguished:

1. Randomly oriented crystals: for randomly oriented crystals the new orientation is determined by drawing random numbers ϱ between [0, 1) for each of the three Euler angles

   $$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{Euler}}=\mathrm{2}\mathit{\pi }\mathit{\varrho },\text{(2)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{Euler}}=\arccos (\mathrm{1}-\mathrm{2}\mathit{\varrho }),\text{(3)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{Euler}}=\mathrm{2}\mathit{\pi }\mathit{\varrho },\end{array}$$

   which corresponds to a DOF=3. This configuration of ice crystal orientation is responsible for formation of the 22 and 46° halo.

2. Singly oriented ice crystals: Singly oriented ice crystals cause the upper and lower tangent arcs in case of columns and sundogs in case of plates. These ice crystals with preferred orientation of their c-axis are implemented by sampling the Euler angles from a Gaussian distribution with mean value μ_β,Euler in the direction of orientation and a standard deviation σ_β,Euler to represent oscillations of the c-axis around the preferred orientation:

   $$\begin{array}{}\text{(4)}& {\mathit{\beta }}_{\mathrm{Euler}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}{\mathit{\sigma }}_{\mathit{\beta },\mathrm{Euler}}^{\mathrm{2}}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{{\left(\mathit{\varrho }-{\mathit{\mu }}_{\mathit{\beta },\mathrm{Euler}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\beta },\mathrm{Euler}}^{\mathrm{2}}}}\right).\end{array}$$

   The standard deviation of the Gaussian orientation distribution is referred to as orientation parameter. Ice crystal plates have aerodynamically stable orientations with their c-axis vertical (similar to leaves falling from trees) which implies ${\mathit{\mu }}_{\mathit{\beta },\mathrm{Euler}}=\mathrm{0}$°. Ice crystal columns instead tend to orient their c-axis horizontally and thus ${\mathit{\mu }}_{\mathit{\beta },\mathrm{Euler}}=\mathrm{90}$°. This configuration has DOF=2 and is referred to as singly oriented since only one crystal axis is oriented and the other two Euler angles α_Euler and γ_Euler are chosen randomly as in Eq. (1).

3. Parry oriented ice crystals: in rare cases horizontally oriented ice crystal columns can also be oriented in their rotation around z^′′′ (cf. Fig. 3c) which reduces their DOF to 1. This means that both β_Euler and γ_Euler are parameterized by a Gaussian distribution in analogy to Eq. (4). Only α_Euler is sampled from a uniform distribution between 0 and 2π. These oriented ice crystal columns form the rare Parry arc which is visible above the upper tangent arc in Fig. 1.

Step 3 – Intersect Light Ray and Crystal Surface. Each ray is described by the intersection of two plane equations with normal vectors corresponding to the oscillating plane of the electric field vectors. The intersection of these two planes with the plane of the crystal provides the point at which the ray hits the crystal surface. The point of intersection between the direction of the incident ray x₀ and the plane of a crystal face is determined at a randomly chosen point. This procedure is illustrated in Fig. 4: a circle is defined in the plane perpendicular to x₀ and centered around the geometric center of the ice crystal s. The radius is determined by the maximum diameter D of the ice crystal. Within this circle a point is chosen by drawing random numbers for the zenith (θ) and azimuth angle (ϕ). The resulting vector of the incident ray is denoted by a red arrow in Fig. 4 and is calculated by

$$\begin{array}{}\text{(5)}& x=s-x_{\mathrm{0}}\cdot D/\mathrm{2}+(\cos \mathit{\varphi }\mathit{u}+\sin \mathit{\varphi }\mathit{v})\sqrt{\mathit{\theta }}\cdot D/\mathrm{2}.\end{array}$$

Random points are repeatedly selected until an intersection is found between the plane of an ice crystal face and the vector of the incident ray x.

Figure 4 Method to determine the initial point x of the incident ray x₀. The red arrow denotes the incident direction of the ray with starting point x. The intersection point with the ice crystal x is determined by drawing random angles θ and ϕ in a circle centered around the barycenter s with diameter D, the maximum dimension of the ice crystal. The vectors u and v are auxiliary vectors which are perpendicular to x and to each other.

Step 4 – Reflection or Transmission? Once a point of intersection with the crystal surface is found, the ray is either reflected or transmitted. Fresnel's equations quantify the intensity of the reflected and transmitted light. Those equations are derived by decomposing the electric fields into their components parallel E_r∥ and perpendicular to the plane of incidence E_r⟂ and applying the respective continuity conditions to obtain the reflection coefficients of the field components parallel (r_∥) and perpendicular (r_⟂) relative to the plane of incidence:

$$\begin{array}{}\text{(6)}& {\displaystyle }{r}_{⟂}={\displaystyle \frac{E_{\mathrm{r}⟂}}{E_{\mathrm{i}⟂}}}=-{\displaystyle \frac{\sin ({\mathit{\theta }}_{\mathrm{i}}-{\mathit{\theta }}_{\mathrm{t}})}{\sin ({\mathit{\theta }}_{\mathrm{i}}+{\mathit{\theta }}_{\mathrm{t}})}},\text{(7)}& {\displaystyle }{r}_{\parallel }={\displaystyle \frac{E_{\mathrm{r}\parallel }}{E_{\mathrm{i}\parallel }}}={\displaystyle \frac{\tan ({\mathit{\theta }}_{\mathrm{i}}-{\mathit{\theta }}_{\mathrm{t}})}{\tan ({\mathit{\theta }}_{\mathrm{i}}+{\mathit{\theta }}_{\mathrm{t}})}}.\end{array}$$

The reflection coefficient is computed by $\mathcal{R}=(r_{\parallel }^{\mathrm{2}}+r_{⟂}^{\mathrm{2}})/\mathrm{2}$ and the transmission coefficient by $\mathcal{T}=\mathrm{1}-\mathcal{R}$. To determine whether the ray is reflected or transmitted, a random number ϱ is sampled uniformly within [0,1). If ϱ<ℛ, the ray is reflected, otherwise it is transmitted. Total reflection occurs if the angle of refraction is ${\mathit{\theta }}_{\mathrm{t}}>\arcsin (\mathrm{1}/n_{\mathrm{Re}})$. In this case, the probability for reflection is 100 % and the ray is externally reflected. The ray paths for refraction and reflection at a hexagonal ice crystal is illustrated in Fig. 5. The scattering geometry is equivalent to a wedge with opening angle Δ. The direction of the reflected and refracted rays can be calculated using Snell's law

$$\begin{array}{}\text{(8)}& n_{\mathrm{i}}\sin {\mathit{\theta }}_{\mathrm{i}}=n_{\mathrm{t}}\sin {\mathit{\theta }}_{\mathrm{t}}.\end{array}$$

Snell's law describes the relation between the angle of incidence θ_i in a medium with refractive index n_i and the angle of refraction θ_t in a medium with refractive index n_t. The refractive index is used from Warren Wiscombe's REFICE code using ice refractive indices from Warren (1984). With the known incident angle θ_i Eq. (8) can be solved for the transmitted angle θ_t. In case of reflection, the refractive indices are the same and cancel in Eq. (8) which reduces to θ_i=θ_r, i.e. the angle of reflection is equal to the angle of incidence.

Figure 5 Refraction and reflection ray paths illustrating the scattering geometry for Snell's law. Δ is the opening angle of the prism, θ_i is the angle of incidence, θ_r the angle of reflection, and θ_t the transmission angle. θ^′ denotes the angles of the second incident and transmitted ray, and Θ is the scattering angle between the incident ray entering the crystal and the refracted ray exiting the crystal. n denotes the refractive index of the incident medium n_i (here air, n_i≈1) and the transmission medium n_t (here ice, n_t≈1.31).

After the first transmission, ray tracing is continued inside the ice crystal as long as the ray is reflected internally off the crystal surfaces until it gets transmitted again and exits the crystal. Once the ray exits the crystal or in case of initial external reflection (the ray never enters the crystal), the CrystalTrace routine ends and the ray gets counted in the respective (μ, ϕ) bin at the top or bottom of the domain depending on the sign of the vector's z component.

2.2 Integrating CrystalTrace in the Monte Carlo Radiative Transfer Model MYSTIC

MYSTIC solves the radiative transfer equation with the Monte Carlo method by tracing photons through the atmosphere. For radiative simulations of ice clouds, MYSTIC uses a look-up table approach which determines the outgoing angle of the light ray for a given incident angle by randomly sampling the scattering angle from the phase function (cf. Fig. 6). For randomly oriented ice crystals, this is an effective approach since the relative azimuth angle of the outgoing ray can be sampled randomly and is independent of the scattering angle. For crystals with preferred orientation however, the scattering angle depends both on the zenith and azimuth angles relative to the plane of scattering, which causes a challenge on computing memory as well as interpolation techniques for the look-up table approach. With CrystalTrace implemented in MYSTIC, the raytracing simply continues inside the crystal once a scattering event takes place and the direction of the outgoing ray directly results from tracing the ray through the crystal. The scattering properties are effectively computed online during the simulations which allows for radiative transfer simulations without explicit knowledge of the scattering phase function. Another advantage is that no delta-scaling is needed since photons refracted by two parallel crystal faces escape the crystal and are traced further through the atmosphere.

Figure 6 Comparison of the scattering phase functions calculated with YANG (black), CrystalTrace (green) and GO (blue) at wavelength of 550 nm. The GO results here were computed with the bugfixed version of Macke (1994) as described in Forster et al. (2017) (see additional details in the Supplement). (a) Scattering phase functions for a solid ice crystal column with AR=1.4 and maximum dimension D=65 µm and (b) for AR=14.4 with D=10⁴ µm. Both the GO and CrystalTrace phase functions were scaled to allow for a comparison of the results for scattering angles excluding the forward scattering peak. The scaling factors are shown in the legend. The results compare well apart from the forward scattering region for the smaller ice crystal (a) due to diffraction which is not accounted for by CrystalTrace.

2.3 Validation of CrystalTrace

CrystalTrace was validated for both single- and multi-scattering simulations. The single-scattering simulations allow comparison with scattering phase functions from Macke (1994)'s raytracing code (referred to by GO) and the database provided by Yang et al. (2013) (referred to by YANG) for ice crystals of different aspect ratios. Figure 6 compares the phase functions calculated with CrystalTrace (green) and GO (blue) for ice crystals columns with aspect ratio 1.4 (Fig. 6a) and long solid columns with aspect ratio 14.4 (Fig. 6b) at a wavelength of 550 nm. The corresponding phase function of the YANG database (black) is shown as a reference. The results compare well apart from the forward scattering region where CrystalTrace shows smaller values than GO and YANG. This is due to diffraction which is not accounted for by CrystalTrace. The difference is more pronounced for small crystals (Fig. 6a) compared to large crystals (Fig. 6b). Due the coarser scattering angle resolution in the case of the GO phase functions the peaks of the 22 and 46° halos are not as pronounced as for CrystalTrace.

Figure 7 Comparison of simulated radiances for a cirrus cloud of optical thickness 1 with randomly oriented ice crystals using optical properties of the YANG database and CrystalTrace for AR=14.4 with D=10⁴ µm at a wavelength of 550 nm. The simulations were performed with 4×10⁷ photons at a solar zenith angle of θ₀=30° and azimuth angle ϕ₀=270°, without atmospheric background. The top panels show a polar projection of the simulated radiances with MYSTIC using (a) CrystalTrace and (b) the YANG look-up table to solve for the ice crystal single scattering properties. Panel (c) shows the difference between CrystalTrace and YANG and panels (d)–(f) show the corresponding standard deviation. The bottom panels display cross sections through the polar plots along the principal (zenith angle transect along ϕ=ϕ₀) (g) and almucantar plane (azimuth angle transect along θ=θ₀) (h). Here, the MYSTIC radiances computed with CrystalTrace are indicated by the green line and YANG in black.

In addition to single scattering, CrystalTrace was validated for multi-scattering radiance simulations as well. As part of libRadtran's MYSTIC solver, CrystalTrace allows simulations ice clouds with different aspect ratios and preferred or random orientations in a realistic atmosphere with air molecules and aerosol particles. To validate the simulated radiances, the results of MYSTIC with CrystalTrace were compared to MYSTIC with tabulated YANG optical properties, assuming randomly oriented ice crystals with the same size and aspect ratio. Figure 7 shows the radiance distribution for an ice cloud of optical thickness of 1, simulated with MYSTIC using (a) CrystalTrace and (b) YANG to solve for the ice crystal single scattering properties. Figure 7c represents the difference between the two simulations. Figure 7d–f show the standard deviation, representing the Monte Carlo noise within a 1σ confidence interval, which is on average about 9 % relative to the simulated radiances. The mean difference in Fig. 7c is (0.0004 ± 0.0573) W m⁻² nm⁻¹ sr⁻¹, effectively zero within the 1σ confidence range. The bottom panel displays cross sections through the principal (Fig. 7g) and almucantar (Fig. 7h) planes.

Table 1 Radiative flux simulations at top-of-atmosphere (TOA) and bottom-of-atmosphere (BOA) with MYSTIC, comparing CrystalTrace and YANG for a solar zenith angle of θ₀=30° and azimuth angle ϕ₀=270°, columns with AR=14.4, D=10⁴ µm at a wavelength of 550 nm, and optical depths $\in [\mathrm{1},\mathrm{10}]$. The simulations were performed without atmospheric background or ice absorption. The net flux is defined as $E_{\mathrm{net}}=E_{\mathrm{dir},\mathrm{TOA}}-E_{\mathrm{dir},\mathrm{BOA}}-E_{\mathrm{dn},\mathrm{BOA}}-E_{\mathrm{up},\mathrm{TOA}}$. Provided are the average and standard deviation over 10⁷ photons. The uncertainty represents a 1σ-confidence interval.

As in Fig. 6, the largest differences occur near the forward peak due to diffraction. Smaller differences are found at the 22° halo peak, which appears sharper in the CrystalTrace results than in YANG. It is important to note that the tabulated scattering phase functions provided by YANG have a resolution of 1° in the scattering angle range, which is not fine enough to resolve the halo peak and therefore smooths out this scattering feature in the radiance distribution. This highlights an important advantage of computing single-scattering properties online with CrystalTrace: the results are not constrained by the angular resolution of a precomputed look-up table. The lower YANG radiances at the 22° halo peak arise from two factors: (1) the coarse sampling of the scattering angle grid of YANG which limits the maximum resolvable peak height and (2) neglecting diffraction may enhance radiances in the CrystalTrace simulations. In addition to radiances, we compared top-of-atmosphere (TOA) and bottom-of-atmosphere (BOA) fluxes for randomly oriented crystals by integrating radiances over the respective hemispheres, with results listed in Table 1. For all configurations, the net fluxes are effectively zero within the 1σ uncertainty, confirming energy closure.

2.4 Application of CrystalTrace

To demonstrate the capabilities of CrystalTrace implemented in MYSTIC, we showcase simulations of cirrus clouds containing mixtures of crystals with different orientations and aspect ratios under realistic atmospheric and solar illumination conditions. Figure 8 shows the radiance distribution simulated with MYSTIC and CrystalTrace for a cirrus with optical thickness 0.8 in an aerosol-free atmosphere with non-reflective surface. A solar zenith angle of 60° and an azimuth angle of 270° were selected. The cirrus cloud contains 40 % oriented ice crystal columns (c-axis horizontal), 40 % oriented ice crystal plates (c-axis vertical), and 20 % randomly oriented ice crystal columns. For both oriented crystal populations an orientation parameter of ${\mathit{\sigma }}_{\mathit{\beta },\mathrm{Euler}}=\mathrm{1}$° was chosen. The oriented plates are responsible for the sundogs (22° parhelia) visible on each side of the sun, the 120° parhelia, and the circumzenithal arc, which is visible at θ≈10°. All parhelia are located on the parhelic circle which forms at the zenith angle of the sun θ=θ₀ and is produced by reflection at vertical crystal faces of both oriented plates and columns. The oriented hexagonal columns produce the upper and lower tangent arcs visible above and below the sun in Fig. 8 and the 22° halo forms due to the presence of randomly oriented hexagons. This simulation was performed using 10⁷ photons and 360 × 180 (ϕ×θ) bins, resulting in about 154 photons per bin. With CrystalTrace, MYSTIC has to be used in the “forward tracing mode” (from the source to the detector), which has a similar computational time as the backward mode (from the detector to the source) if both hemispheres are simulated. The runtime for 10⁷ photons amounts to about 50 s for τ=1 (140 s for τ=10) using a single core Apple M1 Max processor with 3.2 GHz clock frequency and is the same for CrystalTrace and YANG.

Figure 8 Raytracing simulation with MYSTIC and CrystalTrace showing the radiance distribution in a polar projection for the upper hemisphere using 10⁷ photons. The zenith angle θ ranges from 0 to 90° at the horizon, whereas the azimuth angle ϕ ranges from 0 to 360°. The sun is located at an SZA of 60° and an azimuth angle of ϕ₀=270°. The simulation was performed for a wavelength of 550 nm and an ice cloud of optical thickness of 0.8. The cirrus cloud contains 40 % oriented ice crystal columns (c-axis horizontal), 40 % oriented ice crystal plates (c-axis vertical), and 20 % randomly oriented ice crystal columns. For both oriented crystal populations an orientation parameter of ${\mathit{\sigma }}_{\mathit{\beta },\mathrm{Euler}}=\mathrm{1}$° was chosen.

2.4.1 Oriented Ice Crystal Plates

Due to their preferred orientations, the scattering geometry for ice crystals forming sundogs as well as upper and lower tangent arcs strongly depends on the solar elevation. Figure 9 shows CrystalTrace simulations with oriented ice crystal plates for a range of solar elevations starting close to the horizon at 20 to 60°. The raypath for the sundogs is similar to the 22° halo, but only for a solar elevation of 0°, i.e. at sunrise or sunset, when the plane of incidence is perpendicular to the c-axis. For larger solar elevations the raypath is skewed due to the tilted plane of incidence relative to the crystal c-axis, as shown in Fig. 10. For solar elevations larger than about 60° the light rays cannot follow the path of minimum deviation since the scattering plane is not perpendicular to the c-axis. This is in agreement with the CrystalTrace simulations in Figure 9 at 60° solar elevation when the sundogs start to vanish. This cut-off angle can be computed analytically as shown in Minnaert (1993) and Forster (2017).

Figure 9 Sundogs for increasing solar elevations, calculated with CrystalTrace.

Figure 10 Geometry of a skew ray producing a sundog in an ice crystal plate. The solar elevation angle is denoted with h. The incident and transmitted angles are indicated by θ_i and θ_t, respectively.

The Supplement shows additional CrystalTrace simulations for ice crystals mixtures of randomly and preferred orientations. The transition of the sundogs to a 22° halo by gradually increasing the orientation parameter of the ice crystal plates is shown in Fig. S3 of the Supplement.

2.4.2 Oriented Ice Crystal Columns

Upper and lower tangent arcs are caused by ice oriented crystal columns with their c-axis horizontal. Light is refracted through two crystal side faces, similar as for the 22° halo and the sundogs. These arcs are located tangential to the 22° halo (cf. Fig. 11). The shape of the upper and lower tangent arc depends strongly on the solar elevation, as illustrated in Fig. 11 using raytracing simulations. For solar elevations close to the horizon up to about 30°, the upper and lower tangent arcs are separated. As the solar elevation further increases towards the zenith, the wings of the upper and lower tangent arcs are getting closer until, at solar elevations larger than about 30°, they start to merge and form the circumscribed halo. For high solar elevations the circumscribed halo approaches the 22° halo until the two halos coincide for the sun at the zenith. CrystalTrace simulations for ice crystals mixtures of randomly and preferred orientations are shown in the Supplement.

Figure 11 Upper and lower tangent arc merging into circumscribed halo for increasing solar elevations, calculated with CrystalTrace.

2.4.3 Simulation of a Complex Halo Display during the ML-CIRRUS Field Experiment

On 1 April 2014, during the ML-CIRRUS campaign (Voigt et al., 2017) HaloCam_JPG (Forster et al., 2017) recorded halo displays which formed in a very homogeneous cirrus layer over Munich and lasted for about 4 h. The halo displays showed a 22° halo simultaneously with an upper tangent arc and sundogs. High aerosol concentrations due to transported Saharan dust were reported during this period (Voigt et al., 2017). At 08:25 UTC an additional arc was visible on top of the upper tangent arc in Fig. 1, most likely a Parry arc originating from oriented ice crystal columns which have only 1 degree of freedom (see Sect. 2.1, Step 2c). We use CrystalTrace to reproduce this complex and rare halo display and confirm whether this additional arc stems from so-called Parry-oriented crystals. The results of the radiative transfer simulations is shown in Fig. 12 which are performed for a SZA of 57° at a wavelength of 500 nm. The Saharan dust layer is represented by an aerosol layer using OPAC's desert dust mixture (Hess et al., 1998) with an optical thickness of 0.48 (at 500 nm) as provided by the AERONET daily average. The simulation was performed for a COT of 1 and an ice crystal mixture of 45 % randomly oriented columns which form the 22° halo, 10 % oriented plates which cause the sundogs, 35 % oriented columns which are responsible for the upper and lower tangent arcs, and 10 % Parry oriented columns which form the upper Parry arc. The different ice crystal components were estimated to qualitatively reproduce the different halo displays visible in Fig. 1. As in Fig. 1, the Parry arc below the sun in Fig. 12 is not visible due to the large AOT and the long path through the atmosphere due to the low viewing zenith angles. Compared to Fig. 8, which was simulated for an aerosol-free atmosphere, Fig. 12 shows a strong forward scattering signal in the solar aureole region and the brightness of the halo displays is significantly reduced.

Figure 12 CrystalTrace simulation of 1 April 2014 08:25 UTC at an SZA of 57° and a wavelength of 500 nm. The simulation was performed for 45 % randomly oriented columns which form the 22° halo, 10 % oriented plates which cause the sundogs, 35 % oriented columns which are responsible for the upper and lower tangent arcs, and 10 % Parry-oriented columns which form the upper Parry arc. The simulation was performed for 10⁷ photons. A cirrus optical thickness of 1 was chosen and the aerosol was parameterized by the OPAC desert dust mixture with an AOT of 0.48 at 500 nm, which is the daily mean provided by AERONET.

3 Summary and Conclusions

This study presents CrystalTrace, a Monte Carlo raytracing algorithm that simulates light scattering by non-spherical ice crystals in the visible spectrum, accounting for ice crystal orientation. As stand-alone tool, CrystalTrace can model single scattering properties for random as well as preferred ice crystal orientations with varying degrees of orientation, along different crystal symmetry axes. Integrated into the MYSTIC radiative transfer model, CrystalTrace provides an improved framework for simulating multi-scattering radiative transfer in cirrus clouds, in particular optical phenomena produced by oriented ice crystals, such as sundogs and upper/lower tangent arcs which are frequently observed in the atmosphere. This represents a significant advancement over traditional models that assume randomly oriented crystals and rely on precomputed look-up tables.

The ability to simulate real-time single scattering properties of ice crystals allows solving radiative transfer in an efficient and accurate way, not limited by computer memory or pre-defined angular resolution of look-up tables. CrystalTrace dynamically computes scattering properties during simulations, offering higher resolution and greater flexibility for various crystal orientations.

We validated the single scattering properties simulated by CrystalTrace against precomputed scattering phase functions from the Yang et al. (2013) database (YANG). Moreover, we validated the multi-scattering radiances computed by MYSTIC & CrystalTrace against MYSTIC & YANG assuming randomly oriented crystal. The application of CrystalTrace to HaloCam observations from the ML-CIRRUS campaign demonstrates its effectiveness in reproducing the observed halo displays caused by oriented ice crystals. This real-world validation highlights the practical utility of the model and its potential for improving remote sensing retrievals of ice cloud properties. While CrystalTrace marks a significant advancement of simulating realistic radiances of clouds with oriented ice particles efficiently, future evolutions of the model could benefit from implementing the effects of polarization, diffraction, and ice crystal surface roughness, depending on the application.

CrystalTrace, integrated with the radiative transfer model MYSTIC and libRadtran, has demonstrated its potential for enhancing remote sensing retrievals of ice cloud optical and microphysical properties by accounting for ice crystal orientation. Better representation of ice crystal orientation and scattering properties can lead to more accurate satellite retrievals of cirrus cloud properties and improved estimates of cloud radiative forcing. Future developments in this area will continue to bridge gaps in our understanding of ice cloud microphysics, radiative effects, and their role in the Earth's climate system.