# Two-Dimensional Dynamics of Ice Crystal Parcels in a Cirrus Uncinus

Roland P. H. Berton

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

## To cite this version:

Roland P. H. Berton. Two-Dimensional Dynamics of Ice Crystal Parcels in a Cirrus Uncinus. Tellus A, 2023, 75 (1), pp.231 - 270. 10.16993/tellusa.3227. hal-04182292

HAL Id: hal-04182292

https://hal.science/hal-04182292

Submitted on 17 Aug 2023

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Two-Dimensional Dynamics of Ice Crystal Parcels in a Cirrus Uncinus

ORIGINAL RESEARCH
PAPER

ROLAND P. H. BERTON

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

STOCKHOLM
UNIVERSITY PRESS

## ABSTRACT

A simple 2D non-stationary parcel model has been elaborated by coupling kinematics and microphysics of ice crystals, for the purpose of retrieving realistic trajectories in a cirrus uncinus, with a hooked-shape head or cap, a trail or virga, and observed lifetimes. Prescribed kinematics and ambient physics includes vertical profiles of horizontal winds and updraft, and of air temperature and humidity. Modelled microphysics combines mass variation by deposition/sublimation of water vapour, radiative transfer and ventilation due to friction with ambient air. Upon changing the altitude of the updraft base, other parameters being unchanged, two distinct regimes of motions are found: setting the base at an altitude of 10 km yields a steady, non-oscillatory motion, with a hooked cap and a tail, typical of a cirrus uncinus, whereas fixing the base at a lower altitude (7–9 km), yields a self-sustained damped oscillatory motion with a long time period (7–10 hours) and a long wavelength (70–200 km), depending on supersaturation level and radiative conditions. This latter phenomenon may be related to the turrets occurring with “long-lasting” cirrus or to the Mesoscale Uncinus Complex (MUC) reported in published works. The sensitivity to the ice crystal size is analysed in order to estimate the spreading of trajectories, and the impact of radiative transfer is examined on the oscillatory regime at the light of three contrasting examples. A simple analytic model with kinematics and constant updraft is aimed at an outline of salient features and derivation of characteristic time- and spatial scales of the phenomena.

## CORRESPONDING AUTHOR:

Roland P. H. Berton

Office National d'Études et de Recherches Aérospatiales (ONERA)/DOTA, Université Paris Saclay, 91123 Palaiseau, France

roland.berton@onera.fr

## KEYWORDS:

cirrus clouds; ice crystal; thermodynamics; nonlinearity; damped oscillator

## TO CITE THIS ARTICLE:

Berton, RPH. 2023. Two-Dimensional Dynamics of Ice Crystal Parcels in a Cirrus Uncinus. *Tellus A: Dynamic Meteorology and Oceanography*, 75(1): 231-270.
DOI: https://doi.org/10.16993/tellusa.3227
# 1. INTRODUCTION

The modelling of atmospheric clouds has reached a high degree of accuracy, either for weather forecast or for scene simulation. Among these, cirrus clouds are of particular interest, because of their role in greenhouse and screening effects as well as weather concern (Liou et al. 1990; Fusina et al. 2007; Zhang et al. 2014; Köhler & Seifert 2015), and a considerable amount of observational and theoretical works has been devoted to them during the past 60 years. They occur at heights in the range 5–18 km, depending on the altitude of the tropopause between the Equator and the poles (Eguchi et al. 2007; Nazaryan et al. 2008). Of the two types of cirrus that are now distinguished according to their origin, anvil or in-situ (Lawson et al. 2019), we shall consider the second one.

As an alternative to the most widely used Eulerian methods, a simple two-dimensional Lagrangian approach to the modelling of cirrus uncinus, aiming at retrieving their characteristic hooked shape, is considered here. Thus we focus on the interplay between ambient atmospheric dynamics and the microphysics of ice crystals. The present work is guided by the concept of parcel trajectories (Marshall 1953; Heymsfield 1973; Heymsfield 1975a; Heymsfield 1975b; Heymsfield 1975c; Rogers 1975; Cohen et al. 2014) and the synergetic approach (Palmer 1996; Spreitzer et al. 2017), relevant of Lorentz's model, leading to the concept of strange attractor (Lorenz 1963). In such a perspective, the modelling of free convection was considered as an initial value problem (Saltzman 1962). Nevertheless, we do not address the question of cirrus formation itself, which is a difficult problem coped with by complete three-dimensional numerical models (Barahona & Nenes 2008; Lohmann et al. 2008; Köhler & Seifert 2015).

In a first step (Sect. 2.1), a kinematic model, composed of a set of ordinary differential equations (ODE) subject to initial conditions is described in a Lagrangian formulation, including the main sources of the motion: gravity (along the vertical), entrainment by ambient air motions, and drag due to the friction of ice particles (crystals) with air molecules (mainly oxygen and nitrogen). In next step (Sect. 2.2), the relevant parameterization of the microphysical model is presented, as regards the evolution of the crystal mass and shape during its rise and fall. In particular the modelling of radiative transfer (Sect. 2.2.2), and of crystal habit, namely a hollow hexagonal column (Sect. 2.2.3) is fully described. The atmospheric data needed by the coupled integration of the governing ODE, and suitable for mid-latitude cirrus clouds, is described in App. A; it includes ambient temperature, mass density, humidity, radiation flux, and the atmospheric motions driving crystal motion (horizontal winds and updraft).

In Sect. 3 are presented the solving method and results showing the existence of two distinct regimes or modes, depending on the base level of the updraft. These

are commented and compared with published models and observations. In a typical case (we call it a dual case) a steady mode, describing a well-developed cirrus uncinus, is generated for a “high” updraft (above 10 km), and an oscillatory regime is triggered for a “low” updraft (below 9 km). In contrast, in isolated cases an oscillatory mode is generated besides a steady mode reduced to a stretched hooked cap without developed trail.

In Sect. 4, the sensitivity of the trajectory to crystal size is investigated in order to estimate the spreading of the cirrus structure. The impact of radiative transfer is examined in Sect. 5, especially in connection with ice supersaturation level. These results are discussed in Sect. 6 in the light of various analytical models of time evolution of cirrus clouds. The two regimes found are considered in a mesoscale context (Orlanski 1975). In particular, the oscillatory mode is related to observations of “long-lasting cirrus turrets” (Heymsfield 1975b) and Mesoscale Uncinus Complex (Sassen et al. 1989).

In Sect. 7, a conclusion highlights salient features of the present approach and proposes further developments and extensions of our model. An analytic model of trajectory, described in App. D, will prove useful in order to appreciate the ability of the full model including microphysics to yield a realistic trajectory. An analytic model of trajectory (App. D) is aimed at showing the ability of the full model including microphysics to yield a realistic profile. A nomenclature of the mathematical notations used in the paper is summarized in Table 1.

# 2. MODEL

## 2.1. KINEMATICS

For simplicity, we assume that parcel motions are two-dimensional in a vertical plane (O,x,z), where the axes x and z are respectively horizontal and vertical, and z is directed upwards. Then the rectangular coordinates (x,z) of an average crystal's center of mass satisfy the second-order ODE (Pruppacher & Klett 1997):

$$ \begin{cases} \ddot{x} = -k(\dot{x} - U_o) - 2\Omega_o \dot{z} \cos \theta \\ \ddot{z} = -k(\dot{z} - W_o) - g \left(1 - \frac{\rho_o}{\rho_i}\right) + 2\Omega_o \dot{x} \cos \theta \end{cases} \quad (1a, b) $$

where single and double dots respectively mean first and second time-derivative, supplemented with the initial conditions:

$$ \begin{cases} x = x_0 \\ z = z_0 \end{cases} \quad \begin{cases} \dot{x} = 0 \\ \dot{z} = 0 \end{cases} \quad (2) $$

In Eqs. (1a, b) the first terms represent the entrainment acceleration driven by the motions of ambient air, and the last terms, the Coriolis acceleration (Igel & Biello 2020). In the equation of vertical momentum, Eq. (1b),
we included the gravity force corrected for the buoyancy effect (second term in parentheses in right-hand side), although it will reveal to have a negligible impact on the motion. Likewise, as will be checked in the results section, the damping coefficient *k* is approximately 20 s⁻¹ and the rotation rate Ω₀ is of order 10⁻⁴ rad s⁻¹ (Table 1) so that the Coriolis acceleration is negligible in our problem. Consistently, the flat-Earth assumption will be admitted, since the vertical and horizontal extension of the generated cirrus (respectively 2 km and 30–900 km) are much smaller than the Earth’s radius (≈6400 km). The horizontal dimensions take place in the so-called mesoscale range (Orlanski 1975). Moreover, Eqs. (1a, b) describing a variable-mass system with mass included in the coefficient *k* (see Eq. (3) below), an additional term, due to the momentum brought by the deposited mass at the rate *m* governed by Eq. (10) below, has been neglected because the mass loss or gain is isotropical (Plastino & Muzzio 1992).

<table><thead><tr><th>SYMBOL</th><th>MEANING</th></tr></thead><tbody><tr><td>A</td><td>Crystal area (m²)</td></tr><tr><td>a</td><td>Hexagonal column crystal radius (m)</td></tr><tr><td>B<sub>i</sub></td><td>Total radiance of black body (W m<sup>-2</sup>)</td></tr><tr><td>c</td><td>Hexagonal column crystal half-length (m)</td></tr><tr><td>c<sub>B</sub></td><td>Hexagonal column crystal core half-length (m)</td></tr><tr><td>c<sub>pa</sub></td><td>Specific heat of dry air (J K<sup>-1</sup> kg<sup>-1</sup>)</td></tr><tr><td>C</td><td>Capacitance of the ice crystal (m)</td></tr><tr><td>C<sub>D</sub></td><td>Drag coefficient</td></tr><tr><td>C<sub>p</sub></td><td>Phase speed of oscillations (m s<sup>-1</sup>)</td></tr><tr><td>D<sub>i</sub></td><td>Equivalent diameter of ice particle (m)</td></tr><tr><td>D<sub>v</sub></td><td>Diffusivity of water vapour in air (m<sup>2</sup> s<sup>-1</sup>)</td></tr><tr><td>f<sub>v</sub></td><td>Ventilation coefficient</td></tr><tr><td>F<sub>u</sub>, F<sub>d</sub></td><td>Upward, downward radiative flux densities (W m<sup>-2</sup>)</td></tr><tr><td>g</td><td>Gravity field constant (g = 9.81 m s<sup>-2</sup>)</td></tr><tr><td>H<sub>i</sub></td><td>Relative humidity of water vapour over liquid water</td></tr><tr><td>H<sub>v</sub></td><td>Relative humidity of water vapour over ice</td></tr><tr><td>ħ</td><td>Heating rate (K d<sup>-1</sup>)</td></tr><tr><td>IWP</td><td>Ice Water Path (g m<sup>-3</sup>)</td></tr><tr><td>k</td><td>Damping coefficient (s<sup>-1</sup>)</td></tr><tr><td>K<sub>a</sub></td><td>Heat conductivity of ambient dry air (W K<sup>-1</sup> m<sup>-1</sup>)</td></tr><tr><td>L<sub>s</sub></td><td>Latent heat of sublimation of ice (L<sub>s</sub> = 2.837 × 10<sup>6</sup> J kg<sup>-1</sup>)</td></tr><tr><td>L<sub>v</sub></td><td>Latent heat of vaporization of water (L<sub>v</sub> = 2.525 × 10<sup>6</sup> J kg<sup>-1</sup>)</td></tr><tr><td>m</td><td>Crystal mass (kg)</td></tr><tr><td>P</td><td>Period of oscillatory motion (s)</td></tr></tbody></table>

(Contd.)

<table><thead><tr><th>SYMBOL</th><th>MEANING</th></tr></thead><tbody><tr><td>p</td><td>Total pressure (Pa)</td></tr><tr><td>p<sub>r</sub></td><td>Reference pressure (p<sub>r</sub> = 10<sup>5</sup> Pa)</td></tr><tr><td>p<sub>a</sub></td><td>Partial pressure of dry air (Pa)</td></tr><tr><td>p<sub>v</sub></td><td>Partial pressure of water vapour (Pa)</td></tr><tr><td>p<sub>vsi</sub></td><td>Saturation pressure of water vapour over liquid water (Pa)</td></tr><tr><td>p<sub>vsi</sub></td><td>Saturation pressure of water vapour over ice (Pa)</td></tr><tr><td>q<sub>s</sub></td><td>Saturation mixing ratio (kg kg<sup>-1</sup>)</td></tr><tr><td>Q<sub>abs</sub></td><td>Efficiency factor for absorption of radiation (Q<sub>abs</sub> ≈ 1)</td></tr><tr><td>r<sub>i</sub></td><td>Equivalent radius of ice crystal (m)</td></tr><tr><td>R</td><td>Source term in radiative contribution (W)</td></tr><tr><td>R</td><td>Radiative term in the mass equation (%)</td></tr><tr><td>R<sub>a</sub></td><td>Gas constant of dry air (R<sub>a</sub> = 287.05 J K<sup>-1</sup> kg<sup>-1</sup>)</td></tr><tr><td>R<sub>v</sub></td><td>Gas constant of water vapour (R<sub>v</sub> = 461.00 J K<sup>-1</sup> kg<sup>-1</sup>)</td></tr><tr><td>Re</td><td>Reynolds number</td></tr><tr><td>Ri</td><td>Richardson number</td></tr><tr><td>S<sub>i</sub></td><td>Saturation ratio of water vapour over liquid water</td></tr><tr><td>S<sub>i</sub></td><td>Saturation ratio of water vapour over ice</td></tr><tr><td>T</td><td>Absolute temperature (K)</td></tr><tr><td>T<sub>r</sub></td><td>Reference temperature (T<sub>r</sub> = 273.15 K)</td></tr><tr><td>T<sub>a</sub></td><td>Temperature of ambient air (K)</td></tr><tr><td>T<sub>s</sub></td><td>Temperature of crystal surface (K)</td></tr><tr><td>u, w</td><td>Components of crystal velocity (m s<sup>-1</sup>)</td></tr><tr><td>U<sub>a</sub>, W<sub>a</sub></td><td>Components of ambient air velocity (m s<sup>-1</sup>)</td></tr><tr><td>V</td><td>Crystal volume (m<sup>3</sup>)</td></tr><tr><td>W<sub>f</sub></td><td>Free-fall speed (m s<sup>-1</sup>)</td></tr><tr><td>W<sub>f∞</sub></td><td>Terminal free-fall speed (m s<sup>-1</sup>)</td></tr><tr><td>x, z</td><td>Cartesian coordinates of ice crystal (m)</td></tr><tr><td>z<sub>a0</sub></td><td>Critical altitude of horizontal wind (m)</td></tr><tr><td>z<sub>w0</sub></td><td>Critical altitude of updraft (m)</td></tr><tr><td>α</td><td>Coefficient in Eq. (27)</td></tr><tr><td>ε</td><td>Ratio of gas constants R<sub>a</sub>/R<sub>v</sub> (ε = 0.622)</td></tr><tr><td>ζ</td><td>Damping ratio</td></tr><tr><td>η</td><td>Radiative transfer ratio</td></tr><tr><td>θ</td><td>Latitude</td></tr><tr><td>θ<sub>d</sub>, θ<sub>e</sub></td><td>Potential temperature (dry, equivalent)</td></tr><tr><td>Λ</td><td>Wavelength of oscillatory motion (m)</td></tr><tr><td>μ<sub>a</sub></td><td>Dynamic viscosity of air (Pa s)</td></tr><tr><td>ρ<sub>a</sub></td><td>Mass density of dry air (kg m<sup>-3</sup>)</td></tr><tr><td>ρ</td><td>Mass density of wet air (kg m<sup>-3</sup>)</td></tr><tr><td>ρ<sub>i</sub></td><td>Mass density of ice (ρ<sub>i</sub> = 920 kg m<sup>-3</sup>)</td></tr></tbody></table>

(Contd.)
<table><caption>**Table 1** Nomenclature.</caption><thead><tr><th>SYMBOL</th><th>MEANING</th></tr></thead><tbody><tr><td>σ</td><td>Stefan-Boltzmann constant<br>(σ = 5.67 × 10<sup>-8</sup> W m<sup>-2</sup> K<sup>-4</sup>)</td></tr><tr><td>σ<sub>s</sub></td><td>Vertical gradient of ice supersaturation (m<sup>-1</sup>)</td></tr><tr><td>τ</td><td>Non-dimensional time of analytic model</td></tr><tr><td>τ<sub>D</sub></td><td>Drag relaxation time (s)</td></tr><tr><td>φ</td><td>Aspect ratio (φ = c/a)</td></tr><tr><td>ψ</td><td>Hollowness factor (ψ = 1 - c<sub>B</sub>/c)</td></tr><tr><td>ω</td><td>Wind shear (m s<sup>-1</sup> km<sup>-1</sup> or s<sup>-1</sup>)</td></tr><tr><td>Ω</td><td>Angular frequency of oscillations (= 2π/P) (rad s<sup>-1</sup>)</td></tr><tr><td>Ω<sub>o</sub></td><td>Rotation rate of the Earth (Ω<sub>o</sub> = 7.27 × 10<sup>-5</sup> rad s<sup>-1</sup>)</td></tr><tr><td>Ω<sub>BV</sub></td><td>Brunt-Väisälä angular frequency (rad s<sup>-1</sup>)</td></tr></tbody></table>

The wind horizontal component $U_a$ is supposed to flow in a vertical plane (x,z). It is detailed analytically in App. A.1. and plotted in Figure 1. The wind vertical component consists in an updraft $W_a$ that is detailed analytically in App. A.1 and plotted in Figure 2. Concerning the trail formation, we consider the scenario in which the updraft ceases under a critical altitude $z_{w0}$ (Heymsfield 1975b).

<figure><img src="image_5.png" /><figcaption>**Figure 1** Atmospheric profile of horizontal wind.</figcaption></figure>

The damping coefficient $k$ is such that (Pruppacher & Klett 1978/1997):

$$k = 6\pi C_D \frac{\text{Re } \mu_a r_i}{24 m} \qquad (3)$$

Physically speaking, $k$ is the reciprocal of a relaxation time $\tau_D$ for the drag force (Paoli & Shariff 2016). It involves the Reynolds number $Re$, defined from the velocity of the crystal relative to ambient air by:

$$\text{Re} = \frac{\rho_a r_i}{\mu_a} \sqrt{(u - U_o)^2 + (w - W_o)^2} \qquad (4)$$

with $u = \dot{x}$ and $w = \dot{z}$.

Also needed are the drag coefficient $C_D$, such that (Pruppacher & Klett 1978, 1997; Wang and Ji 1997):

$$C_D = \frac{64}{\pi \text{Re}} (1 + 0.078 \text{Re}^{0.945}) \qquad (5)$$

the dynamic viscosity of air $\mu_a$, such that, for $-50 < T - T_r < 50$ ($T_r$ = 273.15 K) (Pruppacher & Klett 1978/1997):

$$\mu_a = 10^{-5} \{1.718 + 0.0049(T - T_r) - 1.2 \times 10^{-5}(T - T_r)^2\} \quad (6)$$

and the mass density $\rho_m$ of moist air (Picard et al. 2008; Nettesheim 2017):

$$\rho_m = \frac{\rho_a}{R_a T_a} \left( 1 - 0.378 \frac{\rho_{vsl}}{\rho} \right) \qquad (7)$$

When the crystal reaches its limit regime defined by vanishing vertical acceleration ($\ddot{z}_\infty = 0$), we obtain, from Eq. (1b), the limit velocity:

$$\dot{z}_\infty = W_a - \frac{g}{k} \left( 1 - \frac{\rho_a}{\rho_i} \right) = W_a + W_{f\infty} \qquad (8)$$

defining the terminal fall velocity $W_{f\infty}$ by:

$$W_{f\infty} = -\frac{g}{k} \left( 1 - \frac{\rho_a}{\rho_i} \right) \qquad (9)$$

Thermodynamic data of atmospheric temperature and pressure are detailed analytically in App. A.2 and plotted in Figure 3. The trajectory is driven by the

<figure><img src="image_6.png" /><figcaption>**Figure 2** Atmospheric profiles of vertical wind.</figcaption></figure>
ambient motion and the variations of the crystal's mass $m$, involved in the damping coefficient $k$, which evolves according to the process described in the following section.

## 2.2. MICROPHYSICS

### 2.2.1. Crystal mass

In this simplified approach, we consider an average ice crystal, once it has been generated by homogeneous nucleation in a highly ice-supersaturated environment ($S_i - 1 > 50\%$), and we follow its subsequent growth only by deposition of water vapour on its surface from an ambient medium slightly or moderately saturated with respect to ice (Köhler & Seifert 2015). Conversely, it loses matter by sublimation of ice into water vapour. In a model of cirrus including radiative transfer (Ramaswamy & Detwiler 1986), the growth of a crystal of 200 µm length is shown to reach an ambient ice supersaturation as small as 5% after 10 minutes. We shall show a posteriori that the oscillatory mode we put into evidence is likely to occur only if the ice supersaturation is of order of the radiative term ($3 < R < 15\%$).

In the capacitance model, crystal's mass growth is governed by the following equation with additional term $R$ and factors $f_v$ and $C$ accounting respectively for radiative transfer ($R$), ventilation ($f_v$), and capacitance ($C$), namely (Houghton 1950; Mason 1953; Hall & Pruppacher 1976; Rogers 1979; Wu et al. 2000; Lebo et al. 2008):

$$ \dot{m} = 4\pi C f_v \frac{S_i - 1 - R}{\frac{R_v T}{p_{vsi} D_v} + \frac{L_s}{K_a T} \left( \frac{L_s}{R_v T} - 1 \right)} \quad (10) $$

subject to the initial condition:

$$ m = m_0 = \rho_i V_0 \quad (11) $$

where the initial volume $V_0$ is estimated from Eq. (31a) below.

The saturation ratios $S_i$ and $S_v$ of water vapour, over liquid water and ice respectively, are defined by (Pruppacher & Klett 1997; Bohren & Albrecht 1998):

