Low and consistent asymmetry parameters in Arctic and mid-latitude cirrus

Cirrus clouds play a critical role in the Earth's radiation budget, yet their shortwave optical properties remain poorly constrained. In particular, the short wave asymmetry parameter (g), which governs the angular distribution of scattered light, is particularly sensitive to ice crystal morphology, a property that varies widely in cirrus. To provide observational constraints on g and to investigate its relationship with ice microphysical properties, we analysed simultaneous in situ measurements of particle morphology and angular light scattering using the Particle Habit Imaging and Polar Scattering (PHIPS) probe. The measurements were obtained during the Cirrus in High Latitudes (CIRRUS-HL) campaign in June and July 2021, sampling both mid-latitude and Arctic cirrus across a range of cloud types and temperatures down to −63 °C. Across both regions, we found consistently low median asymmetry parameters, with a campaign-wide median of 0.738. The observed g values were largely insensitive to variations in temperature, humidity, and crystal aspect ratio, and showed only minor differences between ice habits. In contrast, a systematic decrease in g with increasing particle size was identified, median g ranging from 0.777 for sub-30 µm particles in mid-latitude cirrus to minimum values of 0.719 and 0.713 for 175 µm particles in mid-latitude and Arctic cirrus, respectively. The measured values are substantially lower than those commonly used in current radiative transfer schemes, suggesting that the shortwave warming effect of cirrus clouds may be overestimated in many climate models. These results provide improved observational constraints for the representation of ice cloud optical properties and support efforts to reduce uncertainties in cirrus cloud radiative forcing.

1 Introduction

To improve the representation of cirrus cloud radiative effects in climate models and remote sensing algorithms, it is important to have a good understanding of the ice crystal single-scattering properties (Mishchenko et al., 1996; Zhang et al., 1999; Yang et al., 2001). Within commonly used two-stream radiative transfer frameworks, the key single-scattering properties are the extinction coefficient, single-scattering albedo and the asymmetry parameter (g). Among these, g exhibits a particularly complex dependence on ice crystal morphology in the solar spectrum. Ice crystal morphological properties, like crystal size, habit as well as degree and type of crystal complexity, are known to vary greatly with cloud type, in-cloud temperature, and geographical region (e.g. Heymsfield et al., 2017; Wolf et al., 2018; Woods et al., 2018; Lawson et al., 2019).

As a result of this high morphological variability, constraining the single-scattering properties of ice crystals in numerical simulations remains a major challenge. Numerical ray-tracing studies for the solar spectral range predict a wide variety of g values, typically spanning 0.74–0.94, with higher values associated with pristine hexagonal habits and lower values with complex particle shapes or surface roughness (e.g. Iaquinta et al., 1995; Macke et al., 1996b, 1998; Um and McFarquhar, 2007; Yang et al., 2008). Accurate quantification of g is therefore critical, as biases can lead to a significant uncertainty in estimations of the ice cloud radiative effect in climate models (Järvinen et al., 2018; Yi et al., 2013) and in retrieved cirrus microphysical properties from remote sensing observations (Yang et al., 2008; Yi et al., 2017a). Vogelmann and Ackerman (1995) showed that g must be constrained within approximately 2 %–5 % to limit errors in computed shortwave fluxes to about 5 %.

Observational studies are therefore essential for constraining the magnitude of g in radiative transfer simulations. Unlike numerical estimates, measurements indicate a more moderate variability in g. Early radiometric observations reported a low asymmetry parameter around 0.7 in a case study over the midwestern USA (Stephens et al., 1990), while subsequent satellite-based multi-angle retrievals consistently yielded values for the visible asymmetry parameter below 0.8 (van Diedenhoven et al., 2013, 2014b; Wang et al., 2014; Cole et al., 2014; van Diedenhoven, 2021). Similarly, lidar observations provide another means of retrieving the asymmetry parameter, with recent studies also reporting values around 0.76 (Gil-Díaz et al., 2025).

Estimating g from remote sensing observations is inherently limited by the underlying optical model. The retrieved values are constrained by the model assumptions and cannot fall below the minimum g permitted by the chosen model. Consequently, retrievals often converge near this lower bound (Järvinen et al., 2023b). For instance, studies that use a severely roughened hexagonal aggregate to represent ice crystals consistently report visible-wavelength g values around 0.75 (Cole et al., 2014; Järvinen et al., 2018; Wang et al., 2018; Forster and Mayer, 2022).

A more direct approach for determining g is provided by polar nephelometer measurements. To date, measurements from three airborne polar nephelometers have been reported in the literature: the Polar Nephelometer (PN; Gayet et al., 1998), the Cloud Integrating Nephelometer (CIN; Gerber et al., 2000), and the Particle Habit Imaging and Polar Scattering (PHIPS) probe (Schnaiter et al., 2018). These measurements have yielded median or mean asymmetry parameters typically ranging between 0.73–0.79 in cirrus clouds (Gerber et al., 2000; Auriol et al., 2001; Garrett et al., 2001, 2003; Gayet et al., 2004, 2006, 2012; Febvre et al., 2009; van Diedenhoven et al., 2013).

However, such direct measurements remain sparse and have been largely confined to a limited number of airborne campaigns, predominantly sampling tropical deep convective systems or mid-latitude cirrus. To date, only a single in situ observation has been reported for Arctic cirrus, based on a case study of a cirrus spissatus band with a cloud-top temperature of −46 °C (Garrett et al., 2001). Moreover, previous studies have primarily emphasized bulk optical properties, providing only limited insight into the relationship between g and the concurrently measured microphysical characteristics of the ice crystal population.

The objective of this study is to extend direct in situ observations of g by analyzing a comprehensive, multi-flight dataset obtained in both mid-latitude and Arctic cirrus during the CIRRUS in High Latitudes (CIRRUS-HL) campaign using the PHIPS instrument. This dataset enables a unique dual-regional comparison within the same season and includes single-particle measurements that allow for a direct linkage between ice crystal microphysical properties and their angular light-scattering characteristics. In addition to constraining the magnitude of g, we examine its dependence on key microphysical parameters, including effective radius, ice water content, habit, and aspect ratio. By combining bulk and particle-resolved analyses, the results provide new observational constraints for developing physically consistent optical parameterisations for radiative transfer applications.

2 Methodology

2.1 The CIRRUS-HL airborne campaign

The CIRRUS-HL airborne campaign took place from 6 June to 28 July 2021 using the High Altitude and Long Range (HALO) research aircraft (Jurkat-Witschas et al., 2025). The original aim of the campaign was to investigate the microphysical and optical properties of Arctic cirrus using in situ and remote sensing instrumentation, with operations initially planned from Kiruna, Sweden. However, due to pandemic-related restrictions, the campaign base was relocated to Oberpfaffenhofen, Germany. This logistical change broadened the scientific scope to include both mid-latitude and Arctic cirrus (see Fig. 1), as Arctic flight operations now required long transit legs and refueling stops in Scandinavia or Iceland. The revised campaign strategy focused on contrasting cirrus clouds in different regions and meteorological conditions. Accordingly, flights were designed to sample natural mid-latitude cirrus in diverse weather regimes, including warm conveyor belt systems, high-pressure in situ cirrus, and convective outflow, as well as cirrus potentially influenced by aviation emissions. The time frame of the campaign coincided with pandemic air traffic restrictions, resulting in a 65 %–80 % reduction in daily commercial flight activity compared to the same time frame in 2019 (De La Torre Castro et al., 2023). The campaign period also overlapped with extreme rainfall events in central Europe, resulting in fatal flooding in Germany and Belgium (Lehmkuhl et al., 2022).

Figure 1 Trajectories of CIRRUS-HL measurement flights (F02–F23) included to the analysis. Research flights targeting Arctic cirrus were performed as double flights with a refuelling stop. These were F05/F06, F08/F09, and F10/F11. In total, approximately 5625 km of mid-latitude cirrus and 1060 km of Arctic cirrus were sampled.

