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 (PSCs). Particle backscattering and depolarization coefficients 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 the assumption of a single refractive index and a single morphology to model the scatterers questionable. Here we consider a set of 15 coincident measurements of PSCs above McMurdo Station, Antarctica, using ground-based lidar, a balloon-borne optical particle counter (OPC) and in situ observations taken by a laser backscattersonde and OPC during four balloon stratospheric flights from Kiruna, Sweden. This unique dataset of microphysical and optical observations allows us to test the performances of optical scattering models when both spherical and aspherical scatterers of different composition and, possibly, shapes are present.

We consider particles as STS if their radius is below a certain threshold value R_th and NAT or possibly ice if it is above it. The refractive indices are assumed known from the literature. Mie scattering is used for the STS, assumed spherical. Scattering from NAT particles, considered spheroids of different aspect ratio (AR), is treated with T-matrix results where applicable. The geometric-optics–integral-equation approach is used whenever 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 have been varied between 0.1 and 2 µm and between 0.3 and 3, respectively, and the calculated backscattering coefficient and depolarization were compared with the observed ones. The best agreement was found for R_th between 0.5 and 0.8 µm and for AR less than 0.55 and greater than 1.5. To further constrain the variability of AR within the identified intervals, we have sought an agreement with the experimental data by varying AR on a case-by-case basis and further optimizing the agreement by a proper choice of AR smaller than 0.55 and greater than 1.5 and R_th within the interval 0.5 and 0.8 µm. The ARs identified in this way cluster around the values 0.5 and 2.5. The comparison of the calculations with the measurements is presented and discussed. The results of this work help to set limits to the variability of the dimensions and asphericity of PSC solid particles, within the limits of applicability of our model based on the T-matrix theory of scattering and on assumptions on a common particle shape in a PSD and a common threshold radius for all the PSDs.

1 Introduction

Polar stratospheric clouds (PSCs) appear in the polar stratosphere during winter due to the very low temperatures and the dynamic isolation of air within the polar stratospheric vortex. They are important because of the two roles they play in polar stratospheric ozone depletion: providing surfaces for the heterogeneous reactions that lead to the reactivation of chlorine and decreasing the concentration of HNO₃ in the gaseous phase, thus altering the balance of the chlorine activation–deactivation cycles (Solomon, 1988). A comprehensive review of studies and knowledge acquired can be found in Tritscher et al. (2021).

PSCs can either be formed of liquid droplets composed of supercooled ternary solutions (STS) of sulfuric acid, nitric acid and water or solid nitric acid trihydrates (NATs), the thermodynamically stable form of HNO₃ and H₂O in the polar stratosphere, or possibly – when temperature is low enough – ice. Initially PSCs were classified as three types based on lidar measurements of the intensity of the backscattered light and the amount of depolarization of the returned signals (Browell et al., 1990). With the accumulation of observations, it has been realized that it is not common to observe a PSC of a single, well-defined type (Pitts et al., 2018). More often PSCs appear as external mixtures of liquid STS droplets, NAT, and/or ice, depending on the thermal history that led to their formation. For example, it is believed that the nucleation of NAT could start in droplets of a pre-existing population of STS, but not all liquid droplets may convert into solid NATs, allowing for the coexistence of particles of different composition and phases and thus of intermediate optical characteristics (Peter and Grooß, 2012).

The existence of multi-phase PSCs with particles of different shapes and composition, hence different particle refractive indices, makes the modelling of the scattering characteristics of the cloud problematic. While Mie scattering theory has been used to analyse PSCs consisting of spherical particles (Toon et al., 2000), detailed analysis of observations of different classes of PSCs (Deshler et al., 2000) may be questionable because they may consist of non-spherical solid particles with sizes comparable to the wavelength of the laser, and hence the results may be hampered by biases due to the unverified assumption of spherical scatterers. Because of this, some studies have chosen to limit themselves to liquid clouds only (Jumelet et al., 2008) or to make an effort and use theoretical modelling of light scattering for aspherical scatterers in the analysis. A viable solution which is not so demanding in terms of computational effort is the use of the T-matrix theory (see Liu and Mishchenko, 2001, and references therein). The T-matrix method is an exact technique for the computation of aspherical scattering based on a direct solution of Maxwell's equations assuming homogeneous, rotationally symmetric non-spherical particles or clusters of spheres (Mishchenko et al., 1996), and T-matrix codes are orders of magnitude faster than other approaches used in particle light scattering, like the discrete dipole approximation (Singham and Salzman, 1986) and the finite difference time domain (Yang and Liou, 1996a) techniques. The T-matrix approach to compare microphysical and optical observations has been used in a number of cases (Voigt et al., 2003; Scarchilli et al., 2005), under the assumption of particles as prolate or oblate spheroids. However even this approach could be controversial given that, just as solid particles cannot be modelled as spheres, for similar reasons it is debatable that they can be modelled as spheroids, and biases can arise under that assumption as well. For instance, Reichardt et al. (2002) showed how under the hypothesis of spheroidal particle shapes, the surface area density and volume density of lee wave PSCs are systematically smaller by, respectively, 10 %–30 % and 5 %–25 % than the values found for mixtures of droplets, asymmetric polyhedra and hexagons. Furthermore, there is no clear indication on what kind of aspect ratio can be unequivocally assumed for the spheroidal case (Reichardt et al., 2004; Engel et al., 2013; Woiwode et al., 2016). Finally, T-matrix codes suffer from convergence problems beyond a few dozen size parameters, especially for extreme aspect ratios (Mishchenko and Travis, 1998).