$$ \left\{ \begin{aligned} S_i &= \frac{p_v}{p_{vsi}} \\ S_i &= \frac{p_v}{p_{vsi}} \end{aligned} \right. \quad (12) $$

The relative humidity $H_i$ and $H_v$ of water vapour, with respect to (wrt) liquid water and ice respectively, can be expressed, as recommended by the WMO (World Meteorological Organization), by the following expressions (Pruppacher & Klett 1978/1997; Bohren & Albrecht 1998):

$$ \left\{ \begin{aligned} H_i &= \frac{p_v}{p_{vsi}} \frac{p - p_{vsi}}{p - p_v} \\ H_i &= \frac{p_v}{p_{vsi}} \frac{p - p_{vsi}}{p - p_v} \end{aligned} \right. \quad (13) $$

The atmospheric moisture profiles are detailed in App. A.3, and plotted in Figure 4. All these altitude profiles

<figure><img src="image_7.png" /><figcaption>Figure 3 Atmospheric profiles of temperature and pressure of dry air.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4 Atmospheric profiles of relative humidity and saturation ratios.</figcaption></figure>
of winds and humidity are consistent with radiosonde observations and input data used in realistic simulations with Eulerian numerical codes (Spichtinger 2014). The method of calculation of the partial pressure $p_v$ of water vapour from the given ambient moisture profile $H_i$ is detailed in App. B, and the expressions of the saturation pressures $p_{vsl}$ and $p_{vsi}$ of water vapour over liquid water and ice are summarized in App. C.

The saturation mixing ratio $q_s$, necessary for the calculation of the Kelvin-Helmholtz criterion (Sect. 6.2), is such that:

$$q_s = \frac{\epsilon p_{vsl}}{p - p_{vsl}} \qquad (14)$$

The diffusivity of water vapour $D_v$ in air and the heat conductivity of air $K_a$, involved in Eq. (10), are given by the following expressions (Rogers 1975)

$$\left\{ \begin{aligned} D_v &= 8.28 \times 10^{-3} \left( \frac{T}{p} \right) \frac{293}{T+120} \left( \frac{T}{273} \right)^{3/2} \\ K_a &= 2.42 \times 10^{-2} \frac{293}{T+120} \left( \frac{T}{273} \right)^{3/2} \end{aligned} \right. \qquad (15)$$

Other parameterizations for $D_v$ and $K_a$ exist (see Zeng 2008), that are slightly different, in particular $D_v$ being pressure-dependent, but they would give qualitatively the same results, though with different initial conditions. The ventilation coefficient $f_v$ is calculated using the following formulae, valid for column and plates in a laminar flow (Ji & Wang 1999; Liu et al. 2003a, b):

$$f_v = 1 + 0.039X + 0.1447X^2 \qquad (16)$$

where the variable $X$, such that:

$$X = \sqrt[3]{\frac{Sc}{\sqrt{Re}}} \qquad (17)$$

depends on Reynolds number $Re$ defined by Eq. (4), and Schmidt number $Sc$ defined by:

$$Sc = \frac{\mu_a}{\rho_a D_v} \qquad (18)$$

The factor $S_{i-1-R}$ playing a major role in Eq. (10) since it drives mass growth or decay according to its sign, we shall call it “driving factor” in the following.

### 2.2.2. Radiative transfer

The contribution of heat transfer in Eq. (10), appearing as a correction $R$ to the ice supersaturation term $S_{i-1}$, has been estimated using the formulation (Roach 1976; Wu et al. 2000; Lebo et al. 2008):

$$R = \frac{\mathcal{R}}{4\pi CK_a T} \left( \frac{L_s}{R_v T} - 1 \right) \qquad (19)$$

with the source term:

$$\mathcal{R} = AQ_{\text{abs}} \left( \frac{F_u - F_d}{2} - B_i \right) \qquad (20)$$

where $Q_{\text{abs}}$ denotes the absorption coefficient (Wu et al. 2000), $F_u$ and $F_d$ the upward and downward radiative flux densities respectively, A the crystal area, and $B_i$ the black body radiative flux of crystal surface at temperature $T_s$:

$$B_i = \epsilon \sigma T_s^4 \qquad (21)$$

Eq (19) has been used for fog and droplets (Roach 1976) and extensions to ice crystals in cirrus clouds (Stephens 1983; Wu et al. 2000; Gu & Liou 2000; Zeng 2008; Zeng et al. 2021) have been later widespread.

From the energy budget leading to Eq. (10), it is possible to calculate also the difference of temperature between the crystal's surface ($T_s$), assumed homogeneous, and the ambient air ($T_a$):

$$T_s - T_a = \frac{L_s}{4\pi C K_a} \dot{m} \qquad (22)$$

We shall check a posteriori that this difference is small. Then Wien's displacement law teaches us that the radiative transfer for a black body between - 20°C (253 K) and - 50°C (223 K) is in the range of wavelength 11-13 µm, i.e. infrared. Thus we shall consider infrared fluxes in Eqs. (20) and (23) (Stephens 1983). According to observational measurements (Inoue 1985) and theoretical background (Gierens 1994; Zeng 2008), we also assume that the emissivity $\epsilon$ of the crystal is equal to one.

Moreover we shall consider the transfer in the limit of short wavelength compared to the crystal dimensions, and therefore assume that $Q_{\text{abs}} \approx 1$ (Roach 1976; Liou 2002). It is convenient to express the radiative infrared ratio $\eta$ as (Zeng 2008; Zeng et al. 2022):

$$\eta = \frac{F_u - F_d}{2B_i} \qquad (23)$$

and then recast Eq. (20) as:

$$\mathcal{R} = A(\eta - 1)B_i \qquad (24)$$

Therefore, the crystal is subject to radiative cooling when $\eta < 1$, and to radiative warming when $\eta > 1$ (Zeng et al. 2021). If $F_u$ and $F_d$ are set to the following form at the top of the troposphere (Roach 1976):

$$\left\{ \begin{aligned} F_u &= -B_i \\ F_d &= -0.6F_u \end{aligned} \right. \qquad (25)$$

then the ratio $\eta$ becomes:

$$\eta = 0.9 \qquad (26)$$

This is confirmed by the modelling of the parameter $\eta$ we describe below.

Rather than modelling independently the radiative fluxes $F_u$ and $F_d$, because that would involve too many additional degrees of freedom, it is usual to introduce an atmospheric ratio $n_a$ as in Eq. (23), function of altitude. We shall derive an expression from models (Detwiler & Ramaswamy 1990; Zeng 2008) but consistent with
<figure><img src="image_9.png" /><figcaption>Figure 5 Atmospheric profile of atmospheric infrared ratio $η_a$.</figcaption></figure>

observations (Paltridge and Platt 1981; Philipona et al. 2012), typical of a mid-latitude daytime atmosphere, that we detailed in App. A.4 and plotted in Figure 5.

The crystal-related $η$ in Eq. (24) is related to the atmospheric $η_a$ by a relationship of the form (Zeng 2018):

$$η-1 = α(η_a - 1) \qquad (27)$$

where the coefficient $α$ depends on the shape of the crystal (Zeng 2018). We notice that when $η_a = 1$, then $η = 1$. and reciprocally, therefore radiative transfer cancels, whatever the shape could be. We estimated $α$ from our various test cases in Sect. 5, and we found it to be in the range 0.1–0.4.

For the crystal to grow ($\dot{m} > 0$), Eq. (10) shows that the following relation must be verified:

$$S_i \ge 1 + R \qquad (28)$$

When $η < 1$, i.e. $R < 0$, and therefore $R < 0$, we see from Eq.(28) that the saturation ratio $S_i$ need not be larger than one, and that an ice crystal can develop in a subsaturated environment (Hall and Pruppacher 1976). We shall experiment this situation with the oscillatory mode we found, in Case 3 relevant of uniform nighttime conditions ($η = 0.9$), which is detailed in Sect. 5.4.

Conversely, when $η > 1$, i.e. $R > 0$, and therefore $R > 0$, we see likewise from Eq.(28) that the saturation ratio $S_i$ must be larger than one, and that an ice crystal will develop in a supersaturated environment. We shall experiment this situation for the oscillatory mode, in Case 2 relevant of uniform daytime conditions ($η = 1.1$), which is detailed in Sect. 5.3.

When $η$ is altitude-dependent (Figure 5), we shall show that ice supersaturations as low as 3%, moreover consistent with observational statistics (Ovarlez et al. 2002; Spichtinger et al. 2004), can lead to ice crystal growth (Sect. 3.2.3).

In the present approach, the resulting variation of the ambient saturation ratio $S$ due to crystal growth is not considered, since the crystal trajectory does not go twice through the same region of space, even in the oscillatory regime (Sect. 3.2.3). Moreover Kelvin's effect due to surface tension has been neglected (Gu & Liou 2000).

Of course, involving these effects would be a useful refinement in further developments.

### 2.2.3. Crystal shape

According to theoretical models (Kobayashi, 1961; St-Pierre and Thiérault, 2015) and observations (Heymsfield 1975a; Schmitt & Heymsfield 2007), hollow column and bullet-rosette are thought to be the crystal shapes most likely developing in cirrus clouds, most of them (50 to 80%) having hollow ends (Schmitt & Heymsfield 2007).

A first possible approach to connect crystal mass and dimensions consists in assuming a mass-length relationship, which is the way followed in the reference work (Heymsfield 1973; Heymsfield 1975c) and many others (Rogers 1979; Jensen et al. 2018; Mascio et al. 2017). Because of the crystal's lacunarity, its mass density is then smaller than the actual density of solid ice (Cotton et al. 2013). A considerable amount of investigation has been devoted to the determination of the effective mass density of ice crystals encountered in cirrus clouds, both from in situ observational data (Heymsfield 1972; Heymsfield et al. 2002; Heymsfield et al. 2004; Schmitt & Heymsfield 2007; Cotton et al. 2013) and from laboratory experiments (Jensen & Harrington 2015). This is due to the fact that the cavities and hollow portions that are included in the effective volume make it larger than the actual volume.

An effective mass density $ρ_{ie}$ can be defined from the actual mass $m$ and the maximum particle dimension $D_m$ by the general expression (Cotton et al. 2013):

$$m = \frac{\pi}{6} D_m^3 \rho_{ie} \qquad (29)$$

A more specific definition can be given on the basis of the volume of the solid hexagonal column (Fukuta & Takahashi 1999):

$$m = 3\sqrt{3} a^2 c \rho_{ie} \qquad (30)$$

Dealing with hexagonal hollow columns, we shall use the latter definition.

An alternative approach consists in assuming a prevalent shape of the crystal, with a perfectly known geometry, so that its volume and surface area are easily calculated from its dimensions. Mass and volume are then connected to each other by means of Eq. (11) and (31) at each time, $ρ_i$ being the actual mass density of ice (Fridlind et al. 2016). As far as $ρ_i$ increases when temperature decreases (Pounder 1965), we shall adopt an average value: $ρ_i = 920$ kg m⁻³.

Leaving apart the bullet-rosette as too complex in a first step, we shall focus on the hexagonal hollow column. This kind of crystal has been thoroughly investigated (Chiruta & Wang 2005; Fridlind et al. 2016; Zeng 2018), and it has been shown that the crystal's volume and area can be simply expressed as (Chiruta & Wang 2005; Chen & Wang 2009):
$$ \begin{cases} V = \sqrt{3}(2c + c_B)a^2 \\ A = 6a \left\{ 2c + \sqrt{(c - c_B)^2 + \frac{3}{4}a^2} \right\} \end{cases} \quad (31a,b) $$

The shape of a hollow column with hexagonal cross section, drawn after Chen and Wang (2009), is shown in Figure 6. The core bulk of half-length $c_B$, is limited by two pyramidal hollow ends of depth $c-c_B$.

In their analytical model of crystal growth (Chen & Wang 2009), these authors assume that, during the growth, the parameter $c_B$ is subject to one the conditions:

$$ \begin{align} \text{i) } & c_B = c/2 \notag \\ \text{ii) } & c_B = c \notag \\ \text{iii) } & c_B = 0 \tag{32} \end{align} $$

We shall here extend this assumption and quantify the hollow portion of the crystal by a hollowness factor (Schmitt & Heymsfield 2007) denoted here $\psi$:

$$ \psi = 1 - \frac{c_B}{c} \quad (33) $$

We shall keep $c_B$ constant, equal to its initial value: $c_{B0} = (1-\psi_0)c_0$. For crystals longer than 100 µm, we shall set: $\psi = 0.8$ (Schmitt & Heymsfield 2007).

A further simplification consists in assuming that either the radius $a$ is constant (Harimaya 1968, Chen & Wang 2009) or $a$ and $c$ are connected by some width-length relationship (Chen & Wang 2009; Pruppacher & Klett 1978). Adopting the first assumption, we found that the aspect ratio $\phi$:

$$ \phi = \frac{c}{a} \quad (34) $$

reaches large values ($\phi > 30$) as soon as the saturation exceeds 130%, relevant of thin needles, inconsistent with observations. Therefore, we adopted the second assumption, and assume thereby that, for hollow columns, the aspect ratio $\phi$ is constant during the motion. The value $\phi = 2$ was found for columns with $c < 150$ µm from in situ measurements (Heymsfield 1972; Um et al. 2015). Therefore, using Eq. (31a,b) and the definition of $\phi$ we are led to a third-degree equation for $a$:

$$ 2\phi a^3 + c_B a^2 - \frac{V}{\sqrt{3}} = 0 \quad (35) $$

<figure><img src="image_10.png" /><figcaption>Figure 6 The shape of the mean crystal (meridian and axial cross sections) (after Chen and Wang, 2009).</figcaption></figure>

As the crystal evolves, its mass $m$ changes, and its volume $V$ is calculated by the relation:

$$ V = \frac{m}{\rho_i} \quad (36) $$

The crystal's shape being significantly different from a sphere, it is useful and usual to define the mass equivalent diameter $D_i$ and radius $r_i = D_i/2$ of the crystal as those of the sphere with the same volume $V$ as the crystal, namely:

$$ D_i = \sqrt[3]{\frac{6V}{\pi}} \quad (37) $$

or, using Eq. (36):

$$ D_i = \sqrt[3]{\frac{6m}{\pi \rho_i}} \quad (38) $$

Alternately, by setting a maximum particle size $D_m = \max(2a, 2c)$, as the larger of the two dimensions $2a$ and $2c$ (that is $2c$ if $\phi > 1$ or $2a$ if $\phi < 1$) following the general definition (29) and using Eq. (36), we could also define an effective mass density $\rho_{ie}$ such that:

$$ \rho_{ie} = \frac{6V}{\pi D_m^3} \rho_i \quad (39) $$

Following the more specific definition (30), we obtain the expression (Fukuta & Takahashi 1999):

$$ \rho_{ie} = \frac{V}{3\sqrt{3} a^2 c} \rho_i \quad (40) $$

which will yield values smaller than the actual density $\rho_i$. The first definition (Eq. (39)) overestimates the actual volume, and therefore, it underestimates $\rho_{ie}$, yielding values as low as 200 kg m⁻³ found in observations of complex particles or aggregates (Heymsfield et al. 2004). The second definition (Eq. (40)) yields values in better agreement with observations of simple shapes (Miller & Young 1979; Schmitt & Heymsfield 2007), and we shall follow its variation during the motion of the crystal.

For a hollow column, the capacitance $C$ appearing in Eq. (10) can be written (Chiruta & Wang 2005):

$$ C = 0.751a + 0.491c = (0.751 + 0.491\phi)a \quad (41) $$

using the aspect ratio $\phi$ defined above (Eq. (34)). Moreover, we shall assume that the ventilation coefficient defined by (16) accounts for the different orientations of the crystal during its motion. A detailed account of the flow characteristics around columns and plates with estimation of torques based on Navier-Stokes equations teaches us that the Reynolds number ranges from 2 to 70 (Hashino et al. 2014). Lidar observations show that plane crystals have mainly horizontal orientation during free fall, with a maximum tilt angle of 0.3 degrees (Thomas et al. 1990). The particle radius $r_i$ is set equal to $D_i/2$ in Eqs (3) and (19), and to $D_i$ in Eq. (4) as a possible alternative definition of Reynolds number (Pruppacher and Klett
1997). While it is easy to define an equivalent diameter for liquid spherical particles, this is quite difficult for ice crystals (McFarquhar & Heymsfield 1998). Thus, large discrepancies occur for radii larger than 50 µm, so that we shall keep the above convention. Nevertheless, we shall also display the half-length c in the result of calculations.

# 3. SOLUTIONS

## 3.1. NUMERICAL METHOD

The differential system composed of Eqs. (1a, b) and (10) is a Liénard system (Perko 1990), that we discretized with finite differences, and using an explicit Euler scheme such that, at any time $t_n = n \Delta t$, we write the velocity components of the average crystal's center of mass as:

$$ \begin{cases} u_{n+1} = u_n - k_n(u_n - U_a(z_n))\Delta t \\ w_{n+1} = w_n - \left\{ k_n(w_n - W_a(z_n)) + g\left(1 - \frac{\rho_a(z_n)}{\rho_i}\right) \right\}\Delta t \end{cases} \quad (42) $$

its rectangular coordinates as:

$$ \begin{cases} x_{n+1} = x_n + u_n \Delta t \\ z_{n+1} = z_n + w_n \Delta t \end{cases} \quad (43) $$

and its mass as:

$$ m_{n+1} = m_n + 4\pi C_n f_v \frac{S_n - 1 - R_n}{\frac{R_v T_n}{\rho_{vsi} D_v} + \frac{L_s}{K_a T_n} \left( \frac{L_s}{R_v T_n} - 1 \right)} \Delta t \quad (44) $$

where $T_n = T_a(z_n)$ and $S_n = S_v(z_n)$ are calculated from the data detailed in App. A. Eqs. (42), (43), and (44) constitute a set of strongly coupled nonlinear equations.

Because the radiative correction $R$ and the capacitance $C$ (Eqs. (19) and (41) respectively) are dependent of the axial size $c$, we also need to recalculate $V_{n+1}$ from $m_{n+1}$ using Eq. (36):

$$ V_{n+1} = \frac{m_{n+1}}{\rho_i} \quad (45) $$

and then $c_{n+1}$ from Eq. (31a,b) if $a$ is constant ($a = a_0$):

$$ c_{n+1} = \frac{1}{2} \left( \frac{V_{n+1}}{a_0^2 \sqrt{3}} - c_B \right) \quad (46) $$

or solving Eq. (35) for $a_{n+1}$ if $\phi$ is constant ($\phi = \phi_0$):

$$ 2\phi_0 a_{n+1}^3 + c_B a_{n+1}^2 - \frac{V_{n+1}}{\sqrt{3}} = 0 \quad (47) $$

We neglect the variation of the gravity field $g$ with altitude, between approximately 7 and 10 km height, the Coriolis force, and the influence of the column orientation on the ventilation coefficient $f_v$.

Below is given a detailed account of the two regimes or modes we found with standard initial conditions. Two special sections will be devoted respectively to the sensitivity of various parameters (lifetime, period,

amplitude...) to the crystal's initial size $a_0$ (Sec. 4) and to the impact of radiative transfer (Sect. 5).

## 3.2. RESULTS

### 3.2.1. Data and displays

In this section the simultaneous integration of the discretized differential system composed of Eqs. (42) to (45) and (47) subject to the initial conditions of Table 2 (Case 0) is carried on with the aspect ratio $\phi$ and the core half-length $c_B$ kept constant at each time step: $\phi = \phi_0$ and $c_B = c_{B0}$. We obtain a dual mode composed of a steady well-developed, non-oscillatory, cirrus uncinus (hooked cap + trail) with a high-base updraft, and an oscillatory counterpart associated with a low-base updraft.

For each of the two selected modes are displayed crystal trajectory $z(x)$, hodograph or phase portrait $w(u)$, abscissa and altitude, free fall speed $W_p$, dimensions ($a, c, c_B$), mass equivalent diameter $D_p$, dynamic viscosity of air $\mu_a$, hollowness factor $\psi$, effective mass density of crystal $\rho_{ie}$, drag coefficient $c_D$, Reynolds number $Re$, damping coefficient $k$, temperature difference $T_s - T_a$ between ambient air and crystal surface.

Two kinds of dual solutions are obtained, depending on the altitude $z_{w0}$ of the updraft base, other things being unchanged. Nevertheless, these two solutions are obtained after a minute fitting of the crystal initial half-width $a_0$, the wind velocity amplitude and the updraft intensity, altitude and thickness, in narrow ranges. This behaviour is mostly conditioned by the temperature and humidity profiles we chose (App. A), which are characteristic of a mid-latitude cirrus background. We also show phase relationships in the various time profiles of the damped oscillator solution.

### 3.2.2. Non-oscillatory mode

When the base of the updraft lies at an altitude of 10 km (Figure 2a), the trajectory looks like that obtained from the shape of an average cirrus uncinus, with a hooked head and a trail or virga (Figure 7a), consistently with a conceptual view (Heymsfield 1975b) and the description given by the International Cloud Atlas (WMO 1975). The head is similar to the upper part of a cirrus uncinus as observed during the FIRE IFO II experiment, and reported later (Sassen & Krueger 1993).

The trail or virga, referring to the curved shape with a dot at the origin (similar to the French word "virgule" meaning comma), is seen in the lower part. Its physical nature has been subject to a controversy (Fraser & Bohren 1992; Sassen & Krueger 1993). Actually, in our simplified approach, the crystal sublimation takes place in the virga. Along the trajectory, the maximum height (where $dz/dt = 0$) is reached after approximately 1.5 hours (Figure 8b), and the extremum (where $dx/dt = 0$) after 3.3 hours (Figure 8a).

The hodograph (Figure 7b) has no periodic branch and tends toward free fall with entrainment by horizontal
wind. Our vertical velocities are smaller than reported values (Heymsfield 1975a), consistently with the fact that we assume a weaker updraft (0.6 m s⁻¹ instead of 1 m s⁻¹), but they are in good agreement with FIRE observations (Gultepe et al. 1990).

The scale lengths are compatible with our calculated shape and the integration runs over 5 hours, a duration that represents about 1/6 of a typical cirrus uncinus lifetime (Luo & Rossow 2004). With a constant time step Δt = 0.02 s, 900000 time steps are necessary.

The time profiles of abscissa and altitude (Figure 8) shows that the crystal spends much time in the head where it grows for about three hours, and then decays rapidly, consistently with its mass (Figure 11a). The profile of abscissa (Figure 8a) shows a minimum at 3.3

hours, reached at the left-most point of the trajectory (Figure 7a), where the crystal is pushed by the horizontal wind (Figure 1) directed to the left ($U_0 < 0$). Afterwards, the crystal falls down in the decay phase following the trail.

The profile of free fall speed $W_f$ (Figure 9a) is a good test because it is similar to that obtained with the analytic model developed in App. D.1. The estimate of the theoretical time duration obtained there is consistent with the 5 hours reached here. We find typical values of $W_f$ in the range 25–65 cm s⁻¹ that are in good agreement with observations of downdrafts in cirrus clouds (Heymsfield 1975b), as well as laboratory and field measurements (Heymsfield & Westbrook 2010) and measurements in synoptic cirrus (Mishra et al. 2014).

<figure><img src="image_11.png" /><figcaption>**Figure 7** Profiles of the parcel trajectory and hodograph ($z_0 = 10$ km; $z_{w0} = 10$ km).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8** Time profiles of abscissa and altitude ($z_{w0} = 10$ km).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9** Time profiles of free-fall speed $W_f$ and driving factor $S_j-1-R$ ($z_{w0} = 10$ km).</figcaption></figure>
The variation of ice supersaturation $S_i-1$ (Figure 10a) is roughly stationary, with a local minimum coincident with the height maximum ($t \approx 1.5$ h) and positive between 2 and 3% during the growth phase until 3.4 hours. Then it decreases, becomes negative ($t \approx 3.75$ h) in the decay phase, goes down to -20%, consistent with observational data fitted with a microphysics model (Khvorostyanov & Curry 2008).

The radiative correction term $R$ (Figure 10b) is positive and slightly increasing, amounting to 1%, in the growth phase and becomes negative ($t \approx 3.75$ h), synchronously with $S_i-1$, then decreases to nearly -1% in the decay phase. Therefore this contribution is smaller than the supersaturation ($S_i - 1 \approx 3\%$), and as shown by Eq. (10), its effect counterbalances supersaturation in the heated growth phase ($R > 0$) and contributes in the cooled decay phase ($R < 0$) to lengthen the crystal's life time. Consequently, the driving factor $S_i-1-R$ (Figure 9b) is slightly smaller than the supersaturation $S_i-1$.

The crystal mass $m$ (Figure 11a) is the variable coupled with position $(x, z)$ and its profile shows an inflexion point at the altitude maximum ($t \approx 1.5$ h), then

reaches a maximum ($m \approx 1.2$ µg) at $t \approx 3.75$ hours and decays afterwards down to approximately 0.3 µg. The dimensions $a$ and $c$ (Figure 11b) reach maxima at the same time ($t \approx 3.75$ h), and then decrease to 35 µm and 70 µm respectively, satisfying the initial aspect ratio ($\phi = 2$). According to our assumption, the core half-length $c_B$ remains constant during the motion. The total length $2c$ is of same order as modelled by other researchers (Jensen et al. 2018).

The mass equivalent diameter $D_i$ (Figure 12a) has a maximum ($D_i \approx 137$ µm) at the same time ($t \approx 3.75$ hours), corresponding to a distance of 400 m (Figure 8a), quite consistently with a simple former model (Harimaya 1968). The value $D_i \approx 85$ µm is reached in the virga after 5 hours, in good agreement with observations at similar temperatures (Baum et al. 2000; Kuhn & Heymsfield 2016). Dynamic viscosity (Figure 12b) is, together with the diameter $D_i$ (or radius $r_i = D_i/2$), a key parameter that scales the drag force (Eq. (3)) and Reynolds number (Eq. (4)). It is nearly constant in the head during 3.6 hours, and then it increases in the tail.

<figure><img src="image_14.png" /><figcaption>**Figure 10** Profiles of ice supersaturation $S_i-1$ and radiative correction $R$ ($z_{wO} = 10$ km).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11** Profiles of crystal's mass and dimensions ($z_{wO} = 10$ km).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12** Profiles of mass equivalent diameter and dynamic viscosity ($z_{w0}$ = 10 km).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13** Profiles of hollowness factor and effective mass density of crystal ($z_{w0}$ = 10 km).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14** Profiles of drag coefficient and Reynolds number ($z_{w0}$ = 10 km).</figcaption></figure>

The hollowness factor $\psi$ (Figure 13a) increases during the growth phase, reaching a maximum close to 0.82 at 3.75 h, approximately like mass, and then decreases down to nearly 0.70. This is consistent with the opposite variation of effective mass density (Figure 13b). Actually, $\rho_{ie}$, as calculated from Eq. (40), reaches a minimum when $D_i$ is maximum ($t \approx 3.75$ h) and afterwards increases to large values. The range 670–705 kg m⁻³ covered during crystal motion is slightly smaller than values amounting

to 800 kg m⁻³ published for bullet shapes (Schmitt and Heymsfield 2007).

The drag coefficient $C_D$ (Figure 14a) reaches a minimum at 4.2 h, later than the mass maximum (Figure 11a), but expectingly in phase with the maximum of Reynolds number $Re$ at the same time (Figure 14b). The ranges covered by $C_D$ (9–25) and $Re$ (1–3) are perfectly confirmed by numerical and laboratory experiments (Wang & Ji 1997).
The damping coefficient *k* (Figure 15a) is minimum at 3.75 hours and it increases by 50% in the final phase. It amounts to 15–20 s⁻¹ and therefore the quantity *k* Δ*t* is maintained in the range 0.3–0.6 thus ensuring stability of the algorithm (Eqs. (42)). The temperature difference *T*<sub>s</sub>−*T*<sub>a</sub> (Figure 15b) is very small and positive in the heated growth phase (*T*<sub>s</sub>−*T*<sub>a</sub> < 0.01 K), then it becomes negative and increases in magnitude (*T*<sub>s</sub>−*T*<sub>a</sub> < 0) in the cooled decay phase, though remaining smaller than 0.1 K, consistently with published results (St-Pierre & Thiérault 2015).

We can notice that the angular point appearing on
certain profiles is due to the angular point present on the
theoretical piecewise temperature profile (Figure 3a) with
coordinates: z = 8 km, T_a = -50°C. It affects quantities
depending directly of the temperature (dynamic viscosity,
temperature difference...), and consequently the velocity
components (hodograph). In contrast, the trajectory is
everywhere differentiable because of the smoothing of
coordinates (x,z) due to the time integration of velocity
components (u,w).

3.2.3. Oscillatory mode

When the base of the updraft is lowered down to an
altitude of 9 km (Figure 2b), other parameters being

unchanged, the trajectory is subject to a damped
oscillation (Figure 16a). After the hooked head is formed
with approximately a hundred meter vertical extension,
like in the previous steady, non-oscillatory mode (Sect.
3.3.2), the crystal parcel follows a periodic motion,
whose wavelength *λ* and period *P* are approximately 70
km and 10.5 hours respectively. The amplitude of height
oscillations is about 400 m at the beginning, and it
decreases afterwards (Figure 17b).

Owing to the mass loss of the crystal, which induces
variations of the damping coefficient *k* (Figure 24a)
or equivalently of the viscous drag (Figure 23a), the
updraft acts in the bottom part of the trajectory at an
altitude of 9.2 km, that is 200 m above the base (z<sub>w0</sub> =
9 km), causing the parcel to lift again, like in cumulus
models (Grabowski 1993). The crystal bounces on
the updraft base like a ball on a rigid floor, except
that the crystal equilibrium altitude is not the lower
attained level (z-z<sub>0</sub> ≈ -0.75 km in Figure 16a), but an
intermediate height, unlike the ball that stabilizes at
the floor level.

Let us examine the limit focus of the hodograph
(Figure 16b) which is composed of a straight line followed
by a hook corresponding to the head, and afterwards of

<figure><img src="image_19.png" /><figcaption>**Figure 15** Time profiles of damping coefficient and temperature difference ice-air (z<sub>w0</sub> = 10 km).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16** Profiles of the parcel trajectory and hodograph (z₀ = 10 km; zₘ₀ = 9 km).</figcaption></figure>
a spiral trajectory winding one turn more at each period. The phase velocity $C_p$ corresponding to the period $P$ and wavelength $\Lambda$ is defined by:

$$C_p = \frac{\Lambda}{P} \qquad (48)$$

that is, numerically: $C_p = 70 \times 10^3/(10.5 \times 3600) \approx 1.85 \text{ m s}^{-1}$. Now, let us write the condition for the fixed point:

$$\begin{cases} \dot{u} = 0 \\ \dot{w} = 0 \end{cases} \qquad (49)$$

and using Eqs. (1a, b) and neglecting the buoyancy force, we obtain:

$$\begin{cases} u = U_a(z_\infty) \\ w = W_a(z_\infty) - \frac{g}{k_\infty} \end{cases} \qquad (50)$$

The fixed point satisfies the additional condition:

$$w = 0 \qquad (51)$$

and this finally leads to:

$$\begin{cases} u_\infty = U_a(z_\infty) \\ W_a(z_\infty) = \frac{g}{k_\infty} \end{cases} \qquad (52)$$

Extrapolating the hodograph spiral (Figure 16b) shows that the horizontal velocity tends towards $u_\infty = 1.85 \text{ m s}^{-1}$ approximately. We can then estimate the wind speed at the limit height $z_\infty \approx 9.63 \text{ km}$ (Figure 18a) using Eqs. (68) in layer III: $U_a(z_\infty) = -30 - 5(9.63-16) = 1.85 \text{ m s}^{-1}$. This confirms the value given by the cycle focus ($U_a(z_\infty) = u_\infty$). Moreover, this limit velocity is very close to the phase speed of oscillations Eq. (48), which is approximately 1/10 that of gravity waves.

It is noteworthy that the maxima along the trajectory (Figure 16a) are sharp at the beginning and tend to become blunter and blunter, while the minima have nice round shapes all along the motion, like the extrema of other parameters displayed in the following. This is

due to the fact that the abscissa (Figure 17a) is not a linear function of time, but has a periodic component of decreasing amplitude superimposed on a secular variation (of order $C_p t$; see below).

As can be clearly seen on the time profile of altitude (Figure 17b), the limit height ($z_\infty \approx 9.63 \text{ km}$) is slightly larger than the updraft base (9 km). It is in good agreement with the theoretical model of damped harmonic oscillator built in App. D.2 on the same period and a matched damping ratio (Figure 41a). As it was noticed, the abscissa $x$ has a wavy profile (Figure 17a), and being coupled to altitude $z$ through Eq. (1), it is in quadrature with $z$, since its quasi-horizontal inflexion points are coincident with maxima of $z$, and its oblique inflexions are in phase with minima of $z$. The slope of the straight line passing through these inflexion points is approximately equal to the phase speed $C_p$ and the limit velocity $u_\infty$.

The time profile of fall speed $W_f$ (Figure 18a) oscillates with the same characteristics as the trajectory, and it tends towards the terminal value $W_f \approx -0.6 \text{ m s}^{-1}$, which compensates the updraft ($W_a = 0.6 \text{ m s}^{-1}$), and thus sustains the parcel oscillation. We can notice that, up to 10 hours, its beginning portion quite consistently looks like the theoretical profile relative to the steady mode shown in App. D.1 (Figure 38).

The ice supersaturation $S_i - 1$ (Figure 19a) shows damped oscillations with the same period and it does not tend towards zero, but towards 0.25% with increasing time. Its amplitude is initially small (3%) but larger than the radiative correction (0.75%). Moreover, $S_i - 1$ and $R$ are in phase with altitude (Figure 17b) and in quadrature with mass (Figure 20a). At infinity they both tend towards 0.25%, and consistently the driving factor $S_i - 1 - R$ (Figure 18b) tends towards zero. We shall examine further below (Sect. 5) the relation between $S_i - 1$ and $R$ for the oscillation to develop.

<figure><img src="image_21.png" /><figcaption>Figure 17 Time profiles of abscissa and altitude ($z_{w0} = 9 \text{ km}$).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 18** Time profiles of free-fall speed $W_f$ and driving factor $S_i-1-R$ ($z_{w0}$ = 9 km).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19** Time profiles of ice supersaturation $S_i-1$ and radiative correction $R$ ($z_{w0}$ = 9 km).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20** Time profiles of mass $m$ and dimensions $a, c, c_B$ ($z_{w0}$ = 9 km).</figcaption></figure>

Crystal mass (Figure 20a) is the variable coupled with position (x, z), and careful examination of its profile compared to that of attitude z (Figure 17b) and supersaturation (Figure 19a) shows that they are in phase quadrature. The dimensions *a* and *c* (Figure 20b) naturally show damped oscillations, in phase with mass. The half-length *c* tends towards a limit $c_∞$ ≈ 108 µm at infinity, which is larger than the initial one ($c_0$ = 82.95 µm).

The mass equivalent diameter (Figure 21a) oscillates in phase with the fall speed and half-length *c*, and it tends towards a limit value, approximately 131 µm. Because of its definition, Eq. (38), it is naturally in phase with mass (Figure 20a). Dynamic viscosity (Figure 21b) is an important parameter that scales the viscous drag. We understand from its being in phase opposition with altitude (Figure 17b), that when the crystal reaches the top of a spatial
period, gravity predominates over viscous drag, and the crystal falls down. When it reaches the bottom of a period, the situation is reversed: viscous drag predominates and the crystal is pushed up.

The core half-size $c_B$ being kept constant during the motion (Figure 20b), the hollowness factor $\psi$ (Figure 22a) oscillates with a small amplitude (≈2%) in phase with $c$, and it tends towards a limit value ($\psi_\infty$ ≈ 0.81) slightly larger than its initial value ($\psi_0$ = 0.80). The crystal's effective mass density as calculated from Eq. (40), shows oscillations (Figure 22b) in phase with $W_p$ and phase opposition with $\psi$ and $D_t$. It tends towards the asymptotic approximate value 671.5 kg m⁻³ approximately.

The drag coefficient $C_D$ (Figure 23a) oscillates also at the same frequency and wavelength, in phase opposition with Reynolds number $Re$ (Figure 23b). These two parameters are respectively comprised in the ranges 8.5–11.5 and 2–3 and tend towards 10 and 2.45 at infinity. Like in the non-oscillatory mode (Sect. 3.2.2), the ranges covered by $C_D$ and $Re$ are in good agreement with numerical and laboratory experiments (Wang & Ji 1997). Nevertheless, probably because $C_D$ and $Re$ involve the mass variation through the mass equivalent radius (Eq. (4)), they are phase quadrature with altitude.

The damping coefficient $k$ (Figure 24a) oscillates roughly between 15 and 18, and tends towards the value $k_\infty$ ≈ 16.5 s⁻¹, which is connected to $W_{f\infty}$ (≈ 0.6 m s⁻¹)

<figure><img src="image_25.png" /><figcaption>**Figure 21** Time profiles of mass equivalent diameter and dynamic viscosity ($z_{w0}$ = 9 km).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 22** Time profiles of hollowness factor and effective mass density of crystal ($z_{w0}$ = 9 km).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 23** Time profiles of drag coefficient $C_D$ and Reynolds number $Re$ ($z_{w0}$ = 9 km).</figcaption></figure>
through Eq. (9). The oscillation of the temperature
difference between the crystal's surface and the ambient
air (Figure 24b) tends towards zero as time tends to
infinity, i.e. the surface temperature tends towards the
ambient temperature $T_a$ corresponding to the limit height
$z_\infty \approx 9.63$ km.

In the next sections we shall focus on the impact
of crystal's initial mass or size (Sect. 4) and radiative
transfer (Sect. 5) on the trajectory, as primary causes of
spreading of the cirrus structure.

4. SENSITIVITY TO CRYSTAL SIZE

4.1. NON-OSCILLATORY MODE

Below a critical size ($a_0$ < 53.4 µm), the crystal parcel is
lifted, brought along towards negative x, and does not
fall down below its generating level. On the opposite,
above a size threshold ($a_0$ > 53.4 µm), the parcel is not
lifted and directly falls down below its generating level,
yielding only the virga without hooked head.

Looking in detail the situation of the parcel at
horizontal distances of 10 and 20 km gives some insight
into the crystal's history. The falling heights $z_{10} - z_0$ and

$z_{20} - z_0$ (Figure 25a) are maximum when $a_0 \approx 53.4$ µm
and then decrease quite linearly as function of $a_0$. The
spans of $z_{10} - z_0$ and $z_{20} - z_0$ in the range 55–100 µm reach
600 m and 2 km respectively, and they give an idea of the
trail thickness that is consistent with published results
(Sunilkumar & Parameswaran 2005) for moderately cold
cirrus (–75°C < $T-T_*$ < –50°C). Colder cirrus (–85°C < $T-T_*$ <
–75°C) at higher altitudes are thinner.

Interestingly it appears that the relative maximum
height $z_M - z_0$ reached at the top of the hook (Figure 25b)
plotted as a function of $a_0$ comes to zero at $a_0 \approx 53.3$ µm,
indicating that a hook does not form for larger sizes. In
the opposite, when the crystal is smaller, a head develops
up to more than 100 m above the initial position, and as
expected, the lighter it is, the higher it is lifted.

Simulations of cirrus clouds (Dobbie & Jonas 2001) in
presence of radiative cooling (RC) show that IWP decay is
achieved in approximately 5 hours without RC whatever
the cloud thickness, and this lifetime is increased for thin
clouds up to more 7 hours with RC.

The transit times $t_{10}$, $t_{20}$ at the abscissae x = 10 and 20 km plotted as a function of $a_0$ (Figure 26a) reach several hours for small crystals ($a_0 < 52$ µm) and shorten to less than one hour for large crystals ($a_0 > 55-65$ µm).

<figure><img src="image_28.png" /><figcaption>Figure 24 Time profiles of damping coefficient *k* and temperature difference ice-ambient air ($z_{w0}$ = 9 km).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 25** Profiles of falling height $z_{20}-z_0$ and maximum height in the head $z_M-z_0$ as functions of crystal half-width $a_0$ ($z_{w0}$ = 10 km).</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 26 Profiles of transit times $t_{10}$, $t_{20}$ and half-lengths $c_{10}/c_0$ and $c_{20}/c_0$ as functions of crystal half-width $a_0$ ($z_{w0}$ = 10 km).</figcaption></figure>

The half-lengths $c_{10}$ and $c_{20}$ at the same times, normalized to the initial value $c_0$ and plotted as functions of $a_0$ (Figure 26b) show an interesting behaviour. The profiles have a sharp peak at approximately $a_0 \approx 51.5$ µm, and then decrease with inflexion points corresponding roughly to the maxima of the falling height profiles (Figure 25a). Then, $c_{10}/c_0$ is nearly constant, while $c_{20}/c_0$ decreases nonlinearly down to 80%. Thus, the crystal size changes little up to 10 km and undergoes its major decay at 20 km by more than 10% according to its initial size in the range 55–100 µm. Unlike the falling height, the crystal decay is a strongly nonlinear function of size, reflecting the complexity of the microphysical processes involved in the model despite its many simplifying assumptions.

<figure><img src="image_31.png" /><figcaption>Figure 27 Profiles of pseudo-periods between minima ($P_{min}$) and maxima ($P_{max}$) estimated with Eq. (53), as functions of crystal half-width $a_0$ ($z_{w0}$ = 9 km).</figcaption></figure>

## 4.2. OSCILLATORY MODE

We shall examine the first three (exceptionally four) extrema in the trail. Like in Sect. 4.1, below a critical size ($a_0 < 50.9$ µm), the crystal is lifted, brought along towards negative $x$, and does not fall down below its generating level. On the opposite, above a threshold ($a_0 \approx 56$ µm), a transition occurs at the first maximum (in the $x$-range 25–40 km) where the horizontal velocity $u$ vanishes because of the crystal becoming too heavy and the damping passing through an inflection point. Then horizontal velocity $u$ changes sign and a loop forms, grows until the crystal is trapped in the first period, and the motion loses its spatial periodicity.

From the transit times at the first three minima ($t_{min1}$, $t_{min2}$, $t_{min3}$) and maxima ($t_{max1}$, $t_{max2}$, $t_{max3}$) we can calculate two estimates $P_{min}$ and $P_{max}$ of period by the formulae:

$$ 
\begin{aligned}
P_{\min} & \approx \frac{t_{\min 3} - t_{\min 1}}{2} \\
P_{\max} & \approx \frac{t_{\max 3} - t_{\max 1}}{2}
\end{aligned} 
\quad (53) 
$$

and plot the results as functions of $a_0$ (Figure 27). Ne notice that the average involving the middle minimum or maximum would yield the same result, since: ($t_{min3} - t_{min2}$) + ($t_{min2} - t_{min1}$) = $t_{min3} - t_{min1}$, and: ($t_{max3} - t_{max2}$) + ($t_{max2} - t_{max1}$) = $t_{max3} - t_{max1}$.

$t_{max1}$) = $t_{max3}$ - $t_{max1}$. It is noteworthy that the period is nearly independent of crystal size $a_0$ in the range 51–56 µm and the two curves for $P_{min}$ and $P_{max}$ are nearly coincident within $10^{-3}$. Estimates from these two plots are consistent with the value found in Sect. 3.2.3. The steepening at $a_0 \approx 56.2$ µm reflects the transition explained above.

