Using laboratory and field measurements to constrain a single habit shortwave optical parameterization for cirrus

A single habit parameterization for the shortwave optical properties of cirrus is presented. The parameterization utilizes a hollow particle geometry, with stepped internal cavities as identi ﬁ ed in laboratory and ﬁ eld studies. This particular habit was chosen as both experimental and theoretical results show that the particle exhibits lower asymmetry parameters when compared to solid crystals of the same aspect ratio. The aspect ratio of the particlewasvariedasafunctionofmaximumdimension, D ,inordertoadheretothesamephysicalrelationships assumed in the microphysical scheme in a con ﬁ guration of the Met Of ﬁ ce atmosphere-only global model, concerning particle mass, size and effective density. Single scattering properties were then computed using T-Matrix, Ray Tracing with Diffraction on Facets (RTDF) and Ray Tracing (RT) for small, medium, and large size parameters respectively. The scattering properties were integrated over 28 particle size distributions as used in the microphysical scheme. The ﬁ ts were then parameterized as simple functions of Ice Water Content (IWC) for 6 shortwave bands. The parameterization was implemented into the GA6 con ﬁ guration of the Met Of ﬁ ce Uni ﬁ ed Model along with the current operational long-wave parameterization. The GA6 con ﬁ guration is used to simulate the annual twenty-year short-wave (SW) ﬂ uxes at top-of-atmosphere (TOA) and also the temperature and humidity structure of the atmosphere. The parameterization presented here is compared against the current operational model and a more recent habit mixture model. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY


Introduction
In 2013, the Intergovernmental Panel on Climate Change (IPCC) concluded that the coupling of clouds with the Earth's atmosphere is the largest uncertainty faced in predicting climate change today (Intergovernmental Panel on Climate Change, 2013).One such cloud type that contributes to this uncertainty is cirrus due to their extensive global coverage of about 30%, with coverage reaching 60-80% in the tropics (Sassen et al., 2008).Furthermore, cirrus has diverse microphysical properties, containing a multitude of particle habits which range in size over several orders of magnitude.This variety in size, shape and complexity poses many difficulties for the accurate representation of ice cloud in climate models.The range in size means that currently no single method can be used to calculate the single scattering properties of ice crystals.For smaller sizes, exact solutions can be sought (Mano, 2000;Havemann and Baran, 2004;Groth et al., 2015), but as the size and complexity of the ice crystal increase, more approximate solutions are necessary (Macke et al., 1996b;Hesse et al., 2012).The diversity of particle shape also presents many challenges, and it is well established that the particle habit significantly impacts upon the single scattering properties of ice crystals (Macke et al., 1998;Bacon and Swanson, 2000;Baran, 2012;Baum et al., 2014).
The representation of single ice particles for scattering calculations has improved significantly over the years.Early studies used very simplified shapes such as spheres and cylinders, but these were found to be inadequate approximations for the treatment of ice crystals (Macke and Mishchenko, 1996).More realistic representations of particle habits such as bullet rosettes and aggregates have been constructed (Um and McFarquhar, 2007;Xie et al., 2011;Baran and Labonnote, 2007;Baum et al., 2014;Baran et al., 2014), and are commonly used in habit mixture models to represent cirrus.In addition to particle habit, small scale features such as surface roughness, inclusions and indentations/cavities have gained recognition as potentially important contributors to the scattering behaviour of ice crystals (Schmitt and Heymsfield, 2007;Yang et al., 2008;Schnaiter et al., 2011).Ice particles with deep indentations/cavities (typically on their basal facets) are commonly described as 'hollow'.Hollowness can significantly affect the asymmetry parameter, which is given by: where θ, is the scattering angle, defined as the angle between the incident and the scattered beam.P 11 , is the phase functiona normalised distribution of the intensity of radiation, scattered from a randomly oriented particle (van de Hulst, 1957), given by: The phase function, and by extension, the asymmetry parameter, hold information about the angular distribution of the scattered light, and as such, the asymmetry parameter is commonly used to implement ice crystal optical properties into a GCM.Theoretical results show that hollow hexagonal crystals exhibit a general increasing trend in the asymmetry parameter, suggesting that hollow crystals reflect less than their solid counterparts (Yang et al., 2008).On the contrary, it has been shown that rough particles could reflect up to twice the radiation when compared with pristine crystals (Ulanowski et al., 2006).The asymmetry parameter is typically over predicted in scattering models compared with results estimated from in-situ measurements (Gerber et al., 2000;Gayet et al., 2011).However, these in-situ results could have been affected by particle shattering on the inlets of the microphysical probes, thereby artificially decreasing the size and asymmetry parameter (Korolev et al., 2011).
To mimic naturally occurring surface roughness, many theoretical studies make use of the distortion parameter, which approximates surface roughness by tilting the surface normal for an incoming ray.The tilt angle is defined by a random number from zero up to a maximum, where the maximum is defined by the distortion parameter (Macke et al., 1996b).This method yields smoother phase functions with lower asymmetry parameters than when roughness is not accounted for.Another method used is that of particle inclusions: the incorporation of inclusions into the ice crystal model also smooths peak features such as 'halos' and 'ice bows', and reduces the asymmetry parameter, yielding more realistic values (Macke et al., 1996a;Yang and Liou, 1998;Labonnote et al., 2001;Baran and Labonnote, 2007).For optical parameterizations, many studies make use of a range of particle habits (Baum et al., 2005;Bozzo et al., 2008;Baran and Labonnote, 2007), although the distortion parameter is still widely used as a proxy for surface roughness.Other studies suggest that with the use of distortion, ice clouds can be represented entirely by hexagonal prisms of varying aspect ratio, negating the need to represent the full range of crystal habits (Van Diedenhoven et al., 2014).
Whilst the single scattering properties of cirrus ice particles are affected by micro-scale features, the bulk optical properties of cirrus are also affected by macro-scale properties such as cloud optical depth, ice mass in the cloud column and particle size distribution (PSD).As such, the cirrus net radiative effect is sensitive to the various assumptions made in the microphysical scheme.Typically, in a GCM, the bulk optical properties of clouds are parameterized in terms of the diagnostic variable effective dimension, D e , as a function of the Ice Water Content (IWC) and/or environmental temperature (T) (Bozzo et al., 2008;Fu et al., 1999).Effective dimension of the PSD is given by (McFarquhar and Heymsfield, 1998): where: density of ice or liquid water A i = the mean cross sectional area in bin i n i = number concentration in bin i.
The use of effective dimension assumes that the bulk-optical properties are uniquely defined by D e and either IWC and/or T, but it has been shown that they are also dependent upon the shape of the PSD (Mitchell et al., 2011;Baran, 2005).This dependency is not accounted for in D e based schemes as they tend to be physically inconsistent with the microphysical scheme (Baran, 2012), consequently the microphysical and radiative parameterizations may assume different PSDs.Whilst the use of D e is generally valid for water clouds, the relationship becomes unreliable for ice clouds and for more absorbing wavelengths (Mitchell, 2002;Baran, 2005;Mitchell et al., 2011;Baran et al., 2014;Baran, 2012).Recent parameterizations have bypassed the need for D e by coupling the optical parameterization directly to the GGM prognostic variable IWC (Baran et al., 2014).In order to make the microphysical and optical schemes physically consistent, it has also been argued that particles used in the optical parameterization should adhere to the same mass-dimensional and area ratiodimensional power laws as assumed in the microphysical scheme (Baran, 2012;Baran et al., 2014).
The parameterization presented in this paper (referred to as hex_cav from this point forward) tests the theory that ice clouds can be represented entirely by a single particle habit, as long as the mass-dimensional and area ratio-dimensional relationships are consistent with the those assumed in the microphysical scheme.As asymmetry parameters are typically over-predicted by using pristine particle models, the particle chosen for this parameterization is a stepped hollow column, as observed in recent laboratory studies (Smith et al., 2015(Smith et al., , 2016)).This particle was chosen because it yields lower asymmetry parameters than it's solid counterpart (Smith et al., 2015).By varying the aspect ratio (i.e.ratio of length to radius), as a function of maximum dimension (defined as the distance between opposite corners of the two basal facets, see Fig. 2), the hollow column was fitted to observed ranges of mass and area-ratio relationships.By doing so, the particles satisfy the same power laws that are assumed in the cloud physics scheme of the Met Office 6.0 configuration.The construction of the particle model is discussed in Section 2.1.The single scattering properties are determined using various scattering models for 26 particle sizes across 54 wavelengths (given in Appendix B), in the shortwave only.Single scattering properties were calculated with and without the use of the distortion parameter (which is used to simulate particle roughness) and therefore hex_cav is split into two parameterizations: hex_cav1 (without distortion) and hex_cav2 (with distortion), this is further discussed in Section 2.3.Bulk optical properties were then found by integrating the single scattering properties over 28 particle size distributions (PSDs) (Section 2.4).The bulk properties were used in the GA6 configuration of the Met Office Unified Model along with the current operational longwave parameterization (Section 2.5).The GCM is used to simulate the annual twenty year shortwave (SW) fluxes at the top of the atmosphere (TOA) and the corresponding zonal mean temperatures and specific humidities.In total, four model runs are completed.Each of these runs assumes the same microphysics, but a different optical parameterization.These optical parameterizations are: the current operational model by Edwards et al. (2007), henceforth referred to as Edwards2007; a more recent optical parameterization by Baran et al. (2014), henceforth referred to as Baran2014; hex_cav1 and hex_cav2.Results from the hex_cav model runs are compared against results assuming the Edwards, 2007 & Baran, 2014 optical schemes, and further comparison is made to observations from CERES (Stephens et al., 2012;Loeb et al., 2009).The hex_cav2 predicted zonal mean temperatures and specific humidities are compared against the ERA-Interim re-analysis product (Dee et al., 2011).The parameterizations Edwards2007, Baran2014, hex_cav1 and hex_cav2 are summarised in Table 1.