Given that there is still no completely satisfactory solution to tackle scattering from solid PSC particles due to the ambiguities that still persist about their shape, the use of Mie theory continues to be attractive and is used widely despite the unverified hypothesis of spherical scatterers. The speed of computation it offers is advantageous when used for statistical or climatological studies (Pitts et al., 2009, 2013) or when optical properties less dependent on the assumption of sphericity have to be calculated, as is 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 (Snels et al., 2021; 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 do a very poor job in reproducing the measured depolarization (Cairo et al., 2022).

In our study we employ concomitant microphysical (i.e. particle size distributions, PSDs) and optical (i.e particle backscattering and depolarization coefficients) measurements of PSCs when both liquid and solid particles are present and compare these with optical scattering computations done with codes capable of reproducing depolarization in backscattering. We consider particles as spherical STS if their radius is below a certain threshold value R_th and NAT or possibly ice if above it, the respective refractive indices being known from the literature. Mie scattering is used for the spherical part of the PSD. Scattering from NAT particles, considered spheroids of different aspect ratio (AR), is treated with T-matrix results where applicable. For particle size parameters too large to allow for T-matrix convergence, the geometric-optics–integral-equation approach is used. The parameters R_th and AR of our model have been varied between 0.1 and 2 µm and between 0.3 and 3, respectively, and the calculated backscattering coefficient and depolarization were compared with the observed ones in order to find the R_th and AR providing the best agreement between the computation and the measurements.

The aims of this effort are both to verify the ability of the T-matrix approach to reproduce the observations from lidar/backscattersonde, once the PSDs are supposed known, and to provide a contribution to the estimation of the shape and size limits of the NAT PSC particles.

The question of the shape of NAT particles is in fact far from being clarified and has important implications for the denitrification mechanisms of the polar stratosphere, an important step in the process that leads to the destruction of stratospheric ozone. In fact, the sedimentation of large NAT particles is considered one of the main causes of denitrification of the polar winter stratosphere (Di Liberto et al., 2015). Their settling time influences this process, and it is in turn dependent on the NAT particles' shape and size, both of which determine their settling speed and lifetime and hence their denitrification efficacy. Woiwode et al. (2014) assumed significantly non-spherical NAT particles to simulate the NAT settling speed leading to a the vertical redistribution of HNO₃ observed between two companion flights during the RECONCILE airborne field campaign in the Arctic (von Hobe et al., 2013). Woiwode et al. (2016, 2019) have also suggested that NAT particles may be highly aspherical based on the infrared spectrometer MIPAS-STR limb observations, exhibiting a spectral signature around 820 cm⁻¹ and an overall spectral pattern compatible with large highly aspherical NAT particles. T-matrix calculations assuming randomly oriented highly aspherical NAT particles (aspect ratios 0.1 or 10 for elongated or disk-like spheroids, respectively) were able to reproduce the MIPAS-STR observations to a large degree. Molleker et al. (2014) hypothesized strongly aspherical NAT particles to reconcile the amount of the condensed HNO₃ resulting from PSC cloud spectrometer measurements with the expected stratospheric values and to provide consistency between particles settling velocities and growth times with back trajectories. Moreover, Grothe et al. (2006) observed highly aspherical NAT in laboratory experiments. This is in contrast with earlier studies that assumed an AR = 0.9 for the NAT spheroids to match microphysical model simulation with airborne (Carslaw et al., 1998) or satellite-borne (Hoyle et al., 2013; Engel et al., 2013) lidar observations.

Finally, the methodology illustrated in the present work 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.

2 Data and methods

2.1 Dataset

We have analysed microphysical observations acquired by a balloon-borne optical particle counter (OPC) from 15 Antarctic balloon flights, coincident with measurements of PSC backscattering and depolarization coefficients observed by a ground-based lidar. These observations were taken above McMurdo Station, Antarctica, between 1994 and 1999 (Snels et al., 2021). In addition we analysed four in situ balloon-borne observations carried out by a laser backscattersonde and an OPC during four stratospheric balloon flights from Kiruna, Sweden, between 2000 and 2002 (Weisser et al., 2006). The Antarctic lidar and balloon-borne 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 backscattering coefficients. The main characteristics of the instrumentation will be only briefly recalled here.

The lidar observations have been provided by a system detecting 532 nm backscattered light with parallel and perpendicular polarization with respect to the linear polarization of the emitting laser (Di Donfrancesco et al., 2000), thus allowing for the measurement of backscatter ratio (BR), volume depolarization (δ) and aerosol depolarization (δ_A) from 10 to 23 km. These optical parameters follow 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 backscattering coefficient β (Collis and Russell, 1976).

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{BR}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{par}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{par}}}{{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{par}}}}\text{(2)}& {\displaystyle }\mathit{\delta }={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}}{{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{par}}+{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{par}}}}\text{(3)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{A}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}}{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{par}}}}\end{array}$$