Likewise, from the abscissae at the first three minima ($x_{min1}$, $x_{min2}$, $x_{min3}$) and maxima ($x_{max1}$, $x_{max2}$, $x_{max3}$) we can calculate two estimates $\Lambda_{min}$ and $\Lambda_{max}$ of wavelength by the formulae:

$$ 
\begin{aligned}
\Lambda_{\min} & \approx \frac{x_{\min 3} - x_{\min 1}}{2} \\
\Lambda_{\max} & \approx \frac{x_{\max 3} - x_{\max 1}}{2}
\end{aligned} 
\quad (54) 
$$

and plot the results as functions of $a_0$ (Figure 28). The two curves have a bump in the range 51–56 µm with a maximum in the middle at $a_0 = 53.5$ µm and they are separated by a gap of about 143 m, and a more exact value would be obtained from the transit through the asymptotic height $z_{\infty}$ calculated a posteriori from $z_{min}$ and $z_{max}$ (see below). These results are consistent with the value found in Sect. 3.2.3.
From periods and wavelengths, we can calculate corresponding phase speeds by Eq. (48), such that:

$$C_{pmin} = \frac{\Lambda_{min}}{P_{min}} \qquad (55)$$

$$C_{pmax} = \frac{\Lambda_{max}}{P_{max}}$$

and plot the results as functions of $a_0$ (Figure 29). Consistently with the variations of periods and wavelengths, phase speeds show a slow variation in the range 51–56 µm and they are separated by a small gap, $C_{min}$ and $C_{max}$ having both a maximum at the critical size $a_0 \approx 53.5$ µm.

In contrast, plots of the altitudes at the minima (Figure 30a) and maxima (Figure 30b) reveal that their spacing depends on the crystal size and has a minimum at $a_0 \approx 53.3$ µm, corresponding to the size for which the formation of a hooked head is cancelled (Figure 25b).

Then, from the first three minima and maxima we can obtain a rough estimate of the limit height $z_\infty$ by the formula:

$$z_\infty \approx \frac{z_{min1} + z_{max1} + z_{min2} + z_{max2} + z_{min3} + z_{max3}}{6} \qquad (56)$$

then we derive the limit velocity $U_a(z_\infty)$ using Eqs. (68) in layer III: $U_a(z_\infty) = -30 -5(z_\infty-16)$, and finally plot the results as functions of $a_0$ (Figure 31). They show nice parabolic profiles within 1% in the range 51–56 µm of size, and again an extremum at $a_0 \approx 53.3$ µm.

<figure><img src="image_32.png" /><figcaption>Figure 28 Profiles of wavelengths $\Lambda_{min}$ and $\Lambda_{max}$ derived from distances of minima and maxima by Eq. (54) as functions of crystal half-width $a_0$ ($z_{wo} = 9$ km).</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 29 Profiles of phase speeds $C_{pmin}$ and $C_{pmax}$ derived by Eq. (55) as functions of crystal half-width $a_0$ ($z_{wo} = 9$ km).</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 30 Profiles of relative altitudes at minima and maxima in the trail as functions of crystal half-width $a_0$ ($z_{wo} = 9$ km).</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 31 Profiles of relative limit height $z_\infty - z_0$ and velocity $U_a(z_\infty)$ in the trail as functions of crystal size $a_0$ ($z_{wo} = 9$ km).</figcaption></figure>
The damping ratio $\zeta$ appearing in the theoretical model of damped harmonic oscillator (App. D.2) can be estimated from our simulation data using the following procedure. The minimum and maximum altitudes can be written, within the approximation $\Omega \approx \Omega_0$:

$$ 
\begin{aligned}
z_{\min} - z_{\infty} &\approx \exp(-\zeta \Omega t_{\min}) \sin(\Omega t_{\min} + \Phi_{\min}) \\
z_{\max} - z_{\infty} &\approx \exp(-\zeta \Omega t_{\max}) \sin(\Omega t_{\max} + \Phi_{\max})
\end{aligned}
\quad (57) $$

Therefore, extracting $\zeta$ we obtain, for instance between extrema 1 and 2:

$$ 
\begin{aligned}
\zeta_{\min 12} &\approx \frac{\ln \frac{z_{\min 1} - z_{\infty}}{z_{\min 2} - z_{\infty}} \sin(\Omega t_{\min 2} + \Phi_{\min})}{\Omega(t_{\min 2} - t_{\min 1})} \\
\zeta_{\max 12} &\approx \frac{\ln \frac{z_{\max 1} - z_{\infty}}{z_{\max 2} - z_{\infty}} \sin(\Omega t_{\max 2} + \Phi_{\max})}{\Omega(t_{\max 2} - t_{\max 1})}
\end{aligned}
\quad (58) $$

With the above data providing an average angular frequency $\Omega_0 = 2\pi/(10.24 \times 3600) = 1.7 \times 10^{-4}$ rad s$^{-1}$, and using an initial phase shift $\Phi \approx \pi/4$ (see App. D2), we obtain the value $\zeta \approx 0.019$, which is underestimated compared with the theoretical estimate of App. D.2 ($\zeta \approx 0.05$).

# 5. IMPACT OF RADIATIVE TRANSFER

## 5.1. DATA

Although at temperatures as low as -50°C, ice supersaturation in observed cirrus clouds as large as 50% may be measured (Korolev & Isaac 2006), statistical data (Krämer et al. 2009) show average values equally distributed slightly above 0%, a situation that seems necessary to trigger damped oscillations when other conditions are favourable, but this is actually linked to the smallness of the radiative contribution. Thus, we shall show below that on one hand damped oscillations can develop with zero net radiative transfer (Case 1), and on the other hand we exhibit situations in which ice supersaturation is higher than in the standard case (Case 0) of Sect. 3.2.3 but associated with a higher net radiative transfer (Cases 2 and 3). The parameters and initial conditions are summarized in Table 2.

We shall display only the oscillatory mode and, following the same analysis procedure as in Sect. 3, we shall estimate period and wavelength as functions of the nominal water saturation, namely $H_{a3}$ at the latitude $z_{H3} = 10$ km, Eqs. (79) in App. A.3, and the radiative transfer ratio $\eta$ which will be kept constant in this section. For each case, we show the trajectory, hodograph, and altitude time profile which are sufficient to determine period $P$, wavelength $\Lambda$, and also ice supersaturation $S_i-1$, radiative correction $R$, and driving factor $S_i-1-R$. The maximum height of the initial hooked cap, $z_{max0}-z_0$, and

the limit altitude of the oscillatory motion, $z_{\infty}$, are also recorded in Table 2.

## 5.2. CASE 1

We assume that $\eta = 1$, i.e. the net radiative transfer is zero (Table 2). This is therefore a neutral case from a radiative standpoint. We choose the growth mode such that $\phi = 2$ and $c_B$ is constant. With a radius $a_0 = 50.17$ µm, steady and oscillatory modes are obtained just by

<table><caption>Table 2 Parameters of test-cases/numerical experiments.</caption><thead><tr><th>CASE NB</th><th>0</th><th>1</th><th>2</th><th>3</th></tr><tr><th>NET RAD. TRANSF.</th><th>R ≠ 0</th><th>R = 0</th><th>R ≠ 0</th><th>R ≠ 0</th></tr></thead><tbody><tr><th>Parameter</th><td>η(z)<br>φ = φ<sub>0</sub><br>c<sub>B</sub> = c<sub>B0</sub></td><td>η = η<sub>0</sub> = 1<br>φ = φ<sub>0</sub><br>c<sub>B</sub> = c<sub>B0</sub></td><td>η = η<sub>0</sub> ≠ 1<br>φ = φ<sub>0</sub><br>c<sub>B</sub> = c<sub>B0</sub></td><td>η = η<sub>0</sub> ≠ 1<br>φ = φ<sub>0</sub><br>c<sub>B</sub> = c<sub>B0</sub></td></tr><tr><th>Modes</th><td>Steady - Osc</td><td>Steady - Osc</td><td>(Steady) - Osc</td><td>(Steady) - Osc</td></tr><tr><th>Life time</th><td>5 h/40 h</td><td>5 h/40 h</td><td>4 h/40 h</td><td>4 h/40 h</td></tr><tr><th>Δt (s)</th><td>0.02/0.04</td><td>0.02/0.04</td><td>0.02/0.04</td><td>0.02/0.04</td></tr><tr><th>x<sub>0</sub> (km)</th><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><th>z<sub>0</sub> (km)</th><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><th>u<sub>0</sub> (m s<sup>-1</sup>)</th><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>w<sub>0</sub> (m s<sup>-1</sup>)</th><td>0.6</td><td>0.6</td><td>0.8</td><td>0.8</td></tr><tr><th>a<sub>0</sub> (µm)</th><td>51.1</td><td>50.17</td><td>52.0</td><td>50.2</td></tr><tr><th>c<sub>0</sub> (µm)</th><td>102.2</td><td>100.34</td><td>104.0</td><td>100.4</td></tr><tr><th>φ<sub>0</sub></th><td>2.0</td><td>2.0</td><td>2.0</td><td>2.0</td></tr><tr><th>ψ<sub>0</sub></th><td>0.80</td><td>0.80</td><td>0.80</td><td>0.80</td></tr><tr><th>c<sub>B0</sub> (µm)</th><td>20.44</td><td>20.068</td><td>20.80</td><td>20.08</td></tr><tr><th>m<sub>0</sub> (µg)</th><td>0.935</td><td>0.885</td><td>0.986</td><td>0.887</td></tr><tr><th>c<sub>0</sub> (µm)</th><td>88.6</td><td>86.9</td><td>90.1</td><td>87.0</td></tr><tr><th>D<sub>0</sub> (µm)</th><td>125</td><td>122</td><td>127</td><td>123</td></tr><tr><th>ρ<sub>ie0</sub> (kg m<sup>-3</sup>)</th><td>675</td><td>675</td><td>675</td><td></td></tr><tr><th>η</th><td>0.9 ≤ η ≤ 1.1</td><td>1</td><td>1.1</td><td>0.9</td></tr><tr><th>z<sub>η0</sub> (km)</th><td>9 ≤ z ≤ 10</td><td>0 - 11</td><td>0 - 11</td><td>0 - 11</td></tr><tr><th>W<sub>a0</sub> (m/s)</th><td>0.6</td><td>0.6</td><td>0.8</td><td>0.8</td></tr><tr><th>z<sub>w0</sub> (km)</th><td>10/9</td><td>10/9</td><td>10/7</td><td>10/7</td></tr><tr><th>H<sub>a3</sub></th><td>0.61</td><td>0.61</td><td>0.68</td><td>0.68</td></tr><tr><th>z<sub>H3</sub> (km)</th><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><th>z<sub>max0</sub>-z<sub>0</sub> (m)</th><td>102</td><td>144</td><td>468</td><td>542</td></tr><tr><th>z<sub>∞</sub> (km)</th><td>9.63</td><td>9.58</td><td>8.58</td><td>8.3</td></tr><tr><th>P (hour)</th><td>10.2</td><td>8.87</td><td>7.8</td><td>6.8</td></tr><tr><th>Λ (km)</th><td>69.8</td><td>64.8</td><td>200</td><td>200</td></tr><tr><th>C<sub>p</sub> (m s<sup>-1</sup>)</th><td>1.85</td><td>2.03</td><td>7.1</td><td>8.2</td></tr><tr><th>ζ</th><td>0.019</td><td>0.015</td><td>0.021</td><td>0.011</td></tr></tbody></table>
switching the updraft base altitude, with lifetimes of 5
hours and 40 hours respectively. For small saturations (H<sub>a3</sub>
= 61%, S<sub>i</sub> ≈ 103%), an oscillation is triggered and slowly
damped, much alike the standard case (Sect. 3.2.3). It is
noteworthy that the maxima along the trajectory (Figure
32a) are sharp at the beginning and tend to become
blunter and blunter, while the minima have nice round
shapes all along the motion, like other parameters
displayed below (z, S<sub>i</sub>-1). The initial amplitude of height
oscillations is approximately 400 m, like in Case 0, and
the oscillation pattern during 40 hours spans a little more
than 250 km (Figure 32a).

The maximum height reached in the head is 144 m, and then oscillation sets in, going through four maxima and four minima with moderate damping looking like the standard solution (Case 0) of Sect. 3.2.3. The period and wavelength derived from simulation results are approximately P ≈ 8.87 hours and Λ ≈ 64.8 km, which yield a phase speed C<sub>p</sub> = 2.03 m s<sup>-1</sup>. The damping ratio is ζ ≈ 0.015. As expected, it is consistent both with the velocity limit u<sub>∞</sub> ≈ 2.1 m s<sup>-1</sup> in the slowly spiralling hodograph (Figure 32b) and with the wind speed at the limit height z<sub>∞</sub> ≈ 9.58 km (Figure 32c), that is using Eq. (68) in layer III: U<sub>0</sub>(z<sub>∞</sub>) = -30 -5(9.58-16) = 2.1 m s<sup>-1</sup>. The altitude time profile (Figure 32c) is qualitatively well described by the theoretical model of damped harmonic oscillator built in App. D.2 with moderate damping (Figure 41a).

Ice supersaturation (Figure 33d) oscillates at the
same period and in phase with altitude (Figure 32c). The
amplitude is approximately 2.5 % at the beginning and it
decreases afterwards, tending towards zero.

The absence of net radiative transfer results in a smaller period than in the standard case (Sect. 3.2.3), thus reflecting the property that the cloud lifetime is increased by radiative transfer (Dobbie & Jonas 2001). This can be called a dual case because the steady counterpart, not displayed, is characteristic of a well-developed cirrus uncinus.

5.3. CASE 2

We assume here that η = 1.1, i.e. the net radiative transfer is positive and corresponds to higher layer in the upper troposphere where heating is predominant in daytime (Table 2). Therefore larger water and ice saturations (H<sub>a3</sub> = 68%, S<sub>i</sub> ≈ 115%), and updraft (W<sub>a0</sub> = 0.8 m s<sup>-1</sup>) are possible and the initial radiative correction can be larger (R ≈ 1.5 %) than in the standard case (R ≈ 0.8 %). Again choosing the growth mode such that φ = 2 and c<sub>B</sub> is constant, steady and oscillatory mode are obtained just by switching the updraft base altitude, with lifetimes of 4 hours and 40 hours respectively, upon specifying a radius a<sub>0</sub> = 52.0 µm, and setting the updraft base altitude at 7 km, lower than in Cases 0 and 1.

<figure><img src="image_36.png" /><figcaption>**Figure 32** Profiles of the parcel trajectory, hodograph, altitude, supersaturation (η = 1; W_a = 0.6 m/s; H_a3 = 61%; z_w0 = 9 km).</figcaption></figure>
The oscillatory mode is triggered again (Figure 33a), with a maximum height at 468 m in the cap followed by a trough at -2.5 km, then damped with wavelength $Λ ≈ 200$ km, period $P ≈ 7.8$ hours, and damping ratio $ζ ≈ 0.02$. According to Eq. (48), the phase speed is therefore $C_p = 200 × 10^3 / (7.8 × 3600) ≈ 7.1$ m s⁻¹. As expected, it approximately matches both the focus $u_∞ ≈ 7.1$ m s⁻¹ of the quickly spiralling hodograph (Figure 33b) and the wind speed at the limit height $z_∞ ≈ 8.58$ km (Figure 33c), that is using Eq. (68) in layer III: $U_a(z_∞) = -30 - 5(8.58 - 16) = 7.1$ m s⁻¹.

The larger phase speed attained when the crystal parcel crosses the critical level ($z_{U0} = 8$ km) increases considerably the wavelength and consequently stretching beyond 900 km the span of the oscillation pattern during 40 hours (Figure 33a). The altitude time profile (Figure 33c) is in good agreement with the theoretical model of damped harmonic oscillator built in App. D.2 with the same period and a larger damping ratio (Figure 41b). Nevertheless, a significant period increase ($≈ 7.5%$) over 40 hours shows that the oscillator is not rigorously harmonic.

<figure><img src="image_37.png" /><figcaption>**Figure 33** Profiles of the parcel trajectory, hodograph, altitude, supersaturation $S_i - 1$, driving factor $S_i - 1 - R$ and radiative correction $R$ ($η = 1.1$; $W_a = 0.8$ m/s; $H_{a3} = 68\%$; $z_{w0} = 7$ km).</figcaption></figure>
Ice supersaturation $S_i-1$ (Figure 33d) oscillates at the same period and in phase with altitude (Figure 33c) like in the standard case (Case 0, Sect. 3.2.3). The amplitude of oscillations is approximately 15% at the beginning and it decreases quickly afterwards, as $S_i-1$ tends towards the limit at 1.2% approximately. The radiative correction $R$ (Figure 33f) is in phase quadrature with $S_i-1$ whereas it was in phase in Case 0, and its limit is approximately 1.2%, equal to that of the supersaturation. Consistently, the growth driving factor $S_i-1-R$ (Figure 33e) oscillates about its limit, zero, with a large amplitude ($\approx 10\%$) at the beginning, larger than in Case 0 ($\approx 2\%$), then quickly decreasing.

The hodograph (Figure 33b) shows an amazing feature consisting apparently in a folding of the spiral at the abscissa $u = 10 \text{ m s}^{-1}$ acting as a wall, a phenomenon that can be easily explained (App. D.2). This is well illustrated by the time profile of the horizontal velocity $u$ (Figure 34), which clearly shows the first two eastward elongations truncated in the vicinity of $U_{amax}$. The subsequent elongations are not affected, and $u$ eventually converges towards $u_{\infty} \approx 7.1 \text{ m s}^{-1}$.

Like in Cases 0 and 1 before, the present damping ratio ($\zeta \approx 0.02$) is also underestimated compared with the theoretical estimate of App. D.2 ($\zeta \approx 0.07$). This case can be considered as isolated because the steady counterpart (not shown) is short ($\approx 4$ hours) and composed of a stretched hooked cap ($\approx 20$ km), without trail.

### 5.4. CASE 3

We now assume that $\eta = 0.9$, i.e. the net radiative transfer is negative and corresponds to lower layers in the upper troposphere where cooling is predominant at night time (Table 2). Note that this is the particular radiative situation depicted by Eq. (25) in Sect. 2.2.2. Other parameters are set to values of the previous case ($H_{a3} = 68\%$, $S_i \approx 115\%$, $W_{a0} = 0.8 \text{ m s}^{-1}$; $\phi = 2$;

$c_B = c_{B0}$). Steady and oscillatory modes are obtained with lifetimes of 4 hours and 40 hours respectively, upon specifying a radius $a_0 = 53.9 \mu\text{m}$, and setting the updraft base altitude at 7 km, lower than in Cases 0 and 1. The initial amplitude of altitude oscillations is approximately 1 km, and it decreases slowly afterwards (Figure 35a, c).

Thus the oscillatory mode is slowly damped with average $\Lambda \approx 200$ km, $P \approx 6.8$ hours, and $\zeta \approx 0.01$. According to Eq. (48), the phase speed is therefore $C_p = 200 \times 10^3/(6.8 \times 3600) = 8.2 \text{ m s}^{-1}$. As expected, it matches approximately both the limit $u_{\infty} \approx 8.5 \text{ m s}^{-1}$ in the hodograph (Figure 35b) and the wind speed at the limit height $z_{\infty} \approx 8.3 \text{ km}$ (Figure 23c), that is using Eq. (68) in layer III: $U_a(z_{\infty}) = -30 - 5(8.3-16) = 8.5 \text{ m s}^{-1}$. Like in the previous case, the motion crossing the critical level ($z_{u0} = 8 \text{ km}$) of maximum wind speed ($U_{amax} = 10 \text{ m s}^{-1}$) induces larger phase speed and wavelength, which causes the oscillation pattern during 40 hours to span more than 900 km (Figure 35a).

Ice supersaturation $S_i-1$ (Figure 35d) oscillates at the same period and in phase with altitude (Figure 35c) like in other cases. Its amplitude is approximately 15% at the beginning and it decreases afterwards, as $S_i-1$ tends towards the limit $-1.2\%$ approximately, slightly below saturation. The radiative correction $R$ (Figure 35f) is again in phase quadrature with $S_i-1$, whereas it was in phase in Case 0, and its limit is close to $-1.2\%$. Consistently, the driving factor $S_i-1-R$ (Figure 35e) oscillates about zero, with a large amplitude ($\approx 15\%$) at the beginning, and as expected it tends towards zero.

We notice that the common limit of $S_i-1$ and $R$ at infinity are equal in magnitude but opposite in Cases 2 and 3 ($R = \pm 1.2\%$), a situation corresponding to opposite radiative initial conditions ($\eta - 1 = \pm 0.1$). Moreover, the period shape in $R$ oscillation is slightly asymmetric, the ascending half-period being steeper than the descending one in Case 3 (Figure 35f) and reciprocally in Case 2 (Figure 33f).

The hodograph (Figure 35b) shows here again a folding of the spiral at the abscissa $u = 10 \text{ m s}^{-1}$, but for five complete periods. The phenomenon, explained in App. D.2, is illustrated by the time profile of the horizontal velocity $u$ (Figure 36), which clearly shows that all of the eastward elongations are truncated and folded in the vicinity of $U_{amax}$.

This case can be considered as isolated because the steady counterpart (not shown) is shorter (4 hours) and composed of a hooked cap with a horizontal extension of 20 km, without trail. The comparison of estimated damping ratios $\zeta$ (Table 2) reveals that Case 2 is more damped than Case 3.

<figure><img src="image_38.png" /><figcaption>**Figure 34** Time profile of the horizontal velocity $u$ ($\eta = 1.1$; $W_a = 0.8 \text{ m s}^{-1}$; $H_{a3} = 68\%$; $z_{u0} = 7 \text{ km}$).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>**Figure 35** Profiles of the parcel trajectory, hodograph, altitude, supersaturation $S_i-1$, driving factor $S_i-1-R$ and radiative correction $R$ ($\eta = 0.9$; $W_a = 0.8$ m/s; $H_{a3} = 68\%$; $z_{w0} = 7$ km).</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>**Figure 36** Time profile of the horizontal velocity $u$ ($\eta = 0.9$; $W_a = 0.8$ m/s; $H_{a3} = 68\%$; $z_{w0} = 7$ km).</figcaption></figure>

# 6. DISCUSSION

## 6.1. NON-OSCILLATORY MODE

Early simple analytical models of cirrus produced parabolic or quasi-parabolic shapes (Ludlam 1948; Marshall 1953). Among a lot of works on cirrus uncinus clouds supported by field observations, two series of papers using a Lagrangian approach, the first one published in the sixties (Magono et al. 1967; Yagi et al. 1968; Harimaya 1968; Yagi 1969) and the second one in the seventies (Heymsfield 1973; Heymsfield 1975a; Heymsfield 1975b; Heymsfield 1975c) retained our attention.

The first series reports stereoscopic observations that show real motions to be three-dimensional, the cirrus
uncinus being also bent in a horizontal plane (Magono et al. 1967; Yagi et al. 1968; Harimaya 1968; Yagi 1969). The directions of cloud and wind are slightly tilted by 10° or so at the top of the trail and they are practically coincident at the bottom (Yagi et al. 1968).

The simple model based on these observations (Harimaya 1968) solves the mass growth equation Eq. (10) for crystal size and inserts the result in an analytic expression for the falling velocity, which in turn is integrated to yield the height as a function of time. Likewise, abscissa is obtained by integration of the horizontal drift velocity. Six crystal shapes are investigated and the wind shear is chosen in the range 0–7 × 10⁻³ s⁻¹. Consistently with our result, the maximum fall speed of 0.5 m s⁻¹ is reached 500 m below the origin. Nevertheless, temperature and humidity are assumed constant during the motion and, in the absence of updraft, no hooked cap is retrieved at the root of the trail. The cap is supposed to form within a turbulent layer between two stable layers (Yagi 1969).

The second series investigates the trajectory of five classes of crystals of bullet-rosette type, mostly observed in cirrostratus clouds (Heymsfield 1973; Heymsfield 1975a). By numerical integration of ODE analogous to our Eqs. (42), (43), (44) over a shorter time (3000 s), with an updraft velocity of 1 m s⁻¹, characteristic shapes with a hooked head and a trail are retrieved. The major differences with our work lie in the facts that the buoyancy force and effect of latent heat release are included in the entrainment vertical velocity W_d, that radiative effects are not taken into account and that the crystal habit is more complex and it is based on a mass-size relationship. As a consequence, the vertical extent of the head is somewhat larger, reaching 500 m.

## 6.2. OSCILLATORY MODE

The response of atmospheric layers to a Kelvin-Helmholtz instability is an important issue that may produce wavy patterns in the upper troposphere (billow clouds; cirrus fluctus). This phenomenon is usually discussed by means of the Richardson number Ri, defined by the expressions (Lynch et al. 2002):

$$ 
\begin{aligned}
Ri_d &= \frac{\Omega_{Bvd}^2}{\omega^2} \\
Ri_m &= \frac{\Omega_{BVm}^2}{\omega^2}
\end{aligned}
\qquad (59) $$

with the Brunt-Väisälä angular frequencies Ω_Bvd and Ω_Be, defined by (Durran & Klemp 1982):

$$ 
\begin{aligned}
\Omega_{Bvd}^2 &= g \frac{d}{dz} (\ln \theta_d) \\
\Omega_{Bve}^2 &= g \frac{d}{dz} (\ln \theta_e)
\end{aligned}
\qquad (60) $$

and Ω_BVm in a moist atmosphere, defined by (Durran & Klemp 1982):

$$ \Omega_{BVm}^2 = g \left\{ \frac{1 + \frac{L_v}{R_a} q_s}{1 + \frac{\varepsilon L_v^2 q_s}{c_{pa} R_a T^2} \left( \frac{d}{dz}(\ln \theta_d) + \frac{L_v}{c_p T} \frac{dq_s}{dz} \right) - \frac{dq_l}{dz} } \right\} \qquad (61) $$

