Subvisible cirrus clouds – a dynamical system approach

Ice clouds, so-called cirrus clouds, occur very frequently in the tropopause region. A special class are subvisible cirrus clouds with an optical depth lower than 0.03. Obviously, the : , :::::::: associated :::: with :::: very :::: low ice crystal number concentration of these clouds is very low ::::::::::: concentrations. The dominant pathway for the formation of these clouds is not known well. It is often assumed that heterogeneous nucleation on solid aerosol particles is the preferred mechanism although ho5 mogeneous freezing of aqueous solution droplets might be possible, since these clouds occur in the low temperature regime T < 235K. For investigating subvisible cirrus clouds as formed by homogeneous freezing we develop a simple parcel ::::::: reduced cloud model from first principles; the : , ::::: which :: is :::: close ::::::: enough :: to ::::::: complex ::::::: models ::: but : is :::: also :::::: simple :::::: enough ::: for :::::::::::: mathematical ::::::: analysis. :::: The model consists of a three dimensional set of ordinary differential equations, and includes the rele10 vant processes as ice nucleation, diffusional growth and sedimentation. We study the formation and evolution of subvisible cirrus clouds in the low temperature regime as driven by slow vertical updraughts (0<w ≤ 0.05m s−1). The model is integrated numerically and also investigated by means of theory of dynamical systems. We found two qualitatively different states for the long-term behaviour of subvisible cirrus clouds. The first state is a point attractor state with a stable focus, i.e. 15 the solution of the differential equations performs damped oscillations and asymptotically reaches a constant value (equilibrium ) :: as ::: an :::::::::: equilibrium :::: state. The second state is a limit cycle in phase space, i.e. the solution approaches asymptotically a state of undamped oscillations :::::::::::: asymptotically ::::::::: approaches :: a :::::::::::::: one-dimensional ::::::: attractor :::: with :::::: purely ::::::::: oscillatory ::::::::: behaviour. The transition between the states constitutes :: is ::::::::::: characterised ::: by a Hopf bifurcation and is determined by two parameters 20 – vertical updraughts :::::::: updraught ::::::: velocity : and temperature. In both cases, the microphysical properties of the simulated clouds agree reasonably well with simulations using :::: from : a more detailed model, with former analytical studies : , and with observations of subvisible cirrus. In addition, the : ,


Introduction
Clouds consisting exclusively of ice crystals, so-called cirrus clouds, are frequently found in the tropopause region at low temperatures (T < 235 K).Satellite observations show frequencies of occurrence up to 40 % in extra-tropical storm tracks and up to 60 % in regions of tropical convection (Stubenrauch et al., 2010).Cirrus clouds influence the energy budget of the Earth-atmosphere system like other clouds by reflecting and scattering incoming solar radiation (albedo effect) and by absorbing and re-emitting thermal radiation (greenhouse effect).For liquid clouds, the albedo effect usually dominates (Stocker et al., 2013, chap. 7) but for pure ice clouds both effects (albedo vs. greenhouse effect) are of comparable absolute size.Therefore, microphysical properties (e.g.ice crystal size or shape; see Zhang et al., 1999) or macrophysical properties (e.g.optical depth or spatial inhomogeneity) can influence the balance between both radiative effects, leading to a net warming or cooling.Nevertheless, for cirrus clouds a net warming of the Earth-atmosphere system is often assumed (Chen et al., 2000).Since the for-mation of ice crystals requires high supersaturation (see e.g.Koop et al., 2000;Hoose and Möhler, 2012) and diffusional growth of ice crystals is quite slow in the low-temperature regime (T < 235 K), cirrus clouds mostly exist in a thermodynamic state far away from equilibrium.Thus, in contrast to liquid clouds, which approximately coincide with their (super-)saturated environment, for ice clouds there can be a continuous transition from clear air over very low ice crystal number concentrations to thick cirrus clouds with high mass and number concentrations.Cirrus clouds with optical thickness τ < 0.03 constitute a special class, so-called subvisible cirrus clouds (SVCs; Sassen and Dodd, 1989).These clouds are difficult to measure; remote-sensing techniques such as lidar (e.g.Immler et al., 2008b) or occultation observations (e.g.Wang et al., 1996) are used to detect these very thin cirrus clouds.Only few in situ measurements of subvisible cirrus clouds are available, suggesting very low values in ice crystal number concentrations (Froyd et al., 2010;Kübbeler et al., 2011).Global observations from satellites (Wang et al., 1996;Stubenrauch et al., 2010;Hoareau et al., 2013) as well as observations with stationary lidar systems (Sassen and Campbell, 2001;Hoareau et al., 2013) show frequencies of occurrence of about 10-20 % in the extra-tropics; in the tropics the frequency of occurrence is much higher (up to 50 %; see e.g.Wang et al., 1996).For subvisible clouds, a net warming of the Earth-atmosphere system is almost certain, since the albedo effect is almost negligible.Our knowledge of subvisible cirrus clouds is quite limited.Since the ice crystal number concentration in SVCs is very low, the question about the dominant formation mechanism is still pending.At cold temperatures (T < 235 K), where pure ice clouds occur, two different formation mechanisms are generally possible, namely heterogeneous nucleation at solid aerosol particles (e.g.Dufour, 1861;aufm Kampe and Weickmann, 1951;Hosler, 1951) and homogeneous freezing of aqueous solution droplets (Sassen and Dodd, 1989;Koop et al., 2000).For subvisible cirrus, Kärcher and Solomon (1999) stated that both nucleation mechanisms might be possible; in contrast, Jensen et al. (2001) and Froyd et al. (2010) clearly suggested that the dominant mechanism must be heterogeneous nucleation.However, analytical investigations by Kärcher (2002) indicated that also pure homogeneous nucleation might be possible.
In the present study we focus on the formation of SVCs by homogeneous freezing of aqueous solution droplets (hereafter: homogeneous nucleation).We study the formation and evolution of SVCs in an air parcel that is lifted in slow vertical upward motions (0 < w ≤ 0.05 m s −1 ), as typical for synoptic-scale motions in the extra-tropics (e.g.along warm fronts; see Kemppi and Sinclair, 2011) or in slow ascent regions in the tropics, for example driven by Kelvin waves (Immler et al., 2008a).We concentrate on the coldtemperature regime (T < 235 K); thus, we exclude the possibility of liquid origin ice clouds (Krämer et al., 2016;Wernli et al., 2016), i.e. freezing of pre-existing cloud droplets at states close to water saturation.This is not a strong limitation since the microphysical properties of ice clouds stemming from mixed-phase clouds are quite different, with high ice crystal number and mass concentrations and higher optical depths (Luebke et al., 2016).
For the investigation of subvisible cirrus clouds we develop a parcel model to which we apply numerical and analytical tools.The model is developed on the basis of an evolution equation for mass distributions of ice crystals, including a description of microphysical processes based on former work (Spichtinger and Gierens, 2009).We take into account the relevant processes for ice microphysics, i.e. ice nucleation, ice crystal growth due to diffusion of water vapour, and sedimentation of ice crystals.We make use of some appropriate simplifications in order to obtain a reduced model consisting of an autonomous system of ordinary differential equations (ODEs), suitable for the application of analytical tools.The variables of the system are ice crystal mass and number concentration, respectively, as well as relative humidity with respect to ice.Thus, we have to investigate a three-dimensional autonomous system of ODEs.
To study the qualitative behaviour of the model we use concepts from theory of dynamical systems (see e.g.Verhulst, 1996;Argyris et al., 2010;Hirsch et al., 2013).The qualitative properties of the system near equilibrium states are relevant for the overall behaviour of the system.The stability of these equilibrium states can be investigated by applying perturbations after linearisation and is determined by the eigenvalues of the linearised system.Some theorems are available in order to transfer the qualitative behaviour of the linearised systems to the full nonlinear system.For the characterisation of more complex attractors, e.g.limit cycles, more sophisticated approaches must be used.For instance, limit cycles can be determined using Poincaré sections (Argyris et al., 2010).Investigations of cloud models as dynamical systems were carried out for liquid and mixedphase clouds (Hauf, 1993;Wacker, 1992Wacker, , 1995Wacker, , 2006) ) as well as for cloud-aerosol-precipitation systems (Koren and Feingold, 2011;Feingold and Koren, 2013).For pure ice clouds such investigations have not been carried out yet.In contrast to clouds involving liquid phase, which are close to thermodynamic equilibrium (i.e.RH ∼ 100 %), we have to consider relative humidity as a system variable, which adds another equation to the system and makes the analysis more challenging.The mathematical characterisation of the reduced model allows for a better understanding of the interaction of different nonlinear processes and the impact of external forcings such as vertical updraughts.Finally, the qualitative analysis could be used in future work as starting point for developing cloud parameterisations that represent the qualitative structure of subvisible cirrus clouds.
In Sect. 2 we describe the development of the model.The results of the numerical integration and the mathematical analysis are presented in Sect.3, as well as comparisons with observations and more detailed models.In the final section, we summarise the results, draw some conclusions and give an outlook on future work.