Particle model
The hollow column used in this parameterization is based upon particles observed during laboratory experiments conducted in the Manchester Ice Cloud Chamber (MICC) (Smith et al., 2015).Hollow ice crystals have been observed in many lab and field studies (Walden et al., 2003;Heymsfield et al., 2002;Bailey and Hallett, 2009).Images from these experiments show that the hollow columns have cavities which are pyramidal in structure, which have been modelled in theoretical studies (Schmitt et al., 2006;Yang et al., 2008).When using a rigorous improved geometric approach, the general effect of these pyramidal cavities was to increase the asymmetry parameter (Yang et al., 2008).Laboratory experiments in the MICC found that ice crystals grown at −30 °C tended to have stepped hexagonal intrusions as seen in Fig. 1.There was little variation in the geometry of the cavities at this temperature, and no solid columns were observed.Similar structures can be seen from in-situ studies which catalogued photographs of ice crystals collected in-situ (Weickmann, 1949).
In order to create a particle model based on formvar replicas, similar to the one shown in Fig. 1, an optical microscope was used to take measurements of the crystal facets, averaged values were then used to create a particle model for use in scattering simulations (Smith et al., 2015).The construction of the particle model is shown in Fig. 2.
This particular particle model was chosen for the parameterization because modelled results show that the stepped hollow column causes a reduction in asymmetry parameter compared with a solid column of the same maximum dimension and aspect ratio, in contrast to the pyramidal hollow column which causes a general increase (Smith et al., 2015).Therefore, the hollow column model offers a way of obtaining smaller asymmetry parameters other than the use of the distortion parameter or by embedding air or aerosol inclusions within the volume of the ice.Fig. 3 shows phase functions and asymmetry parameters for the stepped hollow column model, calculated for a wavelength of 632 nm, using both Ray Tracing (Macke et al., 1996b) and Ray Tracing with Diffraction on Facets (Hesse et al., 2012).This latter model takes into account internal diffraction not accounted for by classical geometric optics (Hesse et al., 2009).Calculations from Ray Tracing and RTDF were conducted for a randomly oriented particle, based on 5 × 10 4 particle orientations, and 5 × 10 7 rays per orientation.181 angular bins were used for scattering angles between 0.25°and 179.75°.The hollow particle model was set up as shown in Fig. 2, with basal and prism facets measuring 50 and 100 μm respectively, and a hollowness, h, of 80%.From Fig. 3 it can be seen that Ray Tracing over predicts the halo peaks relative to RTDF but predicts the same g values as RTDF.However, in this paper we prefer to apply the most physically appropriate model, which is RTDF.
To utilize this model in the optical parameterization, the aspect ratio, α, is varied as a function of maximum particle dimension in order to fit within observed mass-dimensional and area ratio-dimensional relationships.In this paper, we define maximum particle dimension, D, as: where: dimension of the basal facet and the aspect ratio α, is defined as:

Area ratio-dimensional relationships
The area ratio, A r , is defined as the ratio of the particle's projected cross-sectional area to the area of a circumscribed circle having a diameter equal to the maximum dimension of the particle.
2.1.1.1.Observed A r (D) relationships.Area ratio is a shape sensitive parameter, and is therefore sensitive to particle habit.Consequently, observed A r (D) laws vary between cloud types.Relationships have been found for various cloud types including mixed-habit cirrus, mixed phase clouds and tropical anvils (Heymsfield and Miloshevich, 2003;McFarquhar et al., 2013;Field et al., 2008).For this parameterization, we fit the particles to A r (D) relationships observed in mixed habit cirrus (Heymsfield and Miloshevich, 2003).This data represents 10 profiles through midlatitude, continental and synoptic cirrus, acquired over 3 field experiments.The combined profile follows the relationship: where D is the maximum particle dimension in cm.This relationship is derived from observations in the size range 0.004-0.320cm, giving area ratios between 0.8 and 0.25, where the area ratio decreases with respect to maximum dimension.
Table 1 Summary of the main features of the optical parameterizations: Edwards2007 (current operational model) and Baran2014 (recent habit mixture model) and hex_cav, which is split into hex_cav1 (no distortion) and hex_cav2 (with distortion).

Parameterization Summary Edwards2007
• An effective dimension D e based scheme, where D e is a function of environmental temperature • Ice particles are represented by the 8 branched hexagonal aggregate with roughened surfaces (Yang and Liou, 1998) • Uses 54 in-situ derived size distributions, from the observational campaign CEPEX (McFarquhar and Heymsfield, 1996).The shattered ice crystal artifacts are not removed from the PSDs (see Section 2.4) Baran2014 • Coupled directly to the microphysics scheme, assumes the same mass-dimensional and area-dimensional relationships • Ice particles are represented by a weighted habit mixture model, using six elements: simple column, six-branched bullet rosette, and three-, five-, eight-, and ten-monomer aggregates Baran and Labonnote (2007) • Uses 28 in-situ derived PSDs (from various field campaigns) parameterized by Field et al. (2007), where shattering artifacts have been removed (see Section 2.4) • Bulk optical properties are parameterized in terms of IWC and wavelength hex_cav1 • Coupled directly to the microphysics scheme, assumes the same mass-dimensional and area-dimensional relationships • Ice particles are represented by a single particle habit: A hollow column with stepped internal cavities (Smith et al., 2015).
• Uses 28 in-situ derived PSDs (from various field campaigns) parameterized by Field et al. (2007), where shattering artifacts have been removed (see Section 2.4) • Bulk optical properties are parameterized in terms of IWC and wavelength hex_cav2 • This parameterization is the same as hex_cav1 but the hollow particle is treated as 'rough' by using distortion in the single scattering calculations