An alternative definition of total volume depolarization δ_T and total aerosol depolarization δ_TA will also be used in the following and is introduced here as

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{T}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{cross}}}{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{cross}}+{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{par}}+{\mathit{\beta }}_{\mathrm{mol}}^{\mathrm{par}}}}\text{(5)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{TA}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}}{{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{par}}+{\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}}}\end{array}$$

being

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{\mathrm{T}}={\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }+\mathrm{1}}}.\end{array}$$

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). Experimental errors in the particle backscatter ratio (R−1) are estimated to be 5 % but not less than 0.05 in absolute value, while the error in volume depolarization is about 10 %–15 %. Additional uncertainty comes from the determination of pressure and temperature from radiosoundings, needed to compute β_A an δ_A (Adriani et al., 2004). Typical measurements are 30–60 min integration over the signal, and the vertical resolution is 75 m in 1994 and 1995 and 225 m in the other years. The OPC, which makes 10 s measurements, corresponding to roughly 50 m vertical resolution, has been averaged to 250 m 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 thorough 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 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 condensation nuclei counter (CNC) 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 bimodal lognormal size distributions, which are the representation of size distribution used in this work. The uncertainties on the determination of the parameters of the mono- or bimodal lognormal distributions were determined by Deshler et al. (2003b) with Monte Carlo simulations. These were 20 % for distribution width, 30 % for median radii and 10 % for modal particle concentrations.

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 an OPC and a CNC, a backscattersonde capable of measuring in situ backscattering and depolarization at 532 nm with 10 s time resolution. Details of these instruments 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 2001 (Deshler et al., 2003b), and 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, corresponding to a time grid spacing of 60 or 75 s depending on flight.

In our study each data point includes the values of BR, δ_A and the PSD defined by the three or six parameters of a mono- or bimodal lognormal distribution. Altitude, pressure and temperature are also present as ancillary data. We identify a data point as a PSC observation when the BR is greater than 1.2 and the temperature at the observation is below 200 K. Moreover, to select the presence of mixed-phase clouds, we require that ${\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}$ is greater than 5 × 10⁻⁶ km⁻¹ sr⁻¹. Under these conditions, a total of 141 data points from the Antarctic and 332 data points from the Arctic flights have been selected for the study.

2.2 Optical model

While the nucleation of NAT and ice is a threshold process, STS particles can grow upon cooling from the ubiquitous liquid stratospheric sulfate 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 a micrometre (Carslaw et al., 1997; Grooß et al., 2014) due to the smaller saturation vapour pressure of the nitric acid with respect to them. Deshler et al. (2003b) provide direct observations of this separation between STS and NAT. Due to the larger availability of water vapour, ice particles can grow even larger, often with linear dimensions exceeding 4–5 µm (Tritscher et al., 2019).