Model
In this section we describe the development of a reduced ice cloud model, which is later used for analytical and numerical investigations.We include the relevant processes for formation and evolution of ice clouds into the model but we try to avoid too much complexity, which makes analysis too complicated (i.e.reducing the complexity paradox; see e.g.Oreskes et al., 1994;Oreskes, 2003).Since we investigate subvisible cirrus clouds in the temperature regime T < 235 K and at low vertical updraughts 0 < w ≤ 0.05 m s −1 , the relevant processes are ice nucleation, diffusional growth and sedimentation, respectively.

Basic equations
An ice cloud is represented by an ensemble of ice particles, which can be described by a mass distribution f (m, x, t) with mass of particles m as internal coordinate and space, x, and time, t, as external coordinates.Notation follows the convention in population dynamics (see e.g.Ramkrishna, 2000).We investigate a test volume with a certain fixed mass of dry air; f is normalised by the total number concentration per unit mass of dry air, N c , with units The evolution of this mass distribution in time and space is determined by a partial differential equation (see e.g.Hulburt and Katz, 1964;Seifert and Beheng, 2006;Beheng, 2010): Here, ρ denotes density of air, u and g are the advection velocities in physical space and phase space of the internal coordinate, and h represents sources and sinks for particles.The divergence in physical space is denoted by ∇ x = (∂/∂x, ∂/∂y, ∂/∂z) T .Note that all functions (u, g, h) generally depend on the full set of variables (m, x, t).The fluid velocity v = v(x, t) describes the motion of the air; cloud particles may experience a velocity v = v (m, x, t) relative to v; thus, the total u is given by u(m, x, t) = v(x, t)+v (m, x, t).
In our study, the only relevant relative velocity of cloud particles is gravitational settling (hereafter: sedimentation), given by a terminal velocity due to balance between gravitational force and drag.The terminal velocity depends on ice crystal mass, i.e. v = (0, 0, −v t (m)).Note the direction towards the Earth's surface, indicated by the minus sign.Instead of solving Eq. (1) for the entire mass distribution, we derive equations for the general moments of f (m, x, t), defined as A bounded mass distribution is uniquely determined by all its integer moments (see e.g.Feller, 1971).The evolution equations for the general moments are derived by multiplication of Eq. ( 1) by m k and integration by parts, using f (0, x, t) = 0 and lim m→∞ f (m, x, t) = 0 as physically meaningful assumptions.Using v = (0, 0, −v t (m)), and the mass continuity equation ∂ρ ∂t + ∇ x • (ρv) = 0, yields particle formation/elimination (3) Since we cannot (and do not want to) treat an infinite number of moment equations, we make the usual ansatz (see e.g.Seifert and Beheng, 2006) for a double-moment scheme (k = 0, 1); i.e. we derive two equations for number concentration (N c = µ 0 ) and mass concentration (q c = µ 1 ) of ice crystals from Eq. ( 1).Note that the units of N c and q c relative to the mass of dry air are [N c ] = kg −1 and [q c ] = kg kg −1 , respectively.For closing the system of equations mathematically, we prescribe a fixed type of mass distribution for the ice crystals.As in the study by Spichtinger and Gierens (2009), we use a log-normal distribution of the following form: with geometric mean mass m m and non-dimensional geometric standard deviation σ m , determining the width of the distribution; log denotes the natural logarithm.The general moments can be described by using the mean mass m = q c /N c = µ 1 /µ 0 .Here, we introduced the dimensionless parameter r 0 is set to a constant; thus, the geometric standard deviation representing the distribution's width is assumed to be constant.Spichtinger and Gierens (2009) suggest a value of r 0 = 3, corresponding to a geometric standard deviation σ m ≈ 2.85.

Parameterisation of relevant processes
In the following the representation of relevant processes is described briefly.For more details we refer the reader to Appendix A. Furthermore, we describe additional assumptions for simplification and present the final equations of the model.

Nucleation
Particle formation in terms of ice nucleation is described by the last term on the right-hand side of Eq. (3).For the formation of ice crystals we exclusively consider homogeneous freezing of aqueous solution droplets (short: homogeneous nucleation; Koop, 2004).We describe the ensemble of solution droplets by a size distribution f a = f a (r), where r denotes the radius.f a is normalised by the total number concentration of solution droplets per unit mass of dry air, We describe homogeneous nucleation as a stochastic process with a nucleation rate J (for details see Appendix A).For the change in the size distribution f a (r) we can formulate the following equation (according to Seifert and Beheng, 2006) assuming J as a volume rate (i.e.[J ] = m −3 s −1 ): Integration of the equation over all radii r leads to an equation for the total loss of solution droplets Assuming a bijective relation between ice crystals and solution droplets, we combine the total gain of ice particles as where µ 3,a [r] denotes the third moment of the size distribution of solution droplets.Here, we assume that ∂J /∂r = 0. Since the ice crystal number concentration in SVCs is very low, we assume that only a minor fraction of solution droplets is converted to ice and the size distribution remains constant in time.Thus, the third moment can be calculated once and is then used as a constant in the resulting equations.We assume f a (r) as a log-normal distribution with a modal radius of r m = 100 nm, a dimensionless geometric standard deviation σ r = 1.5 and a total number concentration ρN a = 3 × 10 8 m −3 , similar to the settings by Spichtinger and Gierens (2009), which are motivated by observations (Minikin et al., 2003).This leads to a formulation of and using a typical droplet mean mass m 0 = 10 −15 kg (size ∼ 1 µm) in the spirit of the mean value theorem.The nucleation rate J is parameterised according to Koop et al. (2000) and can be expressed as a function of relative humidity with respect to ice and temperature.For further details see Appendix A.

