On the relationship between the scattering phase function of cirrus and the atmospheric state

This is the first paper to investigate the relationship between the shape of the scattering phase function of cirrus and the relative humidity with respect to ice (RHi), using space-based solar radiometric angle-dependent measurements. The relationship between RHi and the complexity of ice crystals has been previously studied using data from aircraft field campaigns and laboratory cloud chambers. However, to the best of our knowledge, there have been no studies to date that explore this relationship through the use of remotely sensed space-based angle-dependent solar radiometric measurements. In this paper, one case study of semitransparent cirrus, which occurred on 25 January 2010 off the north-east coast of Scotland, is used to explore the possibility of such a relationship. Moreover, for the first time, RHi fields predicted by a high-resolution numerical weather prediction (NWP) model are combined with satellite retrievals of ice crystal complexity. The NWP model was initialised at midnight, on 25 January 2010, and the mid-latitude RHi field was extracted from the NWP model at 13:00 UTC. At about the same time, there was a PARASOL (Polarization and Anisotropy of Reflectance for Atmospheric science coupled with Observations from a Lidar) overpass, and the PARASOL swath covered the NWP-model-predicted RHi field. The cirrus case was located over Scotland and the North Sea. From the satellite channel based at 0.865 μm, the directionally averaged and directional spherical albedos were retrieved between the scattering angles of about 80 and 130. An ensemble model of cirrus ice crystals is used to predict phase functions that vary between phase functions that exhibit optical features (referred to as pristine) and featureless phase functions. For each of the PARASOL pixels, the phase function that best minimised differences between the spherical albedos was selected. This paper reports, for this one case study, an association between the most featureless phase function model and the highest values of NWP-predicted RHi (i.e. when RHi > 1.0). For pixels associated with NWPmodel-predicted RHi < 1, it was impossible to generally discriminate between phase function models at the 5 % significance level. It is also shown that the NWP model prediction of the vertical profile of RHi is in good agreement with dropsonde, in situ measurements and independent aircraft-based physical retrievals of RHi . Furthermore, the NWP model prediction of the cirrus cloud-top height and its vertical extent is also found to be in good agreement with aircraft-based lidar measurements.