In our study we take advantage of the fact that in a mixed-phase PSC, the large particles are likely NAT or ice, i.e. solid particles which depolarize the backscattered light, while the small ones are liquid STS, that is spherical, and do not depolarize the backscattering. Figure 1 shows our dataset mapped in terms of the BR – as $\mathrm{1}-\mathrm{1}/\mathrm{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 connected with higher fractions of large particles. It is noteworthy that at high BR, high fractions of large particles are related with depolarizations, which, though they are high, are smaller than those observed at medium–low BR. In fact those depolarizations at medium–low BR are associated with lower 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 is higher at medium–low BR. It has also to be noted in the plot that, in the high BR range at $\mathrm{1}-\mathrm{1}/\mathrm{BR}=\mathrm{0.8}$, there are a few cases where high depolarization corresponds to a low fraction of large particles. In these cases, although the relative abundance of large particles was low, the particles were unusually large, exceeding a few micrometres. These very large particles, although in small concentrations, are causing the high depolarization observed. The temperature in these few cases was however high enough to exclude them as ice particles.

Figure 1 Scatterplot of aerosol depolarization vs. $\mathrm{1}-\mathrm{1}/\mathrm{BR}$, where BR is the backscatter ratio. Data are from McMurdo lidar and Kiruna balloon flights backscattersonde, coincident with balloon-borne OPC PSD measurements, fitted with mono- or bimodal lognormal distributions. The colour codes the fraction of particles with radius larger than 1 µm with respect to the total number of particles in the PSDs. We report data points with BR greater than 1.2, ${\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}$ greater than 5 × 10⁻⁶ km⁻¹ sr⁻¹ and temperature at the observation below 200 K.

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 dataset produced by the CALIPSO observations (Hoyle et al., 2013; Pitts et al., 2018) 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 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 PSD into two branches, namely ${\mathrm{PSD}}_{\mathrm{STS}}=\mathrm{PSD}(r<R_{\mathrm{th}})$ and ${\mathrm{PSD}}_{\mathrm{asph}}=\mathrm{PSD}(r>R_{\mathrm{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 ${\mathrm{PSD}}_{\mathrm{NAT}}={\mathrm{PSD}}_{\mathrm{asph}}(r>R_{\mathrm{th}})$, while if T<T_ICE we limit the presence of NAT particles at radii smaller than 4 µm, i.e. ${\mathrm{PSD}}_{\mathrm{NAT}}={\mathrm{PSD}}_{\mathrm{asph}}(R_{\mathrm{th}}<r<\mathrm{4}\mathrm{\mu }\mathrm{m})$, and consider the particles with bigger radii, ${\mathrm{PSD}}_{\mathrm{ICE}}={\mathrm{PSD}}_{\mathrm{asph}}(r>\mathrm{4}\mathrm{\mu }\mathrm{m})$, as ice.

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 website. For NAT and ice we have used the GRASP (Generalized Retrieval of Aerosol and Surface Properties) SPHEROID package. GRASP is the first 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 for fast, fairly accurate and flexible modelling of single scattering properties of randomly oriented spheroids with different size and shape. It includes a software and kernel database. The details of the scientific concept are described in Dubovik et al. (2006). The kernel lookup tables include results of calculations using T-matrix code where convergence was acquired, and when convergence was not achieved, the geometric-optics–integral-equation approach (Yang and Liou, 1996b, 2006) was used, which 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 kernel database provide the kernel matrices for randomly oriented spheroids with aspect ratios (ARs) from ∼ 0.3 (flattened spheroids) to ∼ 3.0 (elongated spheroids) and covering the size parameter range from ∼ 0.012 to ∼ 625 (when a wavelength of 0.44 µm 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

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\beta }}_{\mathrm{A}}& ={\mathit{\beta }}_{\mathrm{A}}(R_{\mathrm{th}};\mathrm{AR})\\ & ={\mathit{\beta }}_{\mathrm{STS}}\left(R_{\mathrm{th}}\right)+{\mathit{\beta }}_{\mathrm{NAT}}(R_{\mathrm{th}};\mathrm{AR})+{\mathit{\beta }}_{\mathrm{ICE}}(R_{\mathrm{th}};\mathrm{AR}).\end{array}\end{array}$$

Here we have used, for the sake of simplicity, the same AR for NAT and ice. Even if this hypothesis is not fully verified, this should not impact our study severely, as only 10 % of our observations have temperatures below 185 K, and no definite ice 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:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{A}}& ={\mathit{\delta }}_{\mathrm{A}}(R_{\mathrm{th}};\mathrm{AR})\\ & ={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{NAT}}(R_{\mathrm{th}};\mathrm{AR})\cdot {\mathit{\delta }}_{\mathrm{A},\mathrm{NAT}}(R_{\mathrm{th}};\mathrm{AR})+{\mathit{\beta }}_{\mathrm{ICE}}(R_{\mathrm{th}};\mathrm{AR})\cdot {\mathit{\delta }}_{\mathrm{A},\mathrm{ICE}}(R_{\mathrm{th}};\mathrm{AR})}{{\mathit{\beta }}_{\mathrm{A}}(R_{\mathrm{th}};\mathrm{AR})}}.\end{array}\end{array}$$