Mass-dimensional relationships
In addition to A r (D) relationships, the particles used in the parameterization must also adhere to observed mass-dimensional power laws.Cirrus ice crystals are observed to obey the following mass-dimensional relationship (Cotton et al., 2013): in the range DN 70 μm where: M(D) = mass of the ice particle, kg D = maximum dimension, m.
In the size range D ≤ 70 μm, ice particles were found to have a constant effective density, given by: where the effective density of ice, ρ ICE , is defined as the mass of the ice crystal, divided by the volume of a sphere with diameter equal to the maximum particle dimension D.
For the hollow particle model used, equations were derived to characterize the relationships between the aspect ratio, α, area ratio (A r ), mass (M) and effective density (ρ ICE ).These equations were fitted to the observed relationships as given in Eqs. ( 6), ( 7) and ( 8).Particles could not be fitted exactly to both mass and area ratio relationships, therefore weighted averages are taken in order to fit the particle models within observed ranges.A full derivation is given in Appendix A.
The chosen aspect ratios fit within observed area-ratio and effective density values in the range D N 70 μm as shown in Fig. 4. The maximum ice effective density that is achievable with the hollow particle model is 384.9 kg m −3 , which is below the observed range, and subsequently particles b70 μm cannot be fitted within observed values.

Single scattering calculations
The single scattering properties for each of the 26 particles were calculated using either T-Matrix, RT or RTDF for 54 wavelengths in the shortwave between 0.2 μm and 5 μm.The wavelengths and refractive indices can be found in Appendix B. The choice of scattering model is dependent upon the size parameter, x, defined as πD/λ, where D is the maximum dimension of the particle and λ is the wavelength (Hesse et al., 2012).The size parameter can be loosely defined as small (x ≤ 20), intermediate (20 ≤ x ≤ 60), or large (x ≥ 60) (Baran, 2004).
Since no scattering model is applicable across the entire range of size parameters found in cirrus, optical parameterization makes use of a range of models, from exact methods for small size parameters (such as T-Matrix), geometric methods for larger size parameters (such as ray tracing) and improved geometric methods for intermediate sizes (such as RTDF).Since the stepped hollow model contains many small facets, which vary in size with aspect ratio, the applicable limits of each scattering model were not well defined.These limits were found by comparing phase function outputs from the three scattering models.At transitional sizes (sizes between small, intermediate and large size parameters), the phase functions were found to be largely similar but below/above these, they were found to deviate.Therefore the limits were defined where the different scattering models showed good agreement.By doing so, the scattering model for each particle size and wavelength was decided on a case by case basis.At smaller size parameters, the hollow hexagonal model was not used because of the very small facets, and therefore solid hexagonal prisms were used.In this case, differences between RT and RTDF were minimal and so RT was used for the small, solid particles.At even smaller sizes, where RT became none-applicable, T-Matrix for solid hexagonal columns was used (Havemann and Baran, 2004).A chart of the chosen models with respect to particle size and wavelength can be found in Appendix C.
For RT and RTDF, each simulation used 5 × 10 4 particle orientations and 5 × 10 7 incident rays.For each of the 26 particles, the single scattering properties (asymmetry parameter, single scattering albedo, extinction cross section and scattering cross section) were calculated for 54 wavelengths in the short-wave ranging from 0.2 μm to 5.0 μm, using complex refractive indices from Warren and Brandt (2008).These calculations form the basis of the hex_cav1 parameterization.
In order to diminish the 22°halo, the simulations were also done using the distortion parameter.Distortion values of 0.1, 0.2, 0.3 and 0.4 were tested for an example column of maximum dimension 100 μm, aspect ratio 2 and wavelength 632 nm (Fig. 5).A distortion value of 0.4 was found to completely remove the halo feature and therefore the single scattering calculations (as done in hex_cav1) were repeated using a distortion of 0.4, forming the basis of the hex_cav2 parameterization.This distortion value was chosen as the halo peak is completely removed, therefore producing a featureless phase function similar to those observed in situ (Labonnote et al., 2001).For hex_cav2, the large values of distortion used caused the outgoing ray paths to be significantly deviated.Near the particle edges, this bending of the outgoing ray path can cause outgoing rays to re-enter the space occupied by the crystal.This can cause errors where the ray is not correctly defined as being either in the scattering particle or the host medium, and the particle can no longer be considered a closed system.This issue limits the applicable size range of RTDF, therefore the applicable size range of RTDF varies between hex_cav1 and hex_cav2 (see Appendix C).