Introduction
Cirrus or pure ice crystal cloud usually forms at temperatures of less than about −40 • C, and at altitudes greater than about 6 km (Wylie et al., 1999;Baran 2012;Guignard et al., 2012).The extent to which ice crystals can grow and form complex shapes is dependent on the environmental temperature, pressure and RH i (Marshall and Langleben, 1954;Nakaya, 1954;Hallett and Mason, 1958;Mason, 1971;Heymsfield, 1977;Liou, 1986;Lynch, 2002;Bailey and Halett, 2004;Bailey and Hallett, 2009;Pfalzgraff et al., 2010;Ulanowski et al., 2013).In cirrus, the supersaturations with respect to ice can range from about 150 to −50 % (Krämer et al., 2009).However, more recent studies of RH i in mid-latitude cirrus report values of about 60 to 120 % (Gayet et al., 2011;Ulanowski et al., 2013), with the latter values being the more typical.With such a range of possible supersaturation values reported by previous authors, the range of ice crystal complexity within the same cirrus is likely to be significant (Bailey and Hallett, 2009).In this paper, ice crystal complexity can mean polycrystals, which may be compact or highly irregular, and ice aggregates, and these ice aggregates may also have low area ratios (i.e. the ratio between projected area of the ice crystal to the projected area of the circumscribing circle of the same maximum dimension as the ice crystal).The monomers that make up the polycrystal or aggregate may also be surfaceroughened on their facets and/or contain air cavities within their volumes.
The role of ice supersaturation in forming complex ice crystals has most recently been studied by Bailey and Hallett (2009), Pfalzgraff et al. (2010), Ulanowski et al. (2013) and Magee et al. (2014).In the cloud chamber study by Pfalzgraff et al. (2010), it was reported that surface roughness was at its greatest when supersaturations were near zero.Moreover, Walden et al. (2003) observed surface roughness on precipitating ice crystals under conditions of near-zero supersaturation at the South Pole.Other laboratory studies by Bacon et al. (2003), Malkin et al. (2012) and Ulanowski et al. (2013) show that ice crystals, under ice-supersaturated conditions, can become surface-roughened through the development of prismatic grooves or dislocations on the ice crystal surface.However, as pointed out by Bacon et al. (2003) and others, the temperature and RH i variables do not uniquely determine the number of monomers that make up the polycrystal or surface roughness.This is because ice crystal complexity and surface roughness may also depend on how the ice crystals were initiated, and thus may depend on the chemical composition of the initiating ice nuclei.A recent paper by Ulanowski et al. (2013) reported for a few cases of mid-latitude cirrus formed in oceanic air that slightly higher values of ice crystal complexity were found than was the case for mid-latitude cirrus formed in continental air (i.e. a polluted air mass).However, in the same study, no correlation was found between ice crystal complexity and instantaneous measurements of RH i .As pointed out by Ulanowski et al. (2013), this lack of correlation could be due to the instantaneous measurements being obtained at a single point in time, whereas the ice crystals, on which the measurements were based, may have gone through different histories of supersaturation, and so for each measurement the history of exposure to supersaturation can never be known.On the other hand, controlled laboratory studies of ice crystal analogues by Ulanowski et al. (2013) show that, under high levels of ice supersaturation, the ice crystals formed can be very complex relative to the regular ice crystals grown under conditions of low ice supersaturation.This latter laboratory study of Ulanowski et al. (2013) is consistent with the findings of Bailey and Hallett (2009).
Theoretical light-scattering studies by (Schmitt and Heymsfield 2010;Macke et al., 1996a;Macke et al., 1996b;Yang and Liou, 1998;Yang et al., 2008;van Diedenhoven, 2014b) have shown that the processes of surface roughness and air inclusions within ice crystals can profoundly alter their scattering phase functions.As surface roughness increases, the 22 and 46 • halos are reduced or completely removed, resulting in featureless phase functions with a high degree of side scattering.This high degree of side scattering results in surface-roughened ice crystals having lower asymmetry parameter values relative to their smooth counterparts.The asymmetry parameter is one of the parameters of importance in climate models, since it affects how much incident solar irradiance is reflected back to space (Stephens and Webster 1981;Yang and Liou 1998;Yang et al., 2008;Ulanowski et al., 2006;Baran 2012).Therefore, it is important to constrain this parameter through observation using a variety of instruments such as those used in the studies by (Gayet et al., 2002(Gayet et al., , 2011;;Field et al., 2003;Garrett et al., 2003;Mauno et al., 2011;Ulanowski et al., 2013;van Diedenhoven et al., 2014a;Cole et al., 2014).Other methods of removing or diminishing halos involve introducing air concavities from the basal ends of hexagonal ice crystals, or embedding spherical air bubbles within the ice crystal volume.The former method removes the 46 • halo and reduces the 22 • halo (Macke et al., 1996b;Yang et al., 2008), and the latter method produces featureless phase functions through multiple scattering between spherical air inclusions (Labonnote et al., 2001; Baran and Labonnote, 2007;Xie et al., 2009).Although recent cloud chamber and theoretical ray-tracing studies by Neshyba et al. (2013) and Shcherbakov (2013), respectively, have shown that surface roughness may not necessarily completely remove the 22 • halo, it is as yet unclear as to whether the results obtained in the laboratory are scalable to the real atmosphere.Indeed, in situ studies on the occurrence of the 22 • halo show that it is a rare event (Field et al., 2003;Gayet et al., 2011;Ulanowski et al., 2013).Clearly, further laboratory studies of ice crystals are needed, which combines angular scattering measurements at visible wavelengths with a detailed analysis of ice crystal habit, surface roughness and the degree of concavity, all obtained, as functions of time.The dimension of time is important to include, as this would be a useful constraint to apply to theoretical studies of ice crystal growth and complexity (Barrett et al., 2012).
Radiometric angle-dependent observations of the transmitted and reflected intensities from cirrus tend to suggest that featureless phase functions best represent those measurements obtained from below and/or above the cloud (Foot, 1988;Baran et al., 1999Baran et al., , 2001;;Doutriaux-Boucher et al., 2000;Jourdan et al., 2003;Baran, 2012;Cole et al., 2013;Cole et al., 2014;and references contained therein).Aircraftbased instruments such as the Polar Nephelometer (PN), which is described by Gayet et al. (1997), have been used to measure the angular scattered intensity of naturally occurring single ice crystals at scattering angles between about 15 and 162 • and at the wavelength of 0.80 µm.Clearly, the PNmeasured polar angle range encompasses the halo regions of 22 and 46 • .Therefore, obtaining the ratio of the ice crystal scattered energy at the polar angle of 22 • to that at 18.5 • (the Atmos.Chem.Phys., 15, 1105-1127, 2015 www.atmos-chem-phys.net/15/1105/2015/latter being an angle at which no halo is formed by prismatic ice crystals) would be a quantitative measure of the presence of the 22 • halo.Halo ratios greater than 1 are likely to be associated with pristine ice crystals, whilst halo ratios less than 1 are likely to be associated with irregular ice crystals.Using a mid-latitude cirrus case, Gayet et al. (2011) used the halo ratio to relate the occurrence of halos to instantaneous measurements of RH i .The study found that, at a temperature of −55 • C at the trailing edge of the cirrus band, the halo ratio is < 1, but at a temperature of −27 • C at the leading edge of the cirrus-band, the halo ratio is > 1.The study did not find any systematic evidence for a relationship between the halo ratio and ice supersaturation, a finding that is consistent with Ulanowski et al. (2013).However, Gayet et al. (2011) did find that halo ratios > 1 were more likely to be found at supersaturation values of around 100 %, and no halo ratios > 1 were found at the highest supersaturation values, which approached 120 %.Moreover, in the recent study by Ulanowki et al. ( 2013), a negative correlation is reported between the occurrence of halos and estimated ice crystal complexity.The measure of ice crystal complexity was derived from in situ observations of spatial light-scattering patterns from single particles obtained in several cases of mid-latitude cirrus.The in situ findings of Gayet et al. (2011) and laboratory studies of Ulanowski et al. (2013) are consistent with previous studies (i.e.Bailey and Hallett, 2009, and references therein), which tend to show that more complex ice crystals are associated with relatively high values of RH i .The relationship between the scattering properties of atmospheric ice and the physical state in which the ice resides is important to explore, as this may lead to an improvement in the parameterisation of ice optical properties in climate models.Such an improvement can only come about through a deeper understanding of how the growth of ice crystal complexity is related to the atmospheric state.This relationship could then be used to predict appropriate ice-scattering properties for some given atmospheric state, rather than assuming the same scattering properties for all states that are found in a climate model,which is what is generally done in presentday studies.The most recent report of the Intergovernmental Panel on Climate Change (IPCC, 2013) stated that the coupling between clouds and the atmosphere was one of the largest uncertainties in predicting climate change.This uncertainty may well be reduced if appropriate parameterisations could be found between ice crystal scattering properties and the atmospheric state.
In this paper, for one case of mid-latitude cirrus, the relationship between the scattering phase function and RH i is further studied by combining with space-based multiangle spectral albedo retrievals a numerical weather prediction (NWP) model prediction of the RH i field.The paper is split into the following sections.Section 2 describes the NWP model and the aircraft-based instruments used in this study.Section 3 describes the single-scattering properties of ice crystals on which the satellite retrievals are based, and a brief description of the radiative transfer model is also given.The retrieval methodology is described in Sect. 4 and results are discussed in Sect. 5. Section 6 presents the conclusions of this study.

The cirrus case, numerical weather prediction model and aircraft instrumentation
The conditions required for this paper are that the cirrus should be sufficiently optically thick to allow for discrimination between various randomisations of the ensemble model using PARASOL retrievals, the aircraft and satellite should be coincident, and there should be no underlying cloud or broken cloud fields.It is practically very difficult to obtain all these necessary conditions at the same time.The cirrus case occurred on 25 January 2010 off the north-east coast of Scotland, which is shown in Fig. 1.The figure shows a highresolution MODerate Imaging Spectroradiometer (MODIS) composite image (Platnick et al., 2003) of semi-transparent cirrus obtained at 13:30 UTC.The semi-transparent cirrus can be clearly seen around the north-east coast of Scotland, whilst further to the east, lower level water cloud underlying the cirrus can be seen.At about the same time as the image shown in Fig. 1 was taken, the FAAM (Facility for Airborne Atmospheric Measurements) BAe-146 aircraft was measuring the same high-cloud field.
The FAAM aircraft is an atmospheric research facility which is jointly owned by the Met Office and the National Environment Research Council (NERC).The cirrus case shown in Fig. 1 was sampled by the aircraft as part of the "Constrain" field programme (Cotton et al., 2013).In this paper, one case from the Constrain field programme is presented.There were several other Constrain cirrus cases, but these did not meet the conditions necessary for this paper.This is because the other cases were either optically too thick, there was no coincidence between PARASOL and the aircraft, or there was substantial underlying cloud.
For the case presented in this paper, the aircraft sampled the cirrus between the latitudes of about 58 and 59 • N, and between the longitudes of about 2.5 and 4.5 • W. The aircraft in situ instrumentation measured the temperatures of the cloud top and base to be about −55 and −30 • C, respectively.The aircraft began sampling the cloud at 11:55:06 UTC, and finished the sampling at 15:49:44 UTC.During this sampling time, the aircraft flew three straight and level runs above the cloud, each of about 10 min duration, commencing at 13:19:00 UTC, 13:27:42 UTC and 15:21:42 UTC, respectively.From the aircraft, a dropsonde (measures vertical profiles of temperature, pressure and relative humidity with respect to water) was released at 13:30:00 UTC.In this study, use is made of the aircraft data from the earlier two runs as well as the dropsonde measurements which were most closely related to the PARASOL overpass.Note also that there was an 8 min interval between the two earlier straight and level runs, during which time the aircraft was manoeuvring into position.In this paper, use is made of observations from four instruments deployed on the aircraft.The first two instruments were the active Leosphere ALS450 elastic backscatter lidar and the passive Airborne Research Interferometer Evaluation System (ARIES).
The nadir-pointing lidar operates at 0.355 µm with an integration time of 2 s and a vertical resolution of 1.5 m (Marenco et al., 2011).Further averaging of the signals has been done at post-processing, bringing the temporal resolution to 10 s (equivalent at aircraft science speed to a 1.5-2 km footprint) and the vertical resolution to 45 m.The volume extinction coefficient is computed from the lidar returns using the Fernald- Klett method described in Fernald (1984) and Klett (1985), assuming a lidar ratio of 20 sr.
The ARIES instrument is fully described in Wilson et al. (1999), but, briefly, it is a modified Bomem MR200 interferometer that measures infrared radiances between the wave numbers of 3030.303 and 500 cm −1 at a spectral resolution of 1 cm −1 .The interferometer is capable of multiple viewing geometries both up and down as well as across-track.The nadir-pointing ARIES and lidar data used here are from the straight and level runs above the cirrus.The other two instruments deployed on the aircraft were used to measure the in situ vertical profile of RH i , these were the General Eastern GE 1011B Chilled Mirror Hygrometer (GE) and the Fluorescence Water Vapour Sensor (FWVS) fast Lyman-alpha hygrometer (Keramitsoglou et al., 2002;Fahey et al., 2009).
The RH i field of the cirrus case, shown in Fig. 1, has been simulated using the Met Office high-resolution NWP model.
The NWP simulation of the RH i field was obtained using a high-resolution limited-area model nested inside a suite of coarser resolution models.The high-resolution domain had a horizontal grid spacing of approximately 1 km and received hourly lateral boundary conditions from a 4 km model on a larger domain.The 4 km model was nested inside a 12 km domain, which was in turn driven by an N216 global model forecast.The 1 km grid was centred on (58.60 • N, 6.45 • W) with 1024 points east-west and 744 points north-south and a zonal and meridional grid spacing of 0.0135 • .The model time step was 50 s and the vertical level set comprised 70 levels, with a grid spacing of approximately 250 m at the altitudes of interest for this study.
The model is non-hydrostatic and employs the semi-Lagrangian advection scheme.
In terms of model physics, the model is broadly comparable to the version of the Met Office operational UKV forecasting system that was used operationally until the autumn of 2011 (Lean et al., 2008).However, in an attempt to represent the simulated cloud system as well as possible, the following changes were made to the ice cloud microphysics.Firstly, the ice particle size distribution (PSD) parameterisation was changed so as to be consistent with the PARA-SOL radiative retrievals (see Sect. 3 for details).Secondly, the mass-diameter relation of the ice crystals was taken directly from the Constrain in situ measurements (Cotton et al., 2013), and is therefore a "best estimate" of this property for the simulated cloud system.For this paper, the NWP model was initialised at midnight on 25 January 2010, and the RH i forecast field was extracted from the model on the same day but at the 13:00 UTC time step.
At about the same time as the NWP model RH i field was extracted there was a PARASOL overpass at about 12:50 UTC.PARASOL has its central channels located at 0.443, 0.490, 0.565, 0.670, 0.763, 0.765, 0.865, 0.910 and 1.02 µm.In this study, use is made of the channel located at 0.865 µm, due to the sea surface at this wavelength being almost black.PARASOL can view the same nadir pixel at up to 14 viewing directions, at scattering angles between approximately 70 and 180 • , and each pixel has a resolution of 5.3 km × 6.3 km.The range of scattering angles sampled by PARASOL depends on the Sun-satellite geometry, latitude of the pixel, and the position of the pixel on the satellite track (i.e.east or west).Given the latitude and time of the year of the case considered here, the range of scattering angles viewed by PARASOL was between about 80 and 130 • , and the total number of viewing directions for each pixel was between 7 and 8.The solar zenith angle at the time of the PARASOL overpass was 75 • and the solar azimuth angle was 187 • .The PARASOL analysis is performed on a pixelby-pixel basis.Moreover, it should also be noted here that the NWP model field of RH i is averaged over approximately the same area as each of the PARASOL pixels.In the next section, the ice crystal model, single-scattering properties and radiative transfer models are briefly described and defined.

Ice crystal model and the particle size distribution function (PSD)
The model of ice crystals used in this study was developed by Baran and Labonnote (2007), and it is referred to as the ensemble model of cirrus ice crystals.The model has previously been fully described by Baran and Labonnote (2007), but a brief description is given here, and the model is shown in Fig. 2. The model consists of six elements.The first element is the hexagonal ice column of aspect ratio unity, and the second element is the six-branched bullet rosette.Thereafter, hexagonal monomers are arbitrarily attached as a function of ice crystal maximum dimension, forming 3-, 5-, 8and 10-monomer polycrystals.The ensemble model has previously been shown to predict the volume extinction coefficient, ice water content (IWC), and column-integrated IWC, as well as the optical depth of mid-latitude and tropical cirrus to within current experimental uncertainties (Baran et al., 2009(Baran et al., , 2011a(Baran et al., , 2013)).Moreover, the model also replicated 1 day of PARASOL cirrus observations of total reflectance, between the scattering angles of about 60 and 180 • ( Baran and Labonnote, 2007).It was further shown by Baran and Labonnote (2007) that the second randomised member of the ensemble model, randomised in such a way as to produce featureless scattering matrix elements, replicated 1 day of the global linear polarised reflectance measurements at close to cloud top.Heymsfield and Miloshevich (2003).
To demonstrate why the ensemble model can replicate in situ estimates of the volume extinction coefficient to within current experimental uncertainties, here the model predictions of the area ratio are compared against in situ estimates of naturally occurring area ratios.The area ratio, A r , is a useful measure of particle non-sphericity using in situ observations obtained from two-dimensional imaging probes.It is defined as the ratio of the projected area of a non-spherical particle to the area of a circumscribing circle of the same maximum dimension as the non-spherical particle.The area ratio value, A r , predicted by each member of the ensemble is shown in Fig. 3, and in the figure, the predicted values are compared against in situ estimates of A r .The in situ observations shown in Fig. 3 are obtained from a number of aircraft field campaigns that took place in the Arctic, at mid-latitudes and in the tropics (Heymsfield and Miloshevich 2003;Field et al., 2008;McFarquhar et al., 2013).The A r values from McFarquhar et al. (2013) were obtained in the Arctic at ice crystal maximum dimensions between about 35 and 60 µm, and the averaged values of A r were reported to be between 0.65 and 0.58, but with a standard deviation of ±45 %.These averaged values were found for 80 % of the estimates compiled by McFarquhar et al. (2013).The in situ estimates of A r from Heymsfield and Miloshevich (2003) are a synthesis of all the mid-latitude data (Eq. 2 from Heymsfield and Miloshevich 2003) compiled in that paper, and those estimates were prescribed an uncertainty of ±50 % for ice crystals less than 3000 µm, the uncertainty reducing to ±25 % for sizes greater than 3000 µm.The data from Field et al. (2008) were obtained in the tropics at ice crystal maximum dimensions between about 200 and 10 000 µm, and in Fig. 3, the upper and lower bounds on the data from Field et al. (2008) are shown, but on those bounds there is an uncertainty of ±30 %.Estimates of A r for ice crystals within 20 µm < D max < 100 µm were reported by Nousiainen et al. (2011) to be between 0.84 and 0.70 using tropical data, and the A r ratios in that paper were based on models of Gaussian random spheres (Muinonen 2000).
To compare members of the ensemble model against the in situ-derived estimates of A r , the maximum dimension is defined as follows.The maximum dimension of the hexagonal column is literally its maximum dimension (McFarquhar et al., 2013).The maximum dimension of the six-branched bullet rosette, and other members of the ensemble, is defined as the maximum distance across the particle when projected onto a two-dimensional plane (Heymsfield and Miloshevich 2003;Field et al., 2008).The area ratio (A r ) of the second member of the ensemble model shown in Fig. 2 assumes it to be in random orientation, which is a reasonable assumption since for the bullet rosette there is little difference between the projected areas if the particle is in random or preferred orientation due to its symmetry.All other members of the ensemble are assumed to be horizontally oriented along their maximum dimension with respect to the incident radiation.The oriented members are randomly oriented about their other two angles in three-dimensional space, with respect to the polar angle, to obtain an average of the projected areas, and it is these averages that are plotted in Fig. 3.In calculating the averaged A r values, the effect of shadowing is taken into account for each of the aggregated ensemble members.Figure 3 shows that the A r ratio calculated for each member of the ensemble model is within the current experimental range of possible A r values reported in the literature.The area ratio of the five-monomer aggregate at about 2500 µm is larger relative to the other ice aggregates.This is due to the other aggregate members of the ensemble model being more longitudinally elongated.Figure 3 shows why the model, as demonstrated by previous studies, can predict the volume extinction coefficient and the optical depth of naturally occurring cirrus to within current measurement uncertainties.In principle, given appropriate weights applied to the ensemble model, the volume extinction coefficient can be calculated for any type of cirrus (Baran et al., 2009(Baran et al., , 2011a(Baran et al., , 2013)), assuming a representative PSD is applied.Currently, the members of the ensemble model are distributed into six equal intervals of the PSD.However, this distribution of the predicted area throughout the PSD can change, given further information on the most general weights to apply to the model.
In this study, the PSD assumed is the moment estimation parameterisation of the PSD developed by Field et al. (2007), hereinafter referred to as F07.The F07 parameterisation relates the second moment (i.e.IWC) to any other moment via a polynomial fit to the in-cloud temperature.The parameterisation is based on 10 000 in situ measurements of the PSD, and the measurements were filtered using the method of Field et al. (2006) to reduce artefacts of ice crystal shattering at the inlet of the microphysical probes (Korolev et al., 2011), and the PSD was truncated at an ice crystal maximum dimension of 100 µm.The in situ observations were obtained from the mid-latitudes and tropics, at in-cloud temperatures between about −60 and 0 • C. The parameterisation does not ignore ice crystals less than 100 µm, but assumes that these ice crystals follow an exponential PSD.For ice crystal sizes greater than 100 µm, the parameterisation uses a gamma function, which was found to best fit the in situ-measured PSDs.The parameterisation adds together the exponential and gamma function to reconstruct the full PSD, given the IWC and incloud temperature.It has been previously shown that the F07 parameterisation is a good fit to in situ measurements of tropical and mid-latitude PSDs (Baran et al., 2011a;Furtado et al., 2014).Since the parameterisation fundamentally relates the second moment of the PSD to any other moment via the in-cloud temperature, in order to estimate representative PSDs, a mass-diameter relationship is required.In this paper, the ensemble-model-predicted mass-diameter relationship is used to generate the F07 PSDs.The ensemble model mass-diameter relationship was previously derived by Baran et al. (2011b), and in that paper it was shown that the ensemble model predicts the ice crystal mass of each particle to be given by 0.04 D 2 , where D is the maximum dimension of each ice crystal and the mass and D are in units of kg and m, respectively.. The ensemble model mass-diameter relationship is within the upper uncertainty of the Constrain-derived mass-diameter relationship derived by Cotton et al. (2013), and it is therefore representative of naturally occurring ice crystal mass.Furthermore, use is made of the F07 parameterisation, as we wish to be consistent with the PSDs assumed in the NWP model cloud scheme used later in the paper.

The single-scattering properties
Incident unpolarised sunlight is assumed to irradiate a collection of randomly oriented non-spherical particles, each of which possesses a plane of symmetry.The single-scattering properties that are applied to the PARASOL measurements are calculated using the Monte Carlo ray-tracing method, which was developed and made generally available by Macke et al. (1996a).Each member of the ensemble is randomised, using the method of distortion, and maximum randomisations are achieved using distortion and embedding within the volume of the ice crystals spherical air bubbles (Shcherbakov et al., 2006).The method of distortion involves randomly tilting the normal vector to the surface of the ice crystal (by assuming a uniform probability distribution; see Macke et al., 1996a, for further details) at the ice-air interface with respect to its original direction.In this way, at each refraction-reflection event, the directions of the ray paths are changed with respect to their original direction, with the result that featureless scattering matrix elements are predicted.The values of distortion can be between 0 and 1, where 0 represents unperturbed scattering matrix elements, and these retain scattering features such as halo and ice bows.As the distortion is increased to higher values, the optical features are removed in order to produce featureless scattering matrix elements.The distortion method attempts to replicate the complex processes that may occur on and within ice crystals, which could lead to featureless phase functions.Other authors refer to distortion as "microscale surface roughness".However, this description of surface roughness may not be accurate, as surface roughness can take on different forms (Mason 1971;Pfalzgraff et al., 2010;Bacon et al., 2003;Malkin et al., 2012).For instance, a theoretical electromagnetic study by Liu et al. (2013) has shown that the method of distortion does not accurately reproduce the scattering phase function at high values of idealised surface roughness.In this study, the method of distortion is merely used to randomise the ice crystal so that featureless scattering phase functions are produced.
Here the distortion parameter is assumed to have the values of 0, 0.15, 0.25 and 0.4.The distortion value of 0.4 is also combined with embedding the ice crystal with spherical air bubbles in order to achieve the maximum randomisation of the ice crystals so as to produce featureless phase functions.The upper distortion value of 0.4 was chosen as this was found to best fit 1 day of global POLDER-2 retrievals of directional spherical albedo and measurements of the linearly polarised reflectance (Baran and Labonnote, 2006).For the most randomised case of assuming a distortion value of 0.4 and embedding the ice crystal with spherical air bubbles, the phase functions are calculated using the modifications by Shcherbakov et al. (2006) applied to the raytracing code of Macke et al. (1996a).The statistics describing the tilt angles were shown by Shcherbakov et al. (2006) to be best represented by using Weibull statistics, where the Weibull distribution is defined by the scale (i.e. the distortion as described above) and shape parameters.This finding was based on cloud chamber measurements of the angular scattered intensity from a collection of ice crystals at a visible wavelength, and comparisons between measurements and ray-tracing results showed that the Weibull statistics were the better match to the measurements.Moreover, the choice of Weibull statistics is consistent with independent cloud chamber results found by Neshyba et al. (2013).For the most randomised case considered in this paper, the Weibull statistics are assumed to have the following scale and shape parameter values of 0.4 and 0.85, respectively, and for the spherical air bubble inclusions, a mean free path of 200 µm is assumed (Baran and Labonnote, 2007).The chosen values describing the Weibull statistics are also consistent with the values derived from independent cloud chamber measurements reported in Neshyba et al. (2013).To interpret the PARASOL measurements, the scattering phase function is required.The bulk-averaged scattering phase function, < P 11 (θ )>, is given by the following equation: where the vector (q) represents the elements of the ensemble model as a function of maximum dimension, n (q) is the F07 parameterised PSD, and C sca (q) is the scattering cross section of each of the ensemble model members.To generate the F07 PSDs, the IWC and in-cloud temperature are assumed to have the values of 0.01 gm −3 and −50 • C, respectively.The bulk-averaged asymmetry parameter, g , is calculated using the following equation: g = g (q) C sca (q) n (q) dq C sca (q) n (q) dq . (2) Another parameter of importance in calculating the total cloud reflectance is the single-scattering albedo, ω 0 , which is the ratio of the scattered radiation to that completely attenuated in the hemisphere of all directions.Here, the wavelength of 0.865 µm is considered.At such a weakly absorbing wavelength the value of ω 0 will be close to unity.The bulk-averaged volume extinction coefficient, β ext , is calculated from the following equation: where C ext (q) is the extinction cross section of each member of the ensemble model, calculated as a function of its maximum dimension.
Figure 4a shows the bulk-averaged ensemble-modelpredicted scattering phase functions, calculated at the wavelength of 0.865 µm, assuming distortion values of 0, 0.15 (slightly distorted), 0.25 (moderately distorted) and 0.4 with spherical air bubble inclusions.The complex refractive index of ice at 0.865 µm has the value of 1.304 + 2.400 × 10 −7 I, where I is the imaginary part of the refractive index (Warren and Brandt, 2008).The total optical properties are tabulated in Table 1 for each of the assumed ensemble models.Figure 5 shows the maximum contribution to the ice crystal scattering cross section per particle, as a function of maximum dimension, assuming the IWC and in-cloud temperature values given above.The figure shows that the maximum contribution to the scattering cross section occurs at a maximum dimension of about 50 µm.Defining the size parameter, x, as π D/λ, where λ is the incident wavelength, gives a value of x of about 182.This value of x means that the Monte Carlo ray-tracing method is within the range of x where the method is applicable (Yang and Liou 1996).As stated previously, the methods adopted throughout this paper to represent ice crystal complexity have been applied to generate a range of phase functions that retain and remove optical features that may be exhibited by naturally oc-curring cirrus phase functions.It is not as yet possible to simultaneously fully represent actual ice crystal complexity (i.e.surface roughening and internal hexagonal cavities) using electromagnetic methods at the size parameters considered in this paper (Baran, 2012).Therefore, approximations to ice crystal light-scattering properties are required at such size parameters.This is achieved, principally, through the method of geometric optics, and as such, approximations are required to represent surface roughness and ice crystal complexity.Here, both of these complexities are represented through the application of distorting ray paths (Macke et al., 1996) and spherical air bubble inclusions (Macke et al., 1995;Labonnote et al., 2001).When applied to the ice crystal model, both randomisations lead to featureless phase functions, which are the phase functions that are generally observed (Baran, 2012;Cole et al., 2013;2014).However, although the methods applied in this study result in featureless phase functions, this, however, does not necessarily mean that the resulting asymmetry parameter values shown in Table 1 cover the actual range of those values.Recent observations by van Diedenhoven et al. (2014a) of the asymmetry parameter derived from global polarimetric space-based remote sensing suggest median values in the range of 0.76 to 0.78.Ulanowski et al. (2006) reported that laboratory estimates of the asymmetry parameter, assuming highly surfaceroughened laboratory-grown rosette crystal analogues, could be as low as 0.61.On the other hand, the same study reported that smooth aggregate crystal analogues had asymmetry parameter values of 0.81.It is yet to be determined as to whether the asymmetry parameter values of actual cirrus ice crystals are as low as 0.61, and the values tabulated in Table 1 are in the upper range of van Diedenhoven et al. (2014a) and Ulanowski et al. (2006).Figure 4a shows as the distortion parameters gradually increase, the halo and ice bow features gradually diminish, and for the most randomised case, the scattering phase function becomes featureless and almost flat at backscattering angles.At the scattering angles of relevance to this study, the figure shows that discrimination between the model phase functions should be possible using the viewing geometry of PARASOL, especially at the scattering angles between 100 and 130 • .To demonstrate the feasibility of PARASOL to discriminate between the different ensemble model randomisations, Fig. 4b shows the ratio of the randomised to the pristine phase functions plotted against scattering angle.Figure 4b shows that, at scattering angles between about 80 • and 130 • , the ratio between the most randomised and pristine phase functions can reach values of about 1.1 to 1.5.At the distortion value of 0.15 (slightly distorted) and at scattering angles greater than 115 • , the ratio can still reach values of 1.1.However, at scattering angles between about 80 and 100 • , the values of the ratio between the pristine and slightly distorted, and moderately distorted phase functions are only slightly greater than unity, which means that discrimination between those models may not be possible at those particular scattering angles.However, due to the increasing values of the ratio at scattering angles between approximately 100 and 125 • , it might be possible to discriminate between models on a pixel-by-pixel basis at those particular scattering angles.
Of course, the phase functions derived from the ensemble model shown in Fig. 4 may not cover the entire range of possible cirrus phase functions as there are many possible cirrus habits that might occur at particular environmental temperatures (see, for example, Baran, 2012, and references therein).However, in the case of aggregates of hexagonal plates or hexagonal columns, it was shown by Baran (2009), using the ice aggregation model of Westbrook et al. (2004), that after three monomers were attached to the ice aggregate, the asymmetry parameters and phase functions asymptote to their limiting values.This asymptote occurs because the ice aggregation model predicts that the ice monomers making up the ice aggregate are well separated from each other.This separation is sufficient to reduce the effects of multiple scattering on the phase function, resulting in only slight modifications to the scattering angle positions of optical features (see, for example, Fig. 5 and Fig. 6 of Baran, 2009).This aggregation process is fundamental, and the same behaviour would be observed independent of the shape of the initial monomer (Westbrook et al., 2004).Therefore, in the case of pristine aggregates, the position of optical features on the phase functions would not be expected to be fundamentally different to those shown in Fig. 4a.If, on the other hand, the monomers that make up the ice crystal aggregate are sufficiently close to each other, then multiple scattering between monomers becomes important, as the scattered energy is increased and therefore also the phase function.However, the positions of the optical features exhibited by the ice aggregate phase functions do not significantly change position with respect to their scattering angles as these are principally determined by the hexagonal geometry (Um andMcFarquhar 2007, 2009).As discussed in the introduction to this paper, the observational evidence indicates that pristine ice crystals are a rarity in nature; therefore the phase functions of highly complex ice crystals exhibiting inclusions, cavities and surface roughness will produce featureless phase functions and the featureless nature of the phase function is invariant with respect to ice crystal habit.
To retrieve the spectral spherical albedo using PARASOL, a radiative transfer model is required; here the model developed by de Haan et al. (1986) is used and its application to PARASOL has been fully described by Labonnote et al. (2001).The radiative transfer model assumes a planeparallel cloud, but it is fully inclusive of multiple scattering.Also taken into account are layers of aerosol below the cloud and Rayleigh scattering above and below the cloud is also taken into account.The aerosol model assumed in the PARASOL retrieval has been previously described by Buriez et al. (2005), and so a description will not be repeated here.However, the aerosol is principally maritime-based, and so its optical depth will be much smaller than the cirrus optical depth, and as such it will not be of any significance for the purposes of this paper.At the wavelength of 0.865 µm, the PARASOL retrieval algorithm assumes that the sea surface has a reflectance value of 0.000612 (the foam contribution outside of sun glint) and a wind speed of 7 ms −1 .See Appendix A in Buriez et al. (2005) for a detailed derivation of the assumed PARASOL sea surface reflectance value.The scattering by aerosol and ocean glint all contribute to the directional variation of the retrieved cloud optical depth, and these effects are taken into account in the PARASOL retrieval algorithm.In the next section, the retrieval methodology is described.

Methodology
The methodology of retrieving the spectral spherical albedo using PARASOL multi-directional measurements of total reflectance has been previously described by (Doutriaux-Boucher et al., 2000;Buriez et al., 2001;and Labonnote et al., 2001), but a brief description of the retrieval is given here.The total reflectance of the cloud is specified by the vertical volume extinction coefficient, the vertical extent of the cloud and the scattering phase function.The cloud optical depth is therefore given by the integral of the vertical extinction over the vertical depth of the cloud.Since the cloud is essentially over a non-reflecting surface, the only directional information, under the assumption of a plane-parallel homogeneous layer, is provided by the assumed scattering phase function.However, inhomogeneity in the cloud can also affect the directional reflection as shown by Buriez et al. (2001), but this effect is not currently accounted for in the PARASOL retrieval algorithm due to its highly variable nature.It has been previously shown by Doutriaux-Boucher et al. (2000) that there is a one-to-one relationship between the cloud optical depth and the cloud spherical albedo (i.e.integral of the plane albedo over all incoming directions, where the plane albedo is a function of solar zenith angle alone) if the surface below the cloud is black.The cloud optical depth is retrieved by matching the simulated cloud reflectance to the measured cloud reflectance at each scattering angle.If the phase function were a perfect representation of the cloud, then the retrieved cloud optical depth will be the same at each scattering angle.Therefore, the retrieved spherical albedo would also be the same at each scattering angle.If the assumed phase function were a poor representation of the cloud, then this would result in a directional dependence on the spherical albedo, which would be unphysical.This retrieval methodology forms the basis of this paper, and it has been applied by other studies that have utilised PARASOL measurements to test ice cloud scattering phase functions (see, for example, Doutriaux-Boucher et al., 2000;Labonnote et al., 2001;Baran et al., 2001;Knapp et al., 2005;Baran and Labonnote, 2006).
As previously stated, the retrievals of spherical albedo are performed on a pixel-by-pixel basis, and the data products derived from the PARASOL observations are only used if the following four conditions are met: (i) for each 6 km × 6 km pixel the cloud fraction is unity, (ii) the total number of view angles ≥ 7, (iii) the difference between the minimum and maximum sampled scattering angle is greater than 50 • , and (iv) only pixels over the sea are considered.The total number of PARASOL pixels that are within the area of interest shown in Fig. 1 is 297.As previously stated, since the true cloud spherical albedo is, by definition, independent of direction, then for each pixel, the retrieved averaged directional spherical albedo, S , should be identical to the directional spherical albedo, S(θ ), where θ is the scattering angle, if the model phase function were a perfect representation of the cloud.The retrievals of S(θ ) depend on the assumed scattering phase function, the vertical volume extinction coefficient and ω 0 .
The averaged spherical albedo, S , for each pixel is defined by where N is the total number of viewing directions, which for the case considered in this paper is between 7 and 8. To find the phase function that best minimises the spherical albedo differences, the root-mean-square error (RMSE) is found for each pixel, which is given by where in Eq. ( 5), S = S − S(θ ).The RMSE is one general measure for choosing the best-fit model to the observa-tions.However, once a set of RMSE minimisers has been identified, one should assess where the identification is statistically significant.This is needed to rule out the possibility that the differences in RMSE for the different distortions could have resulted from chance.To test this, we apply the Levene (1960) test statistic, as it is less sensitive to the condition that the data must be normally distributed than the usual F statistic, which is generally used to test whether the variances between two samples are equal, provided the data follow a normal distribution.In the Levene test, the samples, k, are tested for homogeneity of variances between the k samples.The total number of data points contained in all samples is given by N. The Levene null hypothesis is that variances between k samples are equal.The Levene null hypothesis is rejected, at some level of significance, α, if the Levene test statistic, W , is greater than F α (k − 1, N − k), which is the upper critical value of some F distribution with k − 1 and N − k degrees of freedom.We consider pixels for which the derived RMSE values do not exceed 100 % to require testing for significance using the W test.For these pixels, the W test statistic is applied to test whether the model variances in S j are different at the 5 % significance level.If the null hypothesis is rejected, then that pixel is assigned a particular model phase function.The 5 % or α = 0.05 significance level is chosen, as this is simply between α = 0.1 (10 %) and α = 0.01 (1 %) significance levels, so that the model test is neither too easy nor unrealistically hard, respectively.
In the sections that follow, the model phase functions shown in Fig. 4 and the total optical properties given in Table 1 are applied to the PARASOL measurements, on a pixelby-pixel basis, to retrieve the phase function that best minimises Eq. ( 5) and satisfies rejection of the Levene null hypothesis.Results of this analysis are then used to explore possible relationships between the shape of the scattering phase function and RH i .

Validating the NWP model field of RH i
Before exploring the possible relationship between RH i and the scattering phase function, it is first necessary to show that the NWP-model-predicted field of RH i is sufficiently accurate for the purposes of this paper.Firstly, Fig. 6 shows the NWP-model-predicted field of the water vapour mixing ratio and the location of the PARASOL pixels, and the position of the aircraft within that field.The aircraft positions were predominantly located around the areas of semi-transparent cirrus, with generally no cloud beneath, as shown by Fig. 1. Figure 6 shows that there is considerable variation in the water vapour field at about the cloud top around the north and east of Scotland.As a consequence of this variation, there will be a sufficient change in the RH i field at the cloud top as a function of position to relate retrieval results to the model field.From Fig. 6 we note that the NWP-model-predicted the cloud top to be in the vicinity of 10 km for the region of interest.Figure 7a shows the aircraft-mounted lidar estimate of the volume extinction coefficient as a function of altitude, and the bottom panel shows the derived lidar optical depth as a function of time.In Fig. 7a, note that only the lidar-derived profile of volume extinction coefficient greater than 6 km is shown.This is because, at altitudes less than this, the lidar equation becomes numerically unstable in clear air, and as such there were no meaningful retrievals of cirrus volume extinction coefficient below this altitude.Figure 7a shows that the lidar-estimated cloud-top altitude was at about 10 km at approximately 13:33:00 UTC, which is in good agreement with the NWP model prediction, and the lidar position at that time is marked by the X symbol in Fig. 6.
To validate the NWP model prediction of the RH i field, use is now made of the aircraft-mounted ARIES measurements, which are applied to obtain retrievals of RH i .The retrieval of RH i from the ARIES measurements is achieved using the Havemann-Taylor Fast Radiative Transfer Code (HTFRTC) (Havemann 2006;Havemann et al., 2009) and the retrieval method of Thelen et al. (2012).Moreover, the ARIES-based retrieval of RH i is validated against the dropsonde measurements of RH i .The HTFRTC is a principal component based radiative transfer model and is fully in- clusive of the atmosphere and exact multiple scattering.The ARIES spectrum was averaged over 10 spectra and was denoised using principal components, which act as a low-pass filter.For this case, European Centre for Medium-Range Weather Forecasting (ECMWF) atmospheric profiles are applied as the background state and the optimal estimation (OE) method of Rodgers (1976) is used, which is a Bayesian method.Here, OE is used to retrieve the most likely atmospheric state, and the method includes a rigorous treatment of error.The errors arise from the ARIES instrument itself and from the forward model, as well as from the ECMWF model background fields.The treatment of error by OE assumes that the errors are described by a Gaussian distribution.Here the background errors in the temperature, RH and IWC are assumed to be typically ±0.4 K, ±10 % and ±50 %, respectively.Given the errors, ARIES measurements and simulated measurements using HTFRTC, OE uses a minimisation procedure to find the most likely atmospheric state that best describes the measurement set, given the retrieved parameters or state vector.Currently, the state vector in the HTFRTC retrieval method (Thelen et al., 2012) is composed of the temperature profile, the relative humidity profile, homogeneous cirrus IWC, surface temperature, and surface emissivity.The temperature and relative humidity profiles are retrieved at all 70 levels of the Met Office operational suite of models.
The NWP model prediction of the vertical profile of RH i is compared against the ARIES retrievals, dropsonde measurements and in situ aircraft measurements from the GE and FWVS instruments.The various comparisons are shown in Fig. 8a and b for two different locations.The Figure 8a and b show that the two in situ RH i measurements are in good agreement with each other, whilst the dropsonde took some time to adjust to the prevailing atmospheric conditions.After this adjustment time, the infrared retrievals of RH i , in the presence of cirrus, are in good agreement with the dropsonde and are within the range of RH i measured by the two in situ instruments.The figure demonstrates that the retrieval of RH i using high-resolution passive infrared measurements is sufficiently accurate and can be obtained, in presence of cirrus, on a global scale using space-based high-resolution instruments such as the Infrared Atmospheric Sounding Interferometer (IASI).Furthermore, below the cloud, the dropsonde and retrievals are in very good agreement in the drier regions of the atmosphere, down to pressures of about 600 hPa.Moreover, the retrievals and dropsonde are in good agreement, down to pressures of about 1000 hPa.The two different retrieval colours represent the retrievals based on the two aircraft runs above the cirrus that were previously described.Each of the runs was 10 min in length.There were approximately eight ARIES retrievals per run.Figure 8 demonstrates that the retrievals, dropsonde measurements and in situ measurements are sufficiently consistent to compare against the NWP model.Figure 8a shows the various comparisons at the latitude of 59.14 • N and longitude 3.85 • W, which corresponds to the upper left of Fig. 6.The figure shows that the NWP model prediction of the vertical profile of RH i is consistent with the retrievals and measurements.Figure 8b  The NWP-model-predicted cloud depth is therefore about 3 km, which is also in good agreement with the lidar-derived maximum cloud depth shown in Fig. 7a at approximately 13:33:00 UTC, when the aircraft was above the cloud top.

Estimating the shape of the scattering phase function and its relationship to RH i
In this section, the methodology described in Sect. 4 is used to estimate the ensemble model phase function which best minimises the RMSE and rejects the Levene null hypothesis at the 5 % significance level.The ensemble model phase functions used here were previously described in Sect.3.1.1,and are shown in Fig. 4. The results from the phase function estimates for each pixel, showing the phase function model that best minimised RMSE, are shown in Fig. 9a.The total number of retrievals, showing only those retrievals over the sea, in Fig. 9a is 292.However, 130 of these retrievals correspond to indeterminate results.The reason for the indeterminate results at those pixels is because the retrieved spherical albedo at each of the scattering angles was the same for all ensemble models.The similarity of retrieved results in the indeterminate cases is because the retrieval conditions stated in Sect. 4 were not met.These indeterminate results are shown as black squares in the figure .A comparison between Figs. 9a and 1 show that the indeterminate results generally occurred in the presence of multi-layer cloud.Figure 9b shows the averaged retrieved PARASOL decadal optical thickness (averaged over all available scattering angles) at each of the pixels shown in Fig. 9a.The figure shows that the retrieved PARASOL optical thickness ranged between less than 1 and up to about 250.The largest optical thick-nesses retrieved by PARASOL are associated with the broken frontal cloud shown in Fig. 1 (right-hand side of the figure), and the positions of the broken frontal cloud fields are also predominantly associated with the positions of the indeterminate results shown in Fig. 9a. Figure 9a and b show that even for PARASOL-retrieved optical thicknesses of between about 10 and 30, discrimination between ensemble models is still possible.The physical reason for this was recently given by Zhang et al. (2009).In their paper, it is physically argued that, even if the optical thickness is increased to large values, the shape of the phase function is still retained at top of the atmosphere.This is because scattering within the cloud is dominated by forward scattering, which results from strong diffraction in the forward direction (Macke et al., 1995), and this single-scattering information is still retained in the presence of strong multiple scattering.However, at the largest retrieved optical thicknesses shown in Fig. 9b, multiple scattering will be so strong that discrimination be-tween models will no longer be possible, and some of these largest optical thicknesses are associated with the indeterminate results.Figure 9c shows the estimated randomisations at each PARASOL pixel, but with the indeterminate results removed, again using only the minimised RMSE value to select the best model phase function.The yellow squares in Fig. 9c correspond to the most randomised phase function (i.e.distortion = 0.4 with spherical air bubble inclusions), and the number of pixels associated with this colour is 150, which, from the figure, is clearly the most common.However, 12 of the pixels shown in the top left of the figure are not associated with the most randomised phase function.Rather, these pixels were found to be associated with the pristine phase function (distortion = 0), the slightly distorted phase function (distortion = 0.15), or the moderately distorted (distortion = 0.25) phase function.The retrievals which best minimised the RMSE assuming the pristine phase function are represented by the purple squares.The pixels represented by the dark-and light-brown squares were found to be associated with the slightly (distortion = 0.15) and moderately distorted (distortion = 0.25) ensemble model phase functions, respectively.Therefore, the results shown in Fig. 9c indicate that, on a pixel-by-pixel basis, the most randomised ice crystal model phase functions may not always be the best fit to multi-angular spectral albedo measurements, at least if only the minimised RMSE value is used to select the bestfit phase function.Note, however, that we have so far disregarded whether or not the discrimination of phase function, based on RMSE, is statistically significant.The reliability of the use of minimised RMSE values only to select the best model phase function is further examined in the paragraphs that follow.
The estimated randomisations for two of the pixels shown at the top left of Fig. 9c are further examined in Fig. 10a and  b.The figure shows the spherical albedo differences plotted as a function of scattering angle for each of the two pixels, and in each of the figures, the RMSE values are shown that were derived from the spherical albedo differences assuming the four models.The first pixel shown in Fig. 10a is located at latitude 59.03 • and longitude −3.62 • , and this pixel has been assigned the fully randomised phase function.It can be seen from the figure that, in this case, the spherical albedo differences predicted by the fully randomised phase function are closer to the zero line than the other models for all scattering angles considered.In this case, the RMSE value found for the fully randomised model is a factor of 4.6 smaller than the value of the RMSE found for the pristine model.In contrast to Fig. 10a, Fig. 10b shows the spherical albedo differences for the pixel located at latitude 59.14 • and longitude −3.84 • , and this pixel has been assigned the pristine model phase function.In this case, the pristine model phase function is closest to the zero line at the scattering angles of about 92, 107 and 123 • .However, at the scattering angles of about 99 and 113 • , the pristine model phase function predicted spherical albedo differences are similar to the predictions obtained assuming the moderately and fully randomised phase functions, respectively.The RMSE values shown in the figure indicate that the pristine model phase function best minimises the RMSE, although this value is only a factor of 1.3 smaller than the RMSE value found assuming the fully randomised phase function.In general, if the spherical albedo differences are visually examined at all scattering angles, then one could conclude that, for this pixel, no one phase function model best describes all seven measurements.The question then arises of whether this is true for all the 11 other pixels that are associated with structure in their scattering phase functions at backscattering angles.To test this quantitatively, the Levene test statistic is now applied to all 12 pixels for which Table 2.The Levene test statistic, W , applied to test homogeneity of variances in spherical albedo differences between two groups of scattering phase function models for each set of pixels.In the table, the two phase function models are represented for each set of pixels by their assumed distortion values, referred to as Model pair; the total number of pixels used in each test is n.The null hypothesis is given by H 0 , which is either rejected or accepted; k is the number of samples; N is the total number of observations in the two samples; and F 0.05 (k, N −k) is the value of the tabulated upper critical value at the 5 % significance level composed of k and N − k degrees of freedom.From Fig. 9c it can be seen that using minimised RMSE test, 5 of the 12 pixels are associated with pristine model phase functions (distortion = 0), whilst 4 pixels are associated with slightly distorted phase functions (distortion = 0.15) and the other 3 pixels are associated with moderately distorted phase functions (distortion = 0.25).All pixels associated with each of the above three distortion values were combined together.For each of the distortions, the W statistic was obtained in groups of two, so that k = 2.The variances in the spherical albedo differences obtained with the RMSEdetermined best-fit phase function were compared against the variances obtained assuming all other model phase functions.For each group of two, the test W statistic was computed and then compared against the tabulated upper critical value of theF α (k − 1, N − k) distribution in order to accept or reject the null hypothesis at the 5 % significance level.The results of this analysis are presented in Table 2.

Model pair
In the case of the five pixels associated with pristine phase functions, it can be seen from Table 2 that the Levene null hypothesis must be accepted.Therefore, the variances in the spherical albedo differences determined using the RMSE best-fit model are not sufficiently different from the variances obtained using all other phase function models.A similar result to the above was found for the three pixels, which were associated with the moderately distorted phase function (distortion = 0.25).For the four pixels associated with the slightly distorted model phase function (distortion = 0.15), Table 2 shows that the null hypothesis can be rejected when its variances are compared against the variances obtained assuming the fully distorted model phase function (distortion = 0.4 with spherical air bubble inclusions).However, for all other assumed models, for this group of four pixels, the  2 show that using minimised RMSE values alone may not be sufficient to select model phase functions on a pixelby-pixel basis and that some other test statistic is required to compliment the RMSE method.The Levene test statistic was also applied to some pixels associated with the most randomised phase function in order to test whether W F for these pixels.The results of this test are presented in Table 3.In this case, seven pixels were selected between latitudes 57.92 and 58.92 • , and between longitudes −3.42 and −3.71 • .As before, the seven pixels were combined, and the resulting variances in spherical albedo differences obtained assuming the most randomised phase function were compared against the variances obtained assuming model distortion values of 0.15, 0.25 and 0. The results from Table 3 show that the null hypothesis can be very strongly rejected at α = 0.05 (5 % significance level).The results of this analysis strongly suggest that the selection of the most randomised phase function using minimised RMSE values is acceptable, as illustrated by the example case in Fig. 10a.
Since no one model phase function can be uniquely assigned to any of the 12 pixels, which show small differences in RMSE between models, suggests that the model phase functions do not correctly describe the backscattering properties of the cirrus located at those pixels and/or there might be underlying water cloud affecting the results.To investigate the possibility that there might be an underlying water cloud beneath the cirrus contaminating the 12 pixels, the rangecorrected lidar images were further investigated under higher resolution to see whether there was any water cloud beneath the cirrus.The aircraft passed over the 12 PARASOL pixels at between about 13:20:00 and 13:30:00 UTC.Between these times, the high-resolution lidar images showed only reflection from the sea surface with no evidence of underlying water cloud (results not shown here for reasons of brevity).It is noted here that the averaged retrieved PARASOL optical thickness (averaged over the 12 pixels) was found to be 1.81 ± 0.32, assuming model distortion values of between 0 and 0.25.Furthermore, the retrieved PARASOL averaged optical thickness over the same 12 pixels, assuming the fully randomised phase function, is 1.52 ± 0.26.Both retrievals are within the range of the lidar estimates shown in Fig. 7b, the average of which is about 1.2.The lidar vertical profiles of optical depth were obtained when the aircraft was located above the cirrus, at an altitude of almost 11 km, which occurred during the times shown in the figure.There is a gap of about 3 min shown in Fig. 7b, which is the time required for the aircraft to turn and commence the second straight and level run.The PARASOL retrievals of cirrus optical thickness and the lidar retrievals are both consistent with one another.If there had been an underlying water cloud beneath the cirrus at the time of the PARASOL overpass, then the retrievals would not have been consistent.To examine the issue of underlying water cloud further, at a time nearer to the PARASOL overpass, generally available space-based cloud products were also examined.
The space-based remotely sensed cloud products are available from http://www-pm.larc.nasa.gov/.The cloud products that were examined were obtained at the time of 13:00:00 UTC, which is the time closest to the PARASOL overpass.The available remotely sensed cloud products from the site include detection of multiple cloud layers and cloudtop pressure.Analysis of these images indicated that the cloud overlying the 12 pixels was of a single layer, and the cloud-top pressure of this single-layer cloud was retrieved to be between 100 and 200 mbar (again not shown here for reasons of brevity).These independent space-based remote sensing results indicate that there was no underlying water cloud covering the 12 pixels, and this is consistent with the analysis of the high-resolution range-corrected lidar images as well as the PARASOL and lidar retrievals of cirrus optical depth.
Since it is unlikely that underlying water cloud affected the results discussed above indicates that there might have been backscattering structure present on the cirrus phase function which is not represented by any of the models.Or more simply, there was insufficient scattering angle information available to distinguish between models.Interestingly, Baran et al. (2012) also found that, for a case of mid-latitude, very high IWC anvil cirrus near to the cloud top, the PN-measured averaged scattering phase function also exhibited unusual backscattering features.Clearly, such findings of optical features on the scattering phase function of naturally occurring ice crystals indicate the need for radiometric or in situ observations to sample the scattered angular intensities over a more complete range of scattering angle than is currently possible.Measuring the forward and backscattering intensities alone is not sufficiently general (Baran et al., 2012).However, the most common retrievals shown in Fig. 9a are representative of the most randomised ice crystals, and these have featureless phase functions.For the purposes of retrieving cirrus properties using global radiometric measurements, it is most likely that featureless phase functions are still generally better at representing cirrus radiative properties than their purely pristine counterparts (Foot, 1988;Baran et al., 1999Baran et al., , 2001;;Baran and Labonnote, 2006;Baum et al., 2011;Cole et al., 2013;Ulanowski et al., 2013;Cole et al., 2014).Here, it is also of interest to note the change in the asymmetry parameter values shown in Table 1.From the pristine ensemble model phase function to the most randomised ensemble model phase function, the change in the asymmetry parameter is about 5 %.A change in the asymmetry parameter of 5 % is radiatively important, as illustrated by the following example.Given that the instantaneous solar irradiance arriving at the top of Earth's atmosphere is about 1370 Wm −2 .Under the assumptions of conservative scattering and a dark ocean below the cirrus, a change of 5 % in the asymmetry parameter results in a difference of about 43 W m −2 in the solar irradiance reflected back to space.However, the difference of 43 W m −2 could be an underestimate if actual values of the asymmetry parameter are lower than reported in Table 1 (Ulanowski et al., 2006).Even so, a difference of 43 W m −2 is very significant with regard to the energetics of the Earth's atmosphere, and indicates why it is important to globally constrain values of the asymmetry parameter (Baran, 2012;Ulanowski et al., 2006Ulanowski et al., , 2013;;van Diedenhoven et al., 2014a).
The PARASOL estimations of the shape of the scattering phase function, based on applying the minimised RMSE and the Levene tests, are shown in Fig. 11a.In the figure, the yellow pixels were assigned the most randomised phase function, the brown pixels are the locations where no one model phase function could be uniquely assigned.The blue pixels show the locations where either phase function model, apart from the most randomised phase function, could be assigned.The results shown in Fig. 11a are now directly compared against the NWP-model-predicted RH i field at the cloud top, which is shown in Fig. 11b.The NWP model results are shown at a cloud-top altitude of 10 km.On comparison with Fig. 11a, it can be seen from Fig. 11b that the most randomised phase functions (i.e.yellow squares) generally correspond to model pixels with RH i > 1.0.Conversely, the 12 pixels for which no one model phase function could be assigned generally correspond to NWP model pixels with RH i < 1.0.The results for RH i > 1 are broadly consistent with the findings of Gayet et al. (2011) and Ulanowski et al. (2013).The results of the former paper suggested that featureless phase functions were generally associated with RH i > 1.0.Whilst the laboratory studies of Ulanowski et al. (2013), on ice crystal analogues, indicate that at higher levels of ice supersaturation, surface roughness on the ice crystal increased.This increase in surface roughness would naturally lead to featureless phase functions (Yang and Liou, 1998;Ulanowski et al., 2006;Baran, 2012;and references contained therein).
The NWP results shown in Fig. 11b are at the cloud top.However, the PARASOL retrievals might be based on reflected solar radiation that comes from the extent of the cloud and not just from the cloud top.In reality, solar radiation at 0.865 µm will penetrate to some depth within the cloud layer, and this depth of penetration needs to be calculated to test whether the assumption of cloud-top penetration is correct.
To calculate the depth of penetrating radiation at 0.865 µm, a Monte Carlo radiative transfer model has been used to represent the cirrus layer of relevance to this study.The Monte Carlo model used here is fully described by Cornet al. (2009).A description of the Monte Carlo model setup and defini-tion of the probability of penetration is given in Appendix A. The percent probability of penetration as a function of cloud depth and optical depth is shown in Fig. 12a and b.In the figures, the percent probability of penetration at 0.865 µm is defined as the last position (distance from the cloud top) of the photon before leaving the cloud to reach the sensor.Results are shown in the figure for the forward-and backwardscattered radiation in the solar plane, respectively.
Figure 12a and b show that, by a depth of 1 km from the cloud top, the probability of penetration has been more than approximately halved for optical depths greater than 0.3.By 1.5 km from the cloud top, the probability of penetration is generally less than 5 %.The percent probability of penetration shown in Fig. 12a and b is similar.This is because the scattering phase function used in the Monte Carlo calculations, at backscattering angles, is largely invariant with respect to the scattering angle.This is simply because the scattering phase functions representing the most randomised ice crystals are flat and featureless at backscattering angles.
From Fig. 12a and b, it can be concluded that the PARA-SOL measurements of the total reflectance are biased towards the cloud top, and therefore comparison between the NWP model at the cloud top and PARASOL estimations of the shape of the scattering phase function is acceptable.
This paper has demonstrated the potential of using spacebased remote sensing to investigate relationships between the scattering properties of ice crystals and atmospheric state parameters.However, one drawback of current space-based multi-angle measurements is the limited range of multiangle samplings: in this paper, only seven measurement angles were available.In regions where NWP model values of RH i were generally less than unity, it was impossible to assign a model phase function to the PARASOL observations.Clearly, if more multi-angle samplings were available, coupled together with a greater range of scattering angle, then discriminating between different phase function models may become more likely.
Climate model parameterisations of the asymmetry parameter are currently assumed to be invariant with respect to atmospheric state variables.It is desirable, as argued by Baran et al. (2009Baran et al. ( , 2014) ) and Baran (2012), to relate general circulation model prognostic variables directly to ice optical properties, so that the prognostic variables can then be directly related to space-based radiometric measurements.Only through directly relating general circulation model prognostic variables to radiative measurements can the possibility of error cancellation be removed from within climate models.
This paper has explored the relationship between RH i and the shape of the scattering phase function for one case of mid-latitude cirrus that occurred on 25 January 2010.This relationship has been explored by combining high-resolution NWP model RH i fields with satellite retrieval of the directional spherical albedo at 0.865 µm.The satellite observations were obtained from the PARASOL spherical albedo product at scattering angles between about 80 and 130 • .The satellite observations were analysed on a pixel-by-pixel basis.It was found that featureless phase functions, representing significant ice crystal randomisation, best minimised differences between the directionally averaged spherical albedo and the directional spherical albedo for about 90 % of the pixels for which discrimination was possible.However, for about 10 % of the pixels, it was found that discrimination between model phase functions, based on spherical albedo differences, was not possible.In general, if multi-angular data are not available, given that over 90 % of the spherical albedo differences contained in this study were best described using featureless phase functions, then featureless phase functions are more likely to be a correct assumption for general application to the remote sensing of cirrus properties, rather than phase functions containing optical features.
It has also been demonstrated in this paper that the Met Office nested high-resolution NWP-model-predicted vertical profiles of RH i are sufficiently accurate to combine with remote sensing data to study relationships between atmospheric state variables and the fundamental scattering properties of cirrus.Independent retrievals of the vertical profile of RH i , using aircraft-based high-resolution infrared data, dropsonde measurements and in situ measurements of RH i , showed excellent agreement with the NWP-model-predicted vertical profile of RH i for two very different locations within the cirrus field.Moreover, the NWP model prediction of cloud top, vertical depth and cloud base were shown to be consistent with lidar measurements of the same quantities.Assuming that the NWP RH i fields are representative of truth, the model fields were directly related to the remote sensing observations of the shape of the cirrus scattering phase function.
For this one cirrus case, it is found that featureless phase function models, representing highly randomised ice crystals, were shown to be generally associated with NWP model RH i values greater than unity.In the cases where the NWP model RH i values were found to be generally less than unity, no one single-scattering phase function model could be assigned to the PARASOL pixel using a quantitative statistical measure.The possibility of these pixels being affected by the issue of underlying water cloud below the cirrus was also investigated.Using high-resolution lidar images, retrievals of cirrus optical depth obtained from PARASOL and the aircraft-mounted lidar, and generally available space-based cloud products, it was found that it is unlikely that these pixels were affected by underlying water cloud.Given this finding, the model phase functions did not have the correct structure in the backscattering part of the phase function or, more simply, there was not enough scattering angle information to be able to discriminate more clearly between the different phase function models.Given the latter reason, it would clearly be more desirable if future spacebased instrumentation could more clearly resolve, and over a greater scattering angle range, the backscattering part of the cirrus phase function.This paper has also demonstrated that high-resolution interferometer data can be used, in the presence of optically thin cirrus, to retrieve the vertical profile of RH i .This interferometric capability, which already exists in space through IASI, could be combined with improved resolution of multiple viewing satellites to explore the relationship between atmospheric state parameters and shape of the scattering phase function on a global scale.This paper has demonstrated the potential for obtaining such global spacebased measurements.There already exist aircraft-based instruments that measure the in situ light-scattering properties of atmospheric ice at particular locations, such as those used by Gayet et al. (2011), Ulanowski et al. (2013) andvan Diedenhoven et al. (2014a).Preferably, new in situ instrumentation should be developed that is capable of measuring the scattered intensities over a larger range of scattering angles than is currently possible (Baran et al., 2012).
Currently, the ice radiation scheme in climate models does not take into account ice crystal complexity as a function of atmospheric state.Further research in this area will prove or disprove whether this climate model assumption needs to change.

Appendix A
The radiative transfer model assumes a plane-parallel layer with a vertical extent of 3 km.This cloud depth is consistent with the lidar result shown in Fig. 7a.The vertical resolution of the cloud layer is assumed to be 0.1 km, and the cloud top is situated at an altitude of 10 km, which is also consistent with the lidar result shown in Fig. 7a.The relevant Sun-satellite geometry for this case has been applied to the Monte Carlo calculations; that is, the solar zenith angle is 75 • .The view angle of PARASOL at the time of the overpass has an average value of 50 • .The standard deviation of the PARASOL view angle is generally less than 2 • .For the purposes of this study, an average view angle will suffice.The view angle used in the Monte Carlo model was set to a value of 50 • .Moreover, the PARASOL azimuth angle did not vary significantly and the standard deviation of this angle was no more than 4.5 • .With little variation in the satellite geometry, the Monte Carlo calculations have been performed in the solar plane to obtain the most general results.This means that ϕ-ϕ 0 = 0 for forward-scattered radiation, and ϕ-ϕ 0 = 180 • for backward-scattered radiation, in the principal plane, where ϕ and ϕ 0 are defined as the satellite and solar azimuth angles, respectively.The cloud microphysical model is assumed to be the most randomised ensemble model phase function, which was the most common retrieval shown in Fig. 9a and c, and the values of the volume extinction coefficient and single-scattering albedo were taken from Table 1.The probability of penetration is calculated by computing the total distance the ray travels within each sub-layer of the cloud.From this analysis, the percent probability of penetration is defined as the last position (distance from the cloud top) of the photon before leaving the cloud to reach the sensor.

Figure 1 .
Figure 1.A high-resolution composite MODIS image of the semitransparent cirrus case that occurred on 25 January 2010 located over north-east Scotland.The latitude and longitude grid is superimposed on the image showing latitude 58 to 60 • (left side) and longitude −8 to 0 • (bottom).The composite image was formed by combining the MODIS red, green and blue channels to obtain the closest "true" colour image.The image is from the NERC Satellite Receiving Station, Dundee University, Scotland (http://www.sat.dundee.ac.uk/).

Figure 2 .
Figure 2. The ensemble model as a function of ice crystal maximum dimension, D max .The first element of the model is the hexagonal ice column of aspect ratio unity (first top left), followed by the 6branched bullet rosette (top centre), the 3-monomer hexagonal ice aggregate (top right), 5-monomer ice aggregate (first bottom left), 8-monomer ice aggregate (bottom centre) and the 10-monomer ice aggregate (bottom right).

Figure 3 .
Figure 3.The ensemble model area ratio, A r , as a function of ice crystal diameter or D max .The key is shown on the upper right-hand side of the figure.The members of the ensemble model are represented by the filled cyan circles.The in situ observations are from Field et al. (2008) (red lines), where the plus and cross signs represent the lower and upper range of those observations and those ranges have an uncertainty of ±30 %.The blue error bar represents the mean and range of observations taken from McFarquhar et al. (2013) and the purple error bars represent the uncertainty in the observations taken fromHeymsfield and Miloshevich (2003).

Figure 4 .
Figure 4. (a) The decadal logarithm of the ensemble model normalised scattering phase function as a function of scattering angle assuming a variety of distortions.The model cases shown are the pristine (black), slightly distorted (red), moderately distorted (dashed blue) and fully distorted with spherical air bubble inclusions (dotted purple).(b) The ratio between the distorted and pristine ensemble model phase functions as a function of scattering angle.Shown here are slight distortion (red), moderate distortion (dotted green) and full with spherical air bubble inclusions (dotted blue).The key is shown in each of the figures.

Figure 5 .
Figure 5.The scattering coefficient per particle (m −2 ) as a function of ice crystal maximum dimension, D max .The PSD was generated assuming IWC and in-cloud temperature values of 0.01 gm −3 and −50 • C, respectively.

Figure 6 .
Figure 6.The field of the water vapour mixing ratio (Q v ) on 25 January 2010 at 13:00 UTC, between latitudes 57.8 and 59.7 • and longitudes −5.3 and −1.8 • .The units of Q v are Kg Kg −1 .The PARASOL pixels are represented by the open circles and the aircraft track is represented by the solid line, and X marks the location where the aircraft was directly above the cloud at about 13:33:00 UTC.
Figure 7. (a) The lidar-derived cloud volume extinction coefficient as a function of altitude (m) and time in units of hours after midnight (UTC).The colour bar on the right-hand side of the figure indicates values of the cloud volume extinction coefficient in units of m −1 , and the solid line represents the aircraft altitude.(b) The lidar-derived cloud optical depth from 300 m below the aircraft to the cloud base as a function of UTC time, and the horizontal solid line shown in the figure indicates an optical depth value of unity.

Figure 8 .
Figure 8.A comparison between the retrievals, dropsonde measurements, in situ measurements and NWP model predictions of RH i plotted against the pressure (hPa) for two different locations.(a) The pixel located at longitude −3.84 and latitude 59.14 • and (b) the pixel located at longitude −3.20 and latitude 57.97 • .In (a) and (b), the retrievals are represented by the purple and green plus signs, dropsonde measurements are the solid grey line and filled grey circles, the General Eastern hygrometer is the solid green line, and the FWVS hygrometer is the solid red line.
is similar to Fig. 8a but for the location 57.97 • N and 3.20 • W, which corresponds to the lower left of Fig. 6.In this figure, the NWP model and retrievals can reach values of RH i of up to about 1.20.Figure 8a and b validate the NWP model prediction of RH i , and thus this model can be used to compare against the PARASOL estimates of ice crystal randomisation.Figure 8a and b show that the NWPmodel-predicted cloud top is at about 200 hPa (∼ 10 km), and the cloud base is at about 400 hPa (∼ 7 km).

Figure 9 .
Figure 9.The PARASOL estimates of ensemble model randomisations (based on minimised RMSE) and retrievals of optical thickness as a function of latitude and longitude.(a) The estimated ice crystal randomisation, where the indeterminate results are shown by the black squares, the most randomised phase functions (distortion = 0.4 with spherical air bubble inclusions) by the yellow squares, and the pristine phase functions (distortion = 0) by the purple squares; dark-and light-brown squares represent the slightly distorted (distortion = 0.15) and moderately distorted (distortion = 0.25) phase functions, respectively.(b) The PARASOL-retrieved averaged optical thickness, averaged over all scattering angles, where the decadal logarithm of the retrieved optical thickness is shown by the colour bar on the right-hand side of the figure.(c) The same as (a) but with the indeterminate results removed.

Figure 10 .
Figure 10.Differences between the directionally averaged (< S >) and directional (S(θ)) spherical albedos as a function of scattering angle at two pixel locations.(a) The spherical albedo differences for the pixel located at 59.03 • and longitude −3.62 • , assuming the pristine ensemble model (dist = 0) (open red circles), the slightly distorted model (dist = 0.15) (open green triangles), the moderately distorted model (dist = 0.25) (open blue diamonds), and the fully randomised model (dist = 0.4 with spherical air bubble inclusions) (open purple pentagons).(b) The same as (a) but for the pixel located at latitude 59.14 • and longitude −3.84 • .The zero difference line is shown by the solid bold line, and the RMSE values calculated for each of the models are shown in each of the figures.

Figure 11 .
Figure 11.Associating the PARASOL estimations of shape of the scattering phase function at each pixel to the NWP-model-predicted field of RH i .(a) The estimated shape of the scattering phase function; the yellow squares are as previously defined in Fig. 9.The brown squares represent those PARASOL pixels where no phase function model could be assigned, and the blue squares represent those pixels where phase function models assuming distortion values of between 0 and 0.25 could be assigned.(b) The NWP-modelpredicted cloud-top RH i field, where the colour bar indicates the range in predicted RH i .

Figure 12 .
Figure 12.The percent (%) probability of the penetration depth of solar radiation at 0.865 µm as a function of distance from the cloud top (km), and cloud optical depth for (a) forward-scattered and (b) backward-scattered solar radiation in the principal plane, respectively.The percent probability of penetration is defined as the last position (distance from the cloud top) of the photon before leaving the cloud to reach the sensor.The cloud optical depth colour scale is defined by the key shown on the upper right-hand side of (a).

Table 3 .
Same definitions as Table2but with the Levene test statistic applied to a group of seven pixels, where the fully randomised model phase function was found to best fit spherical albedo differences using minimised RMSE values.The model pair tests are between all other scattering phase function models and the fully randomised scattering phase function model.