​​​​​​​It is useful here to recap how scattering from aspheric particles changes, compared to Mie theory, when the T-matrix theory is used. Obviously, Mie theory cannot reproduce the depolarized backscattering typical of aspherical scatterers. According to the T-matrix theory, the depolarization is negligible for scatterers with size parameters lower than unity (given the wavelength used in our study, this corresponds to a particle radius approximately below 0.1 µm); it grows to a maximum that is reached for size parameters of the order of 10 (i.e. particle radius around 1 µm) and then decreases towards an asymptotic value for size parameters greater than 100 (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 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 and then is progressively reduced to values that can even be one-third of the Mie value when the particle size parameter is above 10. Again this reduction depends in no simple way on the value of AR.

It is quite possible that in some cases, spheroids are not able to fully replicate the scattering properties of PSC particles. Um and McFarquhar (2011) used geometric ray-tracing codes on ice particles of linear dimension of a few tens of micrometres and showed that differences in the backward scattering between different shape models (Chebyshev particles, Gaussian random spheres and droxtals) are higher than 100 %. Furthermore, it is possible that the assumption of the same AR for each particle in the solid part of the PSD, irrespective of its composition or size, can not hold; on the contrary it is quite possible that particle shape changes with size due to the differential condensation growth along preferred dimensions. Laboratory studies (Grothe et al., 2006) have shown that synthesized aspherical NAT develops different morphologies depending on growth conditions. However, our approach maintains an AR common to all PSD particles to be considered the average AR of the PSD. This is a simplistic choice, but in contrast we are not able to justify a particle-size-dependent AR on physical grounds. An additional uncertainty on the measured particle size is that aspherical particles in the OPC will scatter differently than the spherical particles assumed in the OPC retrieval.

2.3 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 ARs ranging respectively from 0.1 to 2 µm and from 0.3 to 3. To find the values that provide the best match with the measured ones ${\mathit{\beta }}_{\mathrm{A}}^{\mathrm{meas}}$ and ${\mathit{\delta }}_{\mathrm{A}}^{\mathrm{meas}}$, we have calculated the respective root mean squared errors (RMSEs) as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{RMSE}}_{{\mathit{\beta }}_{\mathrm{A}}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\mathit{\beta }}_{\mathrm{A}i}-{\mathit{\beta }}_{\mathrm{A}i}^{\mathrm{meas}}{)}^{\mathrm{2}}}{n}}}\text{(10)}& {\displaystyle }{\mathrm{RMSE}}_{{\mathit{\delta }}_{\mathrm{A}}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\mathit{\delta }}_{\mathrm{A}i}-{\mathit{\delta }}_{\mathrm{A}i}^{\mathrm{meas}}{)}^{\mathrm{2}}}{n}}},\end{array}$$

where the index i runs over our dataset. Covariance has also been computed, everywhere resulting to be close to zero, except for R_th smaller than 0.5 µm, where it was positive for AR between 0.55 and 1.5 (and close to unity for AR = 1.25), and it was slightly negative (correlation between −0.3 and 0) for AR smaller than 0.55 and greater than 1.5. Since the range of variability of β_A is 2 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

$$\begin{array}{}\text{(11)}& {\displaystyle }{\mathrm{ReRMSE}}_{{\mathit{\beta }}_{\mathrm{A}}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{A}i}-{\mathit{\beta }}_{\mathrm{A}i}^{\mathrm{meas}}}{{\mathit{\beta }}_{\mathrm{A}i}^{\mathrm{meas}}}}\right)}^{\mathrm{2}}}\text{(12)}& {\displaystyle }{\mathrm{ReRMSE}}_{{\mathit{\delta }}_{\mathrm{A}}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{{\mathit{\delta }}_{\mathrm{A}i}-{\mathit{\delta }}_{\mathrm{A}i}^{\mathrm{meas}}}{{\mathit{\delta }}_{\mathrm{A}i}^{\mathrm{meas}}}}\right)}^{\mathrm{2}}}.\end{array}$$

The colour-coded RMSEs (upper panel) and ReRMSEs (lower panel), with respect to R_th and AR, are shown in Fig. 2 for backscattering and Fig. 3 for depolarization.

Figure 2 Contour plot of the RMSE (a) and relative RMSE (b) of the measured and modelled aerosol backscatter coefficient β_A.

Figure 3 Contour plot of the RMSE (a) and relative RMSE (b) of the measured and modelled aerosol depolarization δ_A.

