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

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):

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

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^–1 and the rotation rate Ω₀ is of order 10^–4 rad s^–1 (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 $\dot{m}$ governed by Eq. (10) below, has been neglected because the mass loss or gain is isotropical (Plastino & Muzzio 1992).

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 1

Figure 2

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

Physically speaking, k is the reciprocal of a relaxation time τ_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:

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):

the dynamic viscosity of air μ_a, such that, for –50 < T – T_r < 50 (T_r = 273.15 K) (Pruppacher & Klett 1978/1997):

and the mass density ρ_m of moist air (Picard et al. 2008; Nettesheim 2017):

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

defining the terminal fall velocity W_{f ∞} by:

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 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.

Figure 3

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):

10

$\dot{m}=4\pi \hspace{0.17em}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)}$

subject to the initial condition:

11

$m=m_0={\rho }_i{V}_0$

where the initial volume V₀ is estimated from Eq. (31a) below.

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

12

$\{\begin{array}{l}{S}_l=\frac{p_v}{p_{vsl}}\\ S_i=\frac{p_v}{p_{vsi}}\end{array}$

The relative humidity H_l and H_i 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):

13

$\{\begin{array}{l}{H}_l=\frac{p_v}{p_{vsl}}\frac{p-p_{vsl}}{p-p_v}\\ H_i=\frac{p_v}{p_{vsi}}\frac{p-p_{vsi}}{p-p_v}\end{array}$

The atmospheric moisture profiles are detailed in App. A.3, and plotted in Figure 4. All these altitude profiles 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_l 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.

Figure 4

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

14

$q_s=\frac{\epsilon p_{vsl}}{p-p_{vsl}}$

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)

15

$\{\begin{array}{l}{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{array}$

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):

16

$f_v=1+0.039X+0.1447X^2$

where the variable X, such that:

17

$X=\sqrt[3]{\text{Sc}}\sqrt{\mathrm{Re}}$

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

18

$\text{Sc}=\frac{{\mu }_a}{{\rho }_a{D}_v}$

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):

19

$R=\frac{ℛ}{4\pi C{K}_{a}T}\left(\frac{L_s}{R_{v}T}-1\right)$

with the source term:

20

$ℛ=A{Q}_{abs}\left(\frac{F_u-F_d}{2}-B_i\right)$

where Q_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:

21

$B_i=\epsilon \sigma T_s^4$

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):

22

$T_s-T_a=\frac{L_s}{4\pi C{K}_a}\dot{m}$

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 ε 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_abs ≈ 1 (Roach 1976; Liou 2002). It is convenient to express the radiative infrared ratio η as (Zeng 2008; Zeng et al. 2022):

23

$\eta =\frac{F_u-F_d}{2B_i}$

and then recast Eq. (20) as:

24

$ℛ=A(\eta -1)B_i$

Therefore, the crystal is subject to radiative cooling when η < 1, and to radiative warming when η > 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):

25

$\{\begin{array}{l}{F}_u=-B_i\\ F_d=-0.6F_u\end{array}$

then the ratio η becomes:

26

$\eta =0.9$

This is confirmed by the modelling of the parameter η 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 η_a as in Eq. (23), function of altitude. We shall derive an expression from models (Detwiler & Ramaswamy 1990; Zeng 2008) but consistent with 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.

Figure 5

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

27

$\eta -1=\alpha ({\eta }_a-1)$

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:

28

$S_i\ge 1+R$

When η < 1, i.e. $ℛ$ < 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. $ℛ$ > 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):

29

$m=\frac{\pi }{6}{D}_m^3{\rho }_{ie}$

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

30

$m=3\sqrt{3}\hspace{0.17em}{a}^{2}c{\rho }_{ie}$

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^–3.

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):

31a,b

$\{\begin{array}{l}V=\sqrt{3}(2c+c_B)a^2\\ A=6a\left\{2c+\sqrt{{(c-c_B)}^2+\frac{3}{4}{a}^2}\right\}\end{array}$

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.

Figure 6

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:

32

$\begin{array}{l}i) c_B=c/2\\ ii) c_B=c\\ iii) c_B=0\end{array}$

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

33

$\psi =1-\frac{c_B}{c}$

We shall keep c_B constant, equal to its initial value: c_B0 = (1–ψ₀) c₀. For crystals longer than 100 μm, we shall set: ψ = 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 ϕ:

34

$\varphi =\frac{c}{a}$

reaches large values (ϕ > 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 ϕ is constant during the motion. The value ϕ = 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 ϕ we are led to a third-degree equation for a:

35

$2\varphi a^3+c_B{a}^2-\frac{V}{\sqrt{3}}=0$

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

36

$V=\frac{m}{{\rho }_i}$

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:

37

$D_i=\sqrt[3]{\frac{6V}{\pi }}$

or, using Eq. (36):

38

$D_i=\sqrt[3]{\frac{6m}{\pi {\rho }_i}}$

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 ϕ > 1 or 2a if ϕ < 1) following the general definition (29) and using Eq. (36), we could also define an effective mass density ρ_ie such that:

39

${\rho }_{ie}=\frac{6V}{\pi D_m^3}{\rho }_i$

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

40

${\rho }_{ie}=\frac{V}{3\sqrt{3}\hspace{0.17em}{a}^{2}c}{\rho }_i$

which will yield values smaller than the actual density ρ_i. The first definition (Eq. (39)) overestimates the actual volume, and therefore, it underestimates ρ_ie, yielding values as low as 200 kg m^–3 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):

41

$C=0.751a+0.491c=(0.751+0.491\varphi )a$

using the aspect ratio ϕ 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.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 ∆t, we write the velocity components of the average crystal’s center of mass as:

its rectangular coordinates as:

and its mass as:

where T_n = T_a(z_n) and S_n = S_i(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):

and then c_n+1 from Eq. (31a,b) if a is constant (a = a₀):

or solving Eq. (35) for a_n+1 if ϕ is constant (ϕ = ϕ₀):

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₀ (Sec. 4) and to the impact of radiative transfer (Sect. 5).

3.2. Results

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).

Figure 7

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).

Figure 8

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^–1 instead of 1 m s^–1), 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_a < 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^–1 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 9

The variation of ice supersaturation S_i –1 (Figure 10a) is roughly stationary, with a local minimum coincident with the height maximum (t ≈ 1.5 h) and positive between 2 and 3% during the growth phase until 3.4 hours. Then it decreases, becomes negative (t ≈ 3.75 h) in the decay phase, goes down to –20%, consistent with observational data fitted with a microphysics model (Khvorostyanov & Curry 2008).

Figure 10

Figure 11

The radiative correction term R (Figure 10b) is positive and slightly increasing, amounting to 1%, in the growth phase and becomes negative (t ≈ 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 ≈ 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 ≈ 1.5 h), then reaches a maximum (m ≈ 1.2 μg) at t ≈ 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 ≈ 3.75 h), and then decrease to 35 μm and 70 μm respectively, satisfying the initial aspect ratio (ϕ = 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 ≈ 137 μm) at the same time (t ≈ 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 ≈ 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 trail.

Figure 12

The hollowness factor ψ (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, ρ_ie, as calculated from Eq. (40), reaches a minimum when D_i is maximum (t ≈ 3.75 h) and afterwards increases to large values. The range 670–705 kg m^–3 covered during crystal motion is slightly smaller than values amounting to 800 kg m^–3 published for bullet shapes (Schmitt and Heymsfield 2007).

Figure 13

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).

Figure 14

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^–1 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_s–T_a (Figure 15b) is very small and positive in the heated growth phase (T_s–T_a < 0.01 K), then it becomes negative and increases in magnitude (T_s–T_a < 0) in the cooled decay phase, though remaining smaller than 0.1 K, consistently with published results (St-Pierre & Thiérault 2015).

Figure 15

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).

Figure 16

Figure 17

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_w0 = 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₀ ≈ –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 a spiral trajectory winding one turn more at each period. The phase velocity C_p corresponding to the period P and wavelength Λ is defined by:

48

$C_p=\frac{\Lambda }{P}$

that is, numerically: C_p = 70 × 10³/(10.5 × 3600) ≈ 1.85 m s^–1. Now, let us write the condition for the fixed point:

49

$\{\begin{array}{l}\dot{u}=0\\ \dot{w}=0\end{array}$

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

50

$\{\begin{array}{l}u=U_a(z_{\infty })\\ w=W_a(z_{\infty })-\frac{g}{k_{\infty }}\end{array}$

The fixed point satisfies the additional condition:

51

$w=0$

and this finally leads to:

52

$\{\begin{array}{l}{u}_{\infty }=U_a(z_{\infty })\\ W_a(z_{\infty })=\frac{g}{k_{\infty }}\end{array}$

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

Figure 18

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_∞ ≈ 9.63 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_∞.

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 ∞ ≈ – 0.6 m s^–1, which compensates the updraft (W_a = 0.6 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 19

Figure 20

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₀ = 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.

Figure 21

The core half-size c_B being kept constant during the motion (Figure 20b), the hollowness factor ψ (Figure 22a) oscillates with a small amplitude (≈2%) in phase with c, and it tends towards a limit value (ψ_∞ ≈ 0.81) slightly larger than its initial value (ψ₀ = 0.80). The crystal’s effective mass density as calculated from Eq. (40), shows oscillations (Figure 22b) in phase with W_f, and phase opposition with ψ and D_i. It tends towards the asymptotic approximate value 671.5 kg m^–3 approximately.

Figure 22

Figure 23

Figure 24

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_∞ ≈ 16.5 s^–1, which is connected to W_f ∞ (≈ 0.6 m s^–1) 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_∞ ≈ 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.1. Non-oscillatory mode

Below a critical size (a₀ < 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₀ > 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₁₀ – z₀ and z₂₀ – z₀ (Figure 25a) are maximum when a₀ ≈ 53.4 μm and then decrease quite linearly as function of a₀. The spans of z₁₀ – z₀ and z₂₀ – z₀ 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_r < –50°C). Colder cirrus (–85°C < T–T_r < –75°C) at higher altitudes are thinner.

Figure 25

Interestingly it appears that the relative maximum height z_M – z₀ reached at the top of the hook (Figure 25b) plotted as a function of a₀ comes to zero at a₀ ≈ 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₁₀, t₂₀ at the abscissae x = 10 and 20 km plotted as a function of a₀ (Figure 26a) reach several hours for small crystals (a₀ < 52 μm) and shorten to less than one hour for large crystals (a₀ > 55–65 μm).

Figure 26

The half-lengths c₁₀ and c₂₀ at the same times, normalized to the initial value c₀ and plotted as functions of a₀ (Figure 26b) show an interesting behaviour. The profiles have a sharp peak at approximately a₀ ≈ 51.5 μm, and then decrease with inflexion points corresponding roughly to the maxima of the falling height profiles (Figure 25a). Then, c₁₀/c₀ is nearly constant, while c₂₀/c₀ 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 simpliflying assumptions.

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₀ < 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₀ ≈ 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:

and plot the results as functions of a₀ (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. It is noteworthy that the period is nearly independent of crystal size a₀ 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₀ ≈ 56.2 μm reflects the transition explained above.

Figure 27

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 Λ_min and Λ_max of wavelength by the formulae:

and plot the results as functions of a₀ (Figure 28). The two curves have a bump in the range 51–56 μm with a maximum in the middle at a₀ = 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_∞ calculated a posteriori from z_min and z_max (see below). These results are consistent with the value found in Sect. 3.2.3.

Figure 28

From periods and wavelengths, we can calculate corresponding phase speeds by Eq. (48), such that:

and plot the results as functions of a₀ (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₀ ≈ 53.5 μm.

Figure 29

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₀ ≈ 53.3 μm, corresponding to the size for which the formation of a hooked head is cancelled (Figure 25b).

Figure 30

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

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

Figure 31

The damping ratio ζ 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 Ω ≈ Ω₀:

Therefore, extracting ζ we obtain, for instance between extrema 1 and 2:

With the above data providing an average angular frequency Ω₀ = 2π/(10.24 × 3600) = 1.7 × 10^–4 rad s^–1, and using an initial phase shift Φ ≈ π/4 (see App. D2), we obtain the value ζ ≈ 0.019, which is underestimated compared with the theoretical estimate of App. D.2 (ζ ≈ 0.05).

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 at 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 η 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 Λ, 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₀, and the limit altitude of the oscillatory motion, z_∞, are also recorded in Table 2.

5.2. Case 1

We assume that η = 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 ϕ = 2 and c_B is constant. With a radius a₀ = 50.17 μm, steady and oscillatory modes are obtained just by switching the updraft base altitude, with lifetimes of 5 hours and 40 hours respectively. For small saturations (H_a3 = 61%, S_i ≈ 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_i–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).

Figure 32

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_p = 2.03 m s^–1. The damping ratio is ζ ≈ 0.015. As expected, it is consistent both with the velocity limit u_∞ ≈ 2.1 m s^–1 in the slowly spiralling hodograph (Figure 32b) and with the wind speed at the limit height z_∞ ≈ 9.58 km (Figure 32c), that is using Eq. (68) in layer III: U_a(z_∞) = –30 –5(9.58–16) = 2.1 m s^–1. 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.

Figure 33

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_a3 = 68%, S_i ≈ 115%), and updraft (W_a0 = 0.8 m s^–1) 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_B 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₀ = 52.0 μm, and setting the updraft base altitude at 7 km, lower than in Cases 0 and 1.

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³ /(7.8 × 3600) ≈ 7.1 m s^–1. As expected, it approximately matches both the focus u_∞ ≈ 7.1 m s^–1 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^–1.

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.

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 (≈ 10 %) at the beginning, larger than in Case 0 (≈ 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 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_∞ ≈ 7.1 m s^–1.

Figure 34

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

5.4. Case 3

We now assume that η = 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 ≈ 115%, W_a0 = 0.8 m s^–1; ϕ = 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₀ = 53.9 μ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).

Figure 35

Thus the oscillatory mode is slowly damped with average Λ ≈ 200 km, P ≈ 6.8 hours, and ζ ≈ 0.01. According to Eq. (48), the phase speed is therefore C_p = 200 × 10³/(6.8 × 3600) = 8.2 m s^–1. As expected, it matches approximately both the limit u_∞ ≈ 8.5 m s^–1 in the hodograph (Figure 35b) and the wind speed at the limit height z_∞ ≈ 8.3 km (Figure 23c), that is using Eq. (68) in layer III: U_a(z_∞) = –30 –5(8.3–16) = 8.5 m s^–1 Like in the previous case, the motion crossing the critical level (z_U0 = 8 km) of maximum wind speed (U_amax = 10 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 (≈ 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 = ± 1.2%), a situation corresponding to opposite radiative initial conditions (η – 1 = ± 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 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.

Figure 36

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 ζ (Table 2) reveals that Case 2 is more damped than Case 3.

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^–3 s^–1. Consistently with our result, the maximum fall speed of 0.5 m s^–1 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^–1, 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_a, 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):

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

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

and the horizontal wind shear ω:

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

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 R_sn and a criterion for onset of convective instability as: 0 < R_sn < 1. With a typical heating rate $H=0.1$ K d^–1 at the altitude z = 9 km, they show that their criterion could be satisfied.

Figure 37

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 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 at 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 ≈ 2–8 m s^–1) and it is much smaller than the wind speeds (U_a ≈ 30 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 (≈ 0.007 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 s).

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 (≈ 6–10 hours) and wavelength (≈ 60–200 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, η_a), supplemented with the 7 initial scalar conditions: position (x₀, z₀), velocity (u₀, w₀), column half-width (a₀), aspect ratio (ϕ), hollowness factor (ψ). Of the total, we matched only five of them: z_w0, W_a0, H_a3, η_a and a₀.

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 (≈ 5 hours) is of same order as the half-period of the associated oscillatory mode (≈ 8–10 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):

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):

In our problem X ≈ 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, ϕ, ψ) 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 Λ 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_∞:

Thus, increasing or decreasing U_a in the layer under concern (III) will increase or decrease Λ 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:

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 Ω and the damping ratio ζ 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).

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:

with the following 6 constants:

It is noteworthy that, with this choice of constants, U_a vanishes above the ground at the altitude z₀ = 10 km. Moreover, the resulting wind shear ω, respectively equal to 1.0, 1.3 and –5 m s^–1 km^–1 in the ranges 0–5, 5–8 and 8–16 km, is quite consistent with observational data (|ω| < 23 ms^–1 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 (1^st mode or regime) or z_w0 = 7 or 9 km (2^nd mode or regime) and zero below z_w1 = z_w0 – ∆z₀, with a transition layer of thickness ∆z₀ = 0.3 km, the transition profile being modelled by a linear function (Figure 2):

with the following constants:

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

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

that is smaller in magnitude than the vertical shear at the same altitude (ω = –5 m s^–1 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:

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:

with the following 10 constants:

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

with the two constants: p₀ = 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:

with the following 12 constants:

On the same figure is then plotted (dotted line) the profile of H_i derived from H_a, first by calculating the water vapour pressure p_v as solution of Eq. (84) below (H_l = 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_l, as it can be easily proved from Eqs. (12) and (13), since p_vsl > p_vsi. We also plotted the profiles of S_l and S_i (Figure 4b), though they are very close to those of H_l 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 η_a(z) defined by Eq. (23) and can be given in the upper troposphere as a piecewise linear function composed of three segments:

with the following 4 constants:

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

D.1. Non-oscillatory mode

Solving system (1) with analytic functions is a difficult task as soon as k and W_a 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_a, constant and uniform vertical updraft W_a and free fall W_f, governed by the set of ODE:

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:

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