A study of optical scattering modelling for mixed phase Polar Stratospheric Clouds

. Scattering codes are used to study the optical properties of Polar Stratospheric Clouds (PSC). Backscattering and extinction can be computed with available scattering codes once the particle size distribution (PSD) is known and a suitable refractive index is assumed. However, PSCs often appear as external mixtures of Supercooled Ternary Solution (STS) droplets, solid Nitric Acid Trihydrate (NAT) and possibly ice particles, making questionable the assumption of a single refractive index and a single morphology to model the scattarers. Here we consider a set of ﬁfteen coincident measurements of PSC above 5 McMurdo Station, Antarctica, by ground-based lidar and balloon-borne Optical Particle Counters (OPC), and in situ observations taken by a laser backscattersonde and an OPC during four balloon stratospheric ﬂights from Kiruna, Sweden. This unique dataset of microphysical and optical observations allows to test the performances of optical scattering models when both spherical and aspherical scatterers of different composition and, possibly, shapes are present. Here we consider particles as STS if their radius is below a certain threshold value R th and NAT or possibly ice if above 10 it. The refractive indices are assumed known from the literature. Mie scattering is used for the STS, assumed spherical, while scattering from NAT particles, considered as spheroids of different Aspect Ratio (AR), is treated with T-Matrix results where applicable, and of geometric-optics-integral-equation approach where the particle size parameter is too large to allow for a convergence of the T matrix method. The parameters R th and AR of our model are chosen to provide the best match with the observed optical backscattering and depolarization. The comparison of the calculations with the measures is satisfactory for the 15 backscattering but not so for the depolarization, and possible causes are discussed. The results of this work help to understand the limits of the application of these scattering theories in modeling the optical response of particles of different composition and morphology

as is in the case of extinction (David et al., 2012;Daerden et al., 2007). In PSC studies, Mie theory was also sometimes used to simulate aspheric particulate backscattering, employing corrections that take into account its reduction due to the asphericity of the scatterers Cairo et al., 2022). However, the Mie theory cannot produce depolarized backscattering and attempts with simple empirical models to use it to mimic depolarization from aspherical particles does a very poor job in reproducing the measured depolarization (Cairo et al., 2022).
The aim of our study is to employ concomitant microphysical and optical measurements of PSCs when both liquid and 65 solid particles are present and compare these with optical scattering computations done with codes capable of reproducing depolarization in backscattering. This allows us to verify the capacities and limits of these codes, tested on a relevant data set which provides both microphysical and optical measurements of mixed-phase particulate distributions. The methodology illustrated is not restricted to the study of mixed phase PSCs, but can find applications in all those cases in which the aerosol appears as an external mixture of solid and liquid particles, distinguishable on the basis of their different typical sizes.