The plots reported in Figs. 2 and 3 suggest how to select R_th and AR to provide the best match between computations and observations. The comparison between the RMSE and ReRMSE for β_A shows similar features: both plots suggest avoiding 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.55 and above 1.5. In this range of variability for AR, the R_ths that reach the best agreement between model and observations are between 0.3 and 1 µm, with 0.5–0.7 µm the most favourable. The analysis of the RMSE and ReRMSE for δ_A clearly shows regions where the model–experiment agreement is really poor, namely for AR below 0.75 and above 1.25 when R_th is below 0.5 µm. It is noteworthy that these regions only partially coincide with those of disagreement for β_A. The relative error is 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.5 µm and AR values below 0.75 or greater than 1.25. The result of this study allows us to identify only the best R_th, resulting around 0.5–0.8 µm, while the ARs compatible with the measurements are all those between 0.3–0.55 and 1.5–3.

To further constrain AR we have kept R_th at a fixed value, chosen between 0.5 and 0.8 µm, and changed this value with a 0.1 µm step. For each of these fixed R_th, and separately for each PSD, in the intervals (0.3–0.55) and (1.5, 3), we identified the value of AR which best matched the observed δ with its computed value. Finally, for each PSD we selected the R_th and AR pair which provided the best match. Once the ARs and R_th were selected by forcing the agreement between the δ_As, the same ones were used for the calculation of the β_As.

3 Results

Figure 4 reports the scatterplot of measured vs. computed β_A, colour-coded in terms of AR. The figure represents the analogue of Fig. 4 in Snels et al. (2021), where in the present case we have used a larger dataset, including now four Arctic balloon flights, and used the T-matrix approach instead of a factor 0.5 reduction in the Mie backscattering. Figure 5 reports the scatterplot of measured vs. computed δ_A, similarly colour-coded in terms of AR. The uncertainties associated with the measured β_A and δ_A derive from the error analysis for the single lidar data, which can be found in Adriani et al. (2004) or from the standard deviation for the averaged data, depending on which is greater. The uncertainties on the calculated β_A and δ_A are 40 %, as determined by Deshler et al. (2003a) for any moment of a PSD derived from the OPC measurements. Deshler et al. determined this through a Monte Carlo simulation, which used the uncertainties of the OPC size and concentration measurements to quantify the uncertainties in the PSD parameters and their subsequent moments.

Figure 4 Scatterplot of computed vs. measured particle backscattering coefficients β_A. The ARs used for the β computations have been selected, case by case, to produce the best agreement between the δ computed and measured values and are represented here by colour coding. Only ARs in the intervals between 0.3 and 0.55 and between 1.5 and 3, have been considered. R_th was also selected within the interval 0.5–0.8 µm to provide the best match. We report data points with BR greater than 1.2, ${\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}$ greater than 5 × 10⁻⁶ km⁻¹ sr⁻¹ and temperature at the observation below 200 K.

Figure 5 Scatterplot of computed vs. measured particle depolarization δ_A. The ARs used for the computations are those that provided the best match between the δ computed and measured values and are represented here by colour coding. Only ARs in the intervals between 0.3 and 0.55 and between 1.5 and 3 have been considered. R_th was also selected within the interval 0.5–0.8 µm to provide the best match. We report data points with BR greater than 1.2, ${\mathit{\beta }}_{\mathrm{A}}^{\mathrm{cross}}$ greater than 5 × 10⁻⁶ km⁻¹ sr⁻¹ and temperature at the observation below 200 K.

Despite the dispersion in Fig. 4, the points cluster around the straight line β^calc=β^meas, indicating the agreement between computation and measurements can be considered fine for β_A, with the exception of β values below 4 × 10⁻⁵ km⁻¹ sr⁻¹, where the β^calc values underestimate the measurements. Such underestimation seems to be of the order of 10⁻⁵ km⁻¹ sr⁻¹, of the same order as the backscattering from the background atmospheric particulate matter in volcanic quiescent conditions, a magnitude compatible with possible inaccuracies in the calibration of the lidar data. The Pearson correlation coefficient for the entire dataset is 0.56, and this increases if the lower values of β are neglected.

The δ_A scatterplot shows the presence of a good number of points that align along the δ^calc=δ^meas correlation line, with ARs selected mainly around the value 0.5. However, for depolarization values greater than 30 % there is no AR that will reproduce the measurements. These points correspond to those presented in Fig. 1, with low values of BR and high values of the concentration ratio of large to total particles. They mainly come from three single observational periods of about 1 min each, characterized by air temperatures between 184–188 K. Given the magnitude of the depolarization, it is possible that those observations are not referable to clouds in mixed phase but rather to clouds of predominantly solid particles. For that particular set of points, we also explored the possibility that all particles were solid, but even under this assumption the comparison with the experimental data did not improve appreciably.

In Fig. 5, for depolarizations lower than 15 %, the points which deviate, by excess or defect, from the 1 : 1 straight line predominantly have an AR greater than 1.5. So it seems that selected ARs greater than 1.5 generally produce a worse correlation. From Fig. 4 we observe that AR values in the range (0.3–0.55) tend to be associated with medium–low β values, while AR values in the range (1.5–3) are mainly associated with medium–high β.