2.2 Microphysical and optical measurements

During CIRRUS-HL in situ measurements of the ice crystal microphysical and optical properties were conducted with the PHIPS airborne cloud instrument (Abdelmonem et al., 2016; Schnaiter et al., 2018). PHIPS is a combination of a single-particle stereo-microscopic imager and a polar nephelometer, which allow for a combined analysis of optical and microphysical properties at the level of individual ice crystals.

2.2.1 Derivation of microphysical properties from stereo-microscopic images

The PHIPS stereo-microscopic imager consists of two CCD camera and microscope assemblies with an angular viewing distance of 120° each acquiring 1360 by 1024 pixel bright field microscopic images of individual cloud particles at a maximum acquisition rate of 10 Hz. During CIRRUS-HL a magnification setting of 6× was used for both camera assemblies corresponding to an optical resolution of about 6.2 µm and a field of view of approximately 1.1 mm×1.5 mm, limiting the maximum observable particle size to about 1.5 mm. For each image pair a set of microphysical parameters were retrieved using an algorithm discussed in Schön et al. (2011). These microphysical parameters include the maximum dimension (D_max), area-equivalent diameter (i.e. equivalent diameter of a sphere having the same cross-sectional area, D_eq), projected area (A_proj), and the aspect ratio (AR). In this study the AR is defined as the ratio between the maximum dimension that is perpendicular with respect to D_max and D_max, and thus is always equal or below 1.

The stereo-microscopic images were also manually visually inspected and assigned to a habit class. A detailed description of the habit classification scheme can be found in the Supplement of Järvinen et al. (2023a). While the original habit classification scheme includes 18 categories, for the purposes of this study we reduced the habit classes to 9 cirrus-relevant classes. These include: plate (encompassing also sectored and skeleton plates), column, frozen droplet, side plane, bullet rosette, mixed rosette, capped column, and other category, which contains crystals that could not be confidently assigned to a specific habit. Additionally, we introduced an aggregate class to account for particles composed of multiple crystals, even though it does not represent a single-crystal habit in the traditional sense.

To account for measurement artefacts, a separate shattered class was introduced for particles identified as fragments from shattering events. Particles were classified as shattered when multiple small fragments (typically <100 µm) were observed within a single image frame, which has a size of approximately 1.1 mm×1.5 mm. The likelihood of multiple unconnected ice particles naturally occurring within such a small field of view is considered low for cirrus conditions, where particle concentrations and clustering tendencies are limited. Additionally, larger crystals were classified as shattered if they exhibited broken edges, truncated branches, or jagged outlines inconsistent with physically intact ice crystal habits. These particles were excluded from further analysis (see Sect. 2.2.2 and 2.2.4). Representative stereo-images for each habit class and the shattered class are provided in the Supplement.

2.2.2 Calculation of particle size distribution and bulk microphysical properties

Although the stereo-microscopic imager is limited to a maximum acquisition rate of 10 Hz, the PHIPS trigger system detects particles at rates up to 13 kHz. The sensitive area of the instrument is defined by the overlap of the trigger system field of view and the scattering laser intensity distribution and varies between approximately 0.3–0.9 mm², depending on particle size (Waitz et al., 2021). This size-dependent sensitive area, together with the aircraft's true air speed, is used to calculate the sampling volume. Owing to the high velocity of a jet aircraft, compression of air occurs near probe and wing edges. To obtain the true ambient sampling volume, a thermodynamic correction described by Weigel et al. (2016) is applied.

Particle sizing is primarily based on the stereo-microscopic images. Images of the same particle acquired from two viewing angles are used to derive the mean D_eq and A_proj for each particle. For trigger events that did not produce stereo images due to the limited imaging rate, D_eq and A_proj were estimated using the relationship between the particle geometrical size and partial scattering cross section measured by the nephelometer as discussed in Waitz et al. (2021). The corresponding calibration curves for the CIRRUS-HL campaign are shown in Figs. S1 and S2 in the Supplement. During most CIRRUS-HL flights, typically 70 %–80 % of the particle triggers also produced stereo-images, so particle sizing was predominantly image based. In four flights, however, fewer than 50 % of triggers produced images, and sizing relied mainly on nephelometer-derived estimates.

Particle size distributions were derived by grouping particles into logarithmically spaced bins between 15–700 µm in D_eq. Number concentrations were calculated by dividing the particle counts in each bin by the corresponding sampling volume. Using the particle-resolved projected area and size information, bulk microphysical properties including ice water content (IWC), extinction coefficient (β_ext), effective radius (r_eff), and mass-equivalent spherical radius ($r_{\overline{m}}$), were calculated.

The extinction coefficient was determined from the projected area concentration in each size bin. For each size bin i, the projected area concentration dA_proj,i was determined as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathrm{d}{A}_{\text{proj},i}={\displaystyle \frac{{\sum }_k{A}_{\text{proj},i,k}}{V_i}},\end{array}$$

where $A_{\text{proj},i,k}$ is the projected area of the kth particle in bin i, and V_i is the corresponding sampling volume for that bin. Summation is performed over all detected particles k within the bin. The total extinction coefficient was then obtained as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathit{\beta }}_{\text{ext}}=\mathrm{2}\cdot \sum _i\mathrm{d}{A}_{\text{proj},i},\end{array}$$

where the factor of 2 assumes randomly oriented particles in the geometrical optics regime (Bohren and Huffman, 1983).

The IWC was calculated analogously by summing the masses of individual particles within each size bin. Particle mass was estimated from the projected area using the empirical power-law relationship derived for mid-latitude ice crystals of Baker and Lawson (2006):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{m}_{i,k}=\mathrm{0.115}\cdot A_{\text{proj},i,k}^{\mathrm{1.218}}.\end{array}$$

The mass concentration in each bin was then given by

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathrm{d}{m}_i={\displaystyle \frac{{\sum }_k{m}_{i,k}}{V_i}},\end{array}$$

and the total ice water content by

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\text{IWC}=\sum _i\mathrm{d}{m}_i.\end{array}$$