Diffusional growth
The growth and evaporation of ice crystals is dominated by diffusion of water vapour.With several simplifications of the growth equation (for details see Appendix A) we obtain the following equation for diffusional growth of a single crystal: with constants C i = 1.02 m kg −α i , α i = 0.4 and using relative humidity over ice with saturation mixing ratio Here, p si (T ) denotes saturation vapour pressure over ice and ε ≈ 0.622 is the ratio of molecular masses of water vapour and air.We can express the term for diffusional growth in the moment Eq. (3) by integration, i.e.

Sedimentation
Following Spichtinger and Gierens (2009), we describe the weighted terminal velocity vk for the flux of the kth moment as (for details see Appendix A).Here, we use a simple power law for the representation of the terminal velocity with γ = 63 292.36 m s −1 kg −δ , δ = 0.57 and a density correction term corr(T , p) (see Appendix A).We can compose the general terms for sedimentation in the moment Eq. (3):

Simplifications
In order to obtain a consistent but simplified system of ODEs we make the following three assumptions: 1. Change to Lagrangian point of view and purely vertical motion: The Eulerian time evolution and advection of a quantity φ in the fluid motion can be seen as total time derivative representing the Lagrangian description.We will exclusively consider vertical motions of the air parcel as driven by a vertical velocity component w, i.e. v = (0, 0, w(t)), and for simplicity, we assume the mass distribution to be horizontally homogeneous.In order to close the system, we additionally derive equations for temperature and pressure rates assuming hydrostatic balance and adiabatic changes.
Here, g denotes acceleration of gravity, M air is the molar mass of dry air and c p is the molar isobaric heat capacity.We would expect additional temperature changes due to phase changes (latent heat release), when ice crystals grow or evaporate by water vapour diffusion.
Since we investigate ice clouds in the low-temperature regime, temperature changes due to latent heat release can be neglected in good approximation.For low temperatures (T < 235 K) the deviation from the dry adiabatic lapse rate is less than 5 % and is decreasing with decreasing temperature.Therefore, we omit temperature change due to latent heat release, which would appear as an additional nonlinear term in the system of equations.
2. Closure using an equation for relative humidity with respect to ice: In our study, we will exclusively consider very low vertical velocities (0 < w ≤ 0.05 m s −1 ), which are typical for formation of SVCs in large-scale upward motions.We are interested in long-time behaviour of the model.A persistence of such weak updraughts for a long time (e.g. 12 h or even longer, resulting in temperature changes smaller than 10 K) is realistic for warm fronts at mid-latitudes (Kemppi and Sinclair, 2011) or Kelvin waves in the tropics (Immler et al., 2008a).In a simple but quite realistic approximation we assume constant vertical velocity.
As temperature decrease at slow upward motions is only very small, in a zeroth-order approximation we assume constant temperature and pressure.In consequence, the parcel's volume remains constant, too.The resulting error for neglecting density changes is usually of order ∼ 10 % (see e.g.Weigel et al., 2016).Since we are primarily interested in a simple conceptual model of reduced complexity, describing the main properties of SVCs, these assumptions are justified.Thus, in our reduced model w, p and T are assumed to be constant and are treated as control parameters.
To close the systems of differential equations we introduce an evolution equation for relative humidity, starting with the total derivative of RH i : While temperature and pressure remain approximately constant during parcel ascent, the relative humidity should nevertheless be affected by terms involving dT /dt and dp/dt, respectively.Neglecting latent heat release as stated above, the first two terms in Eq. ( 19) read as L ice is the molar heat of sublimation; we use the parameterisation for L ice by Murphy and Koop (2005).
As usual, g denotes gravitational acceleration.Note that we only consider temperature and pressure changes in Eq. ( 19), but leave temperature and pressure constant otherwise and thus obtain a reduced model with only three variables: N c , q c , RH i .This approach will be useful for analytical investigations of the long-term behaviour of the system.
The last term in Eq. ( 19) represents the sink due to diffusional growth of ice particles and can be written as We use relative humidity as a control variable instead of specific humidity, which www.nonlin-processes-geophys.net/24/307/2017/ Nonlin.Processes Geophys., 24, 307-328, 2017 has been used in former studies (e.g.Hauf, 1993;Wacker, 1992) for liquid or mixed-phase clouds close to thermodynamic equilibrium (water saturation).Since pure ice clouds often exist at states far away from equilibrium, relative humidity over ice is the relevant thermodynamic variable, controlling growth and nucleation of ice crystals.

Approximation of sedimentation:
Since we are interested in an analytically treatable model of a single air parcel, we need to eliminate the partial derivatives describing sedimentation, which generally lead to a hyperbolic system of partial differential equations, which is too complicated for theoretical analysis.For simplification of the equations we have to consider terms of the form i.e. vertical gradients in the sedimentation flux, j k = ρ vk µ k .Since the volume does not change, we assume a box with volume V = A • z with constant vertical extension z and constant base area A. The sedimentation flux j k is perpendicular to the surface of the base area.
We approximate the vertical change of the flux by centred differences: We investigate the top layer of a cloud; therefore, by definition

Final system of ODEs
In summary, the full system of the model equations reads as where a, b, c, d, e, f > 0 denote positive real constants as indicated in Appendix A. Note that almost all coefficients also depend on the (fixed) parameter T .This reduced model is an autonomous system of ordinary differential equations; i.e. we can write the system in the following form: and F the right-hand side of Eq. ( 25).Note that the assumption of constant temperature, pressure and vertical velocity ensures that the system Eq.( 25) possesses equilibrium states.

Setup
We examine the system for a range of parameter values 0 < w ≤ 0.05 m s −1 and 190 K ≤ T ≤ 230 K, at a constant pressure of p = 300 hPa, which corresponds to upper tropospheric conditions with moderate vertical motions as in synoptic weather situations or slow upward motions in the tropics (e.g.Kelvin waves).
We investigate the reduced model using analytical tools (see details in Sect.3) and also integrate the model numerically.For this purpose, the air parcel is initialised with no ice particles (N c (0) = 0, q c (0) = 0) and at high supersaturation with respect to ice (RH i (0) = 140 %).The prognostic Eq. ( 25) are integrated numerically with the LSODA algorithm from the FORTRAN library ODEPACK (Hindmarsh, 1983).