To conclude, the choice of R_th in a range between 0.5 and 0.8 µm leads to a reasonably good agreement between the β's, but there seems to be a discrepancy between the calculated value and the measurements in their lower range of variability. From Fig. 4, such a mismatch, which makes the measurements larger than the calculations, seems to be of the order of 10⁻⁵ km⁻¹ sr⁻¹. The selection of the AR that produces the best agreement with the observed δ's leads to three results: (i) the ARs in the range 0.3–0.55 tend to be selected in correspondence with medium–low β's and the ARs in the range 1.5–3 in correspondence with medium–high β's; (ii) ARs in the 0.3–0.5 range reproduce the measurements well, except for some observations where the depolarizations are greater than 30 %; (iii) the ARs in the 1.5–3 range reproduce the measurements less well; (iv) there is no AR that will allow the calculations to reproduce the measurements for depolarizations greater than 30 %.

4 Discussion

The identification of the best R_th in the range 0.5–0.8 µm supports what we already know from the theoretical understanding of NAT particle formation in PSC and from measurements (Deshler et al., 2003b). Concerning particle shape, in our model all solid particles in a single PSD share the same AR, but different PSDs can have different ARs. This approach could suggest that the choice of the AR, which, case by case, optimizes the agreement between calculations and measurements, may be the result of chance rather than physics. There are two facts that counter this criticism.

First, it appears that the selected ARs may be related to the shape of the PSD. Figure 6 shows the 2D histogram by occurrence of ARs and of $N(r>\mathrm{0.7}\mathrm{\mu }\mathrm{m})/N_{\mathrm{tot}}$, the ratio between particles with radius greater than 0.7 µm and total particles, which is a parameter related to the PSD shape. In Fig. 6 the ARs are not distributed randomly. Conversely, there is a tendency for the AR to grow as the percentage of large particles increases. In fact AR values tend to peak around 0.5 in the lower $N(r>\mathrm{0.7}\mathrm{\mu }\mathrm{m})/N_{\mathrm{tot}}$ range, while they tend to cluster around 2.5 when $N(r>\mathrm{0.7}\mathrm{\mu }\mathrm{m})/N_{\mathrm{tot}}$ is higher. The shape of the PSD mirrors particle formation conditions and history is linked to the presence of solid particles, as already highlighted in the discussion of Fig. 1, and is likely linked to the average particle shape as well.

Figure 6 A 2D histogram of occurrence of ARs and of $N(r>\mathrm{0.7}\mathrm{\mu }\mathrm{m})/N_{\mathrm{tot}}$, the ratio between particles with radius greater than 0.7 µm and total particles. Only ARs in the intervals between 0.3 and 0.55 and between 1.5 and 3 have been considered.

Second, if we consider the sequences of measurements acquired in individual balloon flights, the corresponding sequences of selected ARs do not evolve randomly but, conversely, are auto-correlated. An example of this behaviour is shown in Fig. 7, where the time series of β and δ are reported respectively with red and blue dots. The ARs that provide the best agreement between the experiment and simulation are shown with black dots. It can be seen that temporally contiguous observations often result in the selection of the same AR. Contiguous observations of PSD are likely to have similar characteristics in terms of microphysics, and this seems to be correctly reflected in the constancy of AR. We are therefore confident that our method produces results with a physics-based content.

Figure 7 Sequences of β (red dots) and δ (blue dots) measured on a balloon flight on 9 December 2001, from Kiruna, Sweden. Each data point represents an average over 60 s. Black dots represent the ARs providing the best match between the δ and those computed from concomitant measurements of PSD.

In general our model leads to good correlations between measured and modelled βs. For the δs the measurements are well reproduced by the calculations in many instances, as is the case for many of the selected ARs in the range 0.3–0.55. However, there are other cases in which the agreement is worse (when the best ARs have been selected in the range 1.5–3) or does not occur at all, as in the cases of observed depolarizations greater than 30 %. In these latter cases, the impossibility of reproducing the observed values even under the hypothesis of a completely solid particles implies that, for those PSDs, our model is not able to produce the observed depolarizations. In these particular cases in which the model performs particularly badly, there may be problems of inhomogeneities of the cloud. These cases come from Antarctic observations, for which the microphysical observations from the balloon and the optical ones from ground-based lidar are separated geometrically, so that the two instruments sample air masses separated by several tens of kilometres, and it may be the case that some clouds were not homogeneous on such spatial scales.