Bulk scattering properties
In order to calculate the bulk scattering properties, we use 28 PSDs as parameterized in Field et al. (2007), referred to as Field2007 from this point forward.The Field2007 parameterization is based on in-situ measurements from the Tropical Rainfall Measuring Mission/Kwajelein Experiment (TRMM/KWAJEX), Cirrus Regional Study of Tropical Anvils and Cirrus Layers-Florida Area Cirrus Experiment (CRYSTAL-FACE) and the First International Satellite Cloud Climatology Project Research Experiment (FIRE).Together, these field campaigns include more than 10,000 measured PSDs for tropical anvils and midlatitude stratiform cloud, covering temperatures from 0 °C to −60 °C.Field2007 improves upon earlier parameterizations as it covers a larger and therefore more representative temperature range.Furthermore, the Field2007 parameterization filters out shattered particles by analysis of ice crystal interarrival times, thus reducing the bias caused by shattering artifacts which is known to have affected historic PSDs (Field et al., 2006).Generally, bulk optical properties are related to the microphysical scheme Fig. 3. Modelled phase functions of randomly oriented solid and hollow hexagonal columns using Ray Tracing (top) and RTDF (bottom).'Hollow 1' is a particle with a pyramidal cavity, whereas 'hollow 2' is a particle with a stepped internal cavity as shown in Fig. 2. through the use of the diagnosed variable, D e , as discussed in Section 1. Instead, we directly couple the bulk optical properties to the prognostic variable IWC.
In order to calculate the bulk scattering properties for each of these PSDs, firstly, the single scattering properties are interpolated onto size bins in each PSD, where the number of size bins in each PSD was 500 and these ranged in size between about 0.4 μm to 28,000 μm.The single scattering properties at each bin size were then integrated over the PSD, thus finding a weighted average of each property.The scattering and extinction cross sections (β sca and β ext , respectively) are weighted by the mass of cloudy air per unit volume (in units of kg m − 3 ).This yields the mass scattering and mass extinction coefficients (K sca and K ext , respectively), which describe the scattering and extinction cross sections per unit mass of cloudy air.The bulk asymmetry parameter was then found by weighting with respect to scattering cross section.These weightings give bulk optical values consistent with the Met Office Unified Model definitions.In the Met Office global model the bulk scattering and extinction coefficients are represented by the mass scattering and extinction coefficients per unit mass of cloudy air, and so the units are m 2 kg −1 .The values of each of the bulk scattering properties are plotted as a function of wavelength for each of the 28 PSDs, these can be seen in Appendix D.1 and Appendix D.2 for parameterizations hex_cav1 and hex_cav2, respectively.These are used to find parameterized fits for g, K ext , K sca and ω 0 for the 6 shortwave bands used in the Met Office configuration 6 atmosphere only model.A table of these fits can be found in Appendix E.

Implementation in the GCM
The hex_cav1 and hex_cav2 short-wave optical parameterizations are used, assuming the current Edwards2007 parameterization applied to the long-wave.In the climate model runs that follow, the Ed-wards2007 parameterization is used as the control model, and comparison is also made with the more recent Baran2014 parameterization.For all four model runs presented here (using the optical parameterizations: Edwards2007, Baran2014, hex_cav1 & hex_cav2), the same microphysical scheme is used, based on PSDs from Field2007, fall speeds parameterized by Furtado et al. (2014) and mass-dimensional relationships derived by Cotton et al. (2013).This is done so that any changes in the short-wave is entirely attributable to the parameterization presented in this paper.The bulk scattering properties are implemented into the GA6 configuration of the Met Office atmosphere only unified model.This is used to simulate the annual twenty year short-wave fluxes (fluxes averaged over 20 one year intervals) at the top of the atmosphere and the corresponding zonal mean temperatures and specific humidities.Details of the GA4 configuration can be found in Walters et al. (2014), the subsequent GA5 and GA6 configurations include a new dynamical core, described by Wood et al. (2014), and the new spectral files for GA6 can be described in Section 3 of Manners et al. (2015).

Comparison of bulk scattering properties
In this section, we compare the bulk scattering properties predicted by the hex_cav parameterizations with the Edwards2007 parameterization and the recent Baran2014 parameterization.The Edwards2007 model is an effective dimension based scheme, with D e as a function of temperature.Both hex_cav models and Baran2014 have no temperature dependence so instead we compare bulk scattering properties at set temperatures of 200 K, 230 K and 270 K with respect to ice mass mixing ratio between 1.0 × 10 − 7 and 1.0 × 10 −3 kg kg −1 as these ranges are found in the GA6 model.We compare results for shortwave band 1 and band 5 (0.2-0.32 μm and 1.19-2.38μm, respectively).These particular bands are chosen for comparison due to their contrasting absorption properties, therefore we expect results to differ largely between the weakly absorbing band 1 and the strongly absorbing band 5.

Mass extinction coefficient
Figs. 6, 7 and 8 show the mass extinction coefficient for hex_cav1, Baran2014, and Edwards2007.From these figures we see that the hex_cav model has the lowest extinction at all fixed values of temperature for short-wave band 1. Results from short-wave band 5 and from hex_cav2 were found to be similar, these are not shown for reasons of brevity.