Dataset
We have analysed fifteen measurements of PSC optical parameters observed by a groundbased lidar, coincident with microphysical observations acquired by a balloon-borne Optical Particle Counter (OPC). These observations were taken above McMurdo Station, Antarctica, between 1994. In addition we analysed four in situ balloonborne 75 observations carried out by a laser backscattersonde and an OPC during 4 stratospheric balloon flights from Kiruna, Sweden, between 2000(Weisser et al., 2006. The Antarctic lidar and balloonborne OPC dataset has been extensively discussed in Snels et al. (2021) where it has been used to provide empirical relations linking particle Surface Area and Volume densities with the backscatter coefficients. The main characteristics of the instrumentation will be here only briefly recalled.
The lidar observations have been provided by a system detecting 532 nm backscattered light with parallel and perpendicular 80 polarization with respect to the linear polarization of the emitting laser (Di Donfrancesco et al., 2000), thus allowing the measurement of Backscatter Ratio BR, Volume Depolarization δ and Aerosol Depolarization δ A from 10 to 23 km. These optical parameters follows the usual definitions (Cairo et al., 1999); in the following the subscripts mol and A denote respectively the molecular and particle contribution to the optical coefficients, and cross and par denote the perpendicular and parallel polarization of the Backscatter Coefficient β (Collis and Russell, 1976 following and is here introduced as: being: The formulas for switching from one to the other can be found in Cairo et al. (1999).
The BR is retrieved using the Klett algorithm where the attenuation correction follows Gobbi (1995). The δ is calibrated with the method described in Snels et al. (2009). Uncertainties in BR are estimated to be 5% while uncertainty in δ A is about 10-15% (Adriani et al., 2004). Typical measurements are 30-60 min integration over the signal, and the vertical resolution 100 is 75 m in 1994 and 1995 and 225 m in the other years. The OPC, which makes 10 second measurements, corresponding to roughly 50 meter vertical resolution, has been averaged to 250 meter bins. For comparison with the lidar the OPC has been interpolated onto the vertical grid of the lidar.
The OPC is described in Hofmann and Deshler (1991) and Deshler et al. (2003a). A throughful revision of its dataset is presented in Deshler et al. (2019). The instrument uses white light to measure scattering at 40 • in the forward direction 105 from particles passing through a dark field microscope. Mie theory and a model of the OPC response function are used to determine particle size throughout the range from 0.19 to 10.0 µm radius. The OPC provides time series of size resolved particle concentration histograms at 8 to 12 size bins, depending on the instrument used. A measurement of total concentration of particles is simultaneously determined by a CN counter close to the OPC, which grows all particles larger than 0.01 µm to an optically detectable size and counts them (Campbell and Deshler, 2014). Particle size histograms are fitted to unimodal or 110 bimodal lognormal size distributions, which are the representation of size distribution used in this work.
A series of balloon launches were carried out from Kiruna, Sweden in the early 2000s. The payload included, among other instruments and in addition to the OPC and the CN, a backscattersonde capable of measuring in situ backscattering and depolarization at 532 nm with 10 seconds time resolution. Details of such instrument are presented in Adriani et al. (1999) and Buontempo et al. (2009). Here we use data from four flights that took place on 19 January 2000 (Voigt et al., 2003), 9 December 115 2001 (Deshler et al., 2003b), 4 and 6 December 2002 (Larsen et al., 2004;Weisser et al., 2006). As these balloon flights were not simple ascents like the antarctic ones, but were commanded to perform altitude changes, by deflating the balloon or releasing ballast, to maximize the transit time in the detected PSCs, the backscattersonde data have been interpolated to the OPC 250 m average of the data, corrisponding to a time grid spacing of 60 or 75 seconds depending on flight.
In our study each data point includes the values of BR, δ A and of the PSD defined by the three or six parameters of a mono 120 or bimodal lognormal distributions. Altitude, pressure and temperature are also present as ancillary data. We identify a data to select the presence of mixed phase clouds, we require that β cross A is greater than 5 10 −6 km −1 sr −1 . Under these condition, a total of 141 data point from the Antarctic and 332 data point from the Arctic flights have been selected for the study.

125
While the nucleation of NAT and ice is a threshold process, STS particles can grow upon cooling from the the ubiquitous liquid Stratospheric Sulphate Aerosol (SSA), by continuously taking up nitric acid and water from the gas phase. They form droplets with volumes varying with temperatures but nevertheless with larger number density and smaller particle dimensions than NAT. Conversely, NAT particles are expected to be of smaller number density, but with dimension that can easily grow larger than the average STS particle radius of a few tenths of µm ( In our study we make use of the fact that in a mixed phase PSC the large particles are likely NAT or ice, solid particles which depolarize the backscatter light, while the small ones are liquid STS, that is spherical, and do not depolarize the backscattering.
135 Figure 1 shows our dataset mapped in terms of the BR -as 1-1/BR for the sake of plot readability -and δ A and colour coded with respect to the fraction of particles larger than 1 µm in radius, with respect to the total number of particles. As can be seen, it is a general feature that for each BR value, higher depolarization values are associated with higher fractions of large particles. Noteworthy, at high BR, high fractions of large particles are associated with depolarizations which, albeit high, are smaller than those observed at medium-low BR. In fact those depolarizations at medium-low BR are associated with lower 140 ratios of big to small particles. In other words, although at medium-low BR the large particles are proportionally less numerous than at high BR, the depolarization in the former case is higher. It has also to be noted in the plot that, in the high BR range at 1-1/BR=0.8, there are a few cases where high depolarization is associated with a low fracton of large particles. In these cases, although the relative abundance of large particles was low, the particles were unusually large, exceeding few µm. These very large particles, although in small concentrations, are causing the high depolarization observed. The temperature in these few 145 cases was however high enough to exclude them as ice particles.
We plan to consider particles as STS when their radius is below a threshold value R th and NAT above it. We use values of 1.44, and 1.48 for the refractive index of, respectively, STS, and NAT. These values are compatible with the large PSC data set produced by the CALIPSO observations Pitts et al., 2018a) and fall within the estimates presented for STS and NAT (Adriani et al., 1995;Deshler et al., 2000;Scarchilli et al., 2005). For completeness, ice particles are considered 150 when radii are larger than 4 µm and temperatures fall below 185 K. This happens only in 10% of the total dataset. In those few cases a value of 1.31 is used (Kokhanovsky, 2004).
For each data point, we split the P SD into two branches, namely P SD ST S = P SD(r < R th ) and P SD asph = P SD(r > R th ). As stated, the presence of ice particles is taken into consideration by inspecting the temperature T observed at the measurement, so if T > T ICE we pose P SD N AT = P SD asph (r > R th ), while if T < T ICE we limit the presence of NAT particles at radii smaller than 4µm, i.e. P SD N AT = P SD asph (R th < r < 4µm ) and consider as ice the particles with bigger radii, P SD ICE = P SD asph (r > 4µm).
The backscattering coefficients and depolarization ratio for the STS, NAT and ice particles are separately computed. For STS we have used a Mie code (Bohren and Huffman, 2008), available from the NASA's OceanColor Web site. For NAT and ice we have used the GRASP (Generalized Retrieval of Aerosol and Surface Properties) Spheroid Package. GRASP is the first 160 unified algorithm developed for characterizing atmospheric properties gathered from a variety of remote sensing observations (Dubovik et al., 2014), whose software packages are available on the Web. The Spheroid Package allows fast, fairly accurate, and flexible modelling of single scattering properties by randomly oriented spheroids with different size and shape. It includes a software and kernels data base. The details of the scientific concept are described in the paper by Dubovik et al. (2006).
The kernel look-up tables include results of calculations using T-Matrix code where convergence was acquired, and when 165 convergence was not achieved, the geometric-optics-integral-equation approach Liou, 1996b, 2006) was used, that is expected to provide accurate optical characteristics for spheroids with size parameter larger than 30 -40. The two methods have been shown to produce comparable results over the size range in which both are applicable (Dubovik et al., 2006). Thus, the software and kernels data base provide the kernel matrices for randomly oriented spheroids with Aspect Ratios (AR) from ∼ 0.3 (flattened spheroids) to ∼ 3.0 (elongated spheroids) and covering the size parameter range from ∼ 0.012 to ∼ 625 (when 170 a wavelenght of 0.44 is used) for a wide range of particle refractive index.
The total particle backscattering from the particles of the PSD can thus be parametrized with R th and AR and written as: Here we have used for the sake of simplicity the same AR for NAT and ice. Even if this hypothesis were not fully verified, this should not impact severely our study, as only 10% of our observations have temperatures below 185 K, and no definite ice 175 observations could be clearly discerned in our database. Now, the particle depolarization δ A can be parametrized in terms of AR and R th as well and be written as a weighted average of the contributions from the different classes of particles: It is useful here to recap how scattering from aspheric particles changes, compared to Mie theory, when T-Matrix is used.
Obviously, Mie theory cannot reproduce the depolarized backscattering typical of aspherical scatterers. According to T-Matrix, 180 the depolarization is negligible for scatterers size parameters lower than unity (given the waveleght used in our study, this corresponds to particle radius approximately below 0.1 µm), it grows up to a maximum that is reached for size parameters of the order of ten (i.e. particle radius around 1 µm) and then decreases towards an asymptotic value for size parameters greater than one hundred (particle radius greater than 10 µm). Both the maximum value and the asymptotic value vary according to the AR considered. In particular, the asymptotic value of the depolarization can assume values from 10% to 40%, and the 185 maximum value from 30% to 80%, the two variabilities are not connected. The dependence of the single particle depolarization on shape and size has been studied extensively by Liu and Mishchenko (2001). It has to be stressed that there is no simple relationship that binds the peak and asymptotic depolarization values to the AR of the particle, although there is a general tendency for small AR values to give large asymptotic depolarization values, and for large AR values to produce medium asymptotic depolarization values. The backscattering itself reproduces the Mie results for size parameters below unity, then 190 is progressively reduced to values that can even be one third of the Mie value when the particle size parameters is above ten, again this reduction depends in no simple way on the value of AR.

Variability with the threshold radius R th and Aspect Ratio AR
We have computed β A (R th ; AR) and δ A (R th ; AR) for a set of threshold radii R th and Aspect Ratios AR ranging respectively from 0.25 to 2 µm and from 0.3 to 3. To find the values that provide the best match with the measured ones β meas A and δ meas A , 195 we have calculated the respective Root Mean Squared Errors (RMSE), as: where the index i runs over our dataset. Covariance has also been computed resulting to be everywhere close to zero, except for R th smaller than 0.5 µm, where it resulted positive for AR between 0.5 and 1.5 (and close to unity for AR=1.25), while it was 200 slightly negative (correlation between -0.3 and 0) for AR smaller than 0.5 and greater than 1.5. Since the range of variability of β A is two orders of magnitude, to estimate the goodness of the agreement independently of the magnitude of β A , we have also calculated the Root Mean Squared relative Error (ReRMSE), defined as:

Results
The plots reported in figure 2 and 3 provide suggestions to select the choice of R th and AR to provide the best match between computations and observations. The comparison between the RMSE and RreRMSE for β A shows similar features: both plots 210 suggest to avoid AR values between 0.6 and 1.5 when an R th below 1 µm is considered, and show similar minimum differences between model and observations for AR below 0.5 and above 1.5. In this range of variability for AR, the R th s that reach the best agreement between model and observations are between 0.25 and 1 µm. The analysis of the RMSE and ReRMSE for δ A clearly shows regions where the agreement model-experiment is really poor, namely for AR below 0.75 and above 1.25 when R th is below 0.75 µm. Noteworthy, these regions only partially coincide with those of disagreement for β A . However, for δ A in relative terms the agreement is unsatisfactory practically everywhere, the relative error being considerably higher than in the case of the comparison of β A . Regions in which the agreement seems better are those with R th greater than 0.75 µm and AR values below 0.75 or greater than 1.5.
If we look for a compromise that ensures at the same time a good accord between model and data for both β A and δ A , the analysis do suggest to maintain R th between 0.75 and 1.25 µm, and keep the AR above 1.5, or below 0.5. Given those 220 constrains, even if an unambiguous choice is not possible, it seems reasonable to pose R th = 1.0 µm and test the two cases for oblate or prolate spheroids, posing AR = 2 .0 and 0.5 respectively, to test the goodness of the agreement. These values meet what we can also expect from what we know about PSC microphysics (Voigt et al., 2000) and about the comparison between backscattering from spherical and scattering from aspherical vs spherical particles (Mishchenko, 2009). With this choice we can present in figures 4 and 5 the scatterplot β A and δ A measured vs computed. The upper panels represent the oblate case, the 225 lower the prolate one. Figure 4 reporting the scatterplot of measured vs modelled β A , colour coded in terms of δ A , represent the analogue of figure   4 in Snels et al. (2021), where in the present case we have used a larger dataset, including Arctic balloon flights, and used T-Matrix instead of a factor 0.5 reduction, common to all particle sizes, for the Mie backscattering computation. In the present 230 case, both for prolate and oblate spheroids, the model is able to reproduce the observations for β A . While the agreement between calculations and measurements can be considered just fine for β A with the exception β < 4 · 10 −5 km −1 sr −1 where the calculations underestimate the measurements, it is not so for δ A as can be clearly seen from figure 5 where the points are colour coded in terms of β A . In that case, although there is some tendency to align the points along the diagonal, the dispersion of the points is very large. We can certainly invoke a noisiness in the data, which is present in the δ A measurements 235 to a larger extent than in the β A . Moreover, another critical aspect of this analysis is that aspherical particles passing through the OPC will scatter differently than a sphere of the same size, assumed in the OPC retrieval, and this causes an additional uncertainty on the particle true size. However it is likely that the lack of agreement has more profound reasons and that the assumptions made in the construction of our optical model are not fully reflected by reality. It is possible that spheroids are not able to fully replicate the scattering properties of PSC particles. Um and McFarquhar (2011) used geometric ray-tracing 240 codes on ice particles of linear dimension of few tens of µm and showed that differences in the backward scattering between different shape models (Chebyshev particles, Gaussian random spheres, and droxtals) are higher than 100%. Moreover it is possible that the assumption of the same value for the AR for each particle, irrespective of its composition or size, can not hold; on the contrary is quite possible that particle shape changes with the size due to the differential condensational growth along preferred dimensions. In fact, laboratory studies (Grothe et al., 2006), have showed that synthesized aspherical NAT 245 develops different morphologies depending on the growth conditions. The situation is even more complicated for ice particles, of which the extreme variety of shape is known, even in the case of pristine crystals whose morphology is influenced by the air temperature, atmospheric pressure, ventilation and ice supersaturation (Heymsfield et al., 2017). We have seen that different shapes produce different polarization, according to T-Matrix. This has also been proven experimentally since the early works of Sassen and Hsueh (1998) and Freudenthaler et al. (1996) that showed how the lidar 250 depolarization ratio in persisting contrails ranged from 10% to 70%, depending on the stage of their growth and on temperature.

Discussion
The oversimplified assumption of a common AR for all solid particles is one probable reason why the optical model behaves so poorly in the simulation of depolarization.
To further investigate the causes of this mismatch we turn to the study of the climatology of PSC observations collected from 255 McMurdo's lidar. The measurement of a PSC composed exclusively of solid particles is a rare and above all uncertain event, as we can never be completely certain of the absence of liquid aerosols. However, Adachi et al. (2001)

demonstrated that in a plot
showing the Total Volume Depolarization δ T towards 1-1/BR, the experimental points of solid, liquid or variously mixed PSCs are distributed within a triangle whose vertices are 0, δ asph T A and δ mol . These vertexes represent respectively the value of δ T in the case of pure liquid clouds and pure solid clouds for BR = ∞, when the δ T coincides with δ T A , and in the case when no 260 particles are present i.e. when δ T attains its molecular value δ mol (Young, 1980). Hence the extrapolated intercept on the y axis at BR = ∞ is precisely δ asph T A . This procedure allows us to estimate this asimptotic value. This requires the assumption that the experimental points that fill the triangle of vertices defined above represent PSC observations in mixed phase in which all solid particles share the same aerosol depolarization. Alternatively, one can interpret differently the presence of the datapoints filling the triangle. These points may as well represent single phase PSC of solid particles but with different shapes, hence producing 265 various and different depolarizations.  (Adriani et al., 2004). Despite the dispersion of the experimental points, a value close to 40% for δ asph T A can be tentatively assumed -the corrisponding value for δ A being 67% according to (6).
If we assume that the failure of our model to reproduce the observed depolarization is due to the incorrect assumption of a 270 common AR for all solid particles, we are led to interpret figure 6 admitting that the experimental points filling the triangle of vertices 0, δ asph T A and δ mol represent both PSC in various degrees of mixed phase, and PSC in purely solid phase but composed of particles of various and unidentical shapes. These various shapes give rise to different δ T A . Thus the found value of δ asph T A at the vertex of the triangle, it must not be considered as the univocal value of the depolarization for a PSC in a totally solid phase, but as the maximum of the possible depolarizations obtainable from solid phase clouds, depending on the many possible 275 shapes that the particulate can assume.
Despite the failure to correctly reproduce the observed depolarization, the simulation of the backscattering gives a good indication on the size transition at which cloud particles may begin to be considered NATs, confirming the common wisdom on the microphysics of mixed NAT/STS clouds.
We have used an optical model to compute with T-Matrix theory the backscatter and depolarization of mixed phase PSC. The model assumes that: i. PSC particles are solid (NAT or possibly ice) above a threshold radius R th , liquid (STS) below; ii. A single shape is common to all solid partiles irrespective of their size or composition. We have tested the model using a data set of coincident lidar, backscattersonde and OPC measurements from Antarctica and Arctic balloon flights.
The analysis suggested the range of optimal R th and AR parameters that best mach the observations. While the agreement 285 between modeled and measured backscatter is reasonable, it is poor for depolarization. Albeit it cannot be excluded a defective accuracy in depolarization data that prevent to fully demonstrate the validity of our model, the most likely reason is an oversimplification in the hypothesis of a common shape for all the solid particles present in the size distribution. However, the result for the backscatter coefficient simulation allows to constrain the model parameters R th and AR into reasonable ranges and should be regarded as the positive result of the study.

290
The illustrated method can find application in other cases of externally mixed aerosols with coexistence of different phases, where the phases can be distinguished on the basis of the particle size.