Different shapes produce different polarization, according to the T-matrix theory. This has also been proven experimentally since the early work of Sassen and Hsueh (1998) and Freudenthaler et al. (1996) that showed how lidar depolarization ratios in persisting contrails ranged from 10 % to 70 %, depending on the stage of their growth and on temperature. In the T-matrix theory, for fixed AR, the depolarization depends on the particle size and maximizes for particular sizes. There is certainly a way to assume a particle-size-dependent AR in our PSDs so as to reconcile the computations with the observed values. However such an approach would have little physical basis and could only be justified to maximize the agreement of calculations. Therefore we have not explored this possibility further, although it is possible that our simplified hypothesis of a common AR for every particle may be the cause of the bad agreement between data and calculations in some case.

To further investigate the causes of the mismatch, we turn to the study of the climatology of PSC observations collected from McMurdo's lidar. The measurement of a PSC composed exclusively of solid particles is a rare and uncertain event. The absence of liquid aerosols is difficult to determine for certain. However, Adachi et al. (2001) demonstrated that in a plot of the total volume depolarization δ_T versus $\mathrm{1}-\mathrm{1}/\mathrm{BR}$, the experimental points of solid, liquid or variously mixed PSCs are distributed within a triangle whose vertices are (1,0), $\left(\mathrm{1},{\mathit{\delta }}_{\mathrm{TA}}^{\mathrm{asph}}\right)$ and (0,δ_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 δ_TA , and in the case when no particles are present, the δ_T attains its molecular value δ_mol (Young, 1980). Hence the extrapolated intercept on the y axis at BR=∞ is precisely ${\mathit{\delta }}_{\mathrm{TA}}^{\mathrm{asph}}$. This procedure allows us to estimate this asymptotic 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 the presence of the data points filling the triangle differently. These points may as well represent single phase PSC of solid particles but with different shapes, hence producing various depolarizations.

Figure 8 reports a 2D histogram of δ_TA towards $\mathrm{1}-\mathrm{1}/\mathrm{BR}$ from 12 years of lidar observations from 1990 to 2002 in the Antarctic station of McMurdo (Adriani et al., 2004). Despite the dispersion of the experimental points, a value close to 40 %, as the highest vertex of the triangle, on the $\mathrm{1}-\mathrm{1}/\mathrm{BR}=\mathrm{1}$ axis, for ${\mathit{\delta }}_{\mathrm{TA}}^{\mathrm{asph}}$ can be assumed. The corresponding value for δ_A is close to 70 % according to Eqs. (3), (4) and (5).

Figure 8 A 2D histogram of total volume depolarization δ_TA vs. $\mathrm{1}-\mathrm{1}/\mathrm{BR}$, where BR is the backscatter ratio. Data are from McMurdo lidar and cover the winters from 1990 to 2002. The colour codes the number of observations.

If we assume that the difficulty of our model to reproduce the observed depolarization in some case is due to the incorrect assumption of a common AR for all solid particles, we are led to interpret Fig. 8 admitting that the experimental points filling the triangle of vertices (1,0), $\left(\mathrm{1},{\mathit{\delta }}_{\mathrm{TA}}^{\mathrm{asph}}\right)$ and (0,δ_mol) represent both PSC in various degrees of mixed phase, and PSC in purely solid phase but composed of particles of various shapes. These various shapes give rise to different δ_TA between 0 and ${\mathit{\delta }}_{\mathrm{TA}}^{\mathrm{asph}}$ at the vertex of the triangle.

5 Conclusions

We have used an optical model to compute the backscattering and depolarization of mixed-phase PSCs with the T-matrix theory. The model assumes that (i.) PSC particles are solid (NAT or possibly ice) above a threshold radius R_th and liquid (STS) below and that (ii.) a single shape is common to all solid particles in a PSD, irrespective of their size or composition. We have tested the model using a dataset of coincident lidar, backscattersonde and OPC measurements from Antarctica and Arctic balloon flights.

While the agreement between modelled and measured backscattering coefficient is generally reasonable, there are cases in which it is less so for depolarization. The most likely reason is our simplified hypothesis of a common shape for all the solid particles present in the size distribution. However, our analysis has provided the range of optimal R_th and AR parameters that best match the observations. To sum up, (i) in an externally mixed PSC, it is reasonable to place a threshold radius R_th between 0.5 and 0.8 µm, which divides the liquid part from the solid part of the particulate; (ii) it is sensible to expect strongly aspherical shapes for the solid part of the cloud; (iii) there are cases, especially those related to observations of high depolarizations, in which such depolarizations are not reproducible with our T-matrix approach. This is possibly due to our assumptions (i.e. a single form of the solid particles and a fixed threshold radius for all PSDs).