General features of the system
The general cloud formation mechanism works as follows: the adiabatic cooling causes the relative humidity, and thus the nucleation rate, to rise until ice nucleation occurs.Due to the steepness of J with respect to RH i , occurrence of ice nucleation corresponds approximately to a threshold in relative humidity (∼ 140-150 %; see e.g.Ren and Mackenzie, 2005;Kärcher and Lohmann, 2002).The stronger the dynamical forcing w, the stronger the nucleation event and the more ice particles form.Ice particle growth then reduces the relative humidity (see Eq. 19, last term) and hence the nucleation rate is also reduced.Crystals grow to larger sizes and begin to sediment out of the air parcel.Sedimentation reduces ice crystal mass and number concentrations, and thus weakens the growth term.Then relative humidity can increase again allowing the cycle to start over.The sedimentation process allows for oscillations in the system; without sedimentation (the only sink for N c and q c ) RH i would drop to values close to saturation and q c would permanently increase; no equilibrium state would be reached for long integration times (see e.g.Kärcher, 2002;Spichtinger and Gierens, 2009).
From the numerical simulations we found that the system exhibits two qualitatively distinct behaviours, depending on values of w and T .First, we give a qualitative overview: The continuous nucleation as well as similar timescales of nucleation, growth and sedimentation leads to a damped oscillation with an equilibrium state for t > 7 h.In phase space, the asymptotic stability is more obvious (see Fig. 5).
-State 1: At rather high temperatures and slow vertical velocities, the three competing microphysical processes (nucleation, growth, sedimentation) are relatively slow and act on similar timescales, so none of them is dominant.In particular, nucleation rates are rather small in these cases; therefore, only few ice crystals are formed initially, which grow and also sediment quite slowly.
The three processes are more or less in balance, resulting in a damped oscillation in all three variables, N c , q c , RH i , asymptotically reaching an equilibrium state, as shown in Fig. 1.Note that, in this state, nucleation is always present, as strong supersaturation with relative humidity close to the nucleation threshold persists at all times and thus the nucleation rates are high enough to produce considerable amounts of ice crystals continuously.This results in smooth oscillations instead of sharp nucleation events, as usually expected (see e.g.Kärcher and Lohmann, 2002).If the air parcel is not dis-turbed and the vertical updraught remains unchanged in the long-term evolution, the cloud persists and has constant ice crystal mass and number concentrations.The cloud in the steady state typically contains low crystal concentrations.The equilibrium state remains at high supersaturations; i.e. the cloud stays far away from thermodynamic equilibrium.
-State 2: When increasing w or decreasing T , respectively, to a certain level, oscillations in variables N c , q c , RH i are no longer damped (see Fig. 2) and no asymptotic equilibrium can be observed (e.g. a point in phase space).Instead, we obtain pulse-like nucleation with distinct nucleation events followed by phases with almost vanishing nucleation rates at low relative humidities.The amplitude of the oscillation is very large in all variables; due to sedimentation, ice particle concentration is reduced to a small fraction of the maximum value once in a period.At colder temperatures and faster vertical velocities, the nucleation rates are much higher, so nucleation is the dominant process in the beginning, leading to pulse nucleation events.After a while, ice crystal growth becomes dominant and when the crystals have become large, sedimentation sets in and crystal numbers decrease rapidly.Finally, the cycle starts over.In this case, the nucleation events are clearly separated, as opposed to the first case.For the time evolution we find that in the beginning, the amplitudes in the three variables decrease slightly from one event to the next, but after a while, the amplitude stays constant.Therefore, it seems that the system asymptotically approaches a limit cycle (one-dimensional attractor).This kind of scenario was also observed in former studies (e.g.Spichtinger and Cziczo, 2010;Kay et al., 2006) but not in a long-term behaviour.
Obviously, we find two qualitatively different states in the numerical solution of the model, depending on parameters w and T , respectively.Next, we investigate the model by means of qualitative theory of dynamical systems.

Qualitative behaviour of the model
For a first investigation we discuss the different terms in Eq. ( 25).The model is driven by an external source; vertical lifting of the air parcel leads to increase of relative humidity.Since temperature and pressure are kept constant, the term e • w • RH i implies a permanent external water vapour source, which is necessary for studying the long-term behaviour of the model.The source of water vapour leads to particle generation.Nucleation can be interpreted as external source, since it is forming ice crystals via water vapour from an external inexhaustible reservoir of solution droplets.
constitute external sinks for cloud variables.Qualitatively, the external sources of water initiate particle generation by nucleation; diffusional growth terms transform water vapour mass into cloud mass until the mass is lost by the external sinks of sedimentation.The model does not fulfil mass conservation due to sources of water vapour and sinks of cloud mass.Because of the external sinks due to sedimentation the system experiences dissipation.Therefore, the system can be seen as an externally forced dissipative system with nonlinear right-hand side.
For a first analysis of the system we compute the divergence of the system (i.e. the trace of the Jacobian DF): using the mean mass m = q c /N c for cloudy states.For clear air (N c = q c = 0), we obtain ∇•F = e•w > 0, hence the system is expanding in phase space.For cloudy air (m > 0) there is competition between different terms determining the sign of ∇ • F .Sedimentation and change of relative humidity due to diffusional growth are sinks (i.e.negative sign in Eq. 27), while the external source term always has a positive sign.Diffusional growth of ice particles can change its sign depending on the thermodynamic state.Since we always investigate a situation with w > 0, the system stays in a supersaturated state (RH i − 100 % > 0); therefore, the last term in Eq. ( 27) is positive.
The balance of terms in Eq. ( 27), i.e. the sign of ∇ • F for cloudy air, is crucially determined by the mean mass of the cloud.Note that for both exponents we have 0 < α i < δ < 1, and thus −1 < α i − 1 < 0. For large ice crystal mass, the terms of form m δ will dominate, leading to a negative sign of ∇ • F and to contraction of the system, mainly due to sedimentation of ice crystals.This is especially the case at higher temperatures, since then diffusional growth is faster and mean masses m tend to larger values.In such cases, the system tends to state 1.
For very small ice crystals, the term including m α i −1 will dominate leading to a positive sign of ∇ • F .For instance, at nucleation events, the ice crystal mass becomes very small; thus, in this situation the system tends to expand explosively (∇ •F > 0).The same is true if almost all particles have fallen out and only small ice crystals are contained in the air parcel.These scenarios are more prevalent at state 2, i.e. at lower temperatures and higher upward velocities.