Asymmetry parameter
Figs. 9 and 10 show asymmetry parameters for Edwards2007, Baran2014, hex_cav1 and hex_cav2 at T = 200 K for short-wave bands 1 and 5, respectively.We see that for the Edwards2007 control model, the asymmetry parameter is invariant with respect to ice mass mixing ratio as the aspect ratio of the particle does not change with particle size.However, the asymmetry values for Edwards2007 do vary slightly with temperature, whilst Baran2014 and the hex_cav parameterizations remain constant, as a function of temperature.For the more absorbing case (Fig. 10), we see that the hex_cav2 parameterization is closest to the fully randomized Baran2014 model.At this band, the asymmetry parameters predicted by Edwards2007 change significantly as a function of temperature due to the larger (and therefore more absorbing) ice crystals, but still remain invariant with respect to ice mass mixing ratio.

Single scattering albedo
Figs. 11, 12 and 13 show single scattering albedos for Edwards2007, Baran2014, hex_cav1 and hex_cav2 at temperatures of 200, 230 and 270 K respectively.At short-wave band 1, ω 0 ≈1, so instead we concentrate on the more absorbing short-wave band 5.The Edwards2007 ω 0 values are larger than both the Baran2014 and the hex_cav models.Both hex_cav values of ω 0 increase with ice mass mixing ratio due to the decrease in volume absorption.

GCM simulations
This section shows results from hex_cav1 and hex_cav 2 from the GA6 configuration of the Met Office unified model, compared with the Edwards2007 control model and CERES observations.
Figs. 14 and 15 show the twenty-year averaged annual downwelling and up-welling short-wave flux at top-of-atmosphere (TOA) as predicted by the hex_cav2 parameterization, respectively.The TOA downwelling short-wave flux is defined as the short-wave irradiance that reaches the Earth's surface from the model top of atmosphere (80 km).Differences between the two parameterizations in predicting the downwelling and upwelling fluxes at top-of-atmosphere can be seen in Figs.14b and 15b, respectively.Results from hex_cav1 were found to be similar and are therefore not shown for reasons of brevity.We see that differences between the hex_cav parameterization and Edwards2007 are highest around the tropics and the southern ocean.When compared with observations, we see that the control model generally under predicts down-welling flux, except in the southern ocean where it tends to be over predicted.On the contrary, the hex_cav2 parameterization tends to over predict down-welling flux when compared with observations, particularly in the tropics and southern ocean.However, there are regional improvements to be seen in the hex_cav prediction of TOA fluxes.Improvements can be seen over the Atlantic, and parts of the Pacific Ocean.Converse to this, Fig. 15c and d shows that the Edwards2007 and hex_cav2 parameterizations generally over predict and under predict the upwelling short-wave flux at TOA, respectively.
Fig. 16 shows the zonal mean temperatures predicted by the Ed-wards2007 and hex_cav2 parameterizations.From this it can be seen that the under prediction of reflected short-wave flux in the tropics (as seen in Figs. 14 and 15) leads to the warming of the tropical troposphere by about 1 K and cooling of the stratosphere by about 0.5 K.Over the North pole this results in a significant reduction in the warming relative to the control model, and over the South Pole there is a reduction in the cooling relative to the control.This warming over the tropics leads to an increase in the specific humidity relative to the  control, reducing the dry bias in the upper tropical troposphere (shown in Fig. 17).
Overall, the predictions of TOA short wave flux, zonal mean temperature and zonal mean specific humidity differ from observations more so than the current operational model.However, there are regional improvements that can be seen.For upwelling and downwelling flux, improvements on the current model are seen over the North Atlantic, Indian and much of the Pacific Ocean.For both the zonal mean temperature and zonal mean specific humidity, although biases in the tropical tropopause are increased, biases in the polar regions are decreased.Many of these differences may be explained by the largely different mass extinction values predicted by each of the parameterizations.As seen in Figs. 6, 7 and 8, the hex_cav parameterizations consistently predicted lower mass extinction values compared with Edwards2007 and Baran2014 (this is due to large aspect ratios needed to fit the particle to within observed relationships).These regional improvements may correspond to areas containing smaller particles, and therefore the larger particles (with large aspect ratios and low mass extinction) have had little influence on the region.Alternatively, it is known that the orientation of ice crystals in the atmosphere is not fully randomized (Yang et al., 2003).Factors such as gravitational sedimentation can cause preferential orientation of ice crystals, particularly for α ≪ 1 or α ≫ 1 (Hashino et al., 2014).In convective systems, electric fields can also cause preferential alignment (Foster and Hallett, 2008).In these cases, the projected area, and hence mass extinction of the ice crystals would be larger than for randomly oriented particles.Therefore the assumption of random orientation in the GCM may lead to larger biases in regions where orientation is not negligible.In Figs. 15 and 14 we see that the largest biases in hex_cav2 occur in the southern ocean and over tropical Asia.The derivation of the area-ratio dimensional relationship is based on data collected in situ via cloud probes.These data are also orientation dependent, and may be affected by particle orientation in the sample volume.Data from a variety of cirrus are used to generate a globally averaged relationship, which may be more representative of certain regions compared with others.