and the horizontal wind shear ω:

$$ \omega = \frac{dU_a}{dz} \qquad (62) $$

The dry and equivalent potential temperatures θ_d and θ_e are usually defined by (Durran & Klemp 1982):

$$ 
\begin{aligned}
\theta_d &= T \left( \frac{p_r}{p} \right)^{R_i/c_p} \\
\theta_e &= \theta_d \exp \left( \frac{L_v q_s}{c_{pa} T} \right)
\end{aligned}
\qquad (63) $$

Vertical profiles of the local Richardson number (Figure 37) obtained with our data (App. A) show that the minimum value, reached at the inversion altitude of 8 km is approximately 2.95 in the dry case (solid line), and since it is everywhere larger than the critical value 0.25, we can conclude according to the criterion that the dry atmosphere is stable (Lynch et al., 2002). In contrast, the moist and equivalent profiles (dotted and dashed lines), which are not very different from each other, clearly show that a Kelvin-Helmholtz instability is likely to occur locally in such an atmosphere above and below 7 km, since Ri_m < 0.25. Published profiles show the same behaviour (Wada et al. 2005; Spichtinger 2014). Other authors (Dobbie & Jonas 2001) derived a radiative stability number RSn and a criterion for onset of convective instability as: 0 < RSn < 1. With a typical heating rate H = 0.1 K d⁻¹ at the altitude z = 9 km, they show that their criterion could be satisfied.

In contrast, the oscillations we have put into evidence in Sect. 3.2.3 are due to the interplay of the crystal's mass variation by deposition-evaporation and the updraft, and therefore it is neither a buoyancy effect relevant of Kelvin-Helmholtz instability nor a radiative instability, though the above remarks clearly suggests that these

<figure><img src="image_41.png" /><figcaption>Figure 37 Atmospheric profiles of Richardson numbers.</figcaption></figure>
kinds of instability may develop in the free atmosphere at the altitude where our cirrus develops, or slightly below. Moreover, we notice that the two kinds of motions we found are such that the life time of a well-developed cirrus uncinus in steady mode is of same order as the half-period of the associated oscillatory pattern.

These oscillations were not found in the quoted work (Heymsfield 1975b), but the mention of “long-lasting cirrus turrets” originating from the trail and produced either by convection above the stable layer or by evaporation of ice crystals and the reference to generators of pulsating-type strongly suggests that they may be related to such observations. Nevertheless, low to moderate ice supersaturation in the range 3–15% seems to be more favorable to the onset of oscillations in our situation, in agreement with statistical data (Krämer et al. 2009), though values as high as 50% could be expected at such low temperature of –50°C according to some other data (Korolev & Isaac 2006).

The parcels being generated continuously in the head and feeding the trail downstream, the phenomenon may also be connected to the more recent concept of “Mesoscale Uncinus Complex” (MUC) that was proposed to explain the grouping of individual cirrus uncinus so as to form mesoscale structures with dimensions ranging from 15 km to 100 km (Sassen et al. 1989), consistent with our long-wavelength oscillatory pattern. The idea of MUC was confirmed afterwards by a series of observations (Wang 2004; Wang & Sassen 2008). Nevertheless, the phase speed $C_p$ we estimate is related to the wind velocity at the altitude of oscillations ($C_p \approx 2–8 \text{ m s}^{-1}$) and it is much smaller than the wind speeds ($U_a \approx 30 \text{ m s}^{-1}$) involved in MUC and related to jet stream (Wang & Sassen 2008).

A first investigation (Demoz et al. 1998) using wavelets to analyze observations of low cirrus uncinus culminating at 8 km, shows spatial scale features of 30 km in extension associated with relatively high frequencies ($\approx 0.007 \text{ Hz}$), corresponding to periods of 140 s. As shown by our Figure 37, convective instability is likely to develop at the altitude of their observations, so that the authors logically conclude that the high frequency spectrum (slope –5/3) is due to turbulence, and the low frequency component (slope –3) is due to the interaction of convection with gravity waves.

Another wavelet analysis (Wang & Sassen 2008) focusing on a high-altitude MUC (10–13.5 km) also reveals that two superimposed spectral components emerge: a small-period one (10–100 s) in a small scale range (0.4–4 km) with slope –5/3 relevant of turbulence in embedded cells, and a large-period one (100–1000 s) in a large-scale range (5–10 km) with slope –3, not relevant of a uniquely identified dynamical phenomenon, but rather resulting of the complex interaction of propagating

gravity wave through turbulent layers. These efforts show that the phenomenon deserves further investigations in a range of larger periods ($> 1000 \text{ s}$).

# 7. CONCLUSION AND PROSPECTS

## 7.1. RESULTS AND FINDINGS

According to our purpose, using a Lagrangian approach with prescribed wind and updraft, we could retrieve the two-dimensional shape as well as dynamical and microphysical features of crystals trajectories leading to a well-developed cirrus uncinus. Moreover, besides this “standard” steady mode shaping a cirrus cloud, we put into evidence a self-sustained, damped quasi-harmonic oscillatory mode of long period ($\approx 6–10 \text{ hours}$) and wavelength ($\approx 60–200 \text{ km}$), never reported before, as far as we know. Despite its apparent simplicity, our model includes many elementary nonlinear processes, and it shows how the height of the updraft base, and radiative transfer competing with supersaturation, by modulating water vapour-ice phase changes that drive the crystal’s mass variation, can significantly modify the development of a mean ice crystal parcel. Nevertheless, small to moderate ice supersaturations in the range 3–15% seem to be more favorable to the onset of oscillations than values as high as 50% that could be expected at such a low temperature of –50°C (Korolev & Isaac 2006). Of course, this an idealized picture that requires an updraft stable over several tens of kilometers in the upper troposphere.

The complexity of our two-dimensional model can be quantified by 43 degrees of freedom, composed of the 36 constants implied in the ambient vertical profile conditions (App. A): dynamics ($U_a, W_a$) and thermodynamics ($T_a, p_a, H_a, \eta_a$), supplemented with the 7 initial scalar conditions: position ($x_0, z_0$), velocity ($u_0, w_0$), column half-width ($a_0$), aspect ratio ($\phi$), hollowness factor ($\psi$). Of the total, we matched only five of them: $z_{w0}, W_{a0}, H_{a0}, \eta_a$ and $a_0$.

By a detailed sensitivity analysis (Sect. 4) we have shown that the crystal size has a significant impact on the head extension and decay duration of the steady mode, but a smaller impact on the period of the oscillatory mode. Moreover, we examined the influence of radiative transfer (Sect. 5) in three specific situations, showing in particular that heating has a larger damping effect than cooling. This variability is a useful structural input for the modelling of cloud texture, as imagined in the prospective section below.

We notice that the two kinds of motions we found are such that the life time of a well-developed cirrus uncinus in steady mode ($\approx 5 \text{ hours}$) is of same order as the half-period of the associated oscillatory mode ($\approx 8–10 \text{ hours}$).
An analytic model in App. D. recovers the computed shapes of the two modes and expressions of the period and damping ratio of the oscillatory mode are derived. A slight asymmetry of oscillations and a significant increase (≈ 7.5%) of period with time reveal some inharmonicity of the oscillator when radiative transfer is of constant sign, either positive i.e. heating (Sect. 5.3) or negative i.e. cooling (Sect. 5.4).

The theoretical background of the underlying Liénard differential system governing our system, which produces bifurcations due to a particular parameter (level of updraft base), is linked with the notion of strange attractor, actually originated in the astrophysical modelling of the dynamo effect maintaining stellar and planetary magnetic fields (Rikitake 1958), although roots are found in the modelling of population dynamics (Verhulst 1838; Vogels et al. 1975). Thus, the phase shift between the oscillations of mass and supersaturation that we put into evidence (Sect. 3.2.3) is quite similar to that appearing in predator-prey models (Koren & Feingold 2011).

In clouds, conceptual one-dimensional models (Wacker 2006; Spreitzer et al. 2017) have shown the importance of stratification and updraft in the interaction between a two-layer cloud system, and how time oscillations are generated according to the choice of parameters. Nevertheless, the periods are much shorter (15 to 50 minutes in the first model, 1.5 hours in the second) than those exhibited in our problem (6 to 10 hours), and the lifetime of the cirrus uncinus.

The possible connection of our theoretical model with the concept of Mesoscale Uncinus Complex (MUC) (Sassen et al. 1989) and periodic generating cells or ascending turrets (Heymsfield 1973; Heymsfield 1975b, 1975c) that we mentioned in the discussion (Sect. 6.2) would be a valuable application that deserves further investigations, looking for cells of several tens of kilometers pulsating over periods of several hours.

Assuming no motion along y coordinate, we also neglected the three-dimensional development of the cloud (see App. A3). Actually, radar and lidar observations usually provide mappings in vertical planes (Hogan & Kew 2005; Wada et al. 2005). At extremely low temperatures (T - T_r ≤ -50°C) we reasonably neglected liquid water in the cloud (Cziczo et al. 2013).

## 7.2. PROSPECTS

As mentioned in Sect. 2.2.3, apart from single column, bullet-rosette is the most common crystal shape in cirrus uncinus (Heymsfield & Iaquinta 2000; Schmitt & Heymsfield 2007) and therefore we could also consider bullet-rosettes instead of hollow columns. The relevant parameters governing the growth are initial mass, capacitance, ventilation coefficient and drag coefficient. Mass will be roughly multiplied by the number N of bullets:

$m = N m_b$, and, instead of Eq. (41), the capacitance could be modelled as (Chiruta & Wang 2003):

$$C = 0.434 a N^{0.257} \qquad (64)$$

In a first step, a simple 4-bullet rosette of plane type 4-4 (Heymsfield and Iaquinta 2000; Westbrook et al. 2008) could be used. Therefore, assuming bullets similar to our hollow column except for the ends, mass would be approximately multiplied by 4, and we would obtain $C/a = 0.62$, instead of $C/a = 1.7$, as shown from our Table 2. In other words, the crystal capacitance would be divided by 3 approximately.

Instead of Eq. (16) for a hollow column, the ventilation coefficient $f_v$ for broad-branch crystals and $1 < Re < 120$ is given by (Ji & Wang 1999):

$$f_v = 1 + 0.354 \frac{X}{10} + 3.55 \left( \frac{X}{10} \right)^2 \qquad (65)$$

In our problem $X \approx 1.5$ and $f_v = 1.4$. A factor 2 on the size would increase Re likewise and consequently would multiply X by 1.4. Now, using (65), we would obtain: $f_v = 1.2$, and this would probably not have a significant impact on crystal growth.

So far, the drag coefficient for bullet-rosettes seems to be an open concern, subject to constant investigations dealing with non-spherical particles and free-falling snowflakes (Heymsfield 1972; Haider & Levenspiel 1989; Heymsfield & Westbrook 2010; Vazquez-Martin et al. 2021; Aguilar et al. 2022). These approaches introduce the concepts of sphericity, projected area and Best number.

With 3 degrees of freedom ($a, c, c_b$) or ($a, \phi, \psi$) for each bullet, we could assume that all of the bullets grow at the same rate, with the same rules as for a single column. Nevertheless, it is known that the equilibrium shape of the crystal at each time should be calculated so as to minimize the total surface energy for a given volume (Pruppacher & Klett 1978). Moreover, at larger supersaturations, the capacitance model based on Eq. (10) presents limitations, that have been discussed by introducing the concept of impedance (Peter and Baker 1996; Nelson & Baker 1996). Alternative models of growth-evaporation based on this concept such as the TLK (terrace-ledge-kink) model have been proposed and applied to hexagonal ice crystals (Wood, Nelson, & Calhoun, 2001).

The vertical shear of the wind $U_a$ and the base of the updraft $W_a$ are essential input of our dynamical conditions, in order to produce the characteristic hooked shape of cirrus uncinus. The profiles of temperature and humidity are suggested by actual measurements in the troposphere but they may be varied. The trajectory is also very sensitive to the initial conditions, which may cancel the formation of a hooked cap.

Especially the wind field $U_a$ we imposed throughout the present work has been chosen realistically according
to published works (Kew 2003; Mace et al. 2005), but there is no standard profile for it, so that its minimum and maximum velocities could be varied in Eqs. (68) (App. A1), thus impacting the wavelength $\Lambda$ of the crystal's oscillation. Actually, the period $P$ being an intrinsic property of the phenomenon, the wavelength can be derived from it via the phase speed $C_p$ defined by Eq. (48), which is equal to the wind speed at the limit height $z_\infty$:

$$ \Lambda = C_p P \approx U_a(z_\infty)P \qquad (66) $$

Thus, increasing or decreasing $U_a$ in the layer under concern (III) will increase or decrease $\Lambda$ accordingly. This gives some flexibility to the scale of our phenomenon, if it can be related to MUC (Sassen et al. 1989) and pulsating-type generators (Heymsfield 1973, 1975b, 1975c). In our situation, fast winds exist in layers (IV) above the layer of interest (III), but the wind shear we adopted, in agreement with published data, makes wind slower at the altitude where we produce the uncinus head, that could not arise otherwise.

Whereas Coriolis force can be neglected in the non-oscillatory mode, it may be included in a more accurate model of the oscillatory mode, since the extension of the structure over several hundreds of kilometers tends to the synoptic scale and may not be completely negligible compared to the Earth's dimensions. Although the crystal's motion is driven mainly by gravity and steady winds, atmospheric turbulence would add random degrees of freedom, and so contribute to the wispy appearance of cirrus. The governing system Eq. (1) would be modified as follows:

$$ \begin{cases} \ddot{x} = -k(\dot{x} - U_a - u_a') \\ \ddot{z} = -k(\dot{z} - W_a - w_a') \end{cases} \qquad (67) $$

In the assembly of crystals forming a parcel, collisions between particles inevitably produce mechanical interactions, aggregation, and also electric interactions via the electric charges that are formed.

Moreover, the influence of atmospheric waves on the evolution of cirrus (Lin et al. 1998; Podglajen et al. 2018; Prasad et al. 2019; Kärcher et al. 2019) is a matter of concern, since it is invoked in MUC. Though they are not the driving mechanism of the oscillatory mode we found, as shown in Sect. 6.2, Kelvin-Helmholtz waves may develop at the altitude of 7 km or above. In a more elaborate model, these additional effects could be included in the modelling equations (67), and they would probably smear the motions about the two basic behaviours reported in the present work.

From an experimental viewpoint, observations would be welcome for investigating the intermediary spatial range between meso- and synoptic scale (300–1000 km), and longer time periods (>1000s) where the damped oscillations we numerically exhibited could be hopefully

detected, in connection with MUC or ascending turrets. From a theoretical viewpoint, approximate analytical solutions of the differential system of Liénard-type governing the crystal's evolution should be searched in order to derive more accurate expressions of the angular frequency $\Omega$ and the damping ratio $\zeta$ of underdamped oscillatory solutions, for which a demonstration has been sketched in App. D.2.

Extending the sensitivity analysis of Sect. 4 to an assembly of crystals of different sizes in a parcel and assuming particle densities would make possible the modelling of size distributions (Heymsfield, Schmitt, & Bansemer, 2013) and consequently the calculation of ice water content (IWC). Actually, apart from the understanding of natural phenomena, the generation of realistic cloudy scenes with radiative transfer is another valuable extension of the present work. In that purpose, we can imagine applying a textured IWC in the spatial and temporal frame elaborated herein, using spectral or fractal methods and the dynamic-microphysical output of our model. Such an approach has been intensively implemented for years to render various types of clouds, from stratocumulus and cumulus to cirrus (Cianciolo 1993; Evans & Wiscombe 2004; Hogan & Kew 2005; Sölch & Kärcher 2010; Szczap et al. 2014). In the case of liquid water clouds, a model of water content was widely used in cloud generators (Feddes 1974).

We also used such methods for stratocumulus (Berton 2008). Nevertheless, it seems that the approach should be different for stratocumulus and in-situ cirrus clouds: while the former are essentially due to convection, the latter evolve under a combination of advection and free fall. This issue has been considered in recent models (Hogan & Kew, 2005 Sölch & Kärcher 2010; Szczap et al. 2014) from a macrophysical point of view in order to render fall streaks.

Moreover, the very nature of the cloud particles – droplets in the former case, ice crystals in the second one – suggests that the medium, owing to the shape and orientation of crystals (Hashino et al. 2014; Hashino et al. 2016), is strongly anisotropic in the second situation. In that respect, a set of differential equations for the angular momentum could supplement the system of Eqs. (67), including the torques caused by fluid friction and the fact that the column is hollow.

The present approach enables the generation of horizontal inhomogeneities related to the history of ice particles along the trajectory, and this point is most important since these are known to be an important input in the modelling of radiative properties of cirrus clouds (Liou & Rao 1996; Buschmann et al. 2022; Kew 2003; Kokhanovsky 2003; Fauchez et al. 2014). Up to now, a few methods have been devised to remedy this problem (Shonk & Hogan 2008).
# APPENDIX A. ATMOSPHERIC PROFILES

## A.1. ATMOSPHERIC MOTIONS

The imposed atmospheric motions are composed of a uniform horizontal wind $U_a$ with a given vertical profile (Figure 1), and a constant vertical updraft $W_a$, operating only above a critical level $z_{w0}$ equal to either 10 km, 9 km, or 7 km (Figure 2). We built our horizontal wind profiles from published data (Kew 2003; Mace et al. 2005). A two-dimensional mapping in a vertical plane of combined observational data (Wada et al. 2005) clearly shows such patterns.

The horizontal wind profile we chose reflects the horizontal shear and the vertical convection necessary for the formation of the virga. The vertical profile of horizontal wind $U_a(z)$ (Figure 1) is given as a piecewise linear function composed of four segments:

$$ egin{cases} U_a(z) = \frac{U_{a1}}{z_{U1}} z & 0 \le z \le z_{U1} & (I) \\ U_a(z) = U_{a1} + \frac{U_{a0} - U_{a1}}{z_{U0} - z_{U1}} (z - z_{U1}) & z_{U1} \le z \le z_{U0} & (II) \\ U_a(z) = U_{a2} + \frac{U_{a0} - U_{a2}}{z_{U0} - z_{U2}} (z - z_{U2}) & z_{U0} \le z \le z_{U2} & (III) \\ U_a(z) = U_{a2} & z_{U2} \le z & (IV) \end{cases} \quad (68) $$

with the following 6 constants:

$$ egin{cases} z_{U1} = 5 \text{ km} \\ z_{U0} = 8 \text{ km} \\ z_{U2} = 16 \text{ km} \end{cases} \quad \begin{cases} U_{a1} = 5 \text{ m s}^{-1} \\ U_{a0} = 10 \text{ m s}^{-1} \\ U_{a2} = -30 \text{ m s}^{-1} \end{cases} \quad (69) $$

It is noteworthy that, with this choice of constants, $U_a$ vanishes above the ground at the altitude $z_0 = 10$ km. Moreover, the resulting wind shear $\omega$, respectively equal to 1.0, 1.3 and -5 m s⁻¹ km⁻¹ in the ranges 0–5, 5–8 and 8–16 km, is quite consistent with observational data ($|\omega| < 23 \text{ m s}^{-1} \text{ km}^{-1}$) (Heymsfield 1975b; Heymsfield 1975c; Wada et al. 2005). As we show in section 3.2, the results are especially sensitive to $z_{w0}$.

The updraft profile $W_a(z)$ is chosen so as ensure that $W_a$ is constant, non-zero above $z_{w0} = 10$ km (1st mode or regime) or $z_{w0} = 7$ or 9 km (2nd mode or regime) and zero below $z_{w1} = z_{w0} - \Delta z_0$, with a transition layer of thickness $\Delta z_0 = 0.3$ km, the transition profile being modelled by a linear function (Figure 2):

$$ egin{cases} W_a(z) = 0 & 0 \le z \le z_{W1} & (I) \\ W_a(z) = W_{a0} \frac{z - z_{W1}}{z_{W0} - z_{W1}} & z_{W1} \le z \le z_{W0} & (II) \\ W_a(z) = W_{a0} & z_{W0} \le z & (III) \end{cases} \quad (70) $$

with the following constants:

$$ egin{cases} \text{Cases 0 & 1} \quad \begin{cases} z_{W0} = 9 \text{ km} ; z_{W1} = 8.7 \text{ m s}^{-1} \\ z_{W0} = 10 \text{ km} ; z_{W1} = 9.7 \text{ m s}^{-1} \end{cases} & W_{a0} = 0.6 \text{ m s}^{-1} \\ \text{Cases 2 & 3} \quad \begin{cases} z_{W0} = 7 \text{ km} ; z_{W1} = 6.7 \text{ m s}^{-1} \\ z_{W0} = 10 \text{ km} ; z_{W1} = 9.7 \text{ m s}^{-1} \end{cases} & W_{a0} = 0.8 \text{ m s}^{-1} \end{cases} \quad (71) $$

We verify that such a two-dimensional flow satisfies the equation of mass conservation for air considered as an incompressible fluid:

$$ \frac{\partial U_a}{\partial x} + \frac{\partial W_a}{\partial z} = 0 \quad (72) $$

Thus in the 300 m-thick transition layer, the gradient of updraft $\partial W_a / \partial z$ produces a shear of the horizontal wind $\partial U_a / \partial x$ in the x-direction equal to

$$ \frac{\partial U_a}{\partial x} = -\frac{\partial W_a}{\partial z} = -\frac{W_{a0}}{z_{W0} - z_{W1}} = -\frac{0.6}{300} \approx -2.0 \text{ m s}^{-1} \text{ km}^{-1} \quad (73) $$

that is smaller in magnitude than the vertical shear at the same altitude ($\omega = -5 \text{ m s}^{-1} \text{ km}^{-1}$). We notice that in a three-dimensional description of the motion, this shear may affect the perpendicular component of wind, $V_a$, here set to zero, according to the complete equation:

$$ \frac{\partial U_a}{\partial x} + \frac{\partial V_a}{\partial y} + \frac{\partial W_a}{\partial z} = 0 \quad (74) $$

and create a horizontal shear of $V_a$ in the y-direction. Now, things being as they are, since we assume that $U_a$ is a function of z alone we shall neglect this horizontal shear of $U_a$ compared to the vertical one, in solving the system (42).

## A.2. TEMPERATURE AND PRESSURE

As suggested by measurements (Baum et al. 2000; Kew 2003), the vertical profile of temperature $T_a(z)$ (Figure 3a) is given as a piecewise linear function composed of five segments:

$$ egin{cases} T_a(z) = \frac{T_{a1} - T_{a0}}{z_{T1} - z_{T0}} z + \frac{T_{a0}Z_{T1} - T_{a1}Z_{T0}}{z_{T1} - z_{T0}} & 0 \le z \le z_{T1} & (I) \\ T_a(z) = \frac{T_{a2} - T_{a1}}{z_{T2} - z_{T1}} z + \frac{T_{a1}Z_{T2} - T_{a2}Z_{T1}}{z_{T2} - z_{T1}} & z_{T1} \le z \le z_{T2} & (II) \\ T_a(z) = \frac{T_{a3} - T_{a2}}{z_{T3} - z_{T2}} z + \frac{T_{a2}Z_{T3} - T_{a3}Z_{T2}}{z_{T3} - z_{T2}} & z_{T2} \le z \le z_{T3} & (III) \\ T_a(z) = \frac{T_{a4} - T_{a3}}{z_{T4} - z_{T3}} z + \frac{T_{a3}Z_{T4} - T_{a4}Z_{T3}}{z_{T4} - z_{T3}} & z_{T3} \le z \le z_{T4} & (IV) \\ T_a(z) = T_{a4} & z_{T4} \le z & (V) \end{cases} \quad (75) $$

with the following 10 constants:

$$ egin{cases} z_{T0} = 0 \text{ km} \\ z_{T1} = 2 \text{ km} \\ z_{T2} = 8 \text{ km} \\ z_{T3} = 14 \text{ km} \\ z_{T4} = 20 \text{ km} \end{cases} \quad \begin{cases} T_{a0} = T_r + 20 \text{ K} \\ T_{a1} = T_r + 0 \text{ K} \\ T_{a2} = T_r - 50 \text{ K} \\ T_{a3} = T_r - 60 \text{ K} \\ T_{a4} = T_r - 60 \text{ K} \end{cases} \quad (76) $$

while the pressure (Figure 3b) is given as an exponential function of height:

$$ p_a(z) = p_0 \exp\left(-\frac{z}{h}\right) \quad (77) $$

with the two constants: $p_0 = 101493$ Pa and $h = 7.5$ km.
### A.3. HUMIDITY

Likewise, the vertical profile of relative humidity with respect to liquid vapour $H_a(z)$ in the atmospheric clear sky (Figure 4a, solid line), as suggested by measurements in clear sky (Baum et al. 2000) and in an environment favorable to the generation of a cirrus cloud (Fusina and Spichtinger 2010), is given as a piecewise linear function composed of six segments:

$$
\left\{
\begin{aligned}
H_a(z) &= \frac{H_{a1} - H_{a0}}{z_{H1} - z_{H0}} z + \frac{H_{a0}z_{H1} - H_{a1}z_{H0}}{z_{H1} - z_{H0}} & 0 \le z \le z_{H1} & \text{(I)} \\
H_a(z) &= \frac{H_{a2} - H_{a1}}{z_{H2} - z_{H1}} z + \frac{H_{a1}z_{H2} - H_{a2}z_{H1}}{z_{H2} - z_{H1}} & z_{H1} \le z \le z_{H2} & \text{(II)} \\
H_a(z) &= \frac{H_{a3} - H_{a2}}{z_{H3} - z_{H2}} z + \frac{H_{a2}z_{H3} - H_{a3}z_{H2}}{z_{H3} - z_{H2}} & z_{H2} \le z \le z_{H3} & \text{(III)} \\
H_a(z) &= \frac{H_{a4} - H_{a3}}{z_{H4} - z_{H3}} z + \frac{H_{a3}z_{H4} - H_{a4}z_{H3}}{z_{H4} - z_{H3}} & z_{H3} \le z \le z_{H4} & \text{(IV)} \\
H_a(z) &= \frac{H_{a5} - H_{a4}}{z_{H5} - z_{H4}} z + \frac{H_{a4}z_{H5} - H_{a5}z_{H4}}{z_{H5} - z_{H4}} & z_{H4} \le z \le z_{H5} & \text{(V)} \\
H_a(z) &= H_{a5} & z_{H5} \le z & \text{(VI)}
\end{aligned}
\right.
$$

with the following 12 constants:

$$
\left\{
\begin{array}{l}
z_{H0} = 0 \text{ km} \\
z_{H1} = 2 \text{ km} \\
z_{H2} = 4 \text{ km} \\
z_{H3} = 10 \text{ km} \\
z_{H4} = 15 \text{ km} \\
z_{H5} = 20 \text{ km}
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
H_{a0} = 0.20 \\
H_{a1} = 0.30 \\
H_{a2} = 0.40 \\
H_{a3} = 0.61 \\
H_{a4} = 0.20 \\
H_{a5} = 0
\end{array}
\right.
\quad
(79)
$$

On the same figure is then plotted (dotted line) the profile of $H_i$ derived from $H_{a1}$ first by calculating the water vapour pressure $p_v$ as solution of Eq. (84) below ($H_i = H_a$), then by estimating $S_i$ by Eq. (12b) and finally by deriving $H_i$ by Eq. (13b). It is noticeable that $H_i > H_p$ as it can be easily proved from Eqs. (12) and (13), since $p_{vsl} > p_{vsi}$. We also plotted the profiles of $S_i$ and $S_i$ (Figure 4b), though they are very close to those of $H_i$ and $H_i$ respectively. Moreover, the profiles of $H_i$ and $S_i$ are plotted down to 2 km, altitude at which the temperature reaches the frozen point ($T_a = 273.15$ K). We notice that our profiles of $S_i$ are quite consistent with those of other simple models of cirrus uncinus (Harimaya 1968).

### A.4. RADIATIVE FLUX DENSITIES

Adopting Zeng's parameterization (Zeng 2008; Zeng et al. 2021), we assume that the vertical profile of the radiative flux density in the infrared is described by the ratio $\eta_a(z)$ defined by Eq. (23) and can be given in the upper troposphere as a piecewise linear function composed of three segments:

$$
\left\{
\begin{array}{l}
\eta_a(z) = \eta_{a0} \\
\\
\eta_a(z) = \frac{\eta_{a1} - \eta_{a0}}{z_{\eta1} - z_{a0}} z + \frac{\eta_{a0}z_{\eta1} - \eta_{a1}z_{a0}}{z_{\eta1} - z_{a0}} \\
\\
\eta_a(z) = \eta_{a1} \\
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
0 \le z \le z_{a0} \qquad \text{(I)} \\
z_{\eta1} \le z \le z_{\eta2} \qquad \text{(II) (80)} \\
z_{\eta2} \le z \qquad \text{(III)}
\end{array}
\right.
$$

with the following 4 constants:

$$
\left\{
\begin{array}{l}
z_{\eta0} = 9 \text{ km} \\
z_{\eta1} = 10 \text{ km}
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
\eta_{a0} = 0.9 \\
\eta_{a1} = 1.1
\end{array}
\right.
\quad
(81)
$$

We lowered the altitude of the transition $\eta_a = 1$ because the original work (Zeng 2008) deals with tropical cirrus clouds while we are concerned with mid-latitude ones. The profile $\eta_a(z)$ is plotted on Figure 5: the transition ($\eta_a = 1$) occurs at an altitude $z = 9.5$ km.

## APPENDIX B. CALCULATION OF THE AMBIENT PRESSURE OF WATER VAPOUR

Since the ambient relative humidity $H_i$ with respect to (wrt) liquid water is given ($H_i = H_a$), it is necessary to derive the pressure $p_v$ of water vapour by inversion of the relation (13a), with the total pressure such that:

$$p = p_a + p_v \quad (82)$$

we obtain the relation:

$$H_i = \frac{p_v p_a + p_v - p_{vsl}}{p_{vsl}} \quad (83)$$

and after rearranging, a quadratic equation for $p_v$ can be finally written:

$$p_v^2 + (p_a - p_{vsl})p_v - H_i p_a p_{vsl} = 0 \quad (84)$$

This equation has two real solutions because its discriminant $\delta$ is essentially positive:

$$\begin{aligned}
p_v &= \frac{-p_a + p_{vsl} \pm \sqrt{\delta}}{2} \\
\delta &= (p_a - p_{vsl})^2 + 4H_i p_a p_{vsl}
\end{aligned}
\quad (85)$$

Nevertheless, only the solution with + is physical, because it yields $p_v = p_{vsl}$ as $H_i \to 1$.

## APPENDIX C. SATURATION PRESSURES OF WATER VAPOUR

According to Sonntag's formulation, the saturation pressures of water vapour over liquid water and ice are respectively expressed as (Sonntag 1990; Alduchov & Eskridge 1996; Murphy & Koop 2005):

$$
\left\{
\begin{aligned}
p_{vsl}(T) &= \exp\left(-\frac{C_{w-1}}{T} + C_{w0} + C_{w1}T + C_{w2}T^2 + C_{w3}\ln T\right) \\
p_{vsl}(T) &= \exp\left(-\frac{C_{i-1}}{T} + C_{i0} + C_{i1}T + C_{i2}T^2 + C_{i3}\ln T\right)
\end{aligned}
\right.
\quad (86)
$$

with the following coefficients:
$$
\left\{
\begin{array}{l}
C_{w-1} = 6096.9385 \\
C_{w0} = 21.2409642 \\
C_{w1} = -2.711193 \times 10^{-2} \\
C_{w2} = 1.673952 \times 10^{-5} \\
C_{w3} = 2.433502
\end{array}
\right.
\qquad
\left\{
\begin{array}{l}
C_{i-1} = 6024.5282 \\
C_{i0} = 29.32707 \\
C_{i1} = 1.0613868 \times 10^{-2} \\
C_{i2} = -1.3198825 \times 10^{-5} \\
C_{i3} = -0.49382577
\end{array}
\right.
\qquad (87)
$$

APPENDIX D. ANALYTICAL MODEL OF
TRAJECTORY

D.1. NON-OSCILLATORY MODE

Solving system (1) with analytic functions is a difficult
task as soon as *k* and *W*<sub>*a*</sub> are not constant. This is the
object of Sect. 3, including the underlying microphysics.
As an alternative, we consider here the non-oscillatory
motion of an ice parcel subject to entrainment by
horizontal wind *U*<sub>*a*</sub>, constant and uniform vertical updraft
*W*<sub>*a*</sub> and free fall *W*<sub>*f*</sub> governed by the set of ODE:

$$
\begin{cases}
\dot{x} = U_a \\
\dot{z} = W_f + W_a
\end{cases}
\qquad (88)
$$

with initial conditions (2). The single dot means first time-
derivative. The horizontal wind is an approximation of the
sheared flow in the range of altitude 8–12 km (Figure 1),
with shear ω, which is cast into the linear form:

$$
U_a = \omega(z - z_{a0}) + U_{a0} \quad (89)
$$

Therefore, the velocity $U_a$ changes sign at the altitude $z_m$ such that:

$$
z_m = z_{a0} - \frac{U_{a0}}{\omega} \qquad (90)
$$

Since the average crystal radius and mass decrease
in the virga under the action of sublimation, the free
fall speed W<sub>f</sub> must also decrease as t tends towards
infinity. We chose a unique analytical function W<sub>f</sub>(t)
defined in [0, +∞[, and satisfying the following conditions
(Figure 38):

$$
\left\{
\begin{array}{l}
W_f(0) = W_{f0} \\
W_f(t_1) = -W_{a0}
\end{array}
\right.
\qquad
\left\{
\begin{array}{l}
W_f(t_2) = -W_{fm} \\
\lim_{t \to +\infty} W_f(t) = 0
\end{array}
\right.
\qquad (91)
$$

<figure><img src="image_42.png" /><figcaption>**Figure 38** Qualitative time profile of analytical free fall speed.</figcaption></figure>

Looking for a solution of the type:

$$
W_f(t) = W_{f0} + Ae^{-kt} \quad (92)
$$

that is subject to the conditions (91), we obtain the
expression:

$$
W_f(t) = W_{f0} - (W_{a0} + W_{f0}) \frac{t}{t_1} e^{-k(t-t_1)} \quad (93)
$$

with the characteristic times $t_1$ and $t_2$ such that:

$$
t_2 = \frac{1}{k} \qquad (94)
$$

and:

$$
\begin{align}
\tau &= \frac{t_1}{t_2} \notag \\
r &= \frac{W_{a0} + W_{f0}}{W_{fm} - W_{f0}} \tag{95}
\end{align}
$$

The non-dimensional parameter $\tau$ is solution of the
classical transcendental Lambert-type equation (Corless
et al. 1996):

$$
\tau = -r e^{r-1} \tag{96}
$$

Substituting the expressions (89) and (93) in the dynamic
system (88), we obtain at any time $t \in [0, +\infty[$:

$$
\begin{cases}
\dot{x} = \omega(z - z_{a0}) + U_{a0} \\
\dot{z} = (W_{a0} + W_{f0}) \left\{ 1 - \frac{t}{t_1} e^{-k(t-t_1)} \right\}
\end{cases}
\qquad (97)
$$

The integration of the system of first order ODE (97)
subject to the following initial conditions:

$$
\begin{cases}
x = x_0 \\
z = z_0
\end{cases}
\quad
\begin{cases}
\dot{x} = 0 \\
\dot{z} = 0
\end{cases}
\qquad (98)
$$

yields the expressions for x and z:

$$
\left\{
\begin{array}{l}
x = x_0 + \{\omega(z - z_{a0}) + U_{a0}\} t \\
+ \omega(W_{a0} + W_{f0}) \left\{ \frac{t^2}{2} - t_2 t \frac{e^r}{\tau} (1 + e^{-kt}) + 2t_2^2 \frac{e^r}{\tau} (1 - e^{-kt}) \right\} \\
z = z_0 + (W_{a0} + W_{f0}) \left\{ t_2 \frac{e^r}{\tau} (e^{-kt} - 1) + t \left( 1 + \frac{e^r}{\tau} e^{-kt} \right) \right\}
\end{array}
\right.
\qquad (99)
$$

Eqs. (99) are the parametric equations of a
transcendental curve, more complex than the basic
parabolic shape (Marshall 1953; Hogan & Kew 2005).

We can notice that the behaviour at infinite
time, involving the secular terms in Eqs. (99), is not
relevant since the variable-mass crystal eventually
vanishes through sublimation and precipitation.
Moreover, our model considers a flat Earth, an
assumption which implies neglecting Coriolis force
and focusing on the behaviour at local scale (less than
100 km).
A maximum occurs where $\dot{z}=0$, at a time $t_M$ that is again a solution of the transcendental Lambert-type equation (Corless et al., 1996):

$$t e^{kt} = t_1 e^{kt_1} \qquad (100)$$

of which $t_M = t_1$ is the solution in the principal branch. Likewise, an extremum occurs where $\dot{x}=0$, at a time $t_m$ solution of a complicated transcendental equation. In contrast, the altitude $z_m$ of the extremum is simply given by Eq. (90).

An example of motion was sampled at 601 points over a 6000-second life time ($\approx 1.7$ hour) with a time step $\Delta t = 10$ s and the following values of parameters:

<table><tr><td>t<sub>1</sub> = 10<sup>3</sup>s = 16.7 min</td><td>W<sub>f0</sub> = -0.5 m s<sup>-1</sup></td><td>x<sub>0</sub> = 2 km</td></tr><tr><td></td><td>W<sub>f<sub>m</sub></sub> = -1.5 m s<sup>-1</sup></td><td>z<sub>0</sub> = 10 km</td></tr><tr><td>U<sub>a0</sub> = 10 m s<sup>-1</sup></td><td>W<sub>a0</sub> = 1.0 m s<sup>-1</sup></td><td>z<sub>a0</sub> = 8 km</td></tr></table>

We choose the following value of wind shear: $\omega \approx -5 \times 10^{-3} \text{ s}^{-1}$, consistently with published values (Harimaya 1968; Wada et al. 2005) which recommend: $|\omega| < 23 \text{ m s}^{-1} \text{ km}^{-1}$. Estimating the solutions of Eq. (96) and the zeros of the derivatives (97), we obtain the following characteristic times and passing times through extrema:

<table><tr><td>maximum</td><td>t<sub>M</sub> = t<sub>1</sub></td><td>x<sub>M</sub> = 1.22 km</td><td>z<sub>M</sub> = 10.23 km</td></tr><tr><td>extremum</td><td>t<sub>m</sub> = 37.0 min</td><td>x<sub>m</sub> = 0.31 km</td><td>z<sub>m</sub> = 10.00 km</td></tr><tr><td>critical time</td><td>t<sub>2</sub> = 71.9 min</td><td></td><td></td></tr></table>

The non-dimensional parameters $r$ and $\tau$ defined by relations (95) respectively take the values: $r = -0.5$ and $\tau = 0.232$. The trajectory (Figure 39a) has a nice hooked-shape head expanding vertically between 0 and $t_1$, like that obtained with the full model (Figure 7a). Nevertheless, the updraft $W_a$ is constant in the whole atmosphere in the present analytical model, while it is zero below $z_{w0}$ and nonzero above in the full model.

The theoretical hodograph (Figure 39b) is qualitatively similar to that obtained with the fully coupled model (dynamics/microphysics) (Figure 7b). Likewise, the

theoretical profile of free fall speed (Figure 40) is qualitatively similar to the assumed profile (Figure 38) and to that of the full model in non-oscillatory mode (Figure 9a).

Casting Eq. (10) as:

$$\dot{m} = \dot{m}_0 \Delta S \qquad (101)$$

with $\Delta S$ denoting the factor $S_i - 1 - R$, and integrating between $t = 0$ and $t > 0$ at approximately constant rate $\dot{m}_0$ yields:

$$m(t) \approx m_0 + \dot{m}_0 \Delta S t \qquad (102)$$

which yields an estimation of the crystal's life duration $t_{max}$ when $m(t_{max})$ vanishes:

$$t_{max} \approx -\frac{m_0}{\dot{m}_0 \Delta S} \qquad (103)$$

With the following typical values obtained from simulations of Sect 3.2: $m_0 \approx 10^{-9}$ kg, $\dot{m}_0 \approx 2 \times 10^{-12}$ kg s$^{-1}$, $\Delta S \approx -0.02$, we obtain: $t_{max} \approx 7$ hours. This elementary calculation shows that a crystal falling through highly saturated layers with a large variation of the growth factor ($|\Delta S| \gg 0.1$), combining the effect of ice saturation and radiative transfer, cannot be long-lived.

<figure><img src="image_43.png" /><figcaption>Figure 40 Time profile of free fall speed for analytic solution, Eq.(99).</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 39 Trajectory and hodograph for analytic solution, Eq.(99).</figcaption></figure>
## D.2. OSCILLATORY MODE

1) Our purpose here is to derive a theoretical ODE proving the existence of damped oscillations as solution, and tentatively expressions of angular frequency and damping ratio. We first note that the full system:

$$ egin{cases} \ddot{x} = -\frac{\dot{m}_1}{m}(\dot{x} - U_\alpha) \\ \ddot{z} = -\frac{\dot{m}_1}{m}(\dot{z} - W_\alpha) - g \\ \dot{m} = \dot{m}_0(S_i - 1 - R) \end{cases} \quad (104a,b,c) $$

is a system of second order ODE of Liénard-type of the general form:

$$ egin{cases} \ddot{x} + f_1(z)\dot{x} + g_1(z) = 0 \\ \ddot{z} + f_2(z)\dot{z} + g_2(z) = 0 \\ \dot{m} + f_3(z)m = 0 \end{cases} \quad (105) $$

that can be cast more explicitly into the form:

$$ egin{cases} u \frac{du}{dx} = -\frac{\dot{m}_1}{m}(u(z) - U_\alpha(z)) \\ w \frac{dw}{dz} = -\frac{\dot{m}_1}{m}(w(z) - W_\alpha(z)) - g \\ w \frac{dm}{dz} = \dot{m}_0(z)(S_i(z) - 1 - R(z)) \end{cases} \quad (106) $$

that is itself a set of coupled Abel equations of the second kind.

Now, neglecting the coupling with horizontal motion, we can thus write the equations of evolution of vertical motion and mass Eqs. (1b) and (10) in the form:

$$ egin{cases} \dot{m} = \dot{m}_0(S_i - 1 - R) \\ \ddot{z} = -\frac{\dot{m}_1}{m_0}(\dot{z} - W_\alpha) - \frac{m}{m_0}g \end{cases} \quad (107a,b) $$

where single and double dots respectively mean first and second time-derivative, and the mass growth rates are defined by:

$$ egin{cases} \dot{m}_0 = \frac{4\pi C_f v}{\frac{R_v T}{p_{vs} D_v} + \frac{L_s}{K_\alpha T} \left( \frac{L_s}{R_v T} - 1 \right)} \\ \dot{m}_1 = 16(1 + 0.078 \text{Re}^{0.945}) \mu_\alpha r_i \end{cases} \quad (108a,b) $$

Then expanding $S_i - 1$ to first order in the vicinity of the equilibrium height $z_\infty$:

$$ S_i - 1 - R \approx \sigma_5(z - z_\infty) \quad (109) $$

and noting that:

$$ \dot{m} = w \frac{dm}{dz} \quad (110) $$

we derive a relation for the vertical mass gradient:

$$ w \frac{dm}{dz} = \dot{m}_0 \sigma_5 (z - z_\infty) \quad (111) $$

Likewise, let us expand the velocity in the vicinity of $z_\infty$:

$$ w \approx \omega(z - z_\infty) \quad (112) $$

and upon inserting this expression into Eq. (111) we obtain a simple ODE for mass:

$$ \frac{dm}{dz} \approx \frac{\dot{m}_0 \sigma_5}{\omega} \quad (113) $$

which yields by integration over $z$ between $z_0$ and $z$:

$$ m \approx m_\infty \left\{ 1 + \frac{\dot{m}_0 \sigma_5}{m_\infty \omega} (z - z_\infty) \right\} \quad (114) $$

Likewise, let us expand the effective radius $r_i$ in Eq. (108) after Eq. (38):

$$ r_i \approx \frac{1}{2} \sqrt{\frac{6m_\infty}{\pi \rho_i}} \left\{ 1 + \frac{\dot{m}_0 \sigma_5}{m_\infty \omega} (z - z_\infty) \right\} \approx \frac{D_{i\infty}}{2} \left\{ 1 + \frac{\dot{m}_0 \sigma_5}{3m_\infty \omega} (z - z_\infty) \right\} \quad (115) $$

with the asymptotic equivalent diameter:

$$ D_{i\infty} = \sqrt[3]{\frac{6m_\infty}{\pi \rho_i}} \quad (116) $$

Substituting these expressions for mass and radius into Eqs. (107b) and (108b) and rearranging, we finally arrive at the fundamental ODE of a damped harmonic oscillator:

$$ \ddot{z} + \frac{\dot{m}_1}{m_0}\dot{z} + \frac{\dot{m}_0\sigma_5}{m_0\omega}\left\{g - \frac{8}{3}\frac{\mu_\alpha D_{i\infty}}{m_\infty}(1 + 0.078\text{Re}^{0.945})W_\alpha\right\}(z - z_\infty) + \frac{m_\infty}{m_0}g - 8\frac{\mu_\alpha D_{i\infty}}{m_0}(1 + 0.078\text{Re}^{0.945})W_\alpha = 0 \quad (117) $$

The second term represents the damping friction and the third one the restoring force producing the oscillation. Eq. (117) clearly shows that it is necessary that an updraft exists ($W_\alpha \neq 0$) below the generating height ($z < z_0$) for a long-period oscillation to take place.

Satisfying the conditions at initial time ($z = z_0$) and infinity ($z = z_\infty$), a solution of the homogeneous equation in the case of an underdamped oscillator can be written as (Meirovitch 1986):

$$ z(t) = z_0 + (z_\infty - z_0) \left\{ \exp(-\zeta \Omega_0 t) \frac{\sin(\sqrt{1-\zeta^2} \Omega_0 t + \Phi)}{\sin\Phi} - 1 \right\} \quad (118) $$

with phase $\Phi \neq 0$. The angular frequency $\Omega_0$ and the damping ratio $\zeta$ are such that:

$$ \left\{ \begin{aligned} \Omega_0^2 &= \frac{\dot{m}_0}{m_0} \frac{\sigma_5}{\omega} \left\{ g - \frac{8}{3} \frac{\mu_\alpha D_{i\infty}}{m_\infty} (1 + 0.078 \text{Re}^{0.945}) W_\alpha \right\} \\ 2\zeta \Omega_0 &= \frac{\dot{m}_1}{m_0} \end{aligned} \right. \quad (119) $$

After substitution of mass growth rates, Eqs. (108), they can be written:
<figure><img src="image_45.png" /><figcaption>Figure 41 Altitude time profile of theoretical underdamped harmonic oscillator, Eq.(118) for moderate (a) and strong (b) damping.</figcaption></figure>

$$ 
\left\{ 
\begin{aligned}
\Omega_0^2 &= \frac{\sigma_s}{\omega m_0} \frac{4\pi C f_v}{\frac{R_v T}{p_{vsi} D_v} + \frac{L_s}{K_a T} \left( \frac{L_s}{R_v T} - 1 \right)} \\
&\quad \times \left\{ g - \frac{8}{3} \frac{\mu_0 D_{vsi}}{m_\infty} (1 + 0.078 \text{Re}^{0.945}) W_\infty \right\} \\
\zeta \Omega_0 &= 8(1 + 0.078 \text{Re}^{0.945}) \frac{\mu_0 r_f}{m_0}
\end{aligned}
\right. 
\quad (120) 
$$

Using the following values derived from the simulations of Sect. 3.2: $m_\infty \approx 10^{-9}$ kg, $\mu_a \approx 10^{-5}$ Pa s, $C \approx 10^{-4}$ m, $r_f \approx 10^{-4}$ m, $f_v \approx 1$, $\sigma_s \approx 7 \times 10^{-5}$ m$^{-1}$, $\omega \approx 5 \times 10^{-5}$ s$^{-1}$, we can estimate $\Omega_0$ and $\zeta$ as:

$$ 
\left\{ 
\begin{aligned}
\Omega_0 &= \sqrt{\frac{7}{5} \frac{12 \times 10^{-4}}{10^{-9} \times 10^9}} \left( 9.8 - 2.7 \frac{1.4 \times 10^{-5} \times 1.3 \times 10^{-4} \times 1.2}{10^{-9}} \right) 0.6 \\
&\approx 0.01 \text{ rad s}^{-1} \\
\zeta &= \frac{4 \times 1.2 \times 1.4 \times 10^{-5} \times 1.3 \times 10^{-4}}{10^{-9} \times 10^{-2}} \approx 8.7 \times 10^2
\end{aligned}
\right. 
\quad (121) 
$$

The angular frequency thus obtained is larger than the effective value obtained in Sect. 3.2.3 and 5.2, namely $\Omega_0 \approx 1.7 \times 10^{-4}$ rad s$^{-1}$, and the damping ratio, compared to $\zeta = 0.05$, has also too large a value, that would characterize an overdamped oscillation. This large discrepancy is probably due to the fact that most variables of the problem are dependent of temperature, itself being altitude-dependent, especially in the mass equation (Eq. (10)), in the mass rates $\dot{m}_0$ and $\dot{m}_1$ (Eq. (108)) we assume constant. A full treatment would necessitate an expansion of all these variables as functions of temperature, and eventually use of Eqs. (75).

In order to illustrate the original behaviour we found in Sections 3 and 5, we plot time profiles of altitude as modelled by Eq. (118) with constants associated to “moderately” and “strongly” damped harmonic oscillators (Table 3). The resulting profiles (Figure 41) are visually matched with the corresponding profiles of numerical solutions, displayed respectively in Sect. 3.2.3 (Figure 17b) and Sect. 5.2 (Figure 32c) for moderate

<table><caption>Table 3 Parameters of theoretical periodic profiles.</caption><thead><tr><th>DAMPING</th><th>MODERATE</th><th>STRONG</th></tr></thead><tbody><tr><td>z<sub>0</sub> (km)</td><td>10.0</td><td>10.0</td></tr><tr><td>z<sub>∞</sub> (km)</td><td>9.63</td><td>8.58</td></tr><tr><td>P (hour)</td><td>10.2</td><td>7.8</td></tr><tr><td>Ω<sub>0</sub> (rad s<sup>-1</sup>)</td><td>1.71 × 10<sup>-4</sup></td><td>2.24 × 10<sup>-4</sup></td></tr><tr><td>ζ</td><td>0.05</td><td>0.07</td></tr><tr><td>Φ (rad)</td><td>π/4</td><td>π/4</td></tr></tbody></table>

damping and Sect. 5.3 (Figure 33c) for strong damping. The agreement is very good except at the beginning of the motion because the time spent by the crystal in the hooked cap broadens the first period ($t < 5$ hours), so that the subsequent periods are slightly shifted to the right in the actual motion. Inharmonicity detected in numerical simulations probably enhances that effect and suggests to analyse it further with the method of Krylov and Bogoliubov (Bose et al. 1989).

2) This short subsection is devoted to the explanation of foldings produced in the hodograph and horizontal velocity profiles of Cases 2 (Sect. 5.3) and 3 (Sect. 5.4). Let us recast the ODE governing horizontal velocity, Eq. (1a), as:

$$ \dot{u} = -k(u - U_a) \quad (122) $$

Since $k$ is positive, Eq. (122) shows that when $u > U_a$, then $\dot{u} < 0$ and $u$ decreases to $U_a$ in a drag relaxation time $\tau_D = 1/k$ (Paoli & Shariff 2016), and reciprocally when $u < U_a$. Since $\tau_D$ is of order of time step $\Delta t$, the relaxation is quasi instantaneous and therefore $u$ follows the variation of $U_a$. Numerical experiments in Sect. 5 reveal that the difference $u-U_a$ is of order $10^{-4}$ m s$^{-1}$, so that with $k \approx 10^{-1}$, the crystal undergoes horizontal accelerations of order $10^{-3}$ m s$^{-2}$. In the vicinity of the maximum wind speed $U_{amin}$, $u$ is therefore constrained to remain smaller than $U_{amin}$, and thus $u$ profiles (Figures 34 and 36) show two spikes per period, corresponding to the ascending and descending branches of $U_a$ profile in the
vicinity of $U_{amax}$ = 10 m s⁻¹ at the altitude of 8 km (Figure 1). The folded loops in hodographs are a result of that phenomenon.

## ACKNOWLEDGEMENTS

The author is very grateful to the reviewer for important and insightful comments which contributed to greatly improve and enlarge the content and discussions of the article.

## COMPETING INTERESTS

The author has no competing interests to declare.

## AUTHOR AFFILIATION

Roland P. H. Berton (1) orcid.org/0000-0002-7363-6363
Office National d'Études et de Recherches Aérospatiales
(ONERA)/DOTA, Université Paris Saclay, 91123 Palaiseau, France

## REFERENCES

Aguilar, B, Trontin, P, Reitter, L, Köbschall, K, Dezitter, F, Roisman, I and Villedieu, P. 2022. Ice crystal drag model extension to snowflakes: Experimental and numerical investigations. AIAA J., 60(12): 6633-6645. DOI: https://doi.org/10.2514/1.J062122

Alduchov, OA and Eskridge, RE. 1996. Improved Magnus form approximation of saturation vapor pressure. *J. Appl. Meteorol.*, 35: 601-609. DOI: https://doi.org/10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2

Barahona, D and Nenes, A. 2008. Parameterization of cirrus cloud formation in large-scale models: Homogeneous nucleation. *J. Geophys. Res.*, 113: D11211. DOI: https://doi.org/10.1029/2007JD009355

Baum, BA, Kratz, DP, Yang, P, Ou, SC, Hu, Y, Soulen, PF and Tsay, S-C. 2000. Remote sensing of cloud properties using MODIS airborne simulator imagery during SUCCESS: 1. Data and models. *J. Geophys. Res.*, 105(D9): 11767-11780. DOI: https://doi.org/10.1029/1999JD901089

Berton, RPH. 2008. Constrained water cloud generator. *Atmos. Res.*, 89(1-2): 12-29. DOI: https://doi.org/10.1016/j.atmosres.2007.09.007

Bohren, CF and Albrecht, BA. 1998. Atmospheric thermodynamics. Oxford: Oxford University Press.

Bose, SK, Dubey, UB and Varma, N. 1989. Study of certain aspects of anharmonic, time-dependent and damped harmonic oscillator systems. *Fortschr. Phys.*, 37(10): 761-818. DOI: https://doi.org/10.1002/prop.2190371002

Buschmann, N, McFarquhar, GM and Heymsfield, AJ. 2022. Effects of observed horizontal inhomogeneities within cirrus clouds on solar radiative transfer. *J. Geophys. Res.*, 107(D20): 4445. DOI: https://doi.org/10.1029/2001JD001273

Chen, C-J and Wang, PK. 2009. Diffusion growth of solid and hollow hexagonal ice columns. *Il Nuovo Cimento*, 124: 87-97. DOI: https://doi.org/10.1393/ncb/i2009-10746-2

Chiruta, M and Wang, PK. 2003. The capacitance of rosette ice crystals. *J. Atmos. Sci.*, 60(6): 836-846. DOI: https://doi.org/10.1175/1520-0469(2003)060<0836:TCORIC>2.0.CO;2

Chiruta, M and Wang, PK. 2005. The capacitance of solid and hollow hexagonal ice column. *Geophys. Res. Lett.*, 32: L05803. DOI: https://doi.org/10.1029/2004GL021771

Cianciolo, ME. 1993. Cumulus cloud scene simulation modeling using fractals and physics. *Proc. SPIE*, 1819: 70-77. DOI: https://doi.org/10.1117/12.142209

Cohen, C, Darbois-Texier, B, Dupeux, G, Brunel, E, Quéré, D and Clanet, C. 2014. The aerodynamic wall. *Proc. R. Soc. A.*, 470: 20130497. 1-20. DOI: https://doi.org/10.1098/rspa.2013.0497

Corless, RM, Gonnet, GH, Hare, DEG, Jeffrey, DJ and Knuth, DE. 1996. On the Lambert W function. *Adv. Comput. Math.*, 5: 329-359. DOI: https://doi.org/10.1007/BF02124750

Cotton, RJ, Field, PR, Ulanowski, Z, Kaye, PH, Hirst, E, Greenaway, RS, Crawford, I, Crosier, J and Dorsey, J. 2013. The effective density of small ice particles obtained from in situ aircraft observations of mid-latitude cirrus. *Q. J. R. Meteorol. Soc.*, 139: 1923-1934. DOI: https://doi.org/10.1002/qj.2058

Cziczo, DJ, Froyd, KD, Hoose, C, Jensen, EJ, Diao, M, Zondlo, MA, Smith, JB, Twohy, CH and Murphy, DM. 2013. Clarifying the dominant sources and mechanisms of cirrus cloud formation. *Science*, 340(6138): 1320-1324. DOI: https://doi.org/10.1126/science.1234145

Demoz, BB, Starr, DO'C, Chan, KR and Bowen, SW. 1998. Wavelet analysis of dynamical processes in cirrus. *Geophys. Res. Lett.*, 25(9): 1347-1350. DOI: https://doi.org/10.1029/97GL03226

Detwiler, A and Ramaswamy, V. 1990. Radiative heating profiles in simple cirrus cloud systems. *J. Atmos. Sci.*, 47(18): 2167-2176. DOI: https://doi.org/10.1175/1520-0469(1990)047<2167:RHPISC>2.0.CO;2

Dobbie, S and Jonas, P. 2001. Radiative influences on the structure and lifetime of cirrus clouds. *Q. J. R. Meteorol. Soc.* B., 127(578): 2663-2682. DOI: https://doi.org/10.1002/qj.49712757808

Durran, DR and Klemp, JB. 1982. On the effects of moisture on the Brunt-Väisälä frequency. *J. Atmos. Sci.*, 39(10): 2152-2158. DOI: https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2

Eguchi, N, Yokota, T and Inoue, G. 2007. Characteristics of cirrus clouds from ICESat/GLAS observations. *Geophys. Res. Lett.*, 34: L09810. DOI: https://doi.org/10.1029/2007GL029529

Evans, KF and Wiscombe, WJ. 2004. An algorithm for generating stochastic cloud fields from radar profile statistics. *Atmos. Res.*, 72(1-4): 263-289. DOI: https://doi.org/10.1016/j.atmosres.2004.03.016
Fauchez, T, Cornet, C, Szczap, F, Dubuisson, P and Rosambert, T. 2014. Impact of cirrus clouds heterogeneities on top-of-atmosphere thermal infrared radiation. *Atmos. Chem. Phys.*, 14: 5599-5615. DOI: https://doi.org/10.5194/acp-14-5599-2014

Feddes, RG. 1974. A synoptic-scale model for simulating condensed atmospheric moisture. USAFETAC TN-74-4. DOI: https://doi.org/10.21236/ADA002118

Fraser, AB and Bohren, CF. 1992. Is virga rain that evaporates before reaching the ground? *Mon. Wea. Rev.*, 120: 1565-1571. DOI: https://doi.org/10.1175/1520-0493(1992)120<1565:IVRTEB>2.0.CO;2

Fridlind, AM, Atlas, R, van Diedenhoven, B, Um, J, McFarquhar, GM, Ackerman, AS, Moyer, EJ and Lawson, R. 2016. Derivation of physical and optical properties of mid-latitude cirrus ice crystals for a size-resolved cloud microphysics model. *Atmos. Chem. Phys.*, 16: 7251-7283. DOI: https://doi.org/10.5194/acp-16-7251-2016

Fukuta, N and Takahashi, T. 1999. The growth of atmospheric ice crystals: a summary of findings in vertical supercooled cloud tunnel studies. *J. Atmos. Sci.*, 56(12): 1963-1979. DOI: https://doi.org/10.1175/1520-0469(1999)056<1963:TGOAIC>2.0.CO;2

Fusina, F and Spichtinger, P. 2010. Cirrus clouds triggered by radiation, a multiscale phenomenon. *Atmos. Chem. Phys.*, 10: 5179-5190. DOI: https://doi.org/10.5194/acp-10-5179-2010

Fusina, F, Spichtinger, P and Lohmann, U. 2007. Impact of ice supersaturated regions and thin cirrus on radiation in the midlatitudes. *J. Geophys. Res.*, 112: D24S14. DOI: https://doi.org/10.1029/2007JD008449

Gierens, KM. 1994. The influence of radiation on the diffusional growth of ice crystals. *Beitr. Phys. Atmos.*, 67: 181-193. https://elib.dlr.de/31864/

Grabowski, WW. 1993. Cumulus entrainment, fine-scale mixing, and buoyancy reversal. *Q. J. R. Meteorol. Soc.*, 119: 935-956. DOI: https://doi.org/10.1002/qj.49711951305

Gu, Y and Liou, KN. 2000. Interactions of radiation, microphysics, and turbulence in the evolution of cirrus clouds. *J. Atmos. Sci.*, 57(15): 2463-2479. DOI: https://doi.org/10.1175/1520-0469(2000)057<2463:IORMAT>2.0.CO;2

Gultepe, I, Heymsfield, AJ and Lenschow, DH. 1990. A comparison of vertical velocity in cirrus obtained from aircraft and Lidar divergence measurements during FIRE. *J. Appl. Meteorol. Clim.*, 7(1): 58-67. DOI: https://doi.org/10.1175/1520-0426(1990)007<0058:ACOVVI>2.0.CO;2

Haider, A and Levenspiel, O. 1989. Drag coefficient and terminal velocity of spherical and nonspherical particles. *Powder Tech.*, 58(1): 63-70. DOI: https://doi.org/10.1016/0032-5910(89)80008-7

Hall, WD and Pruppacher, HR. 1976. The survival of ice particles falling from cirrus clouds in subsaturated air. *J. Atmos. Sci.*, 33(10): 1995-2006. DOI: https://doi.org/10.1175/1520-0469(1976)033<1995:TSOIPF>2.0.CO;2

Harimaya, T. 1968. On the shape of cirrus uncinus clouds: a numerical computation – Studies of cirrus clouds: Part III. *J. Meteorol. Soc. Japan*, 46(4): 272-279. DOI: https://doi.org/10.2151/jmsj1965.46.4_272

Hashino, T, Cheng, KY, Chueh, C-C and Wang, PK. 2016. Numerical study of motion and stability of falling columnar crystals. *J. Atmos. Sci.*, 73(5): 1923-1942. DOI: https://doi.org/10.1175/JAS-D-15-0219.1

Hashino, T, Chiruta, M, Polzin, D, Kubicek, A and Wang, PK. 2014. Numerical simulation of the flow fields around falling ice crystals with inclined orientation and the hydrodynamic torque. *Atmos. Res.*, 150: 79-96. DOI: https://doi.org/10.1016/j.atmosres.2014.07.003

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

Heymsfield, A. 1973. Cirrus uncinus generating cells and the evolution of cirriform clouds. PhD thesis, the University of Chicago, 286 pp.

Heymsfield, A. 1975a. Cirrus uncinus generating cells and the evolution of cirriform clouds. Part I: Aircraft observations of the growth of the ice phase. *J. Atmos. Sci.*, 32(4): 799-808. DOI: https://doi.org/10.1175/1520-0469(1975)032<0799:CUGCAT>2.0.CO;2

Heymsfield, A. 1975b. Cirrus uncinus generating cells and the evolution of cirriform clouds. Part II: The structure and circulations of the cirrus uncinus generating head. *J. Atmos. Sci.*, 32(4): 809-819. DOI: https://doi.org/10.1175/1520-0469(1975)032<0809:CUGCAT>2.0.CO;2

Heymsfield, A. 1975c. Cirrus uncinus generating cells and the evolution of cirriform clouds. Part III: Numerical computations of the growth of the ice phase. *J. Atmos. Sci.*, 32(4): 820-830. DOI: https://doi.org/10.1175/1520-0469(1975)032<0820:CUGCAT>2.0.CO;2

Heymsfield, AJ, Bansemer, A, Schmitt, C, Twohy, C and Poellot, MR. 2004. Effective ice particle densities derived from aircraft data. *J. Atmos. Sci.*, 61(9): 982-1003. DOI: https://doi.org/10.1175/1520-0469(2004)061<0982:EIPD>2.0.CO;2

Heymsfield, AJ and Iaquinta, J. 2000. Cirrus crystal terminal velocities. *J. Atmos. Sci.*, 57(7): 916-938. DOI: https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2

Heymsfield, AJ, Lewis, S, Bansemer, A, Iaquinta, J, Miloshevich, LM, Kajikawa, M, Twohy, C and Poellot, MR. 2002. A general approach for deriving the properties of cirrus and stratiform ice cloud particles. *J. Atmos. Sci.*, 59(1): 3-29. DOI: https://doi.org/10.1175/1520-0469(2002)059<0003:AGAFDT>2.0.CO;2

Heymsfield, AJ, Schmitt, C and Bansemer, A. 2013. Ice cloud particle size distributions and pressure-dependent terminal velocities from in situ observations at temperatures from 0° to -86°C. *J. Atmos. Sci.*, 70(12): 4123-4154. DOI: https://doi.org/10.1175/JAS-D-12-0124.1

Heymsfield, AJ and Westbrook, CD. 2010. Advances in the estimation of ice particle fall speeds using laboratory and
field measurements. *J. Atmos. Sci.*, 67(8): 2469-2482. DOI:
https://doi.org/10.1175/2010JAS3379.1

Hogan, RJ and Kew, SF. 2005. A 3D stochastic cloud model for investigating the radiative properties of inhomogeneous cirrus clouds. *Q. J. R. Meteorol. Soc.* A, 131(601): 2585-2608. DOI: https://doi.org/10.1256/qj.04.144

Houghton, HG. 1950. A preliminary quantitative analysis of precipitation mechanisms. *J. Meteorol.*, 7(6): 363-369. DOI: https://doi.org/10.1175/1520-0469(1950)007<0363:APQA OP>2.0.CO;2

Igel, MR and Biello, JA. 2020. The nontraditional Coriolis terms and tropical convective clouds. *J. Atmos. Sci.*, 77(12): 3985-3998. DOI: https://doi.org/10.1175/JAS-D-20-0024.1

Inoue, T. 1985. On the temperature and effective emissivity determination of semi-transparent cirrus clouds by bi-spectral measurements in the 10 µm window region. *J. Meteorol. Soc. Japan*, 63(1): 88-99. DOI: https://doi.org/10.2151/jmsj1965.63.1_88

Jensen, AA and Harrington, JY. 2015. Modeling ice crystal aspect ratio evolution during riming: a single-particle growth model. *J. Atmos. Sci.*, 72(7): 2569-2590. DOI: https://doi.org/10.1175/JAS-D-14-0297.1

Jensen, EJ, van den Heever, SC and Grant, LD. 2018. The life cycles of ice crystals detrained from the tops of deep convection. *J. Geophys. Res. Atmospheres*, 123: 9624-9634. DOI: https://doi.org/10.1029/2018JD028832

Ji, W and Wang, PK. 1999. Ventilation coefficients for falling ice crystals in the atmosphere at low-intermediate Reynolds numbers. *J. Atmos. Sci.*, 56(6): 829-836. DOI: https://doi.org/10.1175/1520-0469(1999)056<0829:VCFFIC>2.0.CO;2

Kärcher, B, Jensen, EJ and Lohmann, U. 2019. The impact of mesoscale gravity waves on homogeneous ice nucleation in cirrus clouds. *Geophys. Res. Let.*, 46(10): 5556-5565. DOI: https://doi.org/10.1029/2019GL082437

Kew, S. 2003. Development of a 3D fractal cirrus model and its use in investigating the impact of cirrus inhomogeneity on radiation. Unpublished thesis (Master Sci.). University of Reading (UK).

Khvorostyanov, VI and Curry, JA. 2008. Analytical solutions to the stochastic kinetic equation for liquid and ice particle size spectra. Part I: small-size fraction. *J. Atmos. Sci.*, 65(7): 2025-2043. DOI: https://doi.org/10.1175/2007JAS2484.1

Kobayashi, T. 1961. The growth of snow crystals at low supersaturations. *Phil. Mag.*, 6(71): 1363-1370. DOI: https://doi.org/10.1080/14786436108241231

Köhler, CG and Seifert, A. 2015. Identifying sensitivities for cirrus modelling using a two-moment two-mode bulk microphysics scheme. *Tellus B*, 67: 24494. DOI: https://doi.org/10.3402/tellusb.v67.24494

Kokhanovsky, AA. 2003. The influence of horizontal inhomogeneity on radiative characteristics of clouds: an asymptotic case study. *IEEE Trans. Geosci. Rem. Sens.*, 41(4): 817-825. DOI: https://doi.org/10.1109/TGRS.2003.811075

Koren, I and Feingold, G. 2011. Aerosol-cloud-precipitation system as a predator-prey problem. *PNAS*, 108(30):

12227-12232. DOI: https://doi.org/10.1073/pnas.1101777108

Korolev, A and Isaac, GI. 2006. Relative humidity in liquid, mixed-phase and ice clouds. *J. Atmos. Sci.*, 63(11): 2865-2880. DOI: https://doi.org/10.1175/JAS3784.1

Krämer, M, Schiller, C, Afchine, A, Bauer, R, Gensch, I, Mangold, A, Schlicht, S, Spelten, N, Sitnikov, N, Borrmann, S, de Reus, M and Spichtinger, P. 2009. Ice supersaturations and cirrus cloud crystal numbers. *Atmos. Chem. Phys.*, 9: 3505-3522. DOI: https://doi.org/10.5194/acp-9-3505-2009

Kuhn, T and Heymsfield, AJ. 2016. In situ balloon-borne ice particle imaging in high-latitude cirrus. *Pure Appl. Geophys.*, 173: 3065-3084. DOI: https://doi.org/10.1007/s00024-016-1324-x

Lawson, RP, Woods, S, Jensen, E, Erfani, E, Gurganus, C, Gallagher, M, Connolly, P, Whiteway, J, Baran, AJ, May, P, Heymsfield, A, Schmitt, CG, McFarquhar, G, Um, J, Protat, A, Bailey, M, Lance, S, Muehlbauer, A, Stith, J, Korolev, A, Toon, OB and Krämer, M. 2019. A review of ice particle shapes in cirrus formed in situ and in anvils. *J. Geophys. Res. Atmospheres*, 124: 10,049-10,090. DOI: https://doi.org/10.1029/2018JD030122

Lebo, ZJ, Johnson, NC and Harrington, JY. 2008. Radiative influences on ice crystal and droplet growth within mixed-phase stratus clouds. *J. Geophys. Res.*, 113: D09203. DOI: https://doi.org/10.1029/2007JD009262

Lin, H, Noone, KJ, Ström, J and Heymsfield, AJ. 1998. Dynamical influences on cirrus cloud formation process. *J. Atmos. Sci.*, 55: 1940-1949. DOI: https://doi.org/10.1175/1520-0469(1998)055<1940:DIOCCF>2.0.CO;2

Liou, KN. 2002. An introduction to atmospheric radiation. Second Edition. Academic Press.

Liou, KN and Rao, N. 1996. Radiative transfer in cirrus clouds. Part IV: On cloud geometry, inhomogeneity, and absorption. *J. Atmos. Sci.*, 53(21): 3046-3065. DOI: https://doi.org/10.1175/1520-0469(1996)053<3046:RTICCP>2.0.CO;2

Liou, KN, Takano, Y, Ou, SC, Heymsfield, A and Kreiss, W. 1990. Infrared transmission through cirrus clouds: a radiative model for target detection. *Appl. Opt.*, 29(13): 1886-1896. DOI: https://doi.org/10.1364/AO.29.001886

Liu, H-C, Wang, PK and Schlesinger, RE. 2003a. A numerical study of cirrus clouds. Part I: Model description. *J. Atmos. Sci.*, 60(8): 1075-1084. DOI: https://doi.org/10.1175/1520-0469(2003)60<1075:ANSOCC>2.0.CO;2

Liu, H-C, Wang, PK and Schlesinger, RE. 2003b. A numerical study of cirrus clouds. Part II: Effects of ambient temperature, stability, radiation, ice microphysics, and microdynamics on cirrus evolution. *J. Atmos. Sci.*, 60(9): 1097-1119. DOI: https://doi.org/10.1175/1520-0469(2003)060<1097:ANSOCC>2.0.CO;2

Lohmann, U, Spichtinger, P, Jess, S, Peter, T and Smit, H. 2008. Cirrus cloud formation and ice supersaturated regions in a global climate model. *Environ. Res. Lett.*, 3: 045022. DOI: https://doi.org/10.1088/1748-9326/3/4/045022
Lorenz, EN. 1963. Deterministic nonperiodic flow. *J. Atm. Sci.*, 20(2): 130-141. DOI: https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

Ludlam, FH. 1948. The forms of ice-clouds. *Q. J. R. Meteorol. Soc.*, 74: 39-56. DOI: https://doi.org/10.1002/qj.49707431905

Luo, Z and Rossow, WB. 2004. Characterizing tropical cirrus life cycle, evolution, and interaction with upper-tropospheric water vapor using Lagrangian trajectory analysis of satellite observations. *J. Clim.*, 17: 4541-4563. DOI: https://doi.org/10.1175/3222.1

Lynch, DK, Sassen, K, Starr, DO'C and Stephens, G. 2002. *Cirrus*. Oxford: Oxford University Press. DOI: https://doi.org/10.1093/oso/9780195130720.001.0001

Mace, GG, Zhang, Y, Platnick, S, King, MD, Minnis, P and Yang, P. 2005. Evaluation of cirrus cloud properties derived from MODIS data using cloud properties derived from ground-based observations collected at the ARM SGP site. *J. Appl. Meteor.*, 44(2): 221-240. DOI: https://doi.org/10.1175/JAM2193.1

Magono, C, Harimaya, T and Yagi, T. 1967. Observation of cirrus clouds by trigonometric method – Studies of cirrus clouds: part I. *Geoph. Bull. Hokkaido Univ.*, 17: 83-98.

Marshall, JS. 1953. Precipitation trajectories and patterns. *J. Meteorol.*, 10: 25-29. DOI: https://doi.org/10.1175/1520-0469(1953)010<0025:PTAP>2.0.CO;2

Mascio, J, Xu, Z and Mace, GG. 2017. The mass-dimensional properties of cirrus clouds during TC4. *J. Geophys. Res. Atmospheres.*, 122: 10,402-10,417. DOI: https://doi.org/10.1002/2017JD02678

Mason, BJ. 1953. The growth of ice crystals in a supercooled water cloud. *Q. J. R. Meteorol. Soc.*, 79(339): 104-111. DOI: https://doi.org/10.1002/qj.49707933909

McFarquhar, GM and Heymsfield, AJ. 1998. The definition and significance of an effective radius for ice clouds. *J. Atm. Sci.*, 55(11): 2039-2052. DOI: https://doi.org/10.1175/1520-0469(1998)055<2039:TDASOA>2.0.CO;2

Meirovitch, L. 1986. Elements of vibration analysis. McGraw-Hill.

Miller, TL and Young, KC. 1979. A numerical simulation of ice crystal growth from the vapor phase. *J. Atm. Sci.*, 36(3): 458-469. DOI: https://doi.org/10.1175/1520-0469(1979)036<0458:ANSOIC>2.0.CO;2

Mishra, S, Mitchell, DL, Turner, DD and Lawson, RP. 2014. Parameterization of ice fall speeds in midlatitude cirrus: Results from SPartICus. *J. Geophys. Res. Atmos.*, 119: 3857-3876. DOI: https://doi.org/10.1002/2013JD020602

Murphy, DM and Koop, T. 2005. Review of the vapour pressures of ice and supercooled water for atmospheric applications. *Q. J. R. Meteorol. Soc.*, 131: 1539-1565. DOI: https://doi.org/10.1256/qj.04.94

Nazaryan, H, McCormick, MP and Menzel, WP. 2008. Global characterization of cirrus clouds using CALIPSO data. *J. Geophys. Res.*, 113: D16211. DOI: https://doi.org/10.1029/2007JD009481

Nelson, JT and Baker, MB. 1996. New theoretical framework for studies of vapor growth and sublimation of small ice crystals in the atmosphere. *J. Geophys. Res.*, 101(D3): 7033-7047. DOI: https://doi.org/10.1029/95JD03162

Nettesheim, JJ. 2017. A numerical study on the aerodynamics of freely falling planar ice crystals. Unpublished thesis (Master Sci.). Madison, USA: University of Wisconsin.

Orlanski, I. 1975. A rational subdivision of scales for atmospheric processes. *Bul. Am. Meteorol. Soc.*, 56(5): 527-534. DOI: https://doi.org/10.1175/1520-0477-56.5.527

Ovarlez, J, Gayet, J-F, Gierens, K, Ström, J, Ovarlez, H, Auriol, F, Busen, R and Schumann, U. 2002. Water vapour measurements inside cirrus clouds in Northern and Southern hemispheres during INCA. *Geophys. Res. Let.*, 29(16): 1813. 1-4. DOI: https://doi.org/10.1029/2001GL014440

Palmer, AJ. 1996. A spectral model for turbulence and microphysics dynamics in an ice cloud. *Nonlin. Proc. Geophys.*, 3: 23-28. DOI: https://doi.org/10.5194/npg-3-23-1996

Paltridge, GW and Platt, CMR. 1981. Aircraft measurements of solar and infrared radiation and the microphysics of cirrus cloud. *Q. J. R. Meteorol. Soc.*, 107: 367-380. DOI: https://doi.org/10.1002/qj.49710745208

Paoli, R and Shariff, K. 2016. Contrail modeling and simulation. *Ann. Rev. Fluid Mech.*, 48: 393-427. DOI: https://doi.org/10.1146/annurev-fluid-010814-013619

Perko, L. 1990. Differential equations and dynamical systems. *Texts in Applied Mathematics*. Springer. ISBN: 0387974431, 9780387974439, 3540974431. DOI: https://doi.org/10.1007/978-1-4684-0392-3

Peter, Th and Baker, M. 1996. Lifetimes of ice crystals in the upper troposphere and stratosphere. In: Crutzen, PJ and Ramanathan, V (eds.), *Clouds, Chemistry and Climate. NATO ASI Series*, 35: 57-82. Springer, Berlin, Heidelberg. DOI: https://doi.org/10.1007/978-3-642-61051-6_4

Philipona, R, Kräuchi, A and Brocard, E. 2012. Solar and thermal radiation profiles and radiative forcing measured through the atmosphere. *Geophys. Res. Let.*, 39: L13806. DOI: https://doi.org/10.1029/2012GL052087

Picard, A, Davis, RS, Gläser, M and Fujii, K. 2008. Revised formula for the density of moist air (CIPM-2007). *Metrologia*, 45: 149-155. DOI: https://doi.org/10.1088/0026-1394/45/2/004

Plastino, AR and Muzzio, JC. 1992. On the use of and abuse of Newton's second law for variable mass problems. *Cel. Mech. Dyn. Astron.*, 53: 227-232. DOI: https://doi.org/10.1007/BF00052611

Podglajen, A, Plougonven, R, Hertzog, A and Eric Jensen, E. 2018. Impact of gravity waves on the motion and distribution of atmospheric ice particles. *Atmos. Chem. Phys.*, 18: 10799-10823. DOI: https://doi.org/10.5194/acp-18-10799-2018

Pounder, ER. 1965. The physics of ice. Pergamon Press Ltd.

Prasad, AA, Sherwood, S, Reeder, MJ and Lane, TP. 2019. Rapidly evolving cirrus clouds modulated by
convectively generated gravity waves. *J. Geophys. Res.: Atmospheres*, 124: 7327-7338. DOI: https://doi.org/10.1029/2019JD030538

Pruppacher, HR and Klett, JD. 1978, 1997. *Microphysics of clouds and precipitation*. Dordrecht: Kluwer Academic Publishers. ISBN 978-0-306-48100-0.

Ramaswamy, V and Detwiler, A. 1986. Interdependence of radiation and microphysics in cirrus clouds. *J. Atmos. Sci.*, 43(21): 2289-2301. DOI: https://doi.org/10.1175/1520-0469(1986)043<2289:IORAMI>2.0.CO;2

Rikitake, T. 1958. Oscillations of a system of disk dynamos. *Math. Proc. Camb. Phil. Soc.*, 54(1): 89-105. DOI: https://doi.org/10.1017/S0305004100033223

Roach, WT. 1976. On the effect of radiative exchange on the growth by condensation of a cloud or fog droplet. *Q. J. R. Meteorol. Soc.*, 102: 361-372. DOI: https://doi.org/10.1002/qj.49710243207

Rogers, RR. 1975. An elementary parcel model with explicit condensation and supersaturation. *Atmosphere*, 13(4): 192-204. DOI: https://doi.org/10.1080/00046973.1975.9648397

Rogers, RR. 1979. A short course in cloud physics. New York: Pergamon Press.

Saltzman, B. 1962. Finite amplitude free convection as an initial value problem—I. *J. Atmos. Sci.*, 19(4): 329-341. DOI: https://doi.org/10.1175/1520-0469(1962)019<0329:FAFC AA>2.0.CO;2

Sassen, K and Krueger, SK. 1993. Toward an empirical definition of virga: Comments on “Is virga rain that evaporates before reaching the ground?” *Mon. Wea. Rev.*, 121: 2426-2428. DOI: https://doi.org/10.1175/1520-0493(1993)121<2426:TAEDOV>2.0.CO;2

Sassen, K, Starr, DO'C and Uttal, T. 1989. Mesoscale and microscale structure of cirrus clouds: Three case studies. *J. Atmos. Sci.*, 46(3): 371-396. DOI: https://doi.org/10.1175/1520-0469(1989)046<0371:MAMSOC>2.0.CO;2

Schmitt, CG and Heymsfield, AJ. 2007. On the occurrence of hollow bullet rosette- and column-shaped ice crystals in midlatitude cirrus. *J. Atmos. Sci.*, 64(12): 4514-4519. DOI: https://doi.org/10.1175/2007JAS2317.1

Shonk, JKP and Hogan, RJ. 2008. Tripleclouds: An efficient method for representing horizontal cloud inhomogeneity in 1D radiation schemes by using three regions at each height. *J. Clim.*, 21(11): 2352-2370. DOI: https://doi.org/10.1175/2007JCLI1940.1

Sölch, I and Kärcher, B. 2010. A large-eddy model for cirrus clouds with explicit aerosol and ice microphysics and Lagrangian ice particle tracking. *Q. J. R. Meteorol. Soc.*, 136: 2074-2093. DOI: https://doi.org/10.1002/qj.689

Sonntag, D. 1990. Important new values of the physical constants of 1986 vapor pressure formulations based on ITS-90, and psychometer formulae. *Z. Meteor.*, 70: 340-344.

Spichtinger, P. 2014. Shallow cirrus convection – a source for ice supersaturation. *Tellus A*, 66(1): 1993. DOI: https://doi.org/10.3402/tellusa.v66.19937

Spichtinger, P, Gierens, K, Smit, HGJ, Ovarlez, J and Gayet, J-F. 2004. On the distribution of relative humidity in cirrus clouds. *Atmos. Chem. Phys.*, 4: 639-647. DOI: https://doi.org/10.5194/acp-4-639-2004

Spreitzer, EJ, Marschalik, MP and Spichtinger, P. 2017. Subvisible cirrus clouds – a dynamical system approach. *Nonlin. Proc. Geophys.*, 24: 307-328. DOI: https://doi.org/10.5194/npg-24-307-2017

Stephens, GL. 1983. The influence of radiative transfer on the mass and heat budgets of ice crystals falling in the atmosphere. *J. Atmos. Sci.*, 40(7): 1729-1739. DOI: https://doi.org/10.1175/1520-0469(1983)040<1729:TIORTO>2.0.CO;2

St-Pierre, M and Thériault, JM. 2015. Clarification of the water saturation represented on ice crystal growth diagrams. *J. Atmos. Sci.*, 72(7): 2608-2611. DOI: https://doi.org/10.1175/JAS-D-14-0357.1

Sunilkumar, SV and Parameswaran, K. 2005. Temperature dependence of tropical cirrus properties and radiative effects. *J. Geophys. Res.*, 110: D13205. DOI: https://doi.org/10.1029/2004JD005426

Szczap, F, Gour, Y, Fauchez, T, Cornet, C, Faure, T, Jourdan, O, Penide, G and Dubuisson, P. 2014. A flexible three-dimensional stratocumulus, cumulus and cirrus cloud generator (3DCLOUD) based on drastically simplified atmospheric equations and the Fourier transform framework. *Geosci. Model Dev.*, 7: 1779-1801. DOI: https://doi.org/10.5194/gmd-7-1779-2014

Thomas, L, Cartwright, JC and Wareing, DP. 1990. Lidar observations of the horizontal orientation of ice crystals in cirrus clouds. *Tellus B*, 42(2): 211-216. DOI: https://doi.org/10.3402/tellusb.v42i2.16290

Um, J, McFarquhar, GM, Hong, YP, Lee, S-S, Jung, CH, Lawson, RP and Mo, Q. 2015. Dimensions and aspect ratios of natural ice crystals. *Atmos. Chem. Phys.*, 15: 3933-3956. DOI: https://doi.org/10.5194/acp-15-3933-2015

Vázquez-Martín, S, Kuhn, T and Eliasson, S. 2021. Shape dependence of snow crystal fall speed. *Atmos. Chem. Phys.*, 21(10): 7545-7565. DOI: https://doi.org/10.5194/acp-21-7545-2021

Verhulst, PH. 1838. Notice sur la loi que la population suit dans son accroissement. *Corr. Math. Phys.*, 10: 113-121.

Vogels, M, Zoeckler, R, Stasiw, DM and Cerny, LC. 1975. P. F. Verhulst's "Notice sur la loi que la population suit dans son accroissement". *J. Biol. Phys.*, 3: 183-192. DOI: https://doi.org/10.1007/BF02309004

Wacker, U. 2006. Nonlinear effects in a conceptual multilayer cloud model. *Nonlin. Proc. Geophys.*, 13: 99-107. DOI: https://doi.org/10.5194/npg-13-99-2006

Wada, E, Hashiguchi, H, Yamamoto, MK, Teshiba, M and Fukao, S. 2005. Simultaneous observations of cirrus clouds with a millimeter-wave radar and the MU radar. *J. Appl. Met.*, 44: 313-323. DOI: https://doi.org/10.1175/JAM2191.1

Wang, L. 2004. Midlatitude cirrus cloud structural properties analyzed from the extended facility for atmospheric remote sensing dataset. Unpublished thesis (PhD). University of Fairbanks, Alaska (USA).
**Wang, L** and **Sassen, K.** 2008. Wavelet analysis of cirrus multiscale structures from lidar backscattering: a cirrus uncinus complex case study. *J. Appl. Meteorol. Clim.*, 47(10): 2645-2658. DOI: https://doi.org/10.1175/2008JAMC1788.1

**Wang, PK** and **Ji, W.** 1997. Numerical simulation of three-dimensional unsteady flow past ice crystals. *J. Atmos. Sci.*, 54(18): 2261-2274. DOI: https://doi.org/10.1175/1520-0469(1997)054<2261:NSOTDU>2.0.CO;2

**Westbrook, CD**, **Hogan, RJ** and **Illingworth, AJ.** 2008. The capacitance of pristine ice crystals and aggregate snowflakes. *J. Atmos. Sci.*, 65(1): 206-219. DOI: https://doi.org/10.1175/2007JAS2315.1

**WMO.** 1975. International cloud atlas. Geneva, Switzerland: Secretariat of the World Meteorological Organization.

**Wood, SE**, **Baker, MB** and **Calhoun, D.** 2001. New model for the vapor growth of hexagonal ice crystals in the atmosphere. *J. Geophys. Res.*, 106(D5): 4845-4870. DOI: https://doi.org/10.1029/2000JD900338

**Wu, T**, **Cotton, WR** and **Cheng, WYY.** 2000. Radiative effects on the diffusional growth of ice particles in cirrus clouds, *J. Atmos. Sci.*, 57(17): 2892-2904. DOI: https://doi.org/10.1175/1520-0469(2000)057<2892:REOTDG>2.0.CO;2

**Yagi, T.** 1969. On the Relation between the Shape of Cirrus Clouds and the Static Stability of the Cloud Level - Studies of Cirrus Clouds: Part IV. *J. Meteorol. Soc. Japan*, 47(2): 59-64. DOI: https://doi.org/10.2151/jmsj1965.47.2_59

**Yagi, T**, **Harimaya, T** and **Magono, C.** 1968. On the shape and movement of cirrus uncinus clouds by the trigonometric method utilizing stereophotographs - Studies of cirrus clouds: Part II. *J. Meteorol. Soc. Japan*, 46(4): 266-271. DOI: https://doi.org/10.2151/jmsj1965.46.4_266

**Zeng, X.** 2008. The influence of radiation on ice crystal spectrum in the upper troposphere. *Q.J.R. Meteorol. Soc.*, 134: 609-620. DOI: https://doi.org/10.1002/qj.226

**Zeng, X.** 2018. Radiatively induced precipitation formation in diamond dust. *J. Adv. Mod. Earth Syst.*, 10: 2300-2317. DOI: https://doi.org/10.1029/2018MS001382

**Zeng, X**, **Gong, J**, **Li, X** and **Wu, DL.** 2021. Modeling the radiative effect on microphysics in cirrus clouds against satellite observations. *J. Geophys. Res. Atmosphere*, 126: e2020JD033923. DOI: https://doi.org/10.1029/2020JD033923

**Zeng, X**, **Heymsfield, AJ**, **Ulanowski, Z**, **Neely III, RR**, **Li, X**, **Gong, J** and **Wu, DL.** 2022. The radiative effect on cloud microphysics from the Arctic to the Tropics. *Bul. Am. Met. Soc.*, E2108-E2129. DOI: https://doi.org/10.1175/BAMS-D-21-0039.1

**Zhang, Y**, **Hao, S**, **Cong, M** and **Cao, Y.** 2014. A simulation method of 3D cirrus radiance images for space-based missile warning detectors. In: 7th *International Congress on Image and Signal Processing*, Dalian (China), 14-16 October 2014, 541-546. DOI: https://doi.org/10.1109/CISP.2014.7003839

**TO CITE THIS ARTICLE:**

Berton, RPH. 2023. Two-Dimensional Dynamics of Ice Crystal Parcels in a Cirrus Uncinus. *Tellus A: Dynamic Meteorology and Oceanography*, 75(1): 231-270. DOI: https://doi.org/10.16993/tellusa.3227

Submitted: 22 December 2022

Accepted: 17 June 2023

Published: 18 July 2023

**COPYRIGHT:**

© 2023 The Author(s). This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC-BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. See http://creativecommons.org/licenses/by/4.0/.

*Tellus A: Dynamic Meteorology and Oceanography* is a peer-reviewed open access journal published by Stockholm University Press.

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

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