The effective radius r_eff was derived as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{r}_{\text{eff}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\cdot {\displaystyle \frac{\text{IWC}}{{\mathit{\rho }}_{\text{ice}}\cdot {\mathit{\beta }}_{\text{ext}}}},\end{array}$$

where ρ_ice denotes the bulk density of ice.

The mass-equivalent spherical radius was calculated according to

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{r}_{\overline{m}}={\left({\displaystyle \frac{\mathrm{3}\text{IWC}}{\mathrm{4}\mathit{\pi }N{\mathit{\rho }}_{\text{ice}}}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where N is the total number concentration. The quantity $r_{\overline{m}}$ represents the radius of a sphere with the average mass per particle of the ensemble, and is useful for linking two-moment microphysical schemes to bulk optical properties (Baran et al., 2025).

2.2.3 Removal of shattering events

Two particle size distributions (PSDs) were derived: a standard PSD including all detected trigger events, and a screened PSD including only particles classified as intact based on visual inspection of stereo-microscopic images. For the standard PSD, potential shattering events were filtered by applying an interarrival time threshold of 0.1 ms, based on the appearance of a distinct secondary mode in the interarrival time distribution (see Fig. S3).

We emphasize that a visual classification step is essential, as the inter-arrival time method alone cannot detect all shattering events due to the electronic dead time of the PHIPS detection system (12 µs), during which closely spaced fragments cannot be resolved as separate events. At typical aircraft speeds, this corresponds to a spatial separation of roughly 2 mm. However, these shattering events can be captured by the imager, which has a frame size of approximately 1.5 mm×1.1 mm.

Residual shattering artefacts were therefore removed through manual visual classification following the criteria described in Sect. 2.2.1. The impact of this additional screening is illustrated in Fig. 2, which compares the mean particle size and projected area distributions before and after removal of shattering artefacts for mid-latitude and Arctic cirrus. The largest differences occur for D_eq<100 µm, where the interarrival time method alone overestimates the number concentration due to undetected fragments. After image-based screening, the number concentration decreases toward smaller sizes, producing a more physically consistent distribution.

Figure 2 The mean PHIPS particle size distribution (PSD) and projected area size distribution for mid-latitude and Arctic cirrus. The solid lines show the mean based on 1 Hz PHIPS data in cirrus conditions ($T<-\mathrm{38}\mathrm{°}\mathrm{C}$ and $\text{IWC}>{\mathrm{10}}^{-\mathrm{5}}$ g m⁻³) and the dashed lines the mean for only crystals with stereo-images that were manually classified as intact particles.

This refinement primarily affects microphysical parameters sensitive to small particles, such as the mean mass radius $r_{\overline{m}}$, whereas quantities weighted by projected area or extinction coefficient (e.g. r_eff) remain largely unchanged (see Fig. S5). The manually screened dataset thus provides a more physically realistic representation of the shape of the size distribution of the true ice particle population but is underestimating the number concentration and IWC as only a subset of the triggering events accompanied by a stereo-microscopic image are included in the analysis.

2.2.4 Retrieval of ice asymmetry parameter

The asymmetry parameter g was retrieved from PHIPS measurements following the framework of Xu et al. (2022), with methodological refinements and an updated uncertainty assessment as detailed in Appendix A. In this method, the measured angular scattering function between 18–170° is represented by a truncated series of Legendre polynomials. The underlying assumption of the method is that at scattering angles greater than 18°, the measured signal is dominated by geometrical-optics (GO) contributions – namely refraction and reflection – with only minor contributions from diffraction (Takano and Liou, 1989; Macke et al., 1996b). The unmeasured forward range below 18° is assumed to contain the bulk of the diffraction peak, which is estimated from the measured particle size using a parameterised function derived from scalar diffraction theory.

The GO component of the phase function is reconstructed over the full angular range by fitting the measured data with a truncated Legendre polynomial expansion. The Legendre coefficients are obtained using Gauss–Legendre quadrature, which ensures optimal numerical accuracy. To keep all quadrature nodes within the PHIPS angular range, the maximum degree of the Legendre expansion is restricted to seven. This truncation implies that the method is applicable only to smooth phase functions that can be adequately represented by a low-order polynomial expansion.

We evaluated this requirement for our dataset and found that 98 % of the measured angular scattering functions show no distinct features such as 22° halo peaks. An uncertainty analysis was performed using 1612 ray-tracing simulations of hexagonal ice crystals with varying aspect ratios, sizes, and surface roughness, generated with the Macke et al. (1996b) ray-tracing code. This closure study demonstrates that, under realistic mixtures of ice crystal morphologies, the retrieval uncertainty of g ranges from approximately −0.001 to +0.02 when halo features are absent from the phase functions. A detailed description of the underlying assumptions, validation tests, and uncertainty quantification is provided in Appendix A.

3 Results

To analyse cirrus microphysical and optical properties by geographical region, the CIRRUS-HL observations were classified into mid-latitude and Arctic cirrus using a latitude threshold of 65° N. Cirrus clouds were defined based on a combination of in-cloud temperature below −38 °C and IWC greater than 10⁻⁵ g m⁻³, derived from the 1 Hz PHIPS data. The campaign covered a wide geographical range, from 38.2 to 75.6° N, with the highest cirrus cloud sampled at 14.5 km altitude and the coldest cirrus temperature recorded at −62.8 °C. In total, approximately 5625 km of mid-latitude cirrus and 1060 km of Arctic cirrus were sampled, reflecting the greater spatial extent and coverage of mid-latitude observations during the campaign.

3.1 Ice crystal habits in mid-latitudes and Arctic cirrus

Figure 3 shows the fraction of habits as a function of D_eq for the mid-latitude and Arctic cirrus cases. The most frequently occurring ice crystal habit class in both regions was the “other” category, accounting for 45 % and 66 % of all classified particles in mid-latitudes and Arctic, respectively. The high occurrence is especially pronounced for crystals smaller than 100 µm, for which over half of the measured particles fell into this group. In such cases, the classification to the “other” category can be due to the limited optical resolution available for morphological identification. Beyond this technical limitation, the “other” category also reflects the inherent complexity and variability of cirrus ice crystals, which often do not conform to traditional habit classifications (see Fig. 4 for example PHIPS images).

Figure 3 Ice crystal habit distribution as a function of the area equivalent diameter (D_eq) for mid-latitude and Arctic cirrus. The numbers indicate the number of ice crystal stereo-image pairs included in the habit analysis per size bin.

Figure 4 Example ice crystal microscopic images recorded by PHIPS from mid-latitude and Arctic cirrus as a function of their area equivalent diameter (D_eq). The number fractions of the different habits corresponds to Fig. 3.

Irregular and compact ice crystals, broadly corresponding to our “other” category, have frequently been reported as the dominant habit type in both mid-latitude and Arctic cirrus (Garrett et al., 2001; Lawson et al., 2019; Wolf et al., 2018). For example, Wolf et al. (2018) observed that compact or irregular habits accounted for 50 % and 90 % of the crystal populations in Arctic liquid-origin and in-situ cirrus, respectively. In addition to observations from high-resolution imaging instruments, scanning electron microscope analyses have also confirmed the dominance of irregular shapes in cirrus clouds, even among small ice crystals (Magee et al., 2021).

The most notable difference between mid-latitude and Arctic cirrus was a more frequent occurrence of aggregation in mid-latitudes compared to the Arctic. For mid-latitude cirrus aggregates was the second most frequent habit class, accounting for 32 % of all classified crystals, whereas in Arctic cirrus they accounted only for 7 % of all classified crystals. In mid-latitudes the fraction of aggregates was observed to clearly increase towards larger sizes with the highest fraction of 81 % found for the largest size group included in the analysis, between 400–500 µm. The more frequent occurrence of aggregates in mid-latitude cirrus can be explained through the more dynamic atmosphere at the lower latitudes (Barahona et al., 2017) and higher initial ice number concentration (Sourdeval et al., 2018), potentially promoting aggregation.

In mid-latitude cirrus most of the aggregates were consisting of plates (34 %), suggesting a convective origin, or of crystals classified as “other” (27 %). A large fraction of aggregates consisted of frozen droplets (25 %) but these, together with the simultaneously observed single frozen droplets, originated from research flight 12 targeting a strong overshooting convective system (see Schäfler and Rautenhaus, 2023). In contrast, frozen droplets and their aggregates were rarely detected in other mid-latitude cases, indicating that their occurrence was closely linked to specific convective conditions.

The second notable difference between the two regions was a higher occurrence of bullet and mixed rosettes in the Arctic. Bullet and mixed rosettes comprised 19 % of all classified crystals and dominated the population larger than 200 µm, comprising 56 %–89 % in that size range. This makes them the second most frequent habit type in Arctic cirrus. A similar dominance of bullet rosettes as the second most common habit was also reported in Arctic conditions by Wolf et al. (2018). For mid-latitude cirrus bullet and mixed rosettes were common but less frequent comprising 11 % of all crystals and reaching maximum occurrence frequencies of 25 %–28 % in the size range of 250–400 µm.

Other identified habits in both regions, each contributing between 2 %–9 % to the total population, included side planes, as well as single columns. In contrast, single plates were rare in Arctic cirrus (0.5 %) compared to 3.5 % in mid-latitude cirrus.

Figure 5 shows the ice crystal aspect ratios derived from the stereo-microscopic imager for different crystal sizes. The majority of the observed ice crystals exhibit median aspect ratios of approximately 0.8 for sub-100 µm particles and around 0.7 for larger crystals. Mid-latitude cirrus consistently shows lower aspect ratios than Arctic cirrus across all size bins. It should be noted, however, that the aspect ratios derived from the stereo-microscopic images represent the overall crystal geometry. Individual components of aggregates or polycrystals, such as the bullets forming bullet rosettes, can exhibit substantially higher aspect ratios. For example, Wagner et al. (2025) reported aspect ratios between 2–6 for individual bullets within bullet rosettes.

Figure 5 Ice crystal aspect ratio for mid-latitude and Arctic cirrus using a latitude threshold of 65°. Aspect ratio is defined as the minimum of the two orthogonal aspect ratios derived from particle imaging. Boxes represent the interquartile range and the central line denotes the median.

3.2 Cirrus microphysical and optical properties in mid-latitudes and Arctic

Figure 6 presents the frequency distribution of key bulk microphysical properties of cirrus (N, IWC, β_ext, and r_eff) along with g as a function of latitude, derived from PHIPS observations for the D_eq size range of 15–700 µm. The selected microphysical parameters are commonly used as input for radiative transfer schemes and were derived from 10 s averaged PHIPS particle size distributions. The asymmetry parameter was calculated from populations of 50 consecutive ice crystals and includes only manually inspected and classified intact (non-shattered) particles. Similarly, r_eff was derived using only non-shattered, imaged particles to ensure accurate representation of individual crystal morphology. In contrast, N, β_ext, and IWC were calculated using the full dataset, including particles without images, since applying the same shattering removal would lead to a systematic underestimation of these quantities. The minimum detection limits for the 10 s averages were approximately 0.5 L⁻¹ for N, 0.005 km⁻¹ for β_ext, and 0.05 mg m⁻³ for IWC, based on the instrument's sampling volume and estimated particle mass.

Figure 6 Frequency distribution of total concentration (N), ice water content (IWC), extinction coefficient (β_ext), effective radius (r_eff) and asymmetry parameter (g) observations from the PHIPS as a function of latitude. Microphysical observations are 10 s average values and asymmetry parameter is retrieved for only particles with stereo-images where shattering artefacts are excluded. The black solid lines mark the median line. Latitude bin widths are 1° for microphysical observations and 2° for the asymmetry parameter observations. The histograms show the distribution of the values for all mid-latitude (ML) and Arctic cases using a latitude threshold of 65° N. The dashed line over the histograms indicate the median value.

The median total concentration of ice crystals in the 15–700 µm size range varied between 1–52 L⁻¹ and showed no clear latitudinal trend. However, slightly higher overall median values were observed in mid-latitude cirrus (11 L⁻¹) compared to Arctic cirrus (8 L⁻¹) (Table 1). Mean values differed more noticeably, with 120 L⁻¹ in the mid-latitudes and 8 L⁻¹ in the Arctic. This difference was primarily driven by occasional high concentrations exceeding 100 L⁻¹ in the mid-latitude dataset, which were not observed in the Arctic. However, these elevated concentrations are likely influenced by particle shattering effects.

Table 1 Summary statistics (median, mean, and percentiles) of microphysical and optical properties derived from PHIPS measurements in mid-latitude (ML) and Arctic cirrus clouds. Shown are total number concentration (N), ice water content (IWC), effective radius (r_eff), and asymmetry parameter (g). The analysis includes particles in the PHIPS measurement range of 15–700 µm.

^* Lower detection limits for 10 s averages are approximately 0.5 L⁻¹ for N, 0.05 mg m⁻³ for IWC and 0.005 km⁻¹ for β_ext.

Similar to the total ice number concentration, both the IWC and β_ext values did not exhibit a significant latitudinal trend between 42–75° N, where most observations were conducted. Median IWC values ranged from 0.2 to 10 mgm⁻³, with an overall median of 2.9 mgm⁻³ and 3.0 mgm⁻³ for mid-latitudes and Arctic cirrus, respectively. Also the mean values were similar for mid-latitude and Arctic cirrus, 2.6 and 2.4 mgm⁻³, respectively, which confirms that the higher number concentrations observed in the mid-latitudes are from small ice crystals that do not contribute significantly to the IWC. Median β_ext values ranged from 0.04 to 0.4 km⁻¹, with the same overall median value for both mid-latitude and Arctic of 0.13 km⁻¹.

The only microphysical parameter that shows a clear latitudinal dependence was r_eff, with median values increasing from approximately 24 µm within latitudes from 39 to 42° N to 45 µm within latitudes above 70° N. An exception to this trend was observed near 62° N, where median effective radii around 50 µm were recorded in a liquid-origin system sampled off the coast of Norway during the double flight F18/F19. A similar trend towards higher r_eff values in more pristine environments was observed by Gayet et al. (2004) during the Interhemispheric Differences in Cirrus Properties From Anthropogenic Emissions (INCA) field experiments, where the average northern hemispheric effective radius (36 µm) was smaller than that observed in the Southern Hemisphere (42 µm).

In contrast to the trend in r_eff, g showed no systematic latitudinal variation, with median values ranging from 0.717 to 0.777. Notably, the median asymmetry parameter was identical in both regions (0.738) indicating a remarkable similarity in the bulk optical scattering characteristics.

The mid-latitude cases covered a broader range of cirrus conditions, while Arctic sampling was limited to five predominantly liquid-origin cirrus systems north of 65° N. Due to pandemic-related base relocation to Germany, Arctic flights required long transits and thus targeted synoptic cirrus with predictable development: a warm frontal cirrus (F05), two liquid-origin systems (F08/F09), a jet-stream-associated liquid-origin cirrus (F10/F11), and a mix of liquid-origin and in-situ cirrus over Iceland (F13/F14, Fig. 1). This focus introduces a bias toward liquid-origin cirrus, which have been shown to exhibit higher IWC compared to in situ-formed cirrus (Luebke et al., 2016). This sampling bias should be taken into account when interpreting observed latitudinal trends in cirrus microphysical and optical properties.

To better account for thermodynamic influences on cirrus properties, it is useful to examine the observations as a function of in-cloud temperature rather than geographic location. Figure 7 shows the statistical distribution of N, IWC, β_ext, r_eff, and g, comparing mid-latitude (red) and Arctic (blue) observations. Median N was observed to be higher in the mid-latitudes for temperatures down to approximately −55 °C. At colder temperatures, Arctic values become comparable to, or even exceed, those in the mid-latitudes. However, the observed differences are not statistically significant at the coldest temperatures below −55 °C and between −44 and −48 °C (two-sample t test, p>0.05).

Figure 7 Statistical analysis of the observed total number concentration (N), ice water content (IWC), extinction coefficient (β_ext), effective radius (r_eff), and asymmetry parameter (g) as a function of in-cloud temperature in degree Celsius. The solid line shows the median values and the shaded area the interquartile range for mid-latitude (ML; red) and Arctic (blue) observations using a latitude threshold of 65°.

Both IWC and β_ext values were generally similar in both regions at temperatures above −55 °C and showed a decreasing trend with decreasing temperature. At colder temperatures values observed in Arctic cirrus exceeded those observed in mid-latitude, however, due to the limited number of observations this difference is not statistically significant (p>0.05). The decreasing trend in IWC with decreasing temperature is well documented also in previous studies (Krämer et al., 2016; Jeggle et al., 2025).

The r_eff was mostly larger in the Arctic compared to the mid-latitudes with median values exceeding those in the mid-latitudes by 5 µm and showing a decreasing trend with temperature for temperatures down to −55 °C. At temperatures colder than −55 °C, the r_eff of Arctic cirrus increases again, showing up to twice as large r_eff. The observed differences were found to be statistically significant in the temperature ranges between −38 and −40 °C, and between −44 and −52 °C.

The asymmetry parameter showed a slight increasing trend towards colder temperatures with values around 0.72 measured at −40 °C and values between 0.75–0.76 at −60 °C observed in mid-latitudes. Lower values around 0.73 were observed in Arctic cirrus at temperatures below −55 °C, but as with the other parameters, the observed differences at colder temperatures are not statistically significant.

Inspection of Table 1 shows that the microphysical properties derived from the PHIPS instrument are in general agreement with previously reported values in the mid-latitudes and Arctic (Luebke et al., 2016; Heymsfield et al., 2017; Wolf et al., 2018). The median N values are systematically higher than those reported in the literature based on optical array probes with larger sampling volumes (Luebke et al., 2016; Heymsfield et al., 2017; De La Torre Castro et al., 2023). This discrepancy arises because low concentrations, below approximately 0.5 L⁻¹, are not detected with PHIPS within the 10 s averaging period. Values of N and IWC obtained from optical array probe measurements during the CIRRUS-HL campaign are discussed in detail by De La Torre Castro et al. (2023) and Jurkat-Witschas et al. (2025). Furthermore, the observed r_eff in both regions are significantly lower than the values reported by De La Torre Castro et al. (2023) for the same study. This discrepancy is most likely due to their assumption of a simplified spherical ice particle shape to derive this quantity from the measured particle size distribution.

3.3 Sensitivity of g to ice microphysical and environmental factors

To assess how g varies with environmental conditions and ice microphysics, we analysed coincident measurements of in-cloud temperature, relative humidity over ice (RH_i), and particle-resolved microphysical properties from PHIPS. The analysis was performed over the same averaging windows used for the g retrievals, each comprising 50 consecutive intact ice crystals. These correspond to averaging times of 44 s to 4 min (25th–75th percentile) and sampled air volumes of approximately 3–20 L.

Table 2 summarises the Pearson correlation coefficients between g and the investigated parameters. Among the tested parameters, the strongest correlation was observed with D_eq ($r=-\mathrm{0.42}$), followed by $r_{\overline{m}}$ ($r=-\mathrm{0.35}$) and r_eff ($r=-\mathrm{0.31}$), all indicating moderate negative correlations. These results suggest that g decreases with increasing particle size. Only weak correlations were found with temperature, RH_i, AR, and N.

Table 2 Pearson correlation coefficients (r) and p values for the relationship between g and various atmospheric and microphysical parameters. Moderate correlations are highlighted.

Figure 8 shows the frequency distribution of the measured g as a function of r_eff for ML and Arctic cirrus cases. The same figure for $r_{\overline{m}}$ is shown in Fig. S15. It can be seen that g values vary between 0.65–0.79. In the ML dataset, the median g decreases from 0.749 to 0.730 with increasing r_eff up to approximately r_eff=40 µm, beyond which the trend reverses and the median g increases with particle size. However, this upward trend at larger sizes is not statistically robust due to the limited number of observations in that range. In Arctic cirrus, the data are more sparse overall, but a consistent decrease in the median g from 0.747 to 0.709 with increasing r_eff is observed across the measured range.

Figure 8 Frequency distribution of asymmetry parameter (g) observations as a function of effective radius (r_eff). Both the asymmetry parameter and effective radius are retrieved only for particles with stereo-microscopic image that were not classified as shattering. The black line indicated the median value for the effective radius bins with a bin width of 3 µm. Existing short wave optical parameterisations are overlaid to the measurement data. These are the parameterisation by Fu (1996) (Fu; blue line), general habit mixture (GHM, dashed black line; Cole et al., 2013), the severely roughened column aggregate model used in MODIS Collection 6 retrievals (MC6; Platnick et al., 2016) and the optical parameterisation used in ECHAM6 general circulation model (ECHAM; Pincus and Stevens, 2013).

In radiative transfer simulations, ice optical properties are commonly parameterised as a function of r_eff. Figure 8 shows four commonly used parameterisations of g at a wavelength band centered at 532 nm, corresponding to the PHIPS scattering measurements. The first parameterisation, developed by Fu (1996), is based on columnar crystals with single scattering properties calculated using an improved geometric optics method. The asymmetry parameter is expressed as a function of a generalised effective size, which can be related to r_eff via Eq. (3.12) in Fu (1996). The second parameterisation represents the General Habit Mixture (GHM), which assumes a combination of nine crystal habits, all with severe surface roughness (σ=0.5) (Cole et al., 2013). The third is the MODIS Collection 6 (MC6) parameterisation, which represents ice particles as aggregates of severely roughened columns (Platnick et al., 2016). The fourth is the parameterisation implemented in the ECHAM6 general circulation model (Pincus and Stevens, 2013).

Among the four, only the MC6 parameterisation – assuming complex, severely roughened aggregates with a fixed asymmetry parameter of 0.75 – captures the upper range of the observed values, but does not capture the median values or the lower range. The remaining parameterisations overestimate g across the entire size range and fail to reproduce majority of the measured values. Moreover, three ou of four parameterisations predict an increasing trend in the g with increasing r_eff, which contrasts with the observational data showing a clear decreasing trend.

3.3.1 Analysis based on individual particle classes

A key advantage of PHIPS as a single particle polar nephelometer is its ability to directly link the microphysical and optical properties of individual ice crystals. This allows for the selection of particles with specific microphysical characteristics, such as size or aspect ratio, and enables the construction of “artificial” particle populations consisting of crystals with similar properties. This approach results in a very detailed investigation of how individual microphysical parameters influence g, isolated from the natural microphysical variability present in real cloud segments. This approach was recently used to retrieve the asymmetry parameter of bullet rosettes from PHIPS data captured during CIRRUS-HL (Wagner et al., 2025).

Figure 9 presents a statistical analysis of g for all intact ice crystals measured during CIRRUS-HL. The ice crystals were grouped based on their D_eq (for bin definitions, see Waitz et al., 2021), AR (defined as the smaller of the two aspect ratios from the stereo images) using a bin width of 0.1, and simultaneously measured ice supersaturation (S_ice), also using a bin width of 0.1. Each bin contains a minimum of 100 ice crystals. To retrieve g, the crystals within each bin were further divided into groups of 50 consecutive particles, and one g was calculated per group. The number of these groups is indicated in the figure. Additionally, a single g was retrieved for all crystals within each bin, with the corresponding value and uncertainty shown as markers.

Figure 9 Statistical analysis of the asymmetry parameter (g) for mid-latitude (ML) and Arctic cirrus as a function of (a) particle area equivalent diameter (D_eq), (b) aspect ratio, and (c) ice supersaturation (S_ice). The box indicates the 25th and 75th percentiles (quartiles), the line the 50th percentile (median) and the whiskers the 5th and 95th percentiles of the g values retrieved for samples of 50 crystals. The markers indicate the asymmetry parameter retrieval for all crystals within the bin with error bars denoting the retrieval uncertainty. Sample sizes corresponding the number of asymmetry parameter retrievals for groups of 50 particles for each category are annotated in the plots.

Among all variables considered, only D_eq exhibits a clear decreasing trend with g, consistent with the cloud-scale analysis presented earlier. The highest median values (g=0.777), were measured for mid-latitude ice crystals with D_eq between 15–30 µm, while the lowest median values (g=0.719) were observed for D_eq between 150–200 µm. For larger particles, the g increases again, although this trend is not statistically robust due to fewer observations.

Figure 10 further examines this size dependence of g by separating the data according to ice crystal habit. The analysis was conducted for the seven most frequently occurring habit categories: aggregate, column, side plane, frozen droplet, bullet rosette, plate, and the “other” category. For each habit–size combination, a median g was retrieved for each size bin, provided that at least 50 intact crystals were available in the bin.

Figure 10 Mean asymmetry parameter, g, as a function of particle diameter for various ice crystal habits, including aggregates, plates, columns, frozen droplets, side planes, and bullet rosettes. The shaded areas represent the retrieval uncertainty in asymmetry parameter for each habit.

The general decrease of g with increasing size is apparent across most habits, except for frozen droplets, which occur only in the smallest size bins. A maximum difference of 0.05 in g is found between columns, which show the lowest values, and aggregates, which exhibit the highest values. This difference exceeds the measurement uncertainty and is therefore considered to be significant. Side planes, bullet rosettes, plates, and crystals in the “other” category all show similar g values within the uncertainty range for sizes below 200 µm. For larger bullet rosettes, g increases and becomes comparable to that of aggregates of similar size.

No significant difference was found in the median g between mid-latitude and Arctic cirrus with median g between 0.723–0.758 observed for each AR and S_ice category. In general, the variability in g was greater in the mid-latitudes across all categories (g varied between 0.686–0.816 for groups of 50 crystals) compared to the Arctic (g varied between 0.691–0.787). However, this may be partly attributed to the smaller sample size in the Arctic.

4 Discussion

Despite differences in geographic location and the ice crystal habit distributions, no significant variation in g was observed between mid-latitude and Arctic cirrus, with consistently low median values in both environments. Similar low asymmetry parameters below 0.8 have been systematically reported by other airborne polar nephelometers in the visible spectral range. Measurements with the CIN instrument have shown g values around 0.74±0.03 for Arctic ice clouds (Gerber et al., 2000; Garrett et al., 2001) and values around 0.75±0.01 for tropical anvil cirrus (Garrett et al., 2003). In mid-latitudes, median g values between 0.769–0.79 have been measured using the PN instrument (Gayet et al., 2004, 2012; Febvre et al., 2009; Auriol et al., 2001).

Our results show that g is largely insensitive to variations in atmospheric conditions (temperature and relative humidity), as well as to most microphysical properties of ice crystals. This finding is consistent with several previous observational studies. For example, Febvre et al. (2009) reported comparable optical properties for natural cirrus and aged contrails, despite differences in their size distributions. Likewise, Garrett et al. (2001) found that g values were insensitive to habit variations in Arctic cirrus. During the INCA experiment, comparable g values were derived for the northern and Southern Hemispheres (Gayet et al., 2004, 2006), with only a weak temperature dependence noted by Gayet et al. (2006).

Here, the only parameter for which a moderate correlation with g was observed is particle size: we find that g values tend to decrease with increasing D_eq, r_eff or $r_{\overline{m}}$. This behaviour differs from the results of Gayet et al. (2004), who reported an increase in g with effective diameter, whereas Garrett et al. (2003) observed no clear dependence between g and r_eff for small anvil ice crystals with r_eff<20 µm. Interestingly, the aspect ratio does not appear to influence g in our data, which stands in contradiction to theoretical studies predicting a strong dependence between g and AR, even for roughened crystals (Iaquinta et al., 1995; van Diedenhoven et al., 2014a).

The low median g value of 0.738 cannot be explained by pristine hexagonal shapes (Iaquinta et al., 1995; Schmitt et al., 2006). Rather, a significant degree of surface roughness or other deformations are needed in order to numerically reproduce our observations (Macke et al., 1996a; Yang et al., 2013; Wagner et al., 2025). Given that the 95th percentiles of the g distributions are 0.773 for mid-latitude and 0.760 for Arctic cirrus, we conclude that the majority of the measured ice crystals exhibit a significant degree of morphological complexity. These findings suggest that crystal complexity plays a more dominant role in determining g than any other microphysical property. The observed size dependence of g may also reflect an increasing degree of crystal complexity in larger ice particles.

In-situ measurements have reported g values above 0.8 only for a few cases. Garrett et al. (2001) documented g values of 0.81 for sublimating ice crystals at RH_i around 87 %, and Gayet et al. (2012) found similar higher g values for RH_i below 80 %. Febvre et al. (2009) observed g values of 0.827 for young contrails, with a mean RH_i of 118 %. In our dataset, we did not observe an increase in g at lower RH_i values. However, our observations included only a limited number of ice crystals at RH_i below 80 %. Therefore, the effect of sublimation on g remains inconclusive and needs further investigation.

A limitation of this study is that g retrievals were restricted to ice crystals with D_eq below 300 µm. This constraint arises because larger particles frequently saturated the first few measurement channels of our polar nephelometer, preventing a reliable retrieval of g. However, as shown in the particle size and area distributions in Fig. 2 and the cumulative distributions in Fig. S4, ice crystals smaller than 300 µm account for 97 % (99 %) of the total number and 83 % (94 %) of the total cross-sectional area in mid-latitude (Arctic) cirrus. Based on this, we consider our retrieved g values to be representative of the cirrus cloud populations sampled during the CIRRUS-HL campaign.

Our results show a substantial discrepancy between the parameterisations of g used in climate models and the values measured in this study (see Fig. 8). This difference has significant implications for inferred cloud reflectivity or albedo, which can be demonstrated using a simple approximation of the plane albedo for thin cirrus clouds (Meador and Weaver, 1980)

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }R={\displaystyle \frac{{\mathit{\omega }}_{\mathrm{0}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}}(\mathrm{1}-g)\mathit{\tau },\end{array}$$

where ω₀ is the single-scattering albedo, μ₀ is the cosine of the solar zenith angle, and τ is the optical thickness. For constant optical thickness and solar zenith angle, the plane albedo is directly proportional to 1−g. Therefore, a decrease of g from 0.8 to the measured 0.738 corresponds to an increase in albedo of approximately 31 %, making cirrus clouds more reflective at solar wavelengths. This indicates that an overestimation of g can lead to a significant bias in the estimation of the radiative effect of cirrus clouds. Additionally, such an overestimation can introduce errors in the retrieval of cirrus optical depth from space-borne radiance measurements (Yang et al., 2008; Yi et al., 2017a).

Studies using general circulation models (GCMs) or Moderate Resolution Imaging Spectroradiometer (MODIS) retrievals to simulate top-of-the-atmosphere (TOA) fluxes have shown that incorporating optical parameterisations with lower g values in the visible spectrum leads to a closer agreement between these TOA flux simulations and TOA flux estimates from direct measurements of radiances using space-borne remote sensing (Kristjánsson et al., 2000; Yi et al., 2017b). However, Ren et al. (2023) found that optical parameterisations with low g values around and below 0.75 overestimated both shortwave (SW) and longwave (LW) fluxes at TOA compared to observations from the Clouds and the Earth's Radiant Energy System (CERES). Their study emphasised the need for further optimisation of the scattering phase function to achieve radiative closure. Yi et al. (2017b) also showed that lowering g (by changing the optical parameterisation from MC5 to MC6) significantly reduced the estimated TOA SW flux from MODIS observations. This reduction resulted from lower cloud optical thickness retrieved by the MC6 model, which was subsequently used as input to the radiative transfer simulations (Yi et al., 2017a).

Currently, cirrus and ice containing ice clouds have been estimated to exert a globally and annually averaged net warming effect of 2–5 W m⁻² (Hong et al., 2016; Yi et al., 2017b; L'Ecuyer et al., 2019). Modelling studies using GCMs have shown that lowering the g values in visible from values around 0.8–0.75 will decrease the net ice CRE by 1–2 W m⁻² (Järvinen et al., 2018; Yi et al., 2013; Yi, 2022). Our results suggest that an additional decrease of the g to values below 0.75 might be needed, which would lead to an even higher bias in CRE.

A potential bias in CRE estimates is highly significant in the context of the Earth's energy balance. The current estimate of the global radiative imbalance is on the order of 0.5–1 W m⁻² (Loeb et al., 2021), meaning that even small biases in cirrus CRE can mask or amplify the signals of anthropogenic forcing. Cirrus clouds cover up to 30 % of the globe at any given time (Stubenrauch et al., 2013), and errors in their radiative properties can accumulate spatially and temporally, introducing systematic uncertainty into global climate models. For instance, an overestimation of the asymmetry parameter g leads to an underestimation of cirrus reflectivity, allowing more solar radiation to reach the surface. In global climate models, this can trigger dynamic feedbacks: increased surface heating may destabilize the atmosphere, enhance convection, and promote the formation of low-level clouds. In parallel, GCMs are often tuned to maintain top-of-atmosphere energy closure (Mauritsen et al., 2012). If cirrus reflectivity is systematically too low, the tuning can be made at the expense of physical realism. Improving the representation of cirrus optical properties, particularly g, is therefore essential for accurately quantifying cloud–radiation interactions and reducing uncertainty in climate projections.

5 Summary

We analysed ice crystal microphysical and correlated angular light scattering properties of mid-latitude and Arctic cirrus clouds using in situ measurements performed with the PHIPS airborne cloud probe during the CIRRUS-HL campaign. Ice crystal number concentration and ice water content showed no clear latitudinal trends, with similar median values of 11 L⁻¹ and 2.8 mg m⁻³ for mid-latitude cirrus and 8 L⁻¹ and 3.0 mg m⁻³ for Arctic cirrus. In contrast, the effective radius increased toward higher latitudes. Furthermore, clear regional differences in ice crystal habit distributions were observed: aggregates were the dominant habit in mid-latitude cirrus for crystals larger than 200 µm, while bullet and mixed rosettes dominated Arctic cirrus. Smaller crystals in both regions were primarily complex in shape and unclassifiable using standard habit categories.

Despite the differences in habits, both mid-latitude and Arctic cirrus exhibited nearly identical asymmetry parameter values, with the same median asymmetry parameter of 0.738 for both regions and observed values as low as 0.686 (1th percentile). The asymmetry parameter was largely insensitive to environmental conditions such as temperature and relative humidity, as well as to ice crystal aspect ratio and bulk microphysical parameters including IWC and number concentration. Moderate correlations were found only with particle size.

The relatively low asymmetry parameter values cannot be represented using optical models for idealised pristine crystals. Comparison with four commonly used radiative transfer parameterisations revealed that most overestimated asymmetry parameter by up to 0.1. Only the MODIS Collection 6 parameterisation, which assumes severely roughened columnar crystals, captured the upper range of the observed values. Given the inverse relationship between the asymmetry parameter and cloud reflectivity, overestimating asymmetry parameter from 0.738 to 0.800 can lead to an underestimation of cirrus albedo by up to 31 %, with implications for both radiative forcing estimates and satellite retrievals of optical depth. Such biases are particularly important because cirrus clouds have estimated to exert a net warming effect of 2–5 W m⁻² globally, and reducing asymmetry parameter in climate models by just 0.05 has been shown to lower this warming by 1–2 W m⁻². These discrepancies highlight the need for improved observational constraints and parameterisations of asymmetry parameter to better quantify the radiative role of cirrus clouds and reduce uncertainties in climate projections.

Appendix A: Accuracy and Assumptions of the Asymmetry Parameter Retrieval

A1 Overview of the retrieval framework

Representing the scattering phase function as a series of Legendre polynomials is a common approach in atmospheric sciences and is widely used in radiative-transfer and remote-sensing applications (Wiscombe, 1977; Hu et al., 2000). This formulation is convenient because the Legendre polynomials are harmonic functions of cos θ and, thus, provide an efficient mathematical basis for describing the angular dependence of scattered light. As a general rule, the number of Legendre terms, N, required to accurately represent a phase function increases with the particle size parameter, x. For example, in case of large ice crystals (x>100) the narrow forward peak requires a high number of expansion terms, up to several thousands, to resolve the original phase function.

Figure A1 Asymmetry parameter from external diffraction calculated for orientation averaged hexagonal columns with different aspect ratios (AR). Both 3D geometry and 2D projection were considered. The blue line shows calculation for orientation averaged circular aperture with the same projected area. The black line shows the parameterisation used in Xu et al. (2022) and red line the updated parameterisation.

To reduce some of this complexity, the first step in our retrieval framework is to separate the diffraction contribution from the geometrical-optics (GO) component of scattering. Although this decomposition is not fundamental from the electromagnetic perspective, it is a well-established and practical approximation for large atmospheric ice particles and widely used in conventional GO methods (e.g. Takano and Liou, 1989; Macke et al., 1996b). Studies using improved GO methods that avoid this artificial separation have shown convergence toward conventional ray-tracing results for the large size parameters relevant to our study (Yang and Liou, 1996).

This separation is also convenient because the diffraction contribution within the PHIPS measurement range (18–170°) can be considered negligible. This assumption is supported by theoretical studies showing that for large ice particles, the diffraction intensity decreases by more than three orders of magnitude beyond 18° (e.g. Hesse, 2008; Chunlei and Keya, 1996; Travis et al., 2002). Consequently, diffraction primarily affects the unmeasured forward-scattering region and can be parameterised to reconstruct the full phase function.

In earlier implementation by Xu et al. (2022), the diffraction contribution was approximated using scalar diffraction theory for circular apertures, which provides a computationally efficient first-order description. To further improve the physical consistency of the diffraction treatment, here we update the diffraction parameterisation to account for non-spherical particle shapes.

The diffraction calculations were carried out using a vector Fraunhofer integral. The external diffraction is computed by replacing the illuminated particle surface with a set of complementary apertures, in accordance with Babinet's principle. Three levels of geometric representation were considered: (i) a full three-dimensional hull, in which each illuminated facet is treated as an individual aperture and its diffracted field is computed accounting for relative position and orientation, with the total scattered field obtained from the coherent sum of all facets; (ii) a two-dimensional projection, where the illuminated facets are projected onto a plane normal to the incident light, reducing the problem for convex particles to a single effective aperture; and (iii) an area-equivalent circle, in which the 2D projected particle cross section is represented by a circular aperture of cross section-equivalent area. For all geometries, 100 random orientations were simulated, and the orientation-averaged one-dimensional scattering intensity as a function of angle was used to derive the diffraction asymmetry parameter g_D.

The diffraction calculations shown in Fig. A1 confirm that, for orientation-averaged ice crystals, aperture shape and aspect ratio have only a minor influence on the diffraction phase function and the resulting diffraction asymmetry parameter g_D. For particle sizes (d≥15 µm) and aspect ratios ($\mathrm{0.25}<\text{AR}<\mathrm{1.7}$) relevant to this study, the spread in g_D across the different geometric representations remains below 0.002 and decreases further with increasing size. On the other hand, the diffraction parameterisation adopted by Xu et al. (2022), which is based on a fixed circular aperture, yields systematically lower values of g_D(d).

Motivated by this systematic bias, we derived an updated diffraction parameterisation by fitting a fourth-degree polynomial to the ensemble-mean g_D obtained from all diffraction models. The underlying simulations span hexagonal ice crystals with aspect ratios between 0.25–1.7, consistent with the range observed by PHIPS. On a logarithmic size scale, the diffraction asymmetry parameter g_D(d) is approximated by

$$\begin{array}{ll}{\displaystyle }{g}_{\mathrm{D}}\left(d\right)=& {\displaystyle }-\mathrm{3.31551145}\times {\mathrm{10}}^{-\mathrm{4}}{\ln }^{\mathrm{4}}\left(d\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{5.44285886}\times {\mathrm{10}}^{-\mathrm{3}}{\ln }^{\mathrm{3}}\left(d\right)\\ {\displaystyle }& {\displaystyle }-\mathrm{3.36880804}\times {\mathrm{10}}^{-\mathrm{2}}{\ln }^{\mathrm{2}}\left(d\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{9.35390828}\times {\mathrm{10}}^{-\mathrm{2}}\ln \left(d\right)\\ \text{(A1)}& {\displaystyle }& {\displaystyle }+\mathrm{9.01145674}\times {\mathrm{10}}^{-\mathrm{1}},\end{array}$$

where d denotes the particle area-equivalent diameter in micrometres.

The revised parameterisation reproduces the ensemble-mean diffraction behaviour to within 0.002 over the particle size and aspect ratio range relevant to this study and replaces the previous spherical-aperture-based formulation in the retrieval algorithm. Since diffraction contributes one half of the total asymmetry parameter, the remaining uncertainty associated with the diffraction treatment corresponds to a maximum uncertainty of approximately 0.001 in the retrieved total g.

After the diffraction phase function is parameterised, the measured PHIPS angular range can be treated as purely GO scattering. To obtain the GO contribution to g the measurement needs to be extrapolated to the missing forward and backscattering ranges. This is done by fitting a truncated Legendre polynomial expansion to the measured angular intensities. The GO scattering phase function P_GO(θ) is expressed as a finite series of Legendre polynomials P_l(cos θ):

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }{P}_{\text{GO}}\left(\mathit{\theta }\right)\approx \sum _{l=\mathrm{0}}^N(\mathrm{2}l+\mathrm{1}){\widehat{c}}_l{P}_l\left(\cos \mathit{\theta }\right),\end{array}$$

where ${\widehat{c}}_l$ are the Legendre expansion coefficients of degree l, and N is the maximum degree of the expansion. The factor 2l+1 follows the standard normalization of Legendre polynomials over the interval $[-\mathrm{1},\mathrm{1}]$.

The coefficients $\widehat{c_l}$ are obtained through orthogonal projection of the measured phase function onto the Legendre basis. In practice, this integration is carried out numerically using Gauss–Legendre quadrature, which provides optimal accuracy by evaluating the integrand at specific quadrature nodes x_i corresponding to the roots of P_N+1(x_i):

$$\begin{array}{ll}{\displaystyle }\widehat{c_l}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{P}_{\text{GO}}\left(\arccos \right(x\left)\right)P_l\left(x\right)\mathrm{d}x\\ \text{(A3)}& {\displaystyle }& {\displaystyle }\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^N{w}_{i}I\left(\arccos \right(x_i\left)\right)P_l\left(x_i\right),\end{array}$$

where w_i are the Gauss–Legendre weights associated with each node and I(arccos (x_i)) the measured intensity interpolated to the Gauss–Legendre quadrature nodes x_i.

In our retrieval, the number of quadrature nodes is set to eight, corresponding to a Legendre expansion of eight terms (N=7). This truncation ensures that all quadrature nodes remain within the PHIPS angular coverage (18–170°) and avoids extrapolation into unmeasured regions.

Using eight terms implies that the fitted Legendre expansion can accurately represent only smooth angular scattering functions. Therefore, the retrieval method is not applicable to pristine ice crystals that exhibit sharp features in their phase functions. A preliminary uncertainty analysis of this approach was presented by Xu et al. (2022), where the authors compared the retrievals results to true g values using simulations of hexagonal ice crystals with varying degree of surface roughness and aspect ratio. However, the uncertainty analysis was based on a limited number of crystal geometries. In the present work, we extend this analysis using a substantially larger set of simulated ice crystal morphologies to refine the uncertainty estimate and confirm the validity of the truncation choice.

A2 Updated uncertainty analysis

To extend the analysis of Xu et al. (2022), we performed additional closure tests using the ray-tracing model of Macke et al. (1996b). The new simulations include 1612 orientation-averaged configurations of hexagonal ice crystals with sizes between 20–500 µm and aspect ratios ranging from 0.68 to 0.85, motivated by the PHIPS microphysical observations. Surface roughness was implemented using Gaussian tilted-facet perturbations with roughness parameters (σ) between 0–0.9 (interval of 0.1). In addition, internal air inclusions were simulated using the mean free path (MFP) approach (Macke et al., 1996a), with MFP values between 10–10⁸ (no inclusions).

The purpose of this investigation was to evaluate how the retrieval error in g depends on the smoothness of the scattering phase function. To achieve this, the models were grouped according to their roughness parameter σ, which is directly related to the smoothness of the calculated scattering function (higher σ corresponds to smoother phase functions). Since natural cirrus clouds exhibit a diversity of shapes and aspect ratios, the simulations were averaged over ten different crystal configurations before calculating g. Although this ensemble size represents only a moderate level of morphological variability, sensitivity tests (not shown) confirmed that increasing the number of configurations had no effect on the results.

For each case, two values of g were computed: (i) the “true” value derived directly from the full simulated phase functions and (ii) the “retrieved” value obtained from the “true” phase function truncated to the PHIPS angular range (18–170°), corresponding to the angular coverage of the parabolic mirror detectors used in the PHIPS nephelometer. The truncation was carried out by integrating the simulated scattering intensity over the solid angle of each detector. Each detector covers an angular range around its central viewing angle with a half-angle acceptance of 3.5°, corresponding to a solid angle of 0.0117 sr. The simulated phase functions, computed with a 1° angular resolution, were integrated over each detector's acceptance using trapezoidal numerical integration.

Figure A2 illustrates the validation results for simulations with a low roughness parameter (σ=0.1). The simulated phase functions still contain a modest halo feature that is clearly visible even in the truncated phase function at the PHIPS angular resolution. It should be noted that the PHIPS measurements in panel c do not show a halo peak, even in the 1–99 percentile range, suggesting that higher roughness parameters are needed to model the measured phase functions. Panel b shows that the difference between the “true” and “retrieved” g can reach up to 0.08 but decreases rapidly with increasing complexity parameter C_p. The C_p parameter quantifies the decay rate of the Legendre expansion coefficients (Xu et al., 2022), with larger values corresponding to more featureless phase functions. While Xu et al. (2022) previously suggested that the retrieval is accurate within an absolute error of 0.001 for C_p≥0.4, our extended simulations reveal that for relatively pristine particles (σ=0.1), the retrieved C_p values derived from the truncated phase functions can exceed 0.4 (as shown in panel a), making the C_p parameter an unreliable quality indicator for assessing retrieval accuracy for lightly roughened crystals.

Figure A3 presents the corresponding results for σ=0.2. Increasing the surface roughness from 0.1 to 0.2 smooths out the halo feature in the simulations, yielding featureless phase functions that more closely resemble the measured data (Fig. A3c). The retrieval error in g is now below 0.02, although the retrieved values tend to be slightly biased toward lower g. Further increasing the roughness to σ=0.3 or higher values (not shown) reduces the retrieval error below 0.005 and also changes the bias sign, resulting in a slight overestimation of the true g.

These extended simulations confirm that the retrieval algorithm applied to the PHIPS angular range (18–170°) is reliable only for sufficiently smooth phase functions or sufficiently roughened ice crystals (in our simulations, σ≥0.2). We found that relatively pristine particles (σ=0.1) yield C_p values above the previously proposed C_p threshold of 0.4 (Xu et al., 2022), indicating that C_p cannot be used as a criterion for assessing retrieval accuracy. However, in all cases where the retrieval error exceeded 0.02, the simulated phase functions exhibited distinct halo features that would be clearly detectable at the angular resolution of PHIPS. Since such features are entirely absent from the measured phase functions (see panel c in Figs. A2 and A3), we conclude that the natural ice crystals in our dataset exhibit sufficient surface roughness for the retrieval method to remain valid. Based on these results, we adopt a more conservative uncertainty estimate of −0.001 to +0.02 for asymmetry parameter retrievals derived from PHIPS measurements compared to the value given in Xu et al. (2022).

Figure A2 Validation of the asymmetry parameter retrieval for simulations with a roughness parameter σ=0.1. The dataset includes 180 orientation-averaged crystal configurations, with g values retrieved for sets of 10 orientations each. Panel (a) compares the true complexity parameter C_p,true with the retrieved C_p, with the dashed line indicating the 1:1 relation. Panel (b) shows a heatmap of the retrieval error g_true−g_retrieved as a function of C_p, with red dashed lines marking ±0.001 error bounds. Panel (c) displays measured scattering phase functions for mid-latitude (blue) and Arctic (red) cirrus, including median, 25 %–75 %, and 1 %–99 % percentiles, with gray curves representing ray-tracing simulations of individual crystal configurations truncated to PHIPS measurement range.

Figure A3 Same as Fig. A2, but for simulations with a roughness parameter σ=0.2.