Conclusions
It has been argued that, to properly model the optical properties of cirrus ice clouds, the individual particle models used must adhere to observed mass-dimensional and area ratio dimensional relationships.By maintaining these relationships, the optical parameterization not only becomes physically consistent with the microphysics scheme (in which these relationships are assumed), but should ensure that the predicted ice mass and projected areas are accurate.In this paper, we have investigated the ability of a single particle geometry (in this case a hollow hexagonal column) to fit within these constraints.
In order to fit a hexagonal prism (whether solid or hollow) to observed area ratios, preferentially oriented particles had to be assumed, as described in Appendix A. This resulted in very large aspect ratios of up to 20 being assumed.Despite the assumption of preferential orientation for the selection of particle aspect ratio, single scattering properties were found using randomly oriented particles, as required by the GCM.In comparison to preferential orientation, the projected area in random orientation is reduced, which is particularly significant for larger aspect ratios.Therefore the use of such large, and unrealistic aspect ratios caused much lower predictions of mass extinction coefficient compared with other parameterizations, as shown in Figs. 6, 7 and 8.Although the use of the hollow particle model reduced the asymmetry parameter (causing clouds to become brighter) compared with an equivalent solid model, the effect of reducing the asymmetry parameter was cancelled out by the very small mass extinction values (causing clouds to become darker).Therefore, more short-wave radiation will be transmitted to the Earth, which is evident in Figs. 14 and 15, where we can see that hex_cav2 under-predicts the reflected short-wave radiation at TOA.The effect of this is to warm the tropical troposphere (Fig. 16).Generally speaking, the hex_cav predictions of TOA SW flux, zonal mean temperature and specific humidity differed from observation more so than Edwards2007 and Baran2014, however, some regional improvements  were seen.Areas in the North Atlantic, Indian and much of the Pacific Ocean showed improvements for upwelling and downwelling flux, and temperature and humidity biases in polar regions were decreased.This highlights the sensitivity of the climate to small changes in the microphysical properties of ice clouds and it is therefore pivotal to construct parameterizations that are microphysically consistent.Furthermore, it is  The results suggest that a single hexagonal prism cannot be used to approximate ice of all sizes.As seen in Fig. 4, the particle could not be fitted to the high values of ice effective density as observed for smaller particles.In order to conserve the ice mass for such particles, quasispherical particles might be a better approximation (McFarquhar et al., 2002), allowing for higher area ratios to be achieved.As for large  particles, the use of very elongated hexagonal prisms leads to underpredictions in orientation-averaged projected area.These particles may be better represented by spatial aggregate models, which can achieve the low values of area ratio required, but are less sensitive to particle orientation.The stepped hollow particle has been observed in field studies (Weickmann, 1949) and in laboratory studies (Smith et al., 2015(Smith et al., , 2016)), where clouds below − 25°were found to contain almost exclusively stepped hollow particles.Therefore, it is likely that such structures occur frequently in cirrus ice cloud, however the internal structures are often unseen with current measurement techniques.Due to their predominance in these studies, in conjunction with their particular optical properties (Smith et al., 2015(Smith et al., , 2016)), these stepped hollow columns should be incorporated into future habit mixture models.
In the current habit mixture models, perturbations from the pristine form are often treated with the use of distortion as a proxy for surface roughness, or by the use of inclusions.Whilst these methods may yield values of scattered intensity close to observations, they may overlook other properties of the scattered light.Measurements from the A-train now provide us with polarisation measurements from ice cloud, and it has been shown that particles with similar phase functions may differ significantly with respect to degree of linear polarisation (Mishchenko et al., 2007;Baran and Labonnote, 2006;Stephens et al., 2002).It has also been shown in laboratory studies that hollow particles are more weakly depolarising compared with solid crystals (Schnaiter et al., 2007;Smith et al., 2016).In this case, roughness proxies may not be representative of various micro-scale features such as cavities, inclusions, and real surface roughness.

A.1.1. Randomly oriented particles
For a solid convex particle, the average projected cross section is given by S/4, where S is the particle surface area.Although the hollow particle is concave, the projected area is not influenced by the concavities, and therefore the same equation can be applied.For a randomly oriented hexagonal prism, the average projected area is given by: and the area ratio is given by:

. Preferentially oriented particles
In reality, elongated ice particles tend to fall preferentially with their largest projection perpendicular to the direction of propagation (Platt, 1978;Chepfer et al., 1999), although vertically aligned prism facets have also been observed (Westbrook, 2011).As such, columns fall preferentially with their prism facet parallel to the ground, whereas plates fall preferentially with their basal facet parallel to the ground.
For column-oriented particles, the projected cross section is given by: and the area ratio is given by: For plate-oriented particles, the projected cross section is given by: and the area ratio is given by: The area ratios for preferential and random orientations are plotted against aspect ratio in Fig. A.19.
By assuming a randomly oriented particle, the maximum area ratio for a hexagonal column is 0.7271.However, the observed A r (D) relationships exceed this, with a maximum value of 0.8.In order to achieve this value, we must assume oriented plates.Therefore the two orientation specific relationships are used rather than the randomly oriented one.
If we extend the A r (D) relation to cover the full size range used in this parameterization (0.4-28,127 μm) we get a range of values of A r  from 2.8 to 0.1360.Physically, the area ratio for a hexagonal column cannot exceed 0.8270, and therefore this observational relationship cannot be extrapolated to smaller particles.As D tend to infinity, A r tends asymptotically towards 0. So in theory, the relationship can be extrapolated to larger sizes.

