# An Inverse Problem Approach for the Retrieval of Ice Particle Mass from In Situ Measurements

Pierre Coutris, Delphine Leroy, Emmanuel Fontaine, Alfons Schwarzenboeck

► To cite this version:

Pierre Coutris, Delphine Leroy, Emmanuel Fontaine, Alfons Schwarzenboeck. An Inverse Problem Approach for the Retrieval of Ice Particle Mass from In Situ Measurements. Journal of Atmospheric and Oceanic Technology, 2017, 34 (11), pp.2457-2473. 10.1175/JTECH-D-17-0013.1. http://dx.doi.org/10.1175/JTECH-D-17-0013.1

HAL Id: hal-01893713

https://uca.hal.science/hal-01893713

Submitted on 17 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Inverse Problem Approach for the Retrieval of Ice Particle Mass from In Situ Measurements

PIERRE COUTRIS, DELPHINE LEROY, EMMANUEL FONTAINE, AND ALFONS SCHWARZENBOECK

*Laboratoire de Météorologie Physique, Université Clermont Auvergne, CNRS, Clermont-Ferrand, France*

(Manuscript received 17 January 2017, in final form 19 July 2017)

## ABSTRACT

Mass–dimensional relationships (*m* − *D*) have been published for decades to characterize the microphysical properties of ice cloud particles. Classical *m* − *D* retrieval methods employ a simplifying assumption that restricts the form of the mass–dimensional relationship to a power law, an assumption that was proved inaccurate in recent studies. In this paper, a nonstandard approach that leverages optimal use of in situ measurements to remove the power-law constraint is presented. A model formulated as a linear system of equations relating ice particle mass to particle size distribution (PSD) and ice water content (IWC) is established, and the mass re-trieval process consists of solving the inverse problem with numerical optimization algorithms. First, the method is applied to a synthetic crystal dataset in order to validate the selected algorithms and to tune the regularization strategy. Subsequently, the method is applied to in situ measurements collected during the High Altitude Ice Crystal–High Ice Water Content field campaigns. Preliminary results confirm the method is efficient at retrieving size-dependent masses from real data despite a significant amount of noise: the IWC values calculated from the retrieved masses are in good agreement with reference IWC measurements (errors on the order of 10%–15%). The possibility to retrieve ice particle size–dependent masses combined with the flexibility left for sorting datasets as a function of parameters such as cloud temperature, cloud type, or convective index makes this approach well suited for studying ice cloud microphysical properties.

## 1. Introduction

### a. Motivation

The representation of ice particles, especially the parameterization of their microphysical properties, is an important part in atmospheric models and spaceborne sensors retrieval algorithms. It has long been acknowledged that ice particle mass, size, and related quantities such as density are key modules in ice microphysical schemes (Schmitt and Heymsfield 2009; Erfani and Mitchell 2016; Cotton et al. 2013), hence the acute need for improved representation of ice particles. Furthermore, interest in particle microphysical properties has continued to grow since high-altitude ice crystals were identified as a threat to commercial aviation, causing aircraft engine and air data probe icing. This led aircraft industry leaders, aviation regulatory authorities, and research facilities to join forces to launch the High Altitude Ice Crystal–High Ice Water Content (HAIC–HIWC) projects with the overall objective to enhance flight safety when flying in mixed and glaciated icing conditions. Technical developments of detection and awareness technologies and

modeling of ice accretion mechanisms within engine modules and air data probes require better characterization of ice particle mass, size, and concentration. These are the main rationale for the observation of ice particle physical properties and for the development of an empirical model that accurately describes particle mass as a function of size, temperature, and any other relevant physical parameters.

### b. State of the art

To address the abovementioned research topics, many studies have been conducted in the past with the primary goal to characterize cloud microphysical properties and to derive expressions for ice particle mass. In the following subsections, some state-of-the-art studies providing mass–size relationships (*m* − *D*) from particle measurements are reviewed, with emphasis placed on the hypotheses and analytical tools they use.

#### 1) EARLY STUDIES ON ICE CRYSTALS: COLLECTION ON GROUND

Early studies related to natural ice crystal mass focused on solid precipitating particles collected on ground (Nakaya and Terada 1935; Zikmunda and Vali 1972; Kajikawa 1972;

<footnote>Corresponding author: Pierre Coutris, p coutris@opgc.fr</footnote>
Locatelli and Hobbs 1974; and more recently Mitchell et al. 1990, among others). In these studies, ice particles ranging from hundreds of micrometers to several millimeters in size were collected, photographed, and classified into habit categories. The individual crystal mass was estimated from the droplet produced by crystal melting, and habit-dependent *m*–*D* relationships were derived for each subset. The retrieval process consists of fitting the individual mass against size data points using linear regression tools in log–log space, which produces power-law-type *m*–*D* relationships, noted as *m* = *a*D<sup>*b*</sup>, where the prefactor *a* and the exponent *b* are constant values that are specific to each habit. Since high-altitude ice cloud particles and precipitating particles may have very different microphysical properties, the *m*–*D* relationships derived in these studies may not be applicable to airborne in situ measurements.

## 2) IN SITU MEASUREMENTS TO IMPROVE MASS–SIZE RELATIONSHIPS

The use of aircraft in atmospheric research projects has given access to a large panel of cloud types, and with concomitant improvement of airborne instrumentation and data processing techniques, in situ measurement datasets have laid the groundwork for the study of natural ice cloud particles’ microphysical properties. In the early days, Heymsfield (1972) applied the on-ground mass measurement principle on board a research aircraft in order to analyze natural crystals sampled in cirrus clouds. Since then, technical developments of optical array probes (OAPs), total water content (TWC) probes, and cloud radar have made ice crystal ensemble properties, such as particle size distribution (PSD), ice water content (IWC), and radar reflectivity (*Z*<sub>e</sub>), readily accessible. As a result, the mass of ice cloud particles is now indirectly estimated from PSD and IWC (Brown and Francis 1995; Heymsfield et al. 2004, 2010; Cotton et al. 2013; Leroy et al. 2016) or from PSD and *Z*<sub>e</sub> (McFarquhar et al. 2007; Fontaine et al. 2014) rather than measured as in earlier studies.

One of today’s most popular relationships has been established in Brown and Francis (1995) with data from two flights realized in the frame of the International Cirrus Experiment field campaign (Raschke et al. 1990) in 1992. IWC values were estimated by integrating *m*–*D* relationships given in Locatelli and Hobbs (1974) over 5-s average PSDs (2650 in total). Brown and Francis (1995, p. 411) compared these calculated values to the measured IWC values and found that the *m*–*D* relationship established in Locatelli and Hobbs (1974) for “aggregates of unrimed bullets, columns and side-planes” produced the best fits. They concluded that the Locatelli and Hobbs parameters *a* = 7.38 × 10<sup>-11</sup> g µm<sup>-1.9</sup> and *b* = 1.9 were suitable for populations of quasi-spherical irregular crystal populations of midlatitude cirrus clouds.

More sophisticated techniques for determining the parameters of the power law arise in recent publications. For example, Schmitt and Heymsfield (2010) consider ice particles as fractal objects and developed a self-consistent retrieval method in which both *a* and *b* are inferred from 2D images. They used simulated crystal populations to establish a relationship between ice particle two-dimensional and three-dimensional fractal dimensions. The exponent *b*, identified as the three-dimensional fractal dimension, is inferred from 2D projected images for each 5-s data point, assuming that the relationship between fractal dimensions established for synthetic aggregates holds for natural ice crystals. Then they indirectly estimate *a* from area measurements, thereby making several assumptions: 1) the area–dimensional relationship derived from 2D image analysis is accurate; 2) a particle having an area ratio of 1 is spherical and has a density of 0.91 g cm<sup>-3</sup>; and 3) the prefactor *a* calculated from this particle, having an area ratio of 1, is applicable across the whole size range. Similarly, Leroy et al. (2016) found from simulated crystal populations that *b* can be inferred from a combination of exponents of the area–dimensional (σ) and the perimeter–dimensional (τ) power-law relationships. For each 5-s data point, *b* is computed from σ and τ values retrieved from the 2D imagery. Then the prefactor *a* is iteratively estimated as the value that minimizes the discrepancy between the IWC value calculated from PSDs and an independently measured bulk IWC value. While these two dynamical retrieval methods may capture the natural variability in ice particle ensembles along the flight path (e.g., the change in particle morphology, dominant habits, density), they both implicitly assume that the dependence of mass on size conforms to a single power law irrespective of the particle size. Furthermore, they produce as many (*a*, *b*) pairs as there are points in the dataset, which might be impractical for use in models unless reduced to a single averaged *m*–*D*.

As plots in log–log space of relationships found in literature show evidence that *m*–*D* is not linear over a broad size range, Erfani and Mitchell (2016) recently argue that a single *m*–*D* power law cannot satisfactorily represent the variation of mass over the size range spanned by atmospheric ice particles. To account for the curvature, they search for an *m*–*D* expression of the form ln*m* = *a*<sub>0</sub> + *a*<sub>1</sub> × (ln*D*) + *a*<sub>2</sub> × (ln*D*)<sup>2</sup>, where coefficients are derived from a second-order polynomial curve fit of Small Particles in Cirrus (SPARTICUS) data (Mace et al. 2009). The retrieval method combines Mitchell et al.’s (1990) mass measurements, Baker and Lawson’s (2006) mass retrieval technique, and a method for estimating the mass of small ice crystals (*D* < 100 µm) from Cloud Particle Imager (CPI) measurements in order to derive a self-consistent *m*–*D* expression from in situ data corresponding to three temperature and two
cloud-type (synoptic vs anvil cirrus) categories. While this method compares favorably with the *m*–*D* expressions of Cotton et al. (2013) and Heymsfield et al. (2010), IWC values calculated by integration of this new parameterization over PSDs have not been directly compared to an independently measured and accurate bulk IWC value.

In an attempt to better account for the dependence of ice particles' microphysical properties on size, Jackson et al. (2012) introduced a habit-dependent *m*–*D* retrieval method. Ice particles collected in Arctic mixed-phase stratus clouds were classified into nine habit classes using the CPI imagery. Then a composite *m*–*D* relationship was derived from habit-specific power-law relationships reported in Mitchell (1996) and Brown and Francis (1995), thereby weighting by the number fraction of each crystal habit in each bin. The quality of the mass retrieval method was assessed with mass closure tests in which IWC values calculated from measured PSDs and the habit-dependent scheme were compared to IWC values measured by two independent bulk instruments [the deep-cone Nevzorov probe and the Cloud Spectrometer and Impactor (CSI)]. These tests showed that this habit-dependent scheme 1) gives IWC values that agree reasonably well with the measured values and 2) outperforms the mass–area relationship derived in Baker and Lawson (2006), although the latter is based on a detailed morphological analysis of particles' projected images.

### c. Mass–size relationships: Limitation and opportunities for improvement

In conclusion of this partial literature review, there is a consensus in using power laws to describe the dependence of ice particle mass on size with a few exceptions, such as the Erfani and Mitchell (2016) scheme and the Jackson et al. (2012) habit-dependent approach. In practice, the power-law assumption 1) reduces the search of *m*–*D* to the computation of two scalars, *a* and *b*; 2) produces *m*–*D* curves that are globally matching theoretical expectations (in terms of magnitude and tendency for ice particle density to decrease with size); and 3) enables the calculation of particle ensemble properties (IWC or *Z*<sub>e</sub> for instance) that agree reasonably well with observations. However, most of the published *m*–*D* relationships based on this power-law assumption are subject to a fundamental weakness, as acknowledged by recent authors (Schmitt and Heymsfield 2010; Erfani and Mitchell 2016): they might not accurately reflect the contribution of the different growth processes to the ice crystals' masses. Especially, they cannot capture potential irregular variations of the mass with size for crystals ranging from a few microns to a few centimeters in size. Because of the preponderance of different ice growth mechanisms within distinct size intervals, the smaller size classes up to a few hundred micrometers are

generally populated by small pristine ice crystals formed by vapor deposition, whereas crystal aggregation dominates in the largest particles. It is well known that the properties of dense pristine ice and fluffy aggregates may be significantly different (Pruppacher and Klett 2010). The effect of riming could also change the mass properties in some intermediate bins so that the resulting curve may significantly deviate from a regular power-law relationship. The results reported in Jackson et al. (2012) support these considerations, which suggest that accounting for the size dependence of particles' properties, such as crystal habit distribution, could improve the mass retrieval. The situation is even worse when only one (*a*, *b*) pair is derived from an entire dataset. In this case the *m*–*D* parameterization fails to capture the spatio-temporal variability in the microphysical properties of ice particles within clouds.

Since earlier methods based on individual crystal mass measurements are impractical for use in large datasets of mixed-habit ice particles, there is still a strong need for the development of alternative concepts for *m*–*D* retrieval methods.

This study presents a nonstandard approach that leverages the optimal use of common in situ measurements—namely, PSD and IWC—to waive the power-law constraint and therefore better capture the dependence of particle microphysical properties on size. Section 2 introduces the assumptions made to set up a simple forward model that relates particle-size-dependent mass to PSD and bulk IWC values, as well as the mathematical formalism and tools selected to retrieve the particles’ mass by solving the inverse problem. The method is then applied to simulated crystal populations (section 3) and to real cloud in situ measurements extracted from the HAIC–HIWC dataset (section 4). In this article, the emphasis is placed on the problem statement and technical aspects of this new mass retrieval process rather than on the implications of the preliminary results presented in section 4 for modeling and remote sensing retrieval algorithms.

## 2. Conceptual model and numerical optimization tools

### a. Forward model

#### 1) MODEL INPUT: CLOUD IN SITU MEASUREMENTS

The mass retrieval method developed within this study is based on two coincident and collocated 5-s averaged in situ measurements—PSD and IWC—that are shortly described hereafter.

##### (i) PSD

PSDs characterize the size distribution of particles within a population. Individual crystals are sized into size classes,
or bins, using a size parameter derived from the OAP im- agery: in this study, the equivalent projected area diameter $D_{eq}$, defined as the diameter of a disk having the same area as the particle’s 2D projection, is used. In the following, PSD concentrations are calculated as the number of parti- cles per volume of sampled air per bin ($m^{-3} \text{bin}^{-1}$), al- though in atmospheric sciences, conventional number PSDs are usually expressed in $m^{-3} \mu\text{m}^{-1}$, thereby normalizing the number of particles counted in each bin by the bin width. The in situ measured PSDs are derived from 2D projected images of the 2D Stereo (2D-S) probe (Lawson et al. 2006) and the Precipitation Imaging Probe (PIP; Baumgardner et al. 2011), as further detailed in section 4. PSDs are subject to measurement uncertainties in both sizing and concen- tration as a result of optoelectronic limitations and post- processing algorithms applied to remove artifacts, such as shattering/splashing effects or particle coincidence (more than one particle in one image frame). More details on the measurement uncertainties of OAPs can be found in Baumgardner et al. (2017). Also, the attempt to character- ize ice particles’ complex 3D volume with a single size pa- rameter inferred from 2D shadow images is inherently approximate, as exemplified in section 3b.

### (ii) IWC

The IWC value gives the amount of ice in a sampled cloud volume ($\text{g m}^{-3}$). Various measurement techniques exist but they can be sorted into two categories: the “hot wire probes” that rely on the measurement of the elec- trical power necessary to maintain the probe at a con- stant temperature and the “evaporator probes” that rely on humidity measurements. In section 4, the IWC measurements are provided by the IKP-2 isokinetic evaporator probe (Davison et al. 2009; Strapp et al. 2016b). This probe was developed in the frame of the HAIC-HIWC project and is considered as a reference instrument for IWC measurements. IWC measurement uncertainty is documented in Davison et al. (2016).

## 2) PHYSICAL CONSIDERATIONS AND APPROXIMATIONS

The IWC value is by definition the sum of the individ- ual crystal masses within a sampled cloud volume. The mass $M$ of an ice particle is mostly determined by the air and ice fractions and the respective densities of these two components. Actually, it is known from theoretical con- siderations that the particle mass is not evenly distributed over the particle volume as a result of complex interactions of several growth processes, such as water vapor diffusion, aggregation, and riming (Pruppacher and Klett 2010). An accurate computation of the mass of ice particles would require knowledge about how these different mechanisms

contributed to the particle’s growth throughout its life cycle. In this Lagrangian-like approach, in situ measurements would follow the particle transportation from the nucleation site across different cloud parcels. However, the complex transport of ice particles within a cloud—especially within deep convection, where turbulent mixing processes occur over several kilometers in height—makes such data almost impossible to collect. Therefore, the proposed model uses the Eulerian approach, which is prevalent in literature, whereby the ice particles are characterized at a given loca- tion and time by collocated and coincident in situ mea- surements of the cloud environment. It is postulated that a series of measured parameters can be identified and used to determine the mass of ice particles. Among others, one can think of environmental parameters such as static air temperature (ambient temperature) and cloud type (the dependence of particle properties on cloud type is reported in several studies quoted in section 1), and morphological parameters derived from the analysis of OAP images of individual particles (e.g., dominance of some crystal habits, presence of riming) as potentially characteristic parameters. In the following, these parameters are denoted as $P_1, P_2, \dots, P_n$. In summary, the mass of a particle depends on its volume and on some other characteristic parameters, as formally written in Eq. (1):

$$M = M(V, P_1, \dots, P_n) \quad (1)$$

Two simplifying assumptions are made to establish the forward model.

### (i) Characteristic parameters

Time invariance of the parameters is assumed for each 5-s average data point: all particles sampled during a 5-s time interval are characterized by the same parameters (e.g., the time-averaged value computed from the 1-Hz original measurements). Subsequently, datasets can be sorted into subsets of consistent data points using these parameters as sorting criteria. In such a subset where the characteristic parameters are invariant, the mass de- pends only on volume. Equation (1) becomes

$$M = M_{|P_1, \dots, P_n}(V, P_1, \dots, P_n) = M(V). \quad (2)$$

### (ii) Volume

In the $i$th bin of the PSD, the variability in the particle’s size, $D_{eq}$, depends on the bin width $w_i$. As $w_i$ is generally small, the variability in $D_{eq}$ is small. It is further postulated that the spread in volume is also reasonably low within a bin (discussed in more detail in section 3b). This assumption implies that PSD bins are natural subensembles, where the crystals’ mass is uniform, provided data points have been sorted as per assumption $i$. As a consequence, a statistical
“reference mass” can be defined for each bin: $M_{|V}(V) = M_{\text{ref}}$. We note $m_i$ as the reference mass of the ith bin in the PSD, such that $m_1$ and $m_N$ are the reference masses of the first bin (the smallest particles with $D_1 - w_1/2 < D_{\text{eq}} < D_1 + w_1/2$) and the last bin (largest particles with $D_N - w_N/2 < D_{\text{eq}} < D_N + w_N/2$, where N is the number of bins in the PSD), respectively. There are as many reference masses as bins in the PSD, and in

the following we note $\mathbf{m} = (m_1 m_2 \dots m_N)^T$ the vector of reference masses.

With the two abovementioned assumptions, we now develop the mathematical equations that relate bin reference masses to PSD and IWC measurements. For a data point measured at $t = t_1$, the ice water content is calculated by summation over the PSD bins of the number of crystals in each bin $[\text{PSD}(t_1, D_i)]$ times the bins’ reference mass ($m_i$):

$$ \text{PSD}(t_1, D_1) \times m_1 + \text{PSD}(t_1, D_2) \times m_2 + \dots + \text{PSD}(t_1, D_N) \times m_N = \text{IWC}(t_1). \quad (3) $$

Equation (3) can be written as many times as there are data points in a subset satisfying assumption i, say P

points in this explanation. This results in a system of linear equations with N unknowns and P equations:

$$ \left\{ \begin{array}{l} \text{PSD}(t_1, D_1) \times m_1 + \text{PSD}(t_1, D_2) \times m_2 + \dots + \text{PSD}(t_1, D_N) \times m_N = \text{IWC}(t_1) \\ \text{PSD}(t_2, D_1) \times m_1 + \text{PSD}(t_2, D_2) \times m_2 + \dots + \text{PSD}(t_2, D_N) \times m_N = \text{IWC}(t_2) \\ \vdots \\ \text{PSD}(t_p, D_1) \times m_1 + \text{PSD}(t_p, D_2) \times m_2 + \dots + \text{PSD}(t_p, D_N) \times m_N = \text{IWC}(t_p) \end{array} \right. \quad (4) $$

For the sets of synthetic data (section 3) and experimental data (section 4) presented in the next sections, $N = 93$ and 240 and $P = 186$ and 1051, respectively. This system of linear equations can be written as a matrix equation:

$$ \mathbf{PSD} \cdot \mathbf{m} = \mathbf{iwc}, \quad (5) $$

where

* $\mathbf{PSD} \in \mathbb{R}^{P \times N}$ is a matrix in which PSDs are vertically concatenated,

* $\mathbf{iwc} \in \mathbb{R}^P$ is a vector of P IWC values, and

* $\mathbf{m} = (m_1 m_2 \dots m_N)^T$ is the vector of N unknown reference masses.

Thus, the forward model is formulated as a matrix equation that relates in situ measurements to the unknown bin reference masses. It is worth mentioning that the model does not require prior knowledge of the descriptors $P_1, P_2, \dots, P_n$. Sensitivity studies in which reference masses are retrieved from two distinctive subsets will reveal whether the parameter that differentiates between the two subsets influences the mass of ice particles (see section 4 for an example with sampling temperature as a descriptor).

b. Solving the inverse problem

1) SOLVING AN ILL-CONDITIONED INVERSE PROBLEM

The mass retrieval process now consists of solving the inverse problem. Estimates of the problem condition number—defined as $c = \sqrt{\sigma_{\text{max}}/\sigma_{\text{min}}}$, where $\sigma_{\text{max}}$ and

$\sigma_{\text{min}}$ are the largest and smallest singular values of the PSD matrix, respectively—indicate that measurement uncertainties and noise inherent to the measured data significantly influence the solution: the problem written with synthetic data (case presented in section 3) is severely ill-conditioned ($c > 10^4$) and the problem using HAIC-HIWC real data (case presented in section 4) is ill-posed ($c \to \infty$). To circumvent the difficulty related to numerical instabilities, we follow the regularization approach presented in Idier (2008) and search for a solution that minimizes a composite criterion, noted as $J(\mathbf{m})$:

$$ J(\mathbf{m}) \stackrel{\text{def}}{=} \| \mathbf{PSD} \cdot \mathbf{m} - \mathbf{iwc} \|^2 + \lambda \cdot R(\mathbf{m}), \quad (6) $$

where $J(\mathbf{m})$ is twofold:

* The first term is referred to as the least squares term. Minimizing that term forces the solution vector to fit the measurements and to satisfy Eq. (5).

* The second term, referred to as regularization term, allows the introduction of some prior information that the solution should fulfill in order to qualify from a physical point of view and helps in reducing the noise propagation on the solution.

The retrieved solution will be a trade-off between fidelity to the measured data and compliance with prior information embedded in the regularization term. The compromise between the two sources of information is adjusted via the value of the regularization parameter $\lambda$. Since the retrieval process relies on a series of
measurements, a measurement error covariance matrix $S_e$ is added to weight the contribution of each data point according to their respective uncertainty, as presented in Borsdorff et al. (2014). The cost function $J$ becomes

$$J(\mathbf{m}) \stackrel{\text{def}}{=} \|\mathbf{S}_e^{-1/2}(\text{PSD} \cdot \mathbf{m} - \text{iwc})\|^2 + \lambda \cdot R(\mathbf{m}). \quad (7)$$

In the synthetic case, data points are generated from populations of synthetic crystals, so there is no measurement error associated with the data: they are all equally uncertain and the covariance matrix is an identity matrix. The real data are impacted by measurement errors that may vary from one data point to the other: $S_e$ is a diagonal matrix (measurement errors are assumed uncorrelated between data points) and its coefficients are specified in section 4.

## 2) REGULARIZATION

A regularization term is added to the least squares problem in order to stabilize the inversion and to incorporate some prior information about the wanted solution into the minimization process. The definition of the regularization term is qualitative and depends on what known properties the solution should comply with. In this study, a data-based regularization approach has been chosen: $R(m) = \|L_k \Delta^{-1} m\|^2$, where $L_k$ is the discrete approximation to a $k$th-order derivative operator (Tikhonov regularization matrix) and $\Delta$ is a diagonal weighting matrix. The regularization is defined to promote two properties:

* Curve smoothness: The discrete approximation to the first-order derivative operator is chosen [with matrix $L_1$ as defined in Huckle and Sedlacek (2012)]

* Superlinear increase of mass with size: The weighting matrix $\Delta$ balances the effect of the regularization on the elements of the bin reference mass vector

Elements of $\Delta$ are initially set as the mass of spherical ice particles: $\Delta = \text{diag}[m_{\text{sph}}(D_1 + 5), \dots, m_{\text{sph}}(D_N + 5)]$, where $D_1, \dots, D_N$ are bin centers and $m_{\text{sph}}(D) = (\pi/6)\rho_{\text{ice}}D^3$. The rationale for this definition, particularly the +5 offset added to bin centers in the definition of $\Delta$, is explained in section 3. As suggested in Huckle and Sedlacek (2012), an iterative regularization process may be conducted to improve results whereby the elements of $\Delta$ are updated at each new optimization run with the results from the previous run. The mass vectors presented in section 3 and section 4 result from solely one run of optimization, since a double run did not significantly improve or change the results.

## 3) OTHER ALGEBRAIC OPERATIONS

The particle mass is proportional to volume and increases superlinearly with the size parameter. When the size range of interest spans several orders of magnitude, the reference masses are also expected to span several orders of magnitude, which makes the problem poorly scaled. A diagonal scaling may be applied to improve the problem formulation and to decrease the condition number. The optimization problem is solved in terms of the new variable $\mathbf{m}'$, defined as $\mathbf{m} = \Delta' \cdot \mathbf{m}'$, where $\Delta'$ is a diagonal matrix as defined below:

$$\Delta'_{i,j} = \begin{cases} [(\pi/6)\rho_{\text{ice}}D_i^3]^{-1} & \text{if } i=j \\ 0 & \text{elsewhere} \end{cases}. \quad (8)$$

This operation is applied to HAIC-HIWC data, where the size range of interest is $(10 - 12840) \mu\text{m}$. It is not applied to synthetic data, where the size range of interest is limited to $(10 - 940) \mu\text{m}$.

## 4) OPTIMIZATION ALGORITHMS

All algorithms used in this paper are sourced from Nocedal and Wright (2006). The problem is considered as an unconstrained optimization problem and a line search strategy is adopted to identify the minimizer $m^* \in \mathbb{R}^N$ of the objective function. Briefly, the solution is found by creating a convergent sequence of iterates $\{\mathbf{m}^{(k)}\}_k$ such that $\{J(\mathbf{m}^{(k)})\}_k$ decreases. At each iteration, the new iterate is calculated as per Eq. (9):

$$m^{(k+1)} = m^{(k)} + \alpha_k d_k, \quad (9)$$

where $d_k$ is the search direction and $\alpha_k$ is the step length.

At the $k$th iteration, Newton's method is used to compute $d_k$ from the gradient and the Hessian matrix of the objective function [Eq. (10)]. Since the computation of the Hessian matrix is not cumbersome for small size problems, this method was chosen for its fast rate of convergence and because it is relatively insensitive to poor scaling.

$$d_k = -\nabla^2 J(\mathbf{m}^{(k)})^{-1} \cdot \nabla J(\mathbf{m}^{(k)}). \quad (10)$$

The step length $\alpha_k$ is computed either iteratively with a backtracking line search algorithm given in Nocedal and Wright (2006) (applied to the HAIC-HIWC dataset in section 4) or explicitly from Eq. (11) (applied to synthetic data in section 3):

$$\alpha_k = - \frac{\nabla J(\mathbf{m}^{(k)})^\mathrm{T} d_k}{d_k^T \nabla^2 J(\mathbf{m}^{(k)}) d_k}. \quad (11)$$
The new iterate $\mathbf{m}^{(k+1)}$ is computed from Eq. (9), and this calculation is repeated until convergence of the sequence $\{J(\mathbf{m}^{(k)})\}_k$ is reached, with $\varepsilon = 2.2204 \times 10^{-16}$ as the stopping criterion.

The regularization parameter $\lambda$ controls the degree of regularization applied to the problem. Several techniques are available to identify the optimal $\lambda$, that is, the value yielding the best compromise between the two terms of the objective function. The L-curve technique, a simple graphical way of determining $\lambda_{\text{opt}}$ from plots of the least squares term against the regularization term, is used here (Hansen 1992). An example of such a curve is given in section 3c.

### 3. Validation of the method

In this section the algorithms defined in section 2 are applied to a synthetic dataset presented in section 3a. In section 3b, the data are used to improve the understanding of real aircraft data in order to tune the regularization function accordingly. The validation of the retrieval process is presented in section 3c, where the retrieved masses are compared against reference values computed from the actual mass of simulated 3D particles.

#### a. Generation of the synthetic dataset

The generation of a synthetic dataset starts with the simulation of ice particle models using an Interactive Data Language (IDL) program (Fontaine 2014; Fontaine et al. 2014). A model is characterized by its shape and dimensional properties, which are defined in Fig. 1 and summarized in Table 1. Each model is randomly rotated in space and projected to produce a 2D image. This process reproduces the working principles of OAPs and how 2D images are produced from real 3D crystals. The rotation/projection step is repeated 100 times, such that 100 2D images are created from each 3D model. The images are then processed to derive $D_{\text{eq}}$ using the same image processing software as for OAP images. The resulting size interval of $D_{\text{eq}}$ available for each crystal habit is indicated in the last column of Table 1.

Mass is calculated from the 3D volume by prescribing a constant ice density $\rho_{\text{ice}} = 0.917 \text{ g cm}^{-3}$ for all particles: it has the effect of fixing all physical parameters in the forward model (assumption *i* is satisfied), hence the particle mass depends solely on volume. For each habit, a lookup table is created where a single theoretical mass is associated with $100 D_{\text{eq}}$ resulting from the 100 respective projections.

The next step is the creation of a synthetic data point that replicates a 5-s average in situ measurement data point. It consists of forming a population of thousands of simulated ice particles described by one PSD and its corresponding theoretical IWC value, noted $\text{IWC}_{\text{th}}$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Simulated crystals: habit and size definitions.</figcaption></figure>

Two input requirements (user's input) control the generation of a simulated crystal population:

* PSD shape: The size range considered in the synthetic test cases is (10–940) µm. Simulated crystals are sized into 93 equally spaced 10-µm-wide bins, and the number of particles per bin is prescribed by a normal probability density function of the bin midpoint multiplied by a prefactor so that there are approximately 10000 elements in each population (or data point). The normal distribution parameters—namely, its mean $\mu$ and standard deviation $\sigma$—are changed from one data point to the other when generating the synthetic dataset, as detailed below.

* Habit mixing: A mixing law is defined to describe the fraction of each crystal habit within each bin, as in Table 2, for example.

PSD bins are then populated by randomly picking up crystals based on their $D_{\text{eq}}$ from habit tables with respect to the number and habits prescribed by the user's input. For each bin the theoretical bin reference mass, $m_{\text{ref},i}$, for the $i^{\text{th}}$ bin, is calculated as the average of the theoretical masses of the binned particles. This value is stored for later comparison with the bin reference mass retrieved by solving the inverse problem. $\text{IWC}_{\text{th}}$ is calculated by summing the actual mass of all the elements in the population.

Finally, a synthetic dataset is composed of *P* different data points that are generated by randomly varying $\mu$ ($\mu_0 \le \mu \le \mu_0 + 100$) and $\sigma$ ($100 \le \sigma \le 200$). To test the retrieval method, three synthetic datasets are used:

* One optimization dataset, noted as DS1: this is the dataset to which the optimization process is applied. Twice as many PSDs as unknowns are generated ($P=2N$) to ensure that the **PSD** matrix in Eq. (5) has a full rank. In DS1, $\mu_0 = 300$. The unknown mass vector **m** of Eq. (5) (whose elements are the sought bin reference masses) is found by minimizing Eq. (6) written with $\text{PSD}_{\text{DS1}}$ and $\text{IWC}_{\text{th,DS1}}$ quantities.
<table><caption>TABLE 1. Synthetic dataset: habits and sizes of simulated crystal models.</caption><thead><tr><th rowspan="2">Code</th><th rowspan="2">Habit</th><th colspan="2">Size descriptors</th><th colspan="2">Observed sizes</th></tr><tr><th>L, a, and D (μm)</th><th>No. of models</th><th>D<sub>eq</sub> (μm)</th></tr></thead><tbody><tr><td>Sph</td><td>Sphere</td><td>D ∈ (5 – 1000)</td><td>583</td><td>5–1000</td></tr><tr><td>Cyl_R5</td><td>Cylinder</td><td>R = 5, L ∈ (10 – 300)</td><td>146</td><td>11–65</td></tr><tr><td>Cyl_R10</td><td>Cylinder</td><td>R = 10, L ∈ (10 – 500)</td><td>246</td><td>17–116</td></tr><tr><td>Cyl_R30</td><td>Cylinder</td><td>R = 30, L ∈ (30 – 1500)</td><td>317</td><td>46–353</td></tr><tr><td>Cyl_R60</td><td>Cylinder</td><td>R = 60, L ∈ (30 – 1500)</td><td>257</td><td>68–489</td></tr><tr><td>Cyl_R90</td><td>Cylinder</td><td>R = 90, L ∈ (30 – 1500)</td><td>262</td><td>86–597</td></tr><tr><td>Cyl_R120</td><td>Cylinder</td><td>R = 120, L ∈ (30 – 1500)</td><td>197</td><td>102–688</td></tr><tr><td>Cyl_R150</td><td>Cylinder</td><td>R = 150, L ∈ (30 – 1500)</td><td>197</td><td>111–769</td></tr><tr><td>C1e</td><td>Solid column</td><td>c/a = 1, a ∈ (5 – 900)</td><td>395</td><td>5–1114</td></tr><tr><td>C1g</td><td>Solid thick plate</td><td>c/a = 0.2, a ∈ (10 – 1000)</td><td>260</td><td>5–939</td></tr><tr><td>P1a</td><td>Hexagonal plate</td><td>c/a = 0.1, a ∈ (20 – 1000)</td><td>212</td><td>8–924</td></tr><tr><td>P1c</td><td>Crystal with broad branches</td><td>c = 10, a ∈ (11 – 188)</td><td>236</td><td>10–146</td></tr><tr><td>P1d</td><td>Stellar crystal</td><td>c = 10, a ∈ (8 – 125)</td><td>218</td><td>5–48</td></tr><tr><td>CP1a_1</td><td>Column with plate</td><td>a/a' = 3, c/a = 1, a ∈ (10 – 1000)</td><td>289</td><td>9–1137</td></tr><tr><td>CP1a_2</td><td>Column with plate</td><td>a/a' = 5, c/a = 0.2, a ∈ (20 – 1000)</td><td>213</td><td>13–939</td></tr></tbody></table>

* Two validation datasets, noted DS2 and DS3, are used to test the mass vector retrieved from DS1 on populations dominated by small ($\mu_0$ = 100 µm) and large ($\mu_0$ = 600 µm) particles, respectively. IWC values, noted IWC<sub>cal,DS2</sub> and IWC<sub>cal,DS3</sub>, are calculated from the retrieved mass vector **m** and **PSD<sub>DS2</sub>** and **PSD<sub>DS3</sub>** matrices, respectively. These calculated values are then compared to the corresponding IWC<sub>th</sub> for each data point.

Furthermore, some statistical reference vectors can be defined, since the masses of all individual particles in the population are known:

* Lower and upper bounds: vectors whose elements are the minimum and maximum theoretical bin reference masses, respectively, found across the dataset

* Average bin reference mass [$m_{th} = (m_{th,1} m_{th,2} \dots m_{th,N})$]: vector whose elements are the average of the theoretical bin reference mass ($m_{ref,i}$) calculated over the dataset

Some quantitative criteria are calculated to estimate the quality of retrieved values:

* Mass percent error (PE) is computed for each bin from the retrieved mass value ($m_i$, optimization output) and $m_{th,i}$:

$$PE_i(\%) = 100 \frac{m_i - m_{th,i}}{m_{th,i}}$$

* Mass mean absolute percentage error:

$$\text{Mass MAPE}(\%) = \frac{100}{N} \sum_{i=1}^{N} \frac{|m_i - m_{th,i}|}{m_{th,i}},$$

where N is the number of bins in the PSD, indicates the accuracy of the mass retrieval process

* IWC mean absolute percentage error:

$$\text{IWC MAPE}(\%) = \frac{100}{P} \sum_{i=1}^{P} \frac{|\text{IWC}_{\text{cal},i} - \text{IWC}_{\text{th},i}|}{\text{IWC}_{\text{th},i}},$$

where P is the number of data points in the set, indicates the accuracy of the IWC retrieval process

b. Study of input noise in PSD data

The mass retrieval process uses PSD as input data. In addition to the unavoidable measurement uncertainties, approximations made in the forward model also make the PSD a noisy input for inversion. In the frame of the proposed inverse problem approach, the synthetic dataset is

<table><caption>TABLE 2. Mixed-habit population: Mixing laws for different size intervals.</caption><thead><tr><th>Size interval (bin centers in μm)</th><th colspan="4">Crystal habit (fraction in %)</th></tr></thead><tbody><tr><td>15–45</td><td>P1a (25%)</td><td>C1g (25%)</td><td>P1c (25%)</td><td>P1d (25%)</td></tr><tr><td>55–145</td><td>P1a (33%)</td><td>C1g (33%)</td><td>P1c (33%)</td><td></td></tr><tr><td>155–205</td><td>P1a (50%)</td><td>C1g (50%)</td><td></td><td></td></tr><tr><td>215–795</td><td>CP1a_1 (50%)</td><td>CP1a_2 (50%)</td><td></td><td></td></tr><tr><td>805–935</td><td>CP1a_1 (33%)</td><td>CP1a_1 (33%)</td><td>Sph (33%)</td><td></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Mass variability in PSD bins for C1g population; (b) spread in mass in one bin (105 ± 5 µm) for all simulated habits.</figcaption></figure>

first used as an opportunity to study the noise contained in
PSDs based on OAP images.

The concept of bin reference mass relies on the as-
sumption that PSD bins are natural subensembles where
the crystals' mass is uniform, which in turn holds if: 1)
PSD bins are narrow and 2) the size descriptor used to
generate PSD accurately describes the particles' vol-
ume. These conditions are hardly met in reality. On the
one hand, the width of a PSD size class is driven by the
OAP specifications: it cannot be smaller than the probe
resolution, which is nonnegligible (e.g., 10–100 µm, as
detailed later in section 4). On the other hand, any size
descriptor derived from 2D projected images will fail to
accurately describe particle volume: with exception of
spherical particles, whose volume is precisely described
by a single size parameter (e.g., its diameter); an accu-
rate volume description of real ice particles with com-
plex geometry (e.g., plates, columns, bullet rosettes,
capped columns, or even aggregates of those) would
require several size parameters. In addition, the pro-
jection of 3D objects onto a 2D plane inevitably distorts
the volume’s dimensional properties so that the arbitrary
projected 2D image scarcely represents the original 3D
volume information. This results in a large variability of
particle masses when sized into bins based on D<sub>eq</sub> or
similar size parameters derived from 2D projected images.

Figure 2a illustrates this variability for a population
of simulated hexagonal plates (c/a = 0.2, code = C1g in
Table 1). The particles are sized into equally spaced
10-µm-wide PSD bins using D<sub>eq</sub> and the minimum,
mean, and maximum of particle masses in each bin are
plotted against bin centers. For that specific habit, the
particles’ minimum and maximum mass differ by at

least one order of magnitude in the 15–465-µm size
interval. Above 465 µm, the maximum mass curve
plateau observed up to 935 µm indicates that the larg-
est available simulated particle in this population
(L = 1000 µm; M ≈ 1.244 × 10⁻⁴ g) can end up in 48
different bins depending upon its spatial orientation
before projection.

Similarly, the box plot given in Fig. 2b illustrates the mass spread for different crystal shapes ending up in the 105 ± 5-µm size bin, with habit codes as indicated in Table 1. Each box plot gives the median, 25th, and 75th percentiles, as well as extreme values. One can see that generally the spread increases as the particle shape departs from the sphere. It is worth mentioning that the process of sizing multidimensional 3D particles with a size parameter derived from 2D projected images is not the only source of variability: a nonnegligible mass spread is also observed for spherical particles (code = Sph) as a result of the bin width, that is, sorting particles whose diameters range from 100 to 109 µm into the same 10-µm-wide bin. A preliminary study where the population of simulated hexagonal plates (c/a = 0.2; code = C1g in Table 1) is sized using different size parameters reveals that using the maximum dimension Dmax as a size parameter decreases the spread of mass values within bins, compared to using Deq.

This spread in the particles' mass within a bin is ex-
pected to be even more pronounced in real data: OAP
images reveal that several particle habits are mixed
within the ice cloud particle populations sampled in
airborne field campaigns. These considerations give an
insight into the level of noise included in PSD data,
which is detrimental to the accuracy of the mass retrieval
process. It also underlines the need for regularization
when attempting to inverse the forward model presented in section 2a. With this in mind, we proceed with OAP measurements for lack of better inputs.

### c. Application and validation of the mass retrieval process to synthetic data

#### 1) TEST POPULATIONS

The method has been applied to three test cases with increasing degree of complexity:

1) **Spherical particles:** The populations are composed of spherical particles only (code = Sph in Table 1). It limits the source of noise in input data to the binning process, since the characteristic dimension of a sphere is insensitive to 2D projection.

2) **Hexagonal plates:** The populations are composed of solid thick plates with an aspect ratio $c/a = 0.2$ (code = C1g in Table 1) only. Since the geometry deviates significantly from the ideal sphere case, it allows for a preliminary assessment of projection noise effects on the retrieval process.

3) **Mixed habit populations:** In these populations, the crystals' habit varies with size according to a habit mixing law presented in Table 2 with habit codes given in Table 1 (and shown in Fig. 4). This habit mixture law applies to the three datasets (DS1–DS3).

The results obtained for the first two test cases are not reported in the present article, but they helped in the definition of the regularization term and specifically in the choice of $\Delta$ elements: if masses of spherical ice particles of diameter equal to the bin center were selected as weighting coefficients in the regularization matrix for simplicity reasons, it was found by trial and error that adding the +5-$\mu$m offset to sphere diameters yielded improved estimates of the retrieved masses in small size bins, that would have been otherwise underestimated when particle habits significantly depart from spherical shape. The last test case is quite challenging because abrupt changes in crystal habits between two consecutive bins produce sharp variations of mass with size. Classical power-law-based mass retrieval methods would typically fail to capture these features by smoothing discontinuities into a simplistic power-law curve. Given the size range (10–940 $\mu$m) and the number of bins (93 equally spaced 10-$\mu$m-wide bins) of this synthetic test case, there are 93 unknown reference masses to be solved.

#### 2) RESULTS AND DISCUSSION

Figure 3 presents optimization results for the mixed-habit population. Figures 3a and 3b show the mass vectors retrieved for the nonregularized ($\lambda = 0$) and the optimally

regularized ($\lambda = \lambda_{opt} = 7.5 \times 10^{-5}$) least squares problem, respectively. It highlights the importance of the regularization to extract meaningful information from noisy data.

Qualitatively, one can see that the mass vector retrieved for $\lambda = \lambda_{opt}$ matches expectations: masses are all positive and increase with size in close relation with the average bin reference masses, computed from actual masses of simulated particles (Fig. 3b). The regularization has an obvious effect of smoothing sharp discontinuities that are present in the synthetic data: the pronounced discontinuity between 795 and 805 $\mu$m, where bin mean mass is approximately multiplied by 4 due to the introduction of spherical particles, is spread between 695 and 865 $\mu$m. Smaller discontinuities, like the one between 205 and 215 $\mu$m, are totally smoothed out. Figure 3c illustrates how $\lambda_{opt}$ is selected: the curve, referred to as the L curve, of the least squares term against the regularization term has a characteristic L shape, hence its name, and the value of $\lambda$ that yields the best compromise in amplitude between the least squares term and the regularization term can be found at the L-shaped “corner.” Finally Figs. 3d–g give a quantitative assessment of the solution accuracy:

* The mass PE is plotted in Fig. 3d for each bin. If large disagreements are found in boundary bins where particle habits change abruptly ($D_{eq} = 45, 145, 205$ and $795 \ \mu$m, with PE in excess of 50%), the resulting mass MAPE of only 18.1% indicates that bin reference masses computed by optimization algorithms are generally in good agreement with theoretical values.

* Figures 3e–g illustrate how the calculated IWC$_{cal}$ values compare to the theoretical IWC$_{th}$ values for DS1–DS3 [scatterplots, where the red curve underneath the data points (blue symbols) is the $x = y$ line]. The mass vector retrieved from DS1 produces excellent estimates of the IWC in the validation sets: as shown in Table 3, IWC MAPE values calculated for DS2 and DS3 are less than 1.6% and 0.6%, respectively.

The effective density $\rho_{eff}$, defined in this study as the bin reference mass divided by the volume of a sphere of diameter equal to the bin center, is plotted against size in Fig. 4. The retrieved effective density varies between 0.115 and $0.928 \ \text{g cm}^{-3}$ following a pattern controlled by the crystal habits. The retrieved effective density values (blue cross symbols) generally compare well with the theoretical values

$$ \rho_{\text{eff,th}}(i) = \frac{m_{\text{th},i}}{(\pi/6) \times D_i^3}, $$

(black dots in the figure) especially for particle sizes above 250 $\mu$m, where the plateau values are well reproduced notwithstanding the smoothing effect of the
<figure><img src="image_4.png" /><figcaption>FIG. 3. Mass retrieval from synthetic data: (a),(b) m-D plots for $\lambda = 0$ and $\lambda = \lambda_{opt}$, respectively; (c) L curve; (d) mass PE between calculated and theoretical bin reference masses; (e)-(g) comparison between retrieved and theoretical IWC for DS1-DS3, respectively.</figcaption></figure>

regularization. The results are not as good for smaller particles: this is due to the low contribution of these particles to the total ice mass resulting from 1) the small individual weight of these particles and 2) the relatively small number of particles in these bins as a result of the choice of normal distribution centered around 300 and 400 µm. This demonstrates the potential of this method to capture size-dependent features much better than any power-law-based method, producing a regular polynomial curve by construction:

$$ \rho_{\text{eff}}(i) = \frac{6\alpha}{\pi} D_i^{\beta-3}. $$

## 4. Application to HAIC-HIWC field campaign data

In this section the inversion method is tested on in situ measurements extracted from the HAIC-HIWC dataset in order to evaluate its potential in atmospheric sciences problems.
<table><caption>TABLE 3. Optimization results for the synthetic dataset.</caption><thead><tr><th>Mass<br>MAPE (%)</th><th>IWC<sub>DS1</sub><br>MAPE (%)</th><th>IWC<sub>DS2</sub><br>MAPE (%)</th><th>IWC<sub>DS3</sub><br>MAPE (%)</th></tr></thead><tbody><tr><td>18.1</td><td>0.55</td><td>1.56</td><td>0.57</td></tr></tbody></table>

### a. HAIC-HIWC dataset

In the frame of HAIC-HIWC projects (Dezitter et al. 2013; Strapp et al. 2016a), two airborne field campaigns were held in Darwin, Australia, in 2014 and out of Cayenne, French Guyana, in 2015. With the overarching goal of characterizing ice particles in mesoscale convective systems (MCSs), these campaigns resulted in nearly 15 000 data points each, sampled at different stages of cloud life cycle in both oceanic and continental convection. This makes the HAIC-HIWC dataset one of the most comprehensive datasets on MCSs from the standpoint of the microphysical properties of hydrometeors.

In this section we use data taken from three flights, as reported in Table 4. The optimization is conducted on the data of Darwin flight 16 (D16): all data points were sampled in the same oceanic MCS at a single temperature level (236.1 K) such that the dataset is consistent with respect to two parameters (cloud temperature and cloud type) potentially influencing the particles' mass. The other reasons for choosing D16 as the optimization dataset is the high number of data points (1051 data points) and the large size range covered by the sampled particles.

The validation is carried out on data from two arbitrarily chosen flights:

* Darwin flight 19 (D19) was conducted in an oceanic MCS 2 days after D16. Two temperature levels (237 and 226 K) were sampled but only the 212 data points sampled at 237 K (consistent with the temperature range of D16) are used in this study.

* Cayenne flight 26 (C26) was performed one year later in coastal convection. Three temperature levels (264, 244.2, and 229.4 K) were sampled, and all data points (1298 in total) are conserved so that the potential effects of temperature on the retrieved particle mass can be assessed.

The analysis considers only data points measured during level flight periods in areas where IWC values are

<figure><img src="image_5.png" /><figcaption>FIG. 4. Test case 3: Effective density against size.</figcaption></figure>

larger than 0.1 g m⁻³. As documented in Davison et al. (2016), the uncertainty in IWC measurements mostly depends on the temperature and the IWC value: although relative uncertainty remains below a few percent in most of the operating conditions, it can reach 50% for warm temperatures (263.15 K) and small IWC values (0.1 g m⁻³), hence the need to consider measurement errors in the retrieval process. Given the conditions of D16, the coefficients of the matrix S<sub>e</sub> are derived from IWC relative uncertainty values u<sub>IWC</sub>, calculated using Eq. (12) (adapted from Fig. 1 in Davison et al. 2016):

$$ \begin{cases} u_{\text{IWC}}(\%) = 0.0204 \text{ IWC}^{-0.2344} \\ u_{\text{IWC}}(\%) = 0.023\,579 \text{ IWC}^{-0.025\,532} \end{cases} \quad \begin{array}{l} \text{for IWC} \in (0.1 - 0.5) \text{ g m}^{-3} \\ \text{for IWC} > 0.5 \text{ g m}^{-3} \end{array}. \quad (12) $$

The in situ measured PSD are derived from 2D projected images of the 2D-S probe (Lawson et al. 2006;

10–1280-µm size range; 10-µm resolution) and the PIP (Baumgardner et al. 2011; 100–6400-µm size range;
<table><caption>TABLE 4. HAIC-HIWC dataset: Overview of the data from selected flights.</caption><thead><tr><th rowspan="2">Flight label (identification)</th><th rowspan="2">Date, location, cloud type</th><th>Temperature level</th><th>IWC values</th><th>No. of data points</th></tr><tr><th>Mean (std dev) (K)</th><th>Mean/max (g m<sup>-3</sup>)</th><th>Number</th></tr></thead><tbody><tr><td>Darwin 16 (D16)</td><td>7 Feb 2014, Broome, oceanic MCS</td><td>236.1 (0.363)</td><td>1.07/2.75</td><td>1051</td></tr><tr><td>Darwin 19 (D19)</td><td>9 Feb 2014, Broome, oceanic MCS</td><td>236.6 (0.545)</td><td>0.54/1.46</td><td>212</td></tr><tr><td rowspan="3">Cayenne 26 (C26)</td><td rowspan="3">29 May 2015, Cayenne, coastal MCS</td><td>229.4 (0.191)</td><td>0.47/3.03</td><td>936</td></tr><tr><td>244.2 (0.2)</td><td>1.25/2.52</td><td>186</td></tr><tr><td>264 (0.49)</td><td>1.49/3.04</td><td>176</td></tr></tbody></table>

100-µm resolution). The particles are sorted into size classes according to their $D_{eq}$, and a composite PSD with bin centers ranging from 15 to 12800 µm is computed by combining data from the two imaging probes. The composition technique is adapted from Leroy et al. (2016) with two minor deviations: 1) the overlap limits are 800-1250 µm; and 2) the width of the PSD bins is kept as per the OAPs' original resolution except in the overlap (800-1250) µm, where the width of PIP bins is numerically decreased from 100 to 10 µm by linear interpolation to ease merging with 2D-S data. As a result, composite PSDs have 240 bins with a bin width of 10 µm from 15 to 1245 µm and a bin width of 100 µm from 1300 to 12800 µm. For D16, the 1051 × 240 PSD matrix has a rank of 231 because 9 bins out of 240 are empty throughout the flight. These bins are at the largest end of PSDs (bin centers are 1.17, 1.19, 1.2, and 1.23-1.28 cm) where the concentrations are almost insignificant if any. The uncertainty in concentration measurements $u_{PSD}$ is quite hard to establish because of the paucity of data published in literature for these two probes. In this study, estimates provided by Baumgardner et al. (2017) are used with minor deviations:

$$ \begin{cases} u_{\text{PSD}} = +100\%/-50\% & \text{for } D_{\text{eq}} < 100 \, \mu\text{m} \\ u_{\text{PSD}} = +50\%/-33\% & \text{for } D_{\text{eq}} > 100 \, \mu\text{m} \end{cases} \quad (13) $$

The impact of concentration uncertainties on the retrieved masses is evaluated by solving the Eq. (7) for two extreme conditions: once with concentrations systematically underestimated (PSD = PSD + $u_{PSD} \times$ PSD, with negative values of $u_{PSD}$) and once with concentrations overestimated (PSD = PSD + $u_{PSD} \times$ PSD, with positive values of $u_{PSD}$) which gives the upper and lower mass bounds in Fig. 5, respectively.

b. *Mass retrieval from in situ measurements*

The retrieved $m-D$ relationship is plotted in Fig. 5a (blue cross symbols) along with the upper and lower mass bounds (shaded area) computed from the underestimated and overestimated PSDs. The red and yellow dashed lines represent the masses of spherical ice

particles ($\rho = 917 \text{ kg m}^{-3}$) and spherical air particles ($\rho = 1.225 \text{ kg m}^{-3}$), respectively. Qualitatively, the difference between the retrieved curve and the straight line that any power-law based method would have produced is noticeable. Masses are found to globally increase over the size interval although the function is not monotone, especially in the (500-2000) µm subinterval. The plot of effective density (Fig. 5b), calculated from this mass vector, confirms the tendency of the effective density to decrease with size, as commonly reported in literature (e.g., Pruppacher and Klett 2010), albeit for other density definitions. Qualitatively, the shape of the curve reveals quite irregular variations of the effective density as a function of particle size and two distinctive regimes for small ($D_{eq} \le 300 \mu\text{m}$) and large ($D_{eq} \ge 2000 \mu\text{m}$) particles: the effective density of small particles, lying between 0.66 and $0.28 \text{ g cm}^{-3}$ as the size increases from 15 to 300 µm, is one order of magnitude larger than the effective density found for particles above 2000 µm, where the values lie below $0.015 \text{ g cm}^{-3}$. A preliminary assessment of the impact of the choice of the regularization term (different weighting matrices were tested) on the effective density shows that the effective density calculated for the bins of the PSD ends (for particle sizes below 200 µm and for sizes above 3000 µm to a lesser extent) can be significantly influenced by the definition of the regularization term. In contrast, effective density values are more consistent for intermediate size bins. In the $200 \le D_{eq} \le 3000 \mu\text{m}$ size range, the effective density decreases from 0.3 to $0.025 \text{ g cm}^{-3}$ approximately, regardless of the definition of the tested weighting functions (the results of the preliminary sensitivity study are not presented here). The unusually low values retrieved for particles larger than 2000 µm, especially the trough between 2700 and 7700 µm, are subject to discussion. In this size range, the larger bins are populated in only a few data points [90% of the bins above 2000 µm are empty in more than 50% of the data points (529 out of 1051) and 1022 data points have more than 90% of the bins empty above 4000 µm], and for these data points the measured concentrations are extremely small. Therefore, the contribution of largest particles to the total ice mass is negligible. As a consequence of the small significance of the total mass
<figure><img src="image_6.png" /><figcaption>FIG. 5. Mass retrieval from real HAIC-HIWC data: (a) *m* − *D* plot for λ = λopt and (b) effective density. (c)−(e) Comparison between IWCcal and IWCth for D16, D19, and C26 datasets, respectively.</figcaption></figure>

of the smallest and largest particles in the least squares
term, it could be argued that the masses of these par-
ticles are neglected by the minimization algorithms.
This limitation could probably be alleviated with better
scaling of the data (here, the applied scaling focuses on
balancing the reference masses’ contributions with re-
spect to a supposed superlinear increase with size).
Alternatively, truncated PSDs could be used as input
data. However, in this study, the span of PSDs is pur-
posely kept as defined in other HAIC-HIWC studies
[i.e., up to 1.28 cm as in Leroy et al. (2016)] in order to
demonstrate the method’s ability to handle a large
number of unknowns.

Since the particle’s true masses are unknown, the quality
of the mass retrieval is indirectly estimated by comparing
the IWC values calculated from this mass vector to the
reference IWC values measured by the IKP-2 instrument.
A comparison between the calculated and measured IWC
values is presented in Figs. 5c–e for the three datasets, and
the quantitative results are summarized in Table 5. Un-
surprisingly, the best agreement (IWC MAPE = 10.5%) is
found for D16, which supports the minimization process.
The IWC prediction is equally satisfying for D19
(IWC MAPE = 12.1%), where efforts have been made to
sort data points into a coherent subset with respect to
sampling temperature. Perhaps most interesting are the
<table><caption>TABLE 5. Optimization results for HAIC-HIWC partial dataset.</caption><thead><tr><th rowspan="2">Regularization parameter</th><th rowspan="2">Iterations</th><th colspan="3">IWC MAPE(%)</th></tr><tr><th>D16 (1051 data points)</th><th>D19 (212 data points)</th><th>C26 (1298 data points)</th></tr></thead><tbody><tr><td>λ<sub>opt</sub> = 7.5 × 10<sup>6</sup></td><td>2</td><td>10.5</td><td>12.1</td><td>15.1</td></tr></tbody></table>

results of the IWC retrieval from C26 data. In Fig. 5c, points are color coded with temperature and three distinct groups are clearly visible (dark blue for the 229.4 K level, light blue for 244.2 K, and yellow for 264 K). Data points measured at a sampling temperature close to that of D16 yield the most accurate IWC predictions, as they are well distributed along the x = y line. In the other two temperature subsets, the accuracy of IWC prediction decreases as the sampling temperature departs from 236.1 K. At 244.2 K, IWC<sub>cal</sub> compares well with IWC<sub>th</sub>, although the method produces underestimated IWC values. Finally, IWC<sub>th</sub> values are generally underestimated by a factor of 1.5–4 for the points measured at 264 K. The good overall agreement found for C26 (IWC MAPE = 15.1%) contrasts with the noticeable differences observed between the predicted and measured IWC values for the two warmest temperature levels: it can be explained by the dominance of data points sampled at 229.4 K in the C26 dataset (936 out of 1298; see Table 4). These results suggest that 1) ice particles observed in the D16, D19, and C26 MCSs at comparable temperature levels may have similar mass properties and 2) in the same cloud, particles sampled at different temperature levels may have different mass properties.

Since the good agreement between calculated and measured IWC values observed for D19-244.2 K and C26-229.4 K may hide compensating errors, the next step will be to apply the method on subsets sampled at the same temperature level in different clouds and compare the retrieved quantities. Subsequently, the method will be applied to subsets sampled at different temperature levels in order to evaluate whether a first-order approximation model describing the mass of ice particles as a function of their size and the cloud temperature produces consistent results. If not, we will expand the study to a multiparameter approach and search for potential correlations between the crystal mass and some other yet undefined parameters by applying the method to different subsets with specific environmental conditions. This, as well as any specific findings in ice particle microphysical properties and the implications for modeling and remote sensing retrieval algorithms, is well beyond the scope of this technical paper, which is focusing on the development of the technique and the numerical tool.

Last but not least, the convergence of the optimization algorithms is reached after a very small number of

iterations. Practically, the mass vector plotted in Fig. 5 is computed within few seconds on an ordinary computer. It underlines the efficiency of the present method in retrieving particle mass from large datasets of collocated and simultaneous PSD and IWC in situ measurements.

## 5. Conclusions and perspectives

This article presents a new approach for the retrieval of ice particles’ size-dependent masses from in situ measurements. A conceptual model is formulated as a linear system of equations relating particle mass to PSD and IWC measurements. The mass retrieval process consists of solving the inverse problem with numerical optimization tools. The strength of this technique is to make optimal use of common in situ measurements in order to waive the power-law assumption that is classically employed to constrain $m-D$ relationships, although it has been proved inaccurate by recent studies. The new technique simplifies the computation of ice particle mass from large in situ measurement datasets and proves its worth when applied to aircraft data collected during HAIC-HIWC field campaigns: IWC values estimated from the retrieved particle masses match the measured values with 10%–15% accuracy. The preliminary $m-D$ test calculations reveal some interesting size-dependent features otherwise hindered by the power-law assumption. In future work, this data processing technique will be applied to process the whole HAIC-HIWC dataset with the overall objective to develop a statistical model of ice particle mass, taking into account temperature and any other relevant physical parameters (IWC value, cloud type, etc.) susceptible to influence the mass of ice particles. The method is thought to be directly applicable to any collocated PSD and IWC measurement datasets, provided the number of available data points is compatible with the number of unknown elements of the mass vector to be solved. The authors believe such an approach could be adapted to retrieve ice particle mass from PSD and other bulk measurements, such as radar reflectivity factors, as long as an inverse problem can be formulated.
However, the model and mathematical tools have been developed heuristically, which leaves room for improvements in both the formulation of the forward model and the numerical tools selected to solve the inverse problem. On the one hand, the proposed forward model uses particles' 2D imagery as input. It makes the assumption that the particle volume is correctly described by a size parameter derived from the 2D projected images: this greatly simplifies the problem into a linear system of equations but the assumption is known to be inaccurate. Further studies on the relation between particles' 2D projected images and 3D volume could help in establishing a new forward model relating the mass of ice particles to their volume. On the other hand, the technical choices made to solve the inverse problem are open for discussion. In this article a regularized least squares approach is adopted, although a regularized total least squares approach could probably help in handling uncertainties in the PSD data. Since mass is a positive quantity, positivity or box-constrained optimization approaches and a penalty term, such as a logarithmic barrier function, could be advantageously tried. Our regularization term is defined to minimize the roughness of the curve and to account for the nonlinear increase of mass with size, with use of a weighted first-order difference operator. Since the choice of the prior information about the solution is qualitative, other definitions of the regularization terms, such as such those including density-based properties, should be tried.

**Acknowledgments.** The research leading to these results received funding from (i) the European Union's Seventh Framework Programme in research, technological development, and demonstration under Grant Agreement ACP2-GA-2012-314314; (ii) the European Aviation Safety Agency (EASA) Research Program under Service Contract EASA.2013.FC27; and (iii) the Federal Aviation Administration (FAA), Aviation Research Division and Aviation Weather Division, under Agreement CON-I-1301 with the Centre National de la Recherche Scientifique. Funding to support the Darwin flight project was also provided by the NASA Aviation Safety Program, the Boeing Co., and Transport Canada. Additional support was also provided by Airbus Operations SAS, Science Engineering Associates, the Bureau of Meteorology, Environment Canada (now Environment and Climate Change Canada), the National Research Council of Canada, the University of Utah, and the National Science Foundation under Grant AGS 12-13311. Special thanks to SAFIRE for logistical support and for flight mission performance. We are grateful to J. Walter

Strapp for having processed the IKP-2 data that were used in section 4. We are also thankful to the anonymous reviewers for their constructive comments on uncertainty consideration, which helped to improve the paper.

## REFERENCES

Baker, B., and R. P. Lawson, 2006: Improvement in determination of ice water content from two-dimensional particle imagery. Part I: Image-to-mass relationships. *J. Appl. Meteor. Climatol.*, **45**, 1282–1290, doi:10.1175/JAM2398.1.

Baumgardner, D., and Coauthors, 2011: Airborne instruments to measure atmospheric aerosol particles, clouds and radiation: A cook's tour of mature and emerging technology. *Atmos. Res.*, **102**, 10–29, doi:10.1016/j.atmosres.2011.06.021.

—, and Coauthors, 2017: Cloud ice properties: In situ measurement challenges. *Ice Formation and Evolution in Clouds and Precipitation: Measurement and Modeling Challenges, Meteor. Monogr.*, No. 58, Amer. Meteor. Soc., doi:10.1175/AMSMONOGRAPHHS-D-16-0011.1.

Borsdorff, T., O. Hasekamp, A. Wassmann, and J. Landgraf, 2014: Insights into Tikhonov regularization: Application to trace gas column retrieval and the efficient calculation of total column averaging kernels. *Atmos. Meas. Tech.*, **7**, 523–535, doi:10.5194/amt-7-523-2014.

Brown, P. R. A., and P. N. Francis, 1995: Improved measurements of the ice water content in cirrus using a total-water probe. *J. Atmos. Oceanic Technol.*, **12**, 410–414, doi:10.1175/1520-0426(1995)012<0410:IMOTIW>2.0.CO;2.

Cotton, R. J., and Coauthors, 2013: The effective density of small ice particles obtained from *in situ* aircraft observations of mid-latitude cirrus. *Quart. J. Roy. Meteor. Soc.*, **139**, 1923–1934, doi:10.1002/qj.2058.

Davison, C. R., J. MacLeod, and J. W. Strapp, 2009: Naturally aspirating isokinetic total water content probe: Evaporator design and testing. *Proc. First AIAA Atmospheric and Space Environments Conf.*, San Antonio, TX, American Institute of Aeronautics and Astronautics, AIAA 2009-3861, doi:10.2514/6.2009-3861.

—, J. W. Strapp, L. E. Lilie, T. P. Ratvasky, and C. Dumont, 2016: Isokinetic TWC evaporator probe: Calculations and systemic error analysis. *Proc. Eighth AIAA Atmospheric and Space Environments Conf.*, Washington, DC, American Institute of Aeronautics and Astronautics, AIAA 2016-4060, doi:10.2514/6.2016-4060.

Dezitter, F., A. Grandin, J.-L. Brenguier, F. Hervy, H. Schlager, P. Villedieu, and G. Zalamansky, 2013: HAIC (high altitude ice crystals). *Proc. Fifth AIAA Atmospheric and Space Environments Conf.*, San Diego, CA, American Institute of Aeronautics and Astronautics, AIAA 2013-2674, doi:10.2514/6.2013-2674.

Erfani, E., and D. L. Mitchell, 2016: Developing and bounding ice particle mass- and area-dimension expressions for use in atmospheric models and remote sensing. *Atmos. Chem. Phys.*, **16**, 4379–4400, doi:10.5194/acp-16-4379-2016.

Fontaine, E., 2014: Masse des cristaux de glace et facteurs de réflectivité radar dans les systèmes de nuages convectifs de moyenne échelle formés dans les tropiques et la région de la mer méditerranée. Ph.D. thesis, Université Blaise Pascal, 192 pp., http://www.theses.fr/2014CLF22527/document.
—, A. Schwarzenboeck, J. Delanoë, W. Wobrock, D. Leroy, R. Dupuy, C. Gourbeyre, and A. Protat, 2014: Constraining mass–diameter relations from hydrometeor images and cloud radar reflectivities in tropical continental and oceanic convective anvils. *Atmos. Chem. Phys.*, **14**, 11 367–11 392, doi:10.5194/acp-14-11367-2014.

Hansen, P. C., 1992: Analysis of discrete ill-posed problems by means of the L-curve. *SIAM Rev.*, **34**, 561–580, doi:10.1137/1034115.

Heymsfield, A. J., 1972: Ice crystal terminal velocities. *J. Atmos. Sci.*, **29**, 1348–1357, doi:10.1175/1520-0469(1972)029<1348:ICTV>2.0.CO;2.

—, A. Bansemer, C. Schmitt, C. Twohy, and M. R. Poellot, 2004: Effective ice particle densities derived from aircraft data. *J. Atmos. Sci.*, **61**, 982–1003, doi:10.1175/1520-0469(2004)061<0982:EIPDDF>2.0.CO;2.

—, C. Schmitt, A. Bansemer, and C. Twohy, 2010: Improved representation of ice particle masses based on observations in natural clouds. *J. Atmos. Sci.*, **67**, 3303–3318, doi:10.1175/2010JAS3507.1.

Huckle, T., and M. Sedlacek, 2012: Data based regularization matrices for the Tikhonov-Phillips regularization. *Proc. Appl. Math. Mech.*, **12**, 643–644, doi:10.1002/pamm.201210310.

Idier, J., Ed., 2008: *Bayesian Approach to Inverse Problems*. ISTE Ltd. and John Wiley & Sons Inc., 381 pp., doi:10.1002/9780470611197.

Jackson, R. C., and Coauthors, 2012: The dependence of ice microphysics on aerosol concentration in arctic mixed-phase stratus clouds during ISDAC and M-PACE. *J. Geophys. Res.*, **117**, D15207, doi:10.1029/2012JD017668.

Kajikawa, M., 1972: Measurement of falling velocity of individual snow crystals. *J. Meteor. Soc. Japan*, 50, 577–584, doi:10.2151/jmsj1965.50.6_577.

Lawson, R. P., D. O'Connor, P. Zmarzly, K. Weaver, B. Baker, Q. Mo, and H. Jonsson, 2006: The 2D-S (stereo) probe: Design and preliminary tests of a new airborne, high-speed, high-resolution particle imaging probe. *J. Atmos. Oceanic Technol.*, **23**, 1462–1477, doi:10.1175/JTECH1927.1.

Leroy, D., E. Fontaine, A. Schwarzenboeck, and J. W. Strapp, 2016: Ice crystal sizes in high ice water content clouds. Part I: On the computation of median mass diameter from in situ measurements. *J. Atmos. Oceanic Technol.*, **33**, 2461–2476, doi:10.1175/JTECH-D-15-0151.1.

Locatelli, J. D., and P. V. Hobbs, 1974: Fall speeds and masses of solid precipitation particles. *J. Geophys. Res.*, **79**, 2185–2197, doi:10.1029/JC079i015p02185.

Mace, J., E. Jensen, G. McFarquhar, J. Comstock, T. Ackerman, D. Mitchell, X. Liu, and T. Garrett, 2009: SPARTICUS: Small

Particles in Cirrus—Science and operations plan. DOE Tech. Rep. DOE/SC-ARM-10-003, 15 pp., https://www.arm.gov/publications/programdocs/doi-sc-arm-10-003.pdf.

McFarquhar, G. M., M. S. Timlin, R. S. Rauber, B. F. Jewett, J. A. Grim, and D. P. Jorgensen, 2007: Vertical variability of cloud hydrometeors in the stratiform region of mesoscale convective systems and bow echoes. *Mon. Wea. Rev.*, **135**, 3405–3428, doi:10.1175/MWR3444.1.

Mitchell, D. L., 1996: Use of mass- and area-dimensional power laws for determining precipitation particle terminal velocities. *J. Atmos. Sci.*, **53**, 1710–1723, doi:10.1175/1520-0469(1996)053<1710:UOMAAD>2.0.CO;2.

—, R. Zhang, and R. L. Pitter, 1990: Mass-dimensional relationships for ice particles and the influence of riming on snowfall rates. *J. Appl. Meteor.*, **29**, 153–163, doi:10.1175/1520-0450(1990)029<0153:MDRFIP>2.0.CO;2.

Nakaya, U., and T. J. Terada, 1935: Simultaneous observations of the mass, falling velocity and form of individual snow crystals. *J. Fac. Sci., Hokkaido Imp. Univ., Series 2*, **1**, 191–200.

Nocedal, J., and S. J. Wright, 2006: *Numerical Optimization*. Springer Series in Operations Research, Springer-Verlag, 664 pp., doi:10.1007/978-0-387-40065-5.

Pruppacher, H. R., and J. D. Klett, 2010: *Microphysics of Clouds and Precipitation*. Atmospheric and Oceanographic Sciences Library, Vol. 18, Springer, 954 pp., doi:10.1007/978-0-306-48100-0.

Raschke, E., J. Schmetz, and J. Heintzenberg, 1990: The International Cirrus Experiment (ICE)—A joint European effort. *ESA J.*, **14**, 193–199.

Schmitt, C. G., and A. J. Heymsfield, 2009: The size distribution and mass-weighted terminal velocity of low-latitude tropopause cirrus crystal populations. *J. Atmos. Sci.*, **66**, 2013–2028, doi:10.1175/2009JAS3004.1.

—, —, 2010: The dimensional characteristics of ice crystal aggregates from fractal geometry. *J. Atmos. Sci.*, **67**, 1605–1616, doi:10.1175/2009JAS3187.1.

Strapp, J. W., and Coauthors, 2016a: The high ice water content study of deep convective clouds: Science and technical plan. FAA Final Rep. DOT/FAA/TC-14/31, 105 pp.

—, L. E. Lilie, T. P. Ratvasky, C. R. Davison, and C. Dumont, 2016b: Isokinetic TWC evaporator probe: Development of the IKP2 and performance testing for the HAIC-HIWC Darwin 2014 and Cayenne Field Campaigns. *Proc. Eighth AIAA Atmospheric and Space Environments Conf.*, Washington, DC, American Institute of Aeronautics and Astronautics, AIAA 2016-4059, doi:10.2514/6.2016-4059.

Zikmunda, J., and G. Vali, 1972: Fall patterns and fall velocities of rimed ice crystals. *J. Atmos. Sci.*, **29**, 1334–1347, doi:10.1175/1520-0469(1972)029<1334:FPAFVO>2.0.CO;2.