Linear stability of the system
In a first step, the autonomous dynamical system Eq.( 25) can be characterised by its equilibrium states x 0 , i.e. the points in phase space where F (x 0 ) = 0.The equilibrium states of this system cannot be determined analytically, due to strong nonlinearities.We determine the roots of the right-hand side of system (25) numerically.First, we observe that the mass rate of nucleation dq c dt nucleation = a • m 0 • J (RH i , T ) is negligible compared to other mass rates in the system and can be omitted for simplification.This leads to a new system ẋ = F (x).After setting F (x) = 0, the three resulting equations can be combined to a single equation for RH i as follows: For details of the derivation of this equation see Appendix B.
The roots of Eq. (28) determine the equilibrium values of RH i .Then, the values of N c and q c can be derived analytically.Equation ( 28) has a unique solution, i.e. a single point in phase space, because the left-hand side is a strictly monotonic increasing function of RH i and the right-hand side is strictly monotonic decreasing.Therefore, there exists a unique equilibrium point, x 0 , in the relevant domain of the phase space (RH i > 100 %, N c > 0, q c > 0).The roots of Eq. ( 28) are determined numerically for the relevant domain in the parameter space, i.e. 0 < w ≤ 0.05 m s −1 and 190 ≤ T ≤ 235 K.
In order to examine the qualitative behaviour of the solution in a neighbourhood of the equilibrium state, the ODE system is linearised about the equilibrium state x 0 : where DF| x 0 is the Jacobian of F evaluated at x 0 .Note that F (x 0 ) = 0 by definition.The three eigenvalues of the Jacobian, λ 1 , λ 2 , λ 3 , determine the quality of the equilibrium state (Verhulst, 1996, chap. 3).The eigenvalues must be determined numerically for the relevant parameter values w and T .The Jacobian of the system has two complex conjugate eigenvalues, λ 1,2 ∈ C, whose real part can be positive or negative, depending on the parameters, w and T .In Complex eigenvalues of the linearised system indicate oscillatory behaviour, which is prevalent in all simulations.As can be seen in Fig. 3, the real part of the complex eigenvalues λ 1,2 can change its sign depending on parameters w and T .
For negative values of the real part (Re(λ 1,2 ) < 0) the equilibrium state x 0 is stable; i.e. solutions of the ODE (29) starting in a neighbourhood of this point approach this point asymptotically (Verhulst, 1996, Chapter 2).Thus, this equilibrium point can be characterised as stable focus (e.g.Verhulst, 1996;Argyris et al., 2010).According to the Poincaré-Lyapunov theorem (Verhulst, 1996, theorem 7.1), an asymptotically stable linearised system also ensures asymptotic stability of the full nonlinear system Eq.( 25).Therefore, x 0 is asymptotically stable for the nonlinear system Eq.( 25) and constitutes a stable focus.Since there is a unique equilibrium state, all trajectories in phase space tend to this point asymptotically.
In this case the equilibrium point (stable focus) corresponds to state 1 in the numerical simulations.Solutions of the system Eq.( 25) experience damped oscillations until they asymptotically approach the stable focus in phase space.The imaginary part of the complex eigenvalues determines the oscillation period.Figure 5 shows the trajectory of a solution of the system Eq.( 25) in the three-dimensional phase space, spiralling towards the equilibrium point.
For positive values of the real part (Re(λ 1,2 ) > 0) the equilibrium point x 0 is an unstable focus (i.e.Im(λ 1,2 ) = 0).Solutions starting in a neighbourhood of x 0 run away from the unstable equilibrium point.In this case, the identification of an unstable critical point in the linearised system is not sufficient for a general characterisation of the full nonlinear system, since after short time the solutions are too far away from the equilibrium states and linear stability is no longer www.nonlin-processes-geophys.net/24/307/2017/ Nonlin.Processes Geophys., 24, 307-328, 2017 applicable.Numerical integration shows undamped oscillations for solutions that do not start in the equilibrium point; this behaviour points to the possibility of a limit cycle (onedimensional attractor).The transition from a stable equilibrium point to limit cycle is a so-called Hopf bifurcation (Verhulst, 1996) and is associated with a transition from two conjugate complex eigenvalues with negative real part to two conjugate complex eigenvalues with positive real part, via two purely imaginary eigenvalues.For a vanishing real part of λ 1,2 , the Hopf bifurcation occurs.The existence of a limit cycle cannot be shown analytically for this system; however, we can determine the limit cycle numerically.For starting our calculation close to the limit cycle, we compute the Poincaré map of the system (Argyris et al., 2010;Verhulst, 1996).We choose a two-dimensional plane in phase space, which is transverse to the trajectory of the solution of Eq. ( 26); is called a Poincaré section.The sequence of points in phase space where the trajectory crosses converges numerically to the point on the limit cycle that is in .Once we find one such point on the limit cycle, we can use it as the initial condition in Eq. (26) to compute the complete limit cycle.An example of a Poincaré section for determining the respective limit cycle is shown in Appendix C (Fig. C1).The limit cycle itself constitutes a one-dimensional stable attractor, i.e. solutions starting outside of the limit cycle approach the limit cycle asymptotically.Figure 6 shows the trajectory of a solution of the system Eq.( 25) in the three-dimensional phase space, approaching the limit cycle, which forms a warped circle in phase space.
The transition between the two general states of the system (stable point attractor vs. limit cycle) can be represented in a bifurcation diagram of the w-T space (Fig. 7).The bifurcation point is a function of both w and T .The different states are separated by points with vanishing real part of eigenvalues λ 1,2 , indicated by the thick black line.The bifurcation points were obtained numerically.

Quantitative overview
After discussing the different states of the system qualitatively, we now give an overview of the quantitative cloud properties and relative humidity for the stable focus and the limit cycle, respectively.
In the stable focus regime, i.e. state 1 of the system, the equilibrium state corresponds to the properties of the finally persisting cloud.Hence, in this parameter regime, we describe the properties of the modelled cloud by the values of the system variables at the equilibrium point (stable focus).For the limit cycle regime, i.e. state 2 of the system, the unstable equilibrium point x 0 does not describe the changing properties of the cloud since it is only in the centre of the periodic orbit and the trajectory does not approach it.A more revealing measure for the cloud properties in this regime is a probability density of the values the variables take along the limit cycle, or at least median, maximum and minimum values.
Figure 8 shows ice crystal mass and number concentrations, respectively, at the equilibrium state, x 0 , as a function of vertical velocity (q c (w), N c (w)) for different temperature regimes.The solid lines in both panels correspond to state 1 (stable focus, damped oscillations), whereas the dashed lines indicate the values at the equilibrium point, x 0 , for state 2 (limit cycle regime, undamped oscillations); note that for state 2, the equilibrium point x 0 is an unstable focus.
Ice crystal number concentrations at the equilibrium point x 0 take values in the range 3×10 2 kg −1 ≤ N c ≤ 2×10 5 kg −1 (Fig. 8, top), which corresponds to ice crystal number densities of 0.1 L −1 ≤ n c ≤ 110 L −1 .Ice crystal mass concentration ranges between 4 × 10 −9 ≤ q c ≤ 3 × 10 −6 kg kg −1 (Fig. 8, bottom).This corresponds to an ice water content of 2.2 × 10 −9 ≤ IWC ≤ 1.4 × 10 −6 kg m −3 .As expected from theory (e.g.Kärcher and Lohmann, 2002) and from former numerical investigations (e.g.Spichtinger and Gierens, 2009), the ice crystal number concentrations display a strong increase with rising vertical velocity.Due to increased crystal growth rates at higher temperatures, N c decreases with rising T .In the double-logarithmic representation in Fig. 8, the number concentrations N c (w) at the equilibrium point x 0 appear as straight lines.For different temperature regimes, there seems to be a constant shift between the curves N c (w), leading to parallel lines in the double-logarithmic representation.
For the limit cycle regime (state 2), we can still derive the values of mass and number concentrations at the equilibrium state x 0 .However, since this point is an unstable focus, another representation is needed to describe the range of ice crystal concentrations.As indicated in Figs.7 and 8, the limit cycle behaviour occurs for temperatures T < 230 K for the investigated updraught regime 0 ≤ w ≤ 0.05 m s −1 .In Fig. 9 we present maximum and minimum values (dashed lines) and median values (dot-dashed lines) for ice crystal number concentrations in the limit cycle regime for temperatures T = 190, 200, 210, 220 K.In addition, the ice crystal number concentration at the unstable focus, x 0 , is displayed (solid lines).We observe a large variation in the number concentrations of up to two orders of magnitude relative to the median.This behaviour is reasonable since sedimentation reduces the amount of ice crystals in a dominant manner, while new ice crystals are formed by nucleation in a pulsating way.The absolute values are in the range 0.2 ≤ ρN c ≤ 200 L −1 .
The mass concentration of the ice crystals is largely determined by the efficiency of diffusional growth.As indicated in the model description (Sect.2), this term depends on temperature and also on number concentration, leading again to a power law relationship as represented in Fig. 8 (bottom) and to a constant shift between the different temperatures, represented as parallel lines.
For the stable focus regime, we can directly investigate the mean mass of the ice crystals, m = q c /N c , at x 0 , which is displayed in Fig. 10.The variation of m at the equilibrium state x 0 due to the vertical velocity is marginal, as indicated in the figure.Thus, we can assume that m can be approximated by a function of temperature.The mean mass at x 0 ranges between m ∼ 10 −12 kg and m ∼ 2 × 10 −10 kg, which corresponds to mean sizes between L ∼ 16 µm and L ∼ 134 µm.For the limit cycle regime (state 2), we indicate the variation in the mean mass by box-and-whisker plots, displaying the median value (red markers) as well as 25th/75th percentiles and minimum/maximum values.Note here that variation of mean mass is usually of one order of magnitude.For cold temperatures the variation is larger due to a higher variability in ice crystal number concentration (see Fig. 9), whereas the mass concentration in ice clouds is mainly dominated by available water vapour.= 190, 200, 210, 220 K).The solid line represents values at the equilibrium state x 0 (stable or unstable focus).For the limit cycle regime the range of ice crystal number concentrations is indicated by the shaded area bounded by minimum and maximum values for the updraught range 0.001 ≤ w ≤ 0.05 m s −1 ; the median is indicated by the dot-dashed line.Period in s . Oscillation periods for the stable focus regime at x 0 (solid lines), and for the limit cycle regime (dashed lines).For the stable focus regime, the periods are obtained from the imaginary part of the complex eigenvalues; for the limit cycle regime, the periods are calculated using the Poincaré map.
As indicated in Sect.3.3, the imaginary part of the complex eigenvalues λ 1,2 determines the period of the oscillations in state 1 (stable focus regime) near the equilibrium point x 0 .In Fig. 11 the period τ = 2π Im(λ 1,2 ) as calculated from the imaginary part is shown for the stable focus (solid lines; colours indicate different temperature regimes).For the unstable focus, the imaginary part of the eigenvalues is not meaningful, as the limit cycle is not within the linear regime of x 0 .Instead, the periods of the limit cycle is shown (dashed lines; colours indicate different temperature regimes) as calculated from the Poincaré map.Note that for decreasing temperature the period τ becomes very large.

Comparison with observations
For comparison with observations we first consider in-situ measurements of ice crystals in subvisible cirrus clouds.Since it is very difficult to measure low number concentrations, only few measurement studies are available.We compare our results with measurements by Kübbeler et al. (2011), Lawson et al. (2008) and Davis et al. (2010).Our model results lead to ice crystal number concentrations in the range 0.1 L −1 ≤ ρN c ≤ 200 L −1 and mean ice crystal sizes in the range ∼ 16 µm ≤ L ≤ 134 µm.Note that the variation in number concentrations span over 3 orders of magnitude and the variation in mean sizes is still within 2 orders of magnitude.These simulated values agree quite well with the measurements.Kübbeler et al. (2011) observed quite high number concentrations of the order of ∼ 100 L −1 for small ice crystals (L ∼ 10 µm) but quite low number concentrations 0.1 ≤ ρN c ≤ 10 L −1 for large ice crystals (equivalent radius r > 50 µm).Lawson et al. (2008) reported ice crystal number concentrations in the range 22.5 ≤ ρN c ≤ 188.8 L −1 with mean value and standard deviation 66 ± 30.8 L −1 for ice crystals in the size range 1 ≤ L ≤ 200 µm.Finally, Davis et al. (2010) reported very low ice crystal number concentrations with a mean value of 2 L −1 and mean sizes of 14 µm during the tropical measurement campaign TC4.However, in their study values from former measurement campaigns are reported to be in the range 10 ≤ ρN c ≤ 100 L −1 and for effective radii 10 ≤ r ≤ 20 µm.
In a second step we expand our comparison to observations from remote sensing.Since SVCs are optically very thin, we investigate the extinction coefficient for the visible part of the spectrum.For comparing our results with measurements, we calculate the extinction β in the solar range using parameterisations by Fu and Liou (1993): where IWC = q c • ρ denotes ice water content in g m −3 and D e in µm is the generalised size.Constants are given by a = −6.656× 10 −3 m 2 g −1 and b = 3.686 µm m 2 g −1 .As a useful approximation we set D e = L, where the quantity L is calculated from the mean mass m using the mass-length relation L = C i m α i , as indicated in Appendix A. In Fig. 12 the values for β are shown for different temperature regimes as calculated for the mean values at the (stable and unstable) focus (equilibrium point).Note that there is only marginal difference in the values for different temperatures.The values are within the interval 10 −4 ≤ β ≤ 0.02 km −1 .Seifert et al. (2007) report mean values for extinctions of SVCs in the range 0.015 ≤ β ≤ 0.02 km −1 with standard deviations σ ∼ 0.005-0.009km −1 (see their Table 3).Our results are in the same order of magnitude or even smaller for slow vertical updraughts.Davis et al. (2010) report much smaller values of extinction scattered in the range 0 < β < 0.01 with a mean value of β ∼ 0.001 km −1 . .Extinction coefficient at x 0 different temperatures in the stable focus state 1 (solid lines) and the limit cycle state 2 (dashed lines).
These SVCs were measured in the tropics at high altitudes (z ∼ 16 km), i.e. at low temperatures T < 195 K, where slow large-scale updraughts due to Kelvin waves of the order of w < 0.01 m s −1 dominate (Immler et al., 2008b).This is consistent with our results; see Fig. 12. Overall, we can state that regarding the high spread in the measurements the results from our reduced model agree well with in situ measurements.

Comparison with other models
For comparison with a more detailed model, we carried out simulations with the box model described by Spichtinger and Gierens (2009) and Spichtinger and Cziczo (2010).This model includes more sophisticated treatment of microphysical processes, although it is also a two-moment bulk model.It allows a change in the shape of ice crystals from almost spherical droxtals to columns.Homogeneous nucleation is treated in detail, including deliquescence of sulphuric acid and integration over the full size distribution of solution droplets.For diffusional growth, kinetic and ventilation effects are included.Finally, temperature and pressure changes due to vertical upward motions and latent heat release is added to the air parcel's temperature.
Henceforth this model is termed "complex model".We scan through the T -w parameter space using initial temperatures in the range 190 ≤ T ≤ 235 K with a temperature increment of T = 5 K and vertical velocities in the range 0.005 ≤ w ≤ 0.05 m s −1 with a velocity increment of w = 0.005 m s −1 , leading to 90 simulations.Additionally, we fixed initial conditions p = 300 hPa and RH i = 140 %.Generally, the results of these simulations are in good agreement with the results of the reduced model in this study.We can again identify regimes in the T -w parameter space representing the known two different states, i.e. damped oscillations (stable focus regime, state 1) and limit cycle behaviour (state 2).In Fig. 13 the case of damped oscillation is shown in both model simulations.Here, initial temperature of T = 220 K is used with a vertical velocity of w = 0.01 m s −1 .Green lines indicate the evolution in the complex model simulation, whereas blue lines represent the evolution in the reduced model.For the variables number and mass concentration, both models produce almost the same values.The onset of ice nucleation is shifted between the two models due to differently detailed representation of ice nucleation in both models.This leads to the difference in relative humidity values.Qualitatively, the models agree very well -the oscillation periods and the magnitudes of the damping are very similar.
For the complex model simulations the environmental conditions change; i.e. temperature and pressure are decreasing due to adiabatic expansion.Thus, no steady state can be reached.The values for ice crystal number concentrations and relative humidity are slightly rising with time in the quasi-steady state at the end of the simulation.Ice crystal mass concentration is slightly decreasing.In Fig. 14, a case of limit cycle behaviour is shown.As in Fig. 13, green lines indicate the complex model simulations and the reduced model results are represented by blue lines, respectively.The initial conditions for both models are given by T = 210 K and w = 0.02 m s −1 .Again, we find very good agreement in the cloud variables N c , q c between both model simulations.Qualitatively they also agree very well in terms of the periods of the oscillations.
The bifurcation diagram displayed in Fig. 7 cannot be reproduced accurately by the complex simulations.Since in the complex model the parameter T is changed during the simulations, switching from one regime to the other is possible within one simulation.If, for instance, a simulation starts at a point in parameter space within the stable focus regime (e.g.high temperature at low updraughts), the time evolu-tion initially follows the damped oscillations as expected from the bifurcation diagram of the reduced model.However, the temperature change leads to a (horizontal) path in the phase diagram (Fig. 7) and at some stage the boundary between the two states is crossed, and from then on, the system will perform undamped oscillations.Indeed, we observe this transition in the complex model simulations.An example for this situation is given in Fig. 15, with initial conditions T = 225 K and w = 0.035 m s −1 .Note that in the limit cycle regime the properties of the theoretically expected limit cycle also change with decreasing T .This results in increasing amplitudes of the oscillations in N c , q c , RH i and in increasing periods.Thus, we can conclude that for realistic simulations including changes in environmental conditions there could be transitions between the theoretically determined states.However, the behaviour of the actual states can still be explained by the phase diagram as obtained from our analytical considerations.
We also compare our results with the analytical model by Kärcher (2002), which includes a more sophisticated representation of nucleation and growth.The relevant equations are treated using typical timescales and approximation of the occurring integrals.Comparison with results by Kärcher (2002) shows good agreement as well.Actually, in our investigations with the reduced model we found low ice crystal number concentrations similar to results by Kärcher (2002); the dependence of number concentrations on w and T also agrees very well with analytical considerations by Kärcher (2002).However, our approach goes beyond the results by Kärcher (2002) since we allow for sedimentation of ice crystals.This additional process leads to the oscillatory behaviour in both states, whereas in the study by Kärcher (2002) a quasi-steady state at ice saturation is reached soon.Especially the continuous nucleation in the state 1 scenario (stable focus, damped oscillation) is only possible if we allow for sedimentation of ice crystals.Otherwise, the nucleation event would stop after depositional growth has reduced the supersaturation such that nucleation rates become negligible.Thus, our scenarios might be more realistic, although values of mass and number concentrations in both studies are very similar.

Conclusions
In this study we have developed a reduced model for describing subvisible cirrus clouds formed by homogeneous nucleation in the tropopause region.The model consists of a set of autonomous ordinary differential equations for the variables ice crystal mass and number concentration, and relative humidity with respect to ice.It contains the relevant cloud processes ice nucleation, diffusional growth and sedimentation.The model can be viewed as an externally forced dissipative system.The model is integrated numerically and also investigated using (linear) theory of dynamical systems.Transition between stable focus regime (state 1) and limit cycle regime (state 2): simulation with the complex model by Spichtinger and Gierens (2009) for w = 0.035 m s −1 and start temperature: T = 225 K.During the first 2 h of the simulation, the sink property can be clearly seen.After reaching temperatures of about T ∼ 220 K, the regime changes from state 1 (stable focus) to state 2 (limit cycle); see also phase diagram in Fig. 7.After this transition, the amplitudes of number concentrations and relative humidity with respect to ice increase and at the end of the simulation also a shift in the oscillation period can be seen.Increase in amplitude and shift in oscillation period are due to changes of the limit cycle properties for decreasing temperature (see e.g.Fig. 11) Integration and theoretical analysis show that the system contains two different states, a stable focus state and a limit cycle state.The states depend on the environmental parameters vertical updraught, w, and temperature, T .The transition between the states can be described as Hopf bifurcation.Both states show oscillatory behaviour, either damped (stable focus) or basically undamped (limit cycle).
Ice crystal mass and number concentrations of the cloud in both states depend mostly on the environmental conditions as vertical velocity and temperature.However, for the limit cycle case the spread in ice crystal mass and number concentration is obviously larger than in the case of stable equilibrium.For the stable focus, the mean mass depends only slightly on www.nonlin-processes-geophys.net/24/307/2017/ Nonlin.Processes Geophys., 24, 307-328, 2017 vertical velocity; thus, we can approximate the mean mass as a function of temperature.
Comparisons with a more detailed box model by Spichtinger and Gierens (2009) show very good agreement.The qualitative behaviour as determined for the reduced model can also be found for the complex model simulations.Also, in terms of quantitative results both models agree quite well.Former analytical investigations by Kärcher (2002) show good agreement with our reduced model, too.However, since we include sedimentation in our model, our results go clearly beyond the former investigations; the long-term behaviour is different, since the inclusion of sedimentation crucially leads to the bifurcation, depending on environmental conditions.
Since there are only few in situ measurements of subvisible cirrus available, it is quite difficult to carry out solid comparisons.However, we try to compare with measurements as described by Kübbeler et al. (2011), Lawson et al. (2008), and Davis et al. (2010) and find good agreement with our model results.Also, the extinction coefficient as calculated from model results agrees very well with observations obtained with remote sensing techniques (Seifert et al., 2007;Davis et al., 2010).
The major qualitative results can be summarised as follows: -Homogeneous freezing of aqueous solution droplets at low temperatures (T < 235 K) is a possible pathway for the formation of subvisible cirrus clouds at low vertical updraughts.Thus, the question about the dominance of formation mechanisms for these thin clouds remains open (homogeneous vs. heterogeneous nucleation).
-In unperturbed weak large-scale updraughts subvisible cirrus clouds can exist in two different qualitative states, reaching either a stable equilibrium point (stable focus) in the long-term behaviour or experiencing oscillation behaviour in a limit cycle scenario.The state depends on external parameters as large-scale updraught and temperature, respectively.
-The cloud particle properties in the long-term behaviour are very similar for both states.Therefore, we cannot decide from values of mass and/or number concentrations in a certain range in which state the cloud might be.Even if we had more measurements, we probably would not be able to decide between the two states just using the Eulerian measurements without a Lagrangian point of view.
We might derive a minimal model for SVCs from the bifurcation diagram in the following way.If we assume that SVCs are well approximated by their attractors, we could express cloud variables and relative humidity by a simple damped harmonic oscillator of the form ẍ + κ ẋ + ωx = 0, (31) with x ∈ {N c , q c , RH i } and parameters κ = κ(w, T ) and ω = ω(w, T ), respectively.κ describes damping, whereas ω represents oscillation frequency.κ, ω can be determined using eigenvalues λ i for damping and oscillations in the stable focus case (κ = 0).For the limit cycle case (κ = 0), periods as obtained from the Poincaré section (see Fig. 11) can be used for describing ω.Such a minimal model could be used for representing SVCs in large-scale models and can be seen as a prototype for new-generation cloud parameterisations.These models describe the structure of clouds in terms of cloud variables and environmental conditions.They could be used for describing such structures embedded into a coarse-grid model.However, further research in this direction is necessary in order to proceed from pure model prototypes to useful cloud parameterisations.Finally, we can state that we could develop a meaningful reduced model for describing the main features of subvisible cirrus clouds.Former investigations using box models indicated that there might be different regimes in the behaviour of the clouds for longer simulation times.For instance, in studies by Kay et al. (2006) and Spichtinger and Cziczo (2010) oscillatory behaviours as well as asymptotic stability could be seen.However, only a detailed mathematical analysis could show that there is a bifurcation in the long-term behaviour and that it depends mostly on environmental parameters such as updraught velocity and temperature.This analysis was only possible since we developed a reduced model, which is close enough to complex models but is also simple enough for mathematical analysis.
The observed Hopf bifurcation as a transition between two different states shows that clouds might exhibit inherent structures, which are crucially determined by the microphysical cloud processes themselves in addition to environmental conditions.Similar structure formation was already seen in analytical cloud models for liquid and mixed-phase clouds as developed by Wacker (1992Wacker ( , 1995Wacker ( , 2006) ) or Hauf (1993).Investigation and analysis of the microphysical processes in terms of sets of ordinary differential equations are a first but urgently necessary step in order to investigate structure formation inside clouds.Once we understand the possible structures in clouds as determined by microphysics, we can continue to further investigate structure formation as driven by spatial diffusion processes, mixing and others, leading to spatial structures of clouds.A first possible approach might be to investigate equations with additional spatial diffusion terms regarding possible Turing instabilities (Turing, 1952).However, further research is necessary in order to investigate structure formation of ice clouds.
Data availability.The data used in this work are described in Sect.3.5.0 10 000 20 000 30 000 40 000 50 000 60 000 70 000 80 000 .The red dots converge fast to two accumulation points, which determine approximately the section of the limit cycle with the plane .If we start "outside" of the limit cycle, the section points (indicated by blue dots) again converge fast to the same two accumulation points.

Figure 1 .
Figure1.A scenario in state 1 (stable focus regime, damped oscillation) at w = 0.01 m s −1 and T = 220 K.The continuous nucleation as well as similar timescales of nucleation, growth and sedimentation leads to a damped oscillation with an equilibrium state for t > 7 h.In phase space, the asymptotic stability is more obvious (see Fig.5).

Figure 2 .
Figure2.A scenario in state 2 (limit cycle regime) is shown at w = 0.02 m s −1 and T = 210 K. Nucleation events occur as pulses; thus, an undamped oscillation evolves, which describes a limit cycle in phase space (see Fig.6).

Figure 3 .
Figure 3. Real (upper panel) and imaginary part (lower panel) of the complex eigenvalues λ 1,2 of the Jacobian DF| x 0 at the equilibrium point x 0 .

Figure 4 .
Figure 4. Real eigenvalue λ 3 of the Jacobian DF| x 0 at the equilibrium point x 0 .

Figure 5 .
Figure 5. Stable focus for state 1 at T = 220 K, w = 0.01 m s −1 : orbit in phase space approaching the equilibrium state asymptotically.

Figure 6 .Figure 7 .
Figure6.Limit cycle for state 2: orbit in phase space at T = 210 K, w = 0.02 m s −1 .Note that the solution starts "outside" of the limit cycle and approaches the limit cycle attractor asymptotically.

Figure 8 .
Figure8.Ice particle number concentration N c (upper panel) and ice particle mass concentration q c (lower panel) at the equilibrium state x 0 as a function of vertical velocity for different temperatures.Solid lines indicate parameter combinations (w, T ) in the stable focus regime (state 1); dashed lines represent the limit cycle regime (state 2), i.e. at the unstable focus x 0 .

Figure 9 .
Figure 9. Ice crystal number concentrations for different temperature scenarios (T= 190, 200, 210, 220 K).The solid line represents values at the equilibrium state x 0 (stable or unstable focus).For the limit cycle regime the range of ice crystal number concentrations is indicated by the shaded area bounded by minimum and maximum values for the updraught range 0.001 ≤ w ≤ 0.05 m s −1 ; the median is indicated by the dot-dashed line.

Figure 10 .
Figure 10.Mean ice crystal mass m as a function of temperature.For the equilibrium state x 0 , values of m depends only slightly on the vertical velocity, the curve covers the area that corresponds to vertical velocities 0.001 m s −1 ≤ w ≤ 0.05 m s −1 .Additionally, box-and-whisker plots indicate median, 25th/75th percentiles, and minimum/maximum values, respectively, for the limit cycle regime.
Figure12.Extinction coefficient at x 0 different temperatures in the stable focus state 1 (solid lines) and the limit cycle state 2 (dashed lines).

Figure 13 .
Figure13.Stable focus case (state 1): comparison between reduced model (this study) and the complex box model bySpichtinger and Gierens (2009).Updraught w = 0.01 m s −1 , temperature in the reduced model and start temperature of the complex model is T = 220 K.

Figure 14 .
Figure14.Limit cycle case (state 2): comparison between reduced model (this study) and the complex box model bySpichtinger and Gierens (2009).Updraught w = 0.02 m s −1 , temperature in the reduced model and start temperature of the complex model is T = 210 K.
Figure15.Transition between stable focus regime (state 1) and limit cycle regime (state 2): simulation with the complex model bySpichtinger and Gierens (2009) for w = 0.035 m s −1 and start temperature: T = 225 K.During the first 2 h of the simulation, the sink property can be clearly seen.After reaching temperatures of about T ∼ 220 K, the regime changes from state 1 (stable focus) to state 2 (limit cycle); see also phase diagram in Fig.7.After this transition, the amplitudes of number concentrations and relative humidity with respect to ice increase and at the end of the simulation also a shift in the oscillation period can be seen.Increase in amplitude and shift in oscillation period are due to changes of the limit cycle properties for decreasing temperature (see e.g.Fig.11)

Figure C1 .
Figure C1.Example of a Poincaré section in the limit cycle regime.Blue dots indicate intersection points of the trajectory with when starting "outside" the cycle; red dots indicate intersection points when starting near the (unstable) equilibrium point x 0 (red cross).