A.2. Fitting the particles to M(D) relationships
For the hollow column used in this parameterization, the mass is given by: where: M = particle mass, kg ρ = density of ice, kg m −3 h = hollowness described as the combined length of both cavities, expressed as a percentage of the length of prism facet, p.
Varying the hollowness caused little difference in the particle mass, and therefore a constant hollowness of 80% was assumed, as commonly observed in cloud chamber investigations (Smith et al., 2015).In order to fit the hollow particle to observed mass-dimensional relationships, we equate Eqs. ( 7) and (A.9) to get: The relationship between α and D is approximated by a 10th degree polynomial: where c n are the polynomial coefficients, given in Table A.2. Eqs.(A.7) and (A.8) relate the aspect ratio and the maximum dimension of the hollow column in order to adhere to observed area ratio-dimensional relationships, whilst Eq. (A.11) relates aspect ratio and maximum dimension in order to obey observed massdimensional power laws.These equations are not in agreement and therefore the aspect ratio cannot be fitted exactly to both observed relationships.Instead, we take a weighted average in order to fit the values within observed ranges.It was found that a 50:50 weighting gave the best agreement for sizes N70 μm.In the size range 40-70 μm, a 65:35 weighting was used (M(D):A r (D)).These weightings were chosen as they produced the most amount of crystals within the observed ranges of M(D) and A r (D).For particles below 40 μm, there is no established A r (D) relationship and so particles are fitted using only the M(D) relationship.
These equations were used to find the aspect ratios of 26 particles ranging in size from 0.4 to 28127 μm, given in Table A.3.Parameterized fits of K ext , K sca and g for 6 short-wave bands for the hex_cav1 parameterization; where q i is the ice mass mixing ratio in kg per kg.

Fig. 1 .
Fig. 1.Formvar replicas of typical habits of ice crystals found at −30°in the Manchester Ice Cloud Chamber.Viewed from the prism face (a), and the basal face (b).

Fig. 2 .
Fig. 2. Particle model construction based on average measurements from formvar replicas.Panel a shows a plan view of the particle from the basal facet, and panel b shows a cross sectional view of the particle, taken parallel to the prism facet where b and p are the lengths of the basal and prism facets respectively.d is the depth of each cavity, and h is the total combined length of both cavities expressed as a percentage of p.The maximum dimension, D, is given by D ¼

Fig. 4 .
Fig.4.The top graph shows the chosen aspect ratios used for this parameterization.The second graph shows the corresponding area ratios of the chosen particles, and the shaded region shows the observed range.The bottom graph shows the corresponding ice effective density of the chosen particles, and the shaded region shows observed ranges.

Fig. 5 .
Fig.5.Phase functions of a randomly oriented stepped hollow column, aspect ratio 2, with varying values of distortion.Simulations were run using RTDF, with a wavelength of 635 nm.Halo features are evident for distortion values of 0.1, 0.2 and 0.3, but not for 0.4.Therefore a distortion value of 0.4 was used to remove these features.

Fig. 7 .
Fig. 7. Mass extinction plotted against ice mass mixing ratio as predicted by hex_cav1, Edwards2007 and Baran2014 at a temperature of 230 K.

Fig. 14 .
Fig. 14.Annual short-wave down-welling flux at the top of atmosphere.Clockwise from top left: predictions from hex_cav2, hex_cav2 minus control model, hex_cav2 minus observations, control model minus observations.

Fig. 15 .
Fig. 15.Annual short-wave up-welling flux at the top of atmosphere.Clockwise from top left: predictions from hex_cav2, hex_cav2 minus control model, hex_cav2 minus observations, control model minus observations.

Fig. 17 .
Fig. 17.Zonal mean specific humidity predicted by the hex_cav2 parameterization.Clockwise from top left: predictions from hex_cav2, hex_cav2 minus control model, hex_cav2 minus observations, control model minus observations.

Fig
Fig. A.18. Projected area of a hexagonal prism when oriented like a column (left) and a plate (right).Shaded areas represent projected cross sections.'Hollowness', in the form of basal cavities, does not affect the particles projected area, therefore cavities are omitted from the diagram.

Fig
Fig. A.19. Area ratio plotted against aspect for randomly oriented and preferentially oriented hexagonal columns.

Fig. D. 26 .
Fig. D.26.Bulk K ext values calculated using the hex_cav2 parameterization.Each trace represents a different PSD from Field2007.

Fig. D. 27 .
Fig. D.27.Bulk K sca values calculated using the hex_cav2 parameterization.Each trace represents a different PSD from Field2007.
Table A.2 Coefficients of D n for Eq.(A.11).

Table A .
3 Aspect ratios and maximum dimensions of the 26 particles used in the hex_cav parameterizations.For particles ≤ 80 μm, plate orientation is assumed, for particles ≥80 μm, column orientation is assumed.