Technical note: Bimodal parameterizations of in situ ice cloud particle size distributions

The cloud particle size distribution (PSD) is a key parameter for the retrieval of microphysical and optical properties from remote-sensing instruments, which in turn are necessary for determining the radiative effect of clouds. Current representations of PSDs for ice clouds rely on parameterizations that were largely based on aircraft in situ measurements where the distribution of small ice crystals were uncertain. This makes current parameterizations deficient to simulate remote-sensing observations sensitive to small ice, such as from lidar and thermal infrared instruments. In this study we fit the in situ PSDs of ice crystals from the JULIA (JÜLich In situ Aircraft data set) database, which consists of 11 campaigns covering the tropics, midlatitudes and the Arctic, consistently processed and considered more robust in their measurements of small ice. For the fitting, we implement an established approach to PSD parameterizations, which consists of finding an adequate set of parameters for a modified gamma function after normalization of both PSD axes. These parameters are constrained to match in situ measurements when predicting microphysical properties from the PSDs, via a cost function minimization method. We selected the ice water content and the ice crystal number concentration, which are currently key parameters for modern satellite retrievals and model microphysics schemes. We found that a bimodal parameterization yields better results than a monomodal one. The bimodal parameterization has a lower spread for almost all ice crystal sizes over the entire range of analyzed temperatures and fits better the observations, especially for particles between 20 and about 110 µm at temperatures between −60 and −20 ^∘C. For this temperature range, the root mean square error for the retrieved N_ice is reduced from 0.36 to 0.20. This demonstrates a clear advantage to considering the bimodality of PSDs, e.g., for satellite retrievals.

1 Introduction

Ice clouds play an important role in Earth's radiative budget, as their radiative effects can either contribute to a warming or a cooling of the surface (Liou, 1986; Stephens et al., 1990). The balance between these two effects depends on their macro- and microphysical properties, which stem from an array of very complex processes governing ice crystal formation and growth (Zhang et al., 1999; Krämer et al., 2020). Despite extensive research over the past decades most of these processes remain highly uncertain, making ice clouds major unknowns in current climate studies (Bellouin et al., 2020). Ice clouds are particularly challenging for satellite remote-sensing techniques, in great part due to the complexity and variety of their microphysical and optical properties (Baran, 2009). In turn, this lack of accurate global observational constraints leads to critical shortfalls for evaluation efforts of predictions of ice cloud properties and processes in numerical weather forecast and climate models (Lohmann et al., 2007).

The particle size distribution (PSD), which describes the number concentration of ice crystals as a function of their size, impacts most microphysical and radiative properties of clouds. This makes the PSD a central parameter in remote-sensing retrieval techniques (Yang et al., 2001; Vidot et al., 2015). PSD shapes can greatly vary, as they are influenced by formation and growth mechanisms dictated by the environment in which the cloud developed (Heymsfield and McFarquhar, 2002). For instance, the cloud origin (e.g., in situ or liquid origin;  Krämer et al., 2016) was identified as having a major influence on PSD shapes for cirrus (Luebke et al., 2016). However, the exact drivers of PSD shapes remain poorly understood, and their global variability is only addressed by few modeling studies (e.g., Gasparini and Lohmann, 2016). For this reason, retrieval techniques must make critical simplifications regarding PSD shapes, which are commonly reflected by the use of fixed parameterizations that are universally applied globally and for all cloud types. Such recent parameterizations rely on advanced normalization procedures that aim to make PSDs more representative of a large sample of ice clouds (Delanoë et al., 2005; Field et al., 2007, also see Sect. 3). While this allows for more accurate representations of ice clouds in retrieval methods (e.g., Sourdeval et al., 2015, 2016), current parameterizations still have important issues. For instance, they largely rely on in situ measurements where the distribution of small ice crystals (sizes less than 100 µm diameter) are at best very uncertain (Korolev et al., 2011), which makes them deficient in simulating remote-sensing observations sensitive to small ice and leads to erroneous retrievals, especially for estimations of the ice crystal number concentration (N_ice). Reaching more accurate estimations of N_ice is particularly of great importance to understand ice cloud formation mechanisms and for improving their predictions in models, since this parameter is often used as a prognostic variable to predict the evolution of clouds (e.g., Seifert and Beheng, 2006).

An example of a widely used application of normalized PSD parameterizations in satellite retrievals is found in the DARDAR algorithm ($\mathrm{raDAR}/\mathrm{liDAR}$; Delanoë and Hogan, 2008, 2010). DARDAR retrieves vertical profiles of ice cloud microphysical properties by using a combination of coincident measurements from the Cloud Profiling Radar (CPR) onboard CloudSat and from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP). Although the original algorithm does not include retrievals of N_ice, Sourdeval et al. (2018) explored the capabilities of its framework to retrieve this parameter (DARDAR-Nice). However, comparisons of N_ice retrievals by DARDAR-Nice to recent in situ measurements highlighted strong limitations linked to the parameterizations of PSDs used in the algorithm, which are in specific cases not suited to retrieving N_ice (Sourdeval et al., 2018; Krämer et al., 2020). A main reason is that the current parameterization only represents a single ice crystals size mode, which can result in PSDs being significantly misrepresented when they happen to be strongly bimodal, i.e., consisting of two size modes. The bimodality of the PSDs can often occur in cirrus clouds and is a function of temperature and location within the cloud (Zhao et al., 2011), with warmer clouds more often having bimodal PSDs (Jackson et al., 2015). Associated with the study of bimodality is the analysis of the effect of shattering ice crystals at the inlets of the in situ instruments that can artificially increase the concentration of small particles. This problem has been widely discussed, and as a result, new inlets were designed and correction algorithms were applied to minimize this effect (Korolev et al., 2011; Lawson, 2011; Krämer et al., 2020).

The present study investigates the benefits of considering a second mode in existing PSD parameterizations methods. It follows the framework proposed by Delanoë et al. (2005) and Delanoë et al. (2014) (hereafter D05 and D14), which is applied to an extensive database of in situ ice cloud properties from airborne campaigns (JULIA; Krämer et al., 2016), which is considered less sensitive to shattering effects. This study proposes a new set of parameterizations, based on single- and double-mode PSDs, that will be useful to improve remote-sensing retrieval methods. Section 2 describes the JULIA database and Sect. 3 details the PSD normalization method. Section 4 presents and analyzes the newly developed parameterizations by comparison to the original in situ data. Section 5 concludes this study.

2 The JULIA database

JULIA (JÜLich In situ Airborne database) is a compilation of in situ measurements of cirrus, mixed-phase and liquid clouds, water vapor, and also other trace gases, collected at the Research Center Jülich starting in 2008 (Schiller et al., 2008; Luebke et al., 2013; Costa et al., 2017; Krämer et al., 2016, 2018, 2020) and continued until 2021. The measurements were taken on board different research aircraft and cover the tropics, midlatitudes and the Arctic (Fig. 1). For our study, we focus only on the data of ice particles, which include data from 11 field campaigns with measurements of ice crystals with diameters of 3 to 1000 µm taken every second. The analyzed data, i.e., the N_ice and ice water content (IWC), were sampled with (or computed from) the NIXE-CAPS (New Ice eXpEriment: Cloud and Aerosol Particle Spectrometer), which is a combination of CAS (Cloud and Aerosol Spectrometer) and a CIPg (Cloud Imaging Probe – grayscale), and combinations of the FCDP (Fast Cloud Droplet Probe), 2D-S (two-dimensional stereo), CDP (Cloud Droplet Probe) and 2D-C (two-dimensional cloud) instruments (the respective size ranges used are listed in Table 2). These instruments and their data processing procedures are described in the literature (e.g., Lawson et al., 2006; McFarquhar et al., 2007; Lawson, 2011; Krämer et al., 2016; Luebke et al., 2016; Baumgardner et al., 2017; Afchine et al., 2018; Krämer et al., 2020). Since the widths of the size bins of each pair of instruments are different, Krämer et al. (2022) synchronized all PSDs to the same grid with logarithmically equidistant size bins in order to facilitate inter-comparisons and the interpretation of the results. The data included in JULIA have undergone an extensive quality check process to guarantee their validity. In the following section, a brief description of the campaigns and the instruments is given.

Figure 1 Location of all campaigns of the JULIA database. In different colors, the campaigns used for this study; other campaigns in gray (figure by Nicole Spelten).

2.1 Selected airborne campaigns

The analyzed data cover a wide range of meteorological conditions and extend from the tropics to the polar region. In total there are five campaigns in the tropics (ACRIDICON-CHUVA, ATTREX, CONTRAST, POSIDON and STRATOCLIM), three in midlatitudes (COALESC, ML-CIRRUS and START), one covering both midlatitudes and high latitudes (CIRRUS-HL), and two at high latitudes (RACEPAC and VERDI). ACRIDICON consisted of flights over the Amazonian forest in September 2014 with the aim of studying the interaction between aerosols and deep convective clouds (Wendisch et al., 2016). During ATTREX (February–March 2014), measurements inside cirrus clouds were taken over the western tropical Pacific, around the tropical tropopause layer, avoiding convection (Thornberry et al., 2017). CONTRAST, based on Guam, also took place in the western tropical Pacific (January–February 2014) to study tropical convection, oceanic processes and ozone chemistry in the upper troposphere–lowermost stratosphere (UTLS) (Pan et al., 2017). POSIDON (October 2016) was also based on Guam, and among its objectives were the study of cirrus clouds and dehydration in the tropical UTLS (POSIDON, 2022). During POSIDON three flights were near an active typhoon. The last campaign in the tropics is STRATOCLIM (July–August 2017), based in Nepal and focused on water vapor variations and upper-tropospheric clouds (Krämer et al., 2020; Khaykin et al., 2022). COALESC took place during February and March 2011, and the observations were taken over the southeast coast of England and Wales. Its main focus was on stratocumulus clouds, but there are also plenty of observations within mixed-phase and cirrus clouds (Osborne et al., 2014). During March and April 2014, ML-CIRRUS took place, centering its observations over Europe and the North Atlantic, focusing on processes involving cirrus (Voigt et al., 2017). START was based in Colorado, USA, during April–June 2008 (Pan et al., 2010). CIRRUS-HL is a campaign from July–August 2021, based in Oberpfaffenhofen, Germany, that investigated the microphysical properties and climate impact of ice clouds at high latitudes as well as the effect of aviation (CIRRUS-HL, 2021). VERDI (April–May 2012, Inuvik, Canada) (Costa et al., 2017) aimed to investigate, among other goals, the radiation budget of ice clouds in the Arctic and the influence of convective transport in the upper troposphere. During this campaign a persistent anticyclone was present and favored the formation of persistent stratus (Klingebiel et al., 2015). RACEPAC (Costa et al., 2017) is the followup campaign of VERDI and took place in April–May 2014.

In Table 1 a summary of the instruments used for each campaign is given. The NIXE-CAPS is formed by CAS and CIPg, which combined cover particles with diameters from 0.61 to 930 µm. To avoid overlap between particle sizes, the PSDs from NIXE-CAPS are merged between 20 and 25 µm. Detailed explanations about the data analysis methods, including the position of the instrument in the aircraft and shattering effects, are given by Meyer (2013), Krämer et al. (2016), Luebke et al. (2016), and Costa (2017). Other instruments are the light-scattering sensor CDP that measures concentration for particles with diameters between 2 and 50 µm (McFarquhar et al., 2007) and the FCDP for particles between 1 and 50 µm (Baumgardner et al., 2017). The optical imaging cloud probe 2D-C spectrometer (Baumgardner et al., 2017) and the 2D-S spectrometer are used to reconstruct cloud particle shapes and sizes between 25–800 and 5–1280 µm, respectively. The 2D-S includes tips and software to reduce shattering effects (Lawson et al., 2006; Lawson, 2011). FCDP and 2D-S PSDs are merged between 20 and 25 µm, and the ones from CDP and 2D-C are merged at 55 µm. For more information about the campaigns and the instruments, the reader is also referred to Costa (2017) and Krämer et al. (2020). As mentioned in Sect. 1, shattering of the ice particles during the measurements would increase the number of small particles and cause an artificial bimodality in the PSDs. However, as presented in the above references, major efforts were made in the development of anti-shatter probe tips and particle interarrival time algorithms that have resulted in a successful minimization of the shattering of ice particles (see e.g., Krämer et al., 2020). Therefore, we are confident that the bimodality present in the JULIA database is not due to distorted microphysical properties of the PSDs.

Table 1 Instruments used during the field campaigns of the JULIA database (NIXE-CAPS as a combination of CAS, Cloud and Aerosol Spectrometer, and CIPg, Cloud Imaging Probe greyscale; CDP: Cloud Droplet Probe; FCDP: Fast Cloud Droplet Probe; 2D-C: two-dimensional cloud spectrometer; 2D-S: two-dimensional stereo cloud spectrometer). Range refers to the ice particle diameter.

^* The range of the instruments might be larger.

2.2 Ice crystal number concentration and ice water content computation

N_ice is the sum of the ice crystal number concentration of each size bin. The IWC is obtained as follows:

$$\begin{array}{}\text{(1)}& \mathrm{IWC}=\sum _{i=\mathrm{1}}^n{m}_i\cdot \mathrm{\Delta }{N}_{\mathrm{ice},i},\end{array}$$

where m_i and ΔN_ice,i are the mass and number concentration of the ice crystals in the ith size bin, respectively. The mass m is computed using the mass − dimension (m–D) relation of the form $m=a\cdot D^b$ described in Krämer et al. (2016), which is based on the modified m–D relation of Mitchell et al. (2010). The coefficients a and b depend on the ice crystal size. The used m–D relation was compared in Afchine et al. (2018) with other m–D relations from the literature (covering, depending on the m–D relation, temperatures between −65 and 0 ^∘C) and also with the measurements from total water instruments showing good agreement for cirrus clouds. For temperatures higher than the cirrus range, we are aware that the uncertainties in the derived IWC are larger than at lower temperatures. However, we use the same m–D, since as shown in Afchine et al. (2018), the differences between the compared m–D relations is small, even when considering those derived for higher temperatures.

3 Normalization method

The method used to fit the PSDs is based on the methodology described by D05 and D14, which is an adaptation of the normalization method for raindrop spectra introduced by Testud et al. (2001). It consists of computing several moments of the measured in situ PSDs, using them to normalize the in situ PSDs and then fitting the normalized PSDs to a certain function. This universal function can then be used to obtain microphysical and optical properties of cirrus clouds. Here, we present the key points of the method; for a complete description, the reader is referred to the aforementioned references.

The first step is to compute the equivalent melted diameter (D_eq) as follows:

$$\begin{array}{}\text{(2)}& D_{\mathrm{eq}}={\left({\displaystyle \frac{\mathrm{6}\cdot m\left(D\right)}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where m(D) is obtained with the m–D relation indicated in Sect. 2.2 and D_eq is in units of meters.

Using the measured PSD, the moments of the distribution are computed as follows:

$$\begin{array}{}\text{(3)}& M_n=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}^n{N}_{\mathrm{ice}}\left(D\right)\mathrm{d}D\approx \sum _{D=\min }^{D=\max }{D}^n{N}_{\mathrm{ice}}\left(D\right)\mathrm{\Delta }D,\end{array}$$

with n the moment order, N_ice(D) the ice number concentration for the size bin D and ΔD the width of the corresponding size bin.

D05 and D14 use the volume-weighted diameter D_m to scale the PSD in the size space and the intercept parameter $N_{\mathrm{0}}^{*}$ in the number concentration space. These parameters can be defined in terms of the third and fourth moment of the PSD by

$$\begin{array}{}\text{(4)}& {\displaystyle }{D}_m={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{N}_{\mathrm{ice}}\left(D_{\mathrm{eq}}\right)D_{\mathrm{eq}}^{\mathrm{4}}\mathrm{d}{D}_{\mathrm{eq}}}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{N}_{\mathrm{ice}}\left(D_{\mathrm{eq}}\right)D_{\mathrm{eq}}^{\mathrm{3}}\mathrm{d}{D}_{\mathrm{eq}}}},\text{and}\text{(5)}& {\displaystyle }{N}_{\mathrm{ice}}\left(D_{\mathrm{eq}}\right)=N_{\mathrm{0}}^{*}F\left(D_{\mathrm{eq}}/D_m\right),\end{array}$$

considering that N_ice(D_eq)dD_eq=N_ice(D)dD for a given size bin. F is a modified gamma function that describes the shape of the normalized PSDs and, as first given by D05, is

$$\begin{array}{}\text{(6)}& F(\mathit{\alpha },\mathit{\beta },X)=\mathit{\beta }{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{4}\right)}{{\mathrm{4}}^{\mathrm{4}}}}{\displaystyle \frac{\mathrm{\Gamma }{\left(\frac{\mathit{\alpha }+\mathrm{5}}{\mathit{\beta }}\right)}^{\mathrm{4}+\mathit{\alpha }}}{\mathrm{\Gamma }{\left(\frac{\mathit{\alpha }+\mathrm{4}}{\mathit{\beta }}\right)}^{\mathrm{5}+\mathit{\alpha }}}}{X}^{\mathit{\alpha }}\exp \left[-{\left(X{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathit{\alpha }+\mathrm{5}}{\mathit{\beta }}\right)}{\mathrm{\Gamma }\left(\frac{\mathit{\alpha }+\mathrm{4}}{\mathit{\beta }}\right)}}\right)}^{\mathit{\beta }}\right].\end{array}$$

F is defined by four parameters: $N_{\mathrm{0}}^{*}$ and D_m (through $X=D_{\mathrm{eq}}/D_m$) that change for each PSD and α and β that are fixed and can be found by computing a cost function. The normalization is applied to each individual PSD for all campaigns.

The α and β that best fit the normalized measured PSD are chosen using an in situ database (in our study the JULIA data set) and a least square regression linear fit on moments of the PSD (Field et al., 2005, 2007; Delanoë et al., 2014). Following D14, we use a combination of a low and a high moment of the PSD. However, unlike previous studies, we here aim at improving the direct prediction of physical parameters of the PSDs and therefore minimize N_ice and IWC via the cost function. This cost function, J, is commonly used to quantify the consistency between predicted and in situ parameters. Considering the tendency of both IWC and N_ice to follow a log-normal distribution, we used the logarithmic values of these two parameters when computing J:

$$\begin{array}{}\text{(7)}& J=\sum _{i=\mathrm{1}}^n\left(J_{N_{\mathrm{ice}}}+J_{\mathrm{IWC}}\right),\end{array}$$

with $J_{N_{\mathrm{ice}}}$ and J_IWC as

$$\begin{array}{}\text{(8)}& {\displaystyle }{J}_{N_{\mathrm{ice}}}={\left(\mathrm{1}-{\displaystyle \frac{\log \left(N_{\mathrm{param}}(\mathit{\alpha },\mathit{\beta })\right)}{\log \left(N_{\mathrm{insitu}}\right)}}\right)}^{\mathrm{2}},\text{(9)}& {\displaystyle }{J}_{\mathrm{IWC}}={\left(\mathrm{1}-{\displaystyle \frac{\log \left({\mathrm{IWC}}_{\mathrm{param}}\right(\mathit{\alpha },\mathit{\beta }\left)\right)}{\log \left({\mathrm{IWC}}_{\mathrm{insitu}}\right)}}\right)}^{\mathrm{2}},\end{array}$$

where n is the total number of PSDs, N_insitu is the sum of the ice crystal number concentration from the in situ database in each size bin and IWC_insitu is derived using Eq. (1). N_param is computed by integration of the size distribution $N_{\mathrm{ice}}(\mathit{\alpha },\mathit{\beta },D)$ obtained from the normalized function F and using only the size bins that are present in the in situ data. IWC_param is derived by applying Eq. (1). The objective of the cost function approach is to find the optimal coefficient pair of α and β that will minimize J.

4 Results

4.1 Original and normalized PSDs

We understand all clouds colder than 235 K as cirrus (Krämer et al., 2016). In the temperature range just below 235 K, the clouds may originate as mixed-phase clouds ascending from lower altitudes, undergoing complete glaciation at ≥ 235 K. This physical definition of cirrus is based on the ice formation mechanism, which is on the one hand the just mentioned complete glaciation of liquid clouds (liquid origin cirrus) and on the other hand cirrus that forms directly as ice (in situ origin cirrus). The analyzed in situ data are limited to temperatures lower than 255 K. This choice was made to minimize the selection of water droplets in mixed-phase clouds, since below this temperature most of the water droplets are glaciated. No differentiation between contrail and natural cirrus was made in our analysis. In total, 542 719 PSDs were analyzed (≈ 151 h). Of the total number of analyzed PSDs, about 9.8 % are found at temperatures between 235 and 255 K (i.e.,  mixed-phase regime). The individual contribution of each campaign is registered in Table 2.

Table 2 Contribution in percent of each campaign to the total analyzed data considering only T<255 K. The columns indicate the name of the campaign, the number of hours analyzed, the number of PSDs of each campaign, and their contribution in percent. Total number of PSDs: 542 719. Total number of cumulative hours: 150.76.

Figure 2 (right column) shows how PSDs contain information about the characteristics of clouds and their history and may differ according to dynamical conditions. These differences can be seen in the frequency distributions of the concentrations of the PSDs of each campaign. For example, very high number concentrations (> 1 cm⁻³) of the smallest ice crystals can be due to the presence of contrails (e.g., during COALESC, ML-CIRRUS or CIRRUS-HL). As another example, during several flights of ACRIDICON-CHUVA and CIRRUS-HL strong convection was present (i.e., fast updrafts). During fast updrafts, air masses can be lifted and reach temperatures lower than 235 K, and if the conditions are favorable, ice nucleation will take place and result in the formation of small ice crystals. Also during strong updrafts small supercooled liquid water droplets glaciate, giving as a result a large number of small ice crystals. In the PSDs, that is translated into an increase in number concentration between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ –  1 cm⁻³ for ice crystal sizes smaller than ≈ 20 µm. Another indicator of cirrus clouds that have their origin at T<235 K is the presence of large ice crystals that come from heterogeneous drop freezing (≈ 200 µm; more noticeable for particles > 500 µm that reach a number concentration of $\approx \mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ cm⁻³). A more detailed analysis of the characteristics of the PSDs according to dynamical conditions will be given in Krämer et al. (2023).

Figure 2 (left column) shows the results of applying the normalization method described in Sect. 3 for a selection of campaigns (see figure legend for details). As described in D05 and D14, the normalization approach removes some natural variability in PSDs and therefore narrows the data down to smaller size and concentration areas (mostly around $D_{\mathrm{eq}}/D_m=\mathrm{1}$), in comparison to the original observed PSD. Another visible feature of the normalization is that the normalized spectra look similar for all campaigns, although there are some differences depending on the measurement location and the cloud processes involved (Delanoë et al., 2005; Field et al., 2007). Overall, these features are also visible in our study, including as well the larger variability in the PSD tail, which is linked to the temperature. In D14 it is explained that for very low temperatures (less than −60 ^∘C) the tail vanishes, and for higher temperatures the variability in D_m is larger, since the range of possible ice crystal sizes increases due to ice particles coming from different microphysical processes.

Figure 2 The right column shows size-resolved occurrence frequencies of ice crystals of the observed PSDs and the left column their corresponding normalized PSDs for a selection of campaigns covering the tropics (ACRIDICON-CHUVA, STRATOCLIM), midlatitudes (COALESC, ML-CIRRUS, CIRRUS-HL) and polar region (CIRRUS-HL); the bottom panels show the combination of all campaigns. In the panels of the left column, the black dashed and solid lines corresponds to the proposed parameterization by D05 and D14, respectively.

Figure 3 (a–c) Normalized functions for each campaign after selecting the best pair of coefficients. Panel (a) is the result of computing the cost function J using IWC and N_ice, (b) of using only IWC to compute J, and (c) of using only N_ice. Panel (d) shows the selection of the best coefficient pair for each campaign and the coefficients for the parameterization by D05 (black upwards triangle), the parameterization proposed in D14 (black downwards triangle) and the new one (black circle). For this case, the horizontal line symbol in (d) is coincidental with the symbol of the crosses. JULIA 1M refers to the monomodal parameterization obtained in this study. A detailed overview of the coefficients of each parameterization is given in Table 3. Only PSDs measured at T<255 K are considered.

4.2 Investigation of optimal PSD parameters

As described in Sect. 3, a cost function minimization approach is used in this study to find an optimal pair of α and β coefficients needed to define the normalized PSD parameterizations that will fit in situ measurements from JULIA. We seek to optimize the representation of physical properties, IWC and N_ice (see Eqs. 7 to 9), as opposed to optical properties as for instance in D05 and D14. This is done with the intent to produce a PSD parameterization that can be universally used for retrievals of N_ice and IWC using a wide variety of instruments.

In D05 and D14, whose parameterizations are here used for comparison purposes, the lidar extinction coefficient and radar reflectivity factor are used to find the best parameters for the modified gamma function. The selected α and β pairs in the study by D14 yield very similar normalized PSD shapes for each campaign, with almost a complete overlap of the normalized function for $D_m/D_{\mathrm{eq}}$ between 0.5 and 1.5 (see Fig. 9 of D14). Following Fig. 9 of D14, we also plotted the normalized function for each campaign (i.e., $N_{\mathrm{ice}}/N_{\mathrm{0}}^{*}$ vs. $D_m/D_{\mathrm{eq}}$) (Fig. 3). Figure.  3a shows the normalized function for each campaign when selecting the optimal parameters using both the log(IWC) and the log(N_ice) in the cost function. Figure 3b shows the result of only using log(IWC) for the cost function and Fig. 3c the result of only using log(N_ice). Figure 3d is an overview of the selected optimal parameters α and β for each campaign (in different colors) and for all together (black circle; JULIA 1M in Table 3). In our study, although the normalized functions also cluster around $D_m/D_{\mathrm{eq}}$ values between 0.5 and 1.5, the overlap is not as pronounced as in D14 and the spread in the α and β values is also larger (Fig. 3d). To try to obtain a more compact cluster of selected coefficient pairs, we divided the PSDs into ice crystals smaller and larger than 50 µm (not shown) (JULIA small and JULIA large in Table 3). Splitting the PSDs modifies the spread, but our cluster in any case looks as compact as in D14. This might be due to the selected parameters to compute J. Whereas D14 uses optical parameters, i.e., visible extinction and reflectivity, we use physical parameters, i.e., N_ice (which is very sensitive to temperature) and IWC. Moreover, although we follow the indication of applying one moment sensitive to the small particles and another one to the large particles, the selected moments are not the same. In D14 the parameters are proportional to the second and approximately sixth moment of the distribution and in our study to the zeroth moment and between the second and third moment. To summarize, the parameterizations differ in the data used to compute each of them, the m–D relationship used and how the parameters of the modified gamma function were obtained. Table 3 gives an overview of the selected coefficients for each parameterization, and in Sect. 4.3 we discuss the characteristics of each of them.

Table 3 Best α and β coefficients for each parameterization. D05 is from Delanoë et al. (2005) and D14 from Delanoë et al. (2014). D14 was obtained by Delanoë et al. (2014) using extinction and reflectivity in the cost function. All other parameterizations are based on the JULIA database and the use of Eq. (7). JULIA 1M refers to one mode, JULIA large to D≥50 µm, JULIA small to D<50 µm, and JULIA 2M to the combination of the small and large modes.

Figure 4 Median of the percentage error between the single observed PSDs and their corresponding parameterized PSDs. (a, c) Comparison between one-mode (JULIA 1M: yellow) and two-mode parameterization (JULIA 2M: black); (b, d) comparison between parameterizations from the literature (D05: red; D14: blue) and the two-mode parameterization (JULIA 2M: black) of this study. The shaded region in panels (c) and (d) corresponds to the region between the 25th and 75th percentile. The x axis represents the size bins in micrometers. A detailed overview of the coefficients of each parameterization is given in Table 3. The numbers inside the panels correspond to the median percentage error over the complete size range.

Figure 5 Median of the percentage error between the single observed PSDs of all campaigns and their corresponding parameterized PSDs, as in Fig. 4c but for 10 ^∘C temperature intervals. The parameterizations from the literature are D05 in red and D14 in blue. JULIA 2M (in black) is the proposed new bimodal parameterization. The x axis is the size bins in micrometers. A detailed overview of the coefficients of each parameterization is given in Table 3. The numbers inside the panels correspond to the median percentage error over the specified size range.

Figure 6 Comparison between the in situ N_ice and parameterized N_ice for temperatures between −90 and −60 ^∘C ((a) D05, (d) JULIA 2M), temperatures between −60 and −20 ^∘C ((b) D05, (e) JULIA 2M), and temperatures lower than 255 K ((c) D05, (f) JULIA 2M). The y axis is the logarithm of the parameterized variable, and the x axis the logarithm of the observed one. Color code indicates frequency (%). The black line corresponds to the 1-to-1 line, the solid red line to the regression line and the dashed red lines to the regression line shifted by a factor of ± 2. The black ellipse is the 2σ area. A detailed overview of the coefficients of each parameterization is given in Table 3.

4.3 One-mode vs. two-mode parameterization

Ice crystals smaller than about 50 µm present a greater variability in the PSDs than the larger ice crystals (see Fig. 2, right column) because they correspond to the regime of newly formed ice crystals that quickly grow to larger, sedimenting sizes. In cirrus clouds, riming and secondary ice production do not play any role, and aggregation, at the lowest temperatures, is nearly negligible. These processes are of importance for mixed-phase clouds, which, as mentioned in Sect. 4.1 entail 9.8 % of the analyzed data. In Jackson et al. (2015) it was discussed that at temperatures lower than −45 ^∘C, the growth of the ice crystals is likely due to depositional growth, whereas sedimentation and aggregation are less significant. For higher temperatures, the ice particles also grow by vapor deposition, but sedimentation from above can also be a source of large particles that causes bimodality (Zhao et al., 2011). Another process that can lead to bimodality is a two-step ice nucleation. Initially, the few heterogeneously nucleated ice crystals may grow to larger sizes, followed by homogeneous nucleation of more and smaller ice crystals. However, the main reason for the bimodality of cirrus PSDs is the superposition of in situ origin and liquid origin cirrus. Ice crystals of liquid origin are significantly larger than those of in situ origin because they stem from lower altitudes where there is more water to allow them to grow to large sizes, especially since only very few drops out of a liquid cloud freeze, so the available water vapor is deposited only among them. Examples of measured bimodal PSDs are shown in Appendix A. In our study we compare the results of using one or two size modes, i.e., a parameterization for small particles (D<50 µm) combined with another parameterization for large particles (D≥50 µm). This cutting diameter agrees well with the division between the small and large modes when plotting the median PSD of all data (not shown).

Figure 4 shows a comparison between the median of the percentage error between the single PSDs for different parameterizations. In the upper row, Fig. 4a shows the comparison between the parameterizations JULIA 1M and JULIA 2M proposed in this study. Figure 4b shows the parameterizations from the literature and JULIA 2M. The bottom row (Fig. 4c–d) is the same as the upper row but adding for each parameterization the region between the 25th and 75th percentile in shaded color to illustrate their variability. It is important to note, when considering the large values of the percentage error, that between the number concentration of the largest ice crystals and the number concentration of the smallest crystals, there can be a difference of 6 orders of magnitude (Fig. 2).

From Fig. 4a and c it is clear that the bimodal parameterization (JULIA 2M: black) is closer to the observations than the proposed monomodal one (JULIA 1M: yellow) and that the variability is reduced for the bimodal parameterization. Between the bimodal parameterization and the ones from the literature (Fig. 4b and d), the clearer improvement is for ice crystals between 20 and ≈ 110 µm. For the small particles (smaller than 20 µm), D14 (dark blue) presents the highest deviation. This parameterization underestimates the number concentration of particles in this size range. In the case of D05 (red), the number concentration is also underestimated but to a lesser extent. The bimodal JULIA 2M (black) is the only parameterization that slightly overestimates the number concentration of the smaller particles. Figure 4d shows that the variability is considerable for all parameterizations, and the peak at around 20 µm is most probably caused by the merging of two instruments at this size. A detailed analysis by temperature of the comparison from Fig. 4b, presented in Fig. 5, shows that the lowest temperatures (−90 and −60 ^∘C) have the lowest deviation between observations and parameterizations for the smaller particles. For larger particles, all parameterizations tend to overestimate the concentration (PSD percentage error < 0), with JULIA 2M underestimating the concentration for particles larger than ≈ 600 µm. For temperatures between −60 and −50 ^∘C, where bimodality starts playing an important role, the parameterization from D05 (red) and D14 (blue) in Fig. 5d overestimate the number concentration for particles between 20 and 110 µm. This overestimation is also observed between −50 and −20 ^∘C. However, the tendency for particles smaller than 20 µm for D05 and D14 is to underestimate their concentration. For this range of temperatures, the bimodal parameterization is closer to the measured number concentrations, especially between 20 and 110 µm. As indicated by the median of the percentage error for particles smaller and larger than 50 µm, using a bimodal parameterization improves the representation of both the small and large modes, improving the large mode especially for higher temperatures.

In terms of IWC (not shown), correlation plots show that there are no significant differences between the results of the parameterizations (all parameterizations have a correlation factor of 1.0 and a RMSE of 0.18 or 0.19 when considering T<255 K). There is a slight underestimation (about 2 %) of the IWC for values between about $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ gm⁻³ and an overestimation between about $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ and 1 gm⁻³ (about 7 %). Since all parameterizations have a similar behavior for the large particles and IWC is sensitive to large particles (≳ 300 µm), this result was expected.

In terms of N_ice, correlation plots between the measurements and parameterizations (Fig. 6) show that the JULIA 2M parameterization (lower row) significantly reduces the spread seen in the monomodal D05 parameterization (upper row) and that the highest frequency for D05 is found slightly above the 1-to-1 line. Considering temperatures lower than 255 K, the root mean square error (RMSE) is 0.34 for D05 and 0.21 for the JULIA 2M parameterization (Fig. 6c and f, respectively). The bimodal parameterization presents a slightly higher correlation factor (0.98 vs. 0.96) than D05. Binning this data into 10 ^∘C temperature intervals between −90 and −60 ^∘C, where there is only one mode, JULIA 2M shows the lowest RMSE, with 0.22 vs. 0.30 in the monomodal D05 (Fig. 6a, d). For temperatures between −60 and −20 ^∘C, where two modes are clearly visible, D05 has a RMSE of 0.36 vs. a RMSE of 0.20 in the two-mode parameterization (Fig. 6b, e). A comparison of temperature intervals of 10 ^∘C (not shown) gives the same conclusions as already discussed: similar results for the lowest temperatures, with the two-mode parameterization being slightly closer to the observations, and better results of the two-mode parameterization for the higher temperatures.

These analyses confirm what was suggested by Sourdeval et al. (2018); i.e., defining a two-mode parameterization instead of a one-mode one improves the reconstruction of PSDs and retrieval of N_ice for higher temperatures because it adjusts better to the bimodal shape of the PSDs when it occurs.

5 Conclusions

In the present study, recent airborne in situ measurements of ice PSDs from the JULIA database were used to assess the ability of existing PSD parameterizations to represent N_ice and IWC, and to investigate the added value of considering two-mode PSDs. One of the main advantages of JULIA with respect to other data sets used for constructing current PSD parameterizations is that it includes observations of ice crystals down to 3 µm, which are consistently processed for all field campaigns, in particular to minimize the impact of ice-shattering effects. To find a new parameterization, we have followed the method by Delanoë et al. (2014), which consists of normalizing the PSDs and fitting them to a modified gamma function whose coefficients are chosen by minimizing a cost function J. The variables chosen to define the cost function were IWC and N_ice. We found that considering a possible bimodality of PSDs by combining two parameterizations, one for particles with D<50 µm and one for D≥50 µm, yields better results when compared to the observations than the hitherto used monomodal parameterizations. Considering a second mode improves the PSD prediction of both small and large ice crystals despite the large measurement uncertainties associated with the latter. Also, the variability of the retrieved PSDs when compared with the observations is reduced across all analyzed temperatures, and there is a better fit, especially for ice particles between 20 and ≈ 110 µm and for temperatures between −60 and −20 ^∘C. For this temperature range, the RMSE for the retrieved N_ice is reduced from 0.36 to 0.20. In conclusion, we propose here a new bimodal ice particle PSD parameterization including ice crystals smaller than 50 µm. An important next step would be to test the feasibility of implementing two parameterizations, one for smaller particles and another for larger particles, in the retrieval algorithms of remote-sensing instruments.

Appendix A: Examples of in situ observed bimodal ice particle PSDs

Figure A1 Selection of in situ observed bimodal PSDs belonging to different airborne campaigns included in the JULIA database. In the figures the corresponding campaign, flight, number of PSD within the flight and temperature are specified. The y axis is the concentration per cubic centimeter, and the x axis is the diameter in micrometers.

Figure A2 Selection of in situ observed bimodal PSDs belonging to different airborne campaigns included in the JULIA database. In the figures the corresponding campaign, flight, number of PSD within the flight and temperature are specified. The y axis is the concentration per cubic centimeter, and the x axis is the diameter in micrometers.

Figure A3 Selection of in situ observed bimodal PSDs belonging to different airborne campaigns included in the JULIA database. In the figures the corresponding campaign, flight, number of PSD within the flight and temperature are specified. The y axis is the concentration per cubic centimeter, and the x axis is the diameter in micrometers.

Figure A4 Selection of in situ observed bimodal PSDs belonging to different airborne campaigns included in the JULIA database. In the figures the corresponding campaign, flight, number of PSD within the flight and temperature are specified. The y axis is the concentration per cubic centimeter, and the x axis is the diameter in micrometers.