Pattern formation in clouds via Turing instabilities

Pattern formation in clouds is a well-known feature, which can be observed almost every day. However, the guiding processes for structure formation are mostly unknown, and also theoretical investigations of cloud pattern are quite rare. From many scientific disciplines the occurrence of pattern in non-equilibrium systems due to Turing instabilities is known, i.e. unstable modes grow and form spatial structures. In this study we investigate a generic cloud model for the possibility of Turing instabilities. For this purpose, the model is extended by diffusion terms. We can show that for some cloud models, i.e special cases of the generic model, no Turing instabilities are possible. However, we also present an example of a cloud model, where Turing instabilities can occur. Using numerical simulations of this special case, we provide spatial pattern of clouds in one and two spatial dimensions. From the numerical simulations we can see that the competition between collision terms and sedimentation is an important issue for the existence of pattern formation.


Introduction
Pattern formation is a general feature in nature. We find pattern in many different locations and research fields, e.g. sand ripples at sand dunes or at the beach, stripes on zebras and fishes, convective cells in Rayleigh-Benard convection, spiral states in chemical reaction systems as e.g. the famous Belousov-Zhabotinski system, and many other examples. The generation of structures is a common feature for systems out of thermodynamic equilibrium. In contrast to states at equilibrium, which tend to be homogeneous, an external forcing driving a system out of equilibrium has the potential to form new structures. These structures can have different forms, i.e. homogeneous or inhomogeneous in space and stationary or oscillatory in time (see, e.g. [Cross, M. C., Hohenberg, P. C. (1993)]). Pattern formation is an emergent process, and is usually not predictable a priori from the underlying micro states of the system; the structures on larger scales often appear in a spontaneous way. Research on pattern formation is an important field in many disciplines in natural sciences e.g. mathematical biology ( [Murray, J. D. (2003)]), chemistry ( [Kondepudi, D. and Prigogine, I. (2015)]), fluid dynamics (e.g. Rayleigh-Benard convection, see [Bodenschatz, E. et al. (2000)]) and many other fields.
There are several approaches to represent pattern formation in models. One of the first approaches was presented by [Turing, A. M. (1952)] in his seminal article on morphogenesis. Chemical reactions are represented by a system of ordinary differential equations (ODEs). This set of equations is extended by diffusion terms, i.e. a Laplacian in spatial directions is added to each equation representing the concentration of a chemical species. It can be shown by linear stability analysis that under certain conditions (e.g. different diffusion coefficients) stable fix points of the ODE system can be destabilized, i.e. some Fourier modes become unstable and grow, until they become saturated by nonlinear effects. Since only wave numbers out of a finite interval become unstable, spatial structures become visible. This phenomenon is called Turing instability. There are other attempts to represent structures in models; a whole zoo of structure equations is available (see, e.g., the review by [Cross, M. C., Hohenberg, P. C. (1993)]. However, these approaches are often empirical and the variables are not directly linked to physical quantities. Sometimes, it is possible to reduce or reformulate an underlying physical system of equations to a known structure equation ( [Monroy, D. L. and Naumis, G. G. (2020)]). The approach of using reaction-diffusion equations is more direct, but often ignores other feedback due to the simplistic starting point. Nevertheless, reaction-diffusion equations provide an important class for pattern formation, and are directly linked to the physical variables.
In atmospheric physics, a very prominent example of emerging structures is pattern formation in clouds, which can be seen nicely from surface observation as well as obtained by remote sensing techniques (e.g. from satellites). Surprisingly, the investigation of pattern formation in clouds is not a widespread topic in atmospheric physics. There are only few studies on pattern formation, mainly in the connection with investigations of open and closed cells in marine stratocumulus (see, e.g. [Glassmeier, F. and Feingold, G. (2017)] or [Khouider, B. and Bihlo, A. (2019)]); however, rigorous and theoretic investigations are lacking. This is surprising, since internal structures of clouds constitute a serious uncertainty in terms of radiative feedback. Radiative transfer in homogeneous media is completely different than in inhomogeneous media. For the investigations of Earth's energy budget, clouds play a major role due to scattering and reflection of sunlight as well as trapping infrared radiation by absorption and re-emission. In structured clouds, many assumptions of radiative transfer in homogeneous media do not work anymore; for instance multiple scattering occur frequently, and horizontal radiative transport becomes more important. Thus, in this respect the investigation of structured (i.e. inhomogeneous) clouds and their origin and evolution is quite essential for meaningful estimations of cloud radiative forcings.
There is another difficulty concerning the representation of cloud pattern in models. Clouds constitute an ensemble of many water particles. In cloud physics, we often consider processes on particle scale, which are partly understood until now. The description of the statistical ensemble of cloud particles, forming the macroscopic "object" cloud is not very precise, and is lacking a kind of rigorous description and formulation. There are some attempts based on Boltzmann-type evolution equations (see, e.g., [Morrison, H. et al. (2020)], [Beheng, K. D. (2010)], however there is no general theory of clouds and no basic set of equations as a common ground to start is available, as e.g. something comparable to the Navier-Stokes equations describing the underlying motion of dry air. For the description of clouds, often averaged variables as the number concentration or mass concentration of particles are used. It is possible to relate these quantities to general moments of the underlying size/mass distribution of the particle ensemble ( see, e.g. [Beheng, K. D. (2010)], [Khain, A. P. et al. (2015)]. For these averaged variables, the process rates of cloud processes are formulated as ordinary differential equations (ODEs), often using nonlinear functions for representing the reaction terms. Since the basic theory is lacking, the process rates are different for many available cloud models, and often the formulations are not mathematically consistent. For instance, the uniqueness of solutions is not always guaranteed, and often requires a more rigorous treatment (see, e.g., [Hanke, M, and Porz, N. (2020)]. Nevertheless, these clouds models are often used, and they are useful for scientific investigations as well as for operational weather forecasts or climate predictions.
In this study, we investigate the potential of generic cloud models, as formulated in a former study ( [Rosemeier, J. et al. (2018)]), to form spatial structures. We couple the model equations with diffusion terms, i.e. Laplacians in the spatial directions; this leads to reaction-diffusion equations for cloud physics schemes, which will be investigated in terms of Turing instabilities.
The study is structured as follows: In the next section we will briefly describe the generic cloud model and the represented processes. In section 3 we present the approach of linear stability theory, leading to conditions for Turing instabilities in reaction-diffusion equations. The generic cloud model always allows a trivial equilibrium (no clouds, only rain); in section 4 we show that this equilibrium state cannot form Turing instabilities. Actually, in case of an unstable equilibrium, diffusion leads to stabilization. In section 5 we present a special case of a cloud model, which does not show pattern formation; this case contains standard cloud models. In contrast, in section 6 we show a case, which allows pattern formation from theory; in the following section 7 this cloud model is numerically simulated for a 1D and a 2D scenario, and these results are shown. We end the study with a summary and some conclusions.

Generic cloud model
We present the generic cloud model, as formulated in the former study by [Rosemeier, J. et al. (2018)]. The model represents clouds consisting exclusively of liquid water droplets, so-called warm clouds. The droplet population is divided into two different regimes, namely cloud droplets and rain drops, respectively. Cloud droplets are water droplets of small sizes (radius smaller than ∼ 40 µm), whereas rain drops are larger water particles. This separation can be seen in detailed simulations (see, e.g., [Beheng, K. D. (2010)] his figure 4), and partly in measurements. All water particles fall in vertical direction due to gravity. However, since for small droplets the fall velocities are very small due to friction of air, we can assume that these droplets are stationary, in contrast to large rain drops, which fall out faster. This separation was firstly proposed by [Kessler, E. (1969)] in an ad hoc measure; however, it might be justified by the simulations mentioned above. We consider the mass concentrations of these two populations as variable in the model. The phase transitions (water vapour vs. liquid water) are guided by the saturation ratio S l = pv ps(T ) , i.e. the ratio of partial pressure of water vapour, p v and its temperature dependent saturation vapour pressure, p s (T ). Thermodynamic equilibrium, i.e. coexistence of gaseous and liquid water, is then fulfilled at S l = 1; for simplification of the notation, we also introduce the supersaturation S := S l − 1, i.e. equilibrium is reached for S = 0. For liquid clouds, the following processes must be taken into account.
• Condensation and diffusional growth/evaporation Cloud droplets are formed at thermodynamic conditions slightly beyond thermodynamic equilibrium, i.e. at supersaturation (S > 0); actually, aerosol particles are activated and after passing a critical size, as given by Köhler theory ( [Köhler, H. (1936)]), they constitute cloud droplets. In simple cloud models, this process of condensation is simplified and represented together with diffusional growth. Cloud droplets can grow or shrink by uptake or evaporation of water vapour, which is provided by diffusion; this process is also driven by the supersaturation, which controls the thermodynamic equilibrium. Diffusional growth is quite inefficient for large droplets, thus this process is only relevant for small droplets, i.e. for the cloud droplet category. Both processes, condensation and diffusional growth (or on the contrary evaporation for S < 0) are represented by a rate C = c Sq c with a suitable constant c depending on temperature and pressure only. For simplification, we assume in our investigations a permanent source of supersaturation, e.g. driven by a vertical upward motion. Thus, we can also neglect evaporation (of rain drops) for the system.

• Collision processes
Since water particles fall with different velocities as depending on their masses, there will be collision between neighbouring particles of different size, and these particles will form a single droplet after collision (so-called collision-coalescence). Because of the artificial splitting into two categories, we have to consider two processes: (1) Two cloud droplets collide and form a large rain drop; this process is called autoconversion.
(2) A large rain drop collect a small cloud droplet by collision; this process is called accretion. These processes are usually modelled in the spirit of population dynamics, using nonlinear terms; however, derivations from integrals over size/mass distributions lead to similar descriptions (see, e.g., [Seifert, A. and Beheng, K. D. (2006)], [Beheng, K. D. (2010)]. Autoconversion can be represented by terms A 1 = a 1 q γ c with a suitable constant a 1 and an exponent γ > 0. For accretion, the terms can be formulated as A 2 = a 2 q βc c q βr r with a suitable constant a 2 and exponents β c , β r , mimicking a generalized predator-prey process.

• Sedimentation of particles
For the representation of rain drops falling out of a cloud level, in general we would have to consider a hyperbolic term in the vertical direction. However, for simplicity, we assume just one atmospheric layer with a prescribed vertical extension. Thus, we discretize the hyperbolic term and assume a constant flux of mass from above. Then, the sedimentation term can be approximated by D = B −dq ζ r , with constants B, d > 0 and an exponent ζ > 0. Note, that terminal velocities of cloud particles can be parameterized by power laws (see, e.g., [Seifert, A. et al. (2014)]).
Using the representation of the processes as stated above, we obtain the generic cloud scheme as described in [Rosemeier, J. et al. (2018)]: To simplify the notation we will write c instead of c S in the remaining of the paper. For the analysis of the equations, we assume constant environmental conditions, i.e. constant temperature, pressure, and supersaturation (S > 0), respectively. This assumption leads to an idealized situation, however it could be shown that similar conditions can be encountered in the atmosphere (see, e.g., [Korolev, A.V. and Mazin, I.P. (2003)]) for quite long times. Assuming these constant conditions allows us to investigate the asymptotic states of the system.
The ODE system (1) was discussed in detail in [Rosemeier, J. et al. (2018)]. In the presented work the equations (1) are extended by diffusion terms, so we obtain the following system This is a reaction-diffusion system (or Turing system). Note, that the added diffusion terms do not represent molecular dynamics, as in chemical systems. Actually, these terms can be seen as a representation of unresolved (dynamical) processes, as e.g. small eddies or turbulence. For the representation of turbulence in subgrid scale schemes or entrainment due to unresolved eddies, often gradient terms are used (see, e.g., [Deardorff, J. W. (1972)], [Stull, R. B. (1988)]). This approach leads to diffusion terms in the equations for the mean variables. The different values of the diffusion constants for the two water species can be motivated as follows: Small clouds droplets will mainly follow the small scale motions in the system, thus the diffusion coefficient D 1 for this species should be large. On the other hand, rain droplets are mostly accelerated by gravity, thus they are less affected by small scale motions. For this species, the diffusion coefficient D 2 can be chosen different from the coefficient D 1 , e.g. we would assume D 2 < D 1 .
In the sequel this system (2) is investigated with respect to pattern formation. The occurrence of patterns cannot be guaranteed for the generic model, i.e. for all possible choices of parameters, but in some cases linear stability analysis predicts pattern formation. These findings can be confirmed by numerical simulations.

Linear stability analysis
The ideas of linear stability analysis (e.g. [Turing, A. M. (1952)]) can be used for determination of stable and unstable modes of the system of equations. A classical example for the analysis of reaction diffusion equations using linear stability analysis is the investigation of the Brusselator as a simple system describing chemical reactions (see, e.g., [Cross, M. and Greenside, H. (2009)] pp. 105-108). In this section we mostly follow the exposition given in [Cross, M. and Greenside, H. (2009)] for a 2D system of reaction diffusion equations, as, e.g., given by (2). The subsequent 2D reaction diffusion system is giveṅ In a first step, we determine the stationary and homogeneous equilibrium states, thus we omit the diffusion terms. By neglecting the Laplacians, we obtain a system of ordinary differential equationṡ The right hand side is called the reaction term. We want to derive conditions for a stable equilibrium of (4) which can be destabilized by diffusion terms. First, we consider an equilibrium solution u e1 , u e2 of the system (4). By definition it satisfies the equations Next we compute the Jacobian of (4) evaluated at the equilibrium solution u e1 , u e2 where the entries of the matrix are determined by The (potentially complex) eigenvalues of the Jacobian at the equilibrium state are denoted by σ 1 , σ 2 The equilibrium solution u e1 , u e2 is asymptotically stable if and only if the following relations are fulfilled This is equivalent to the more common condition for asymptotic stability, i.e. Re(σ i ) < 0, i = 1, 2 Now we consider the system (3) including the diffusion terms. For this purpose, we use spatial coordinates x = (x 1 , . . . , x n ) T and a generalised wave number vector k = (k 1 , . . . , k n ) T . For each spatial direction, we consider linear waves with wave lengths The spatial dimension is given by n ≥ 1.
For the linear stability analysis of the reaction diffusion system, we replace the reaction term by its linearization evaluated at u e1 , u e2 , i.e. with the linearization u = u e + u p , u p small perturbation around the constant equilibrium state u e we obtaiṅ We want to derive conditions for the destabilization of u e1 , u e2 due to the diffusion terms. For simplification we assume periodic boundary conditions; therefore a Fourier discretization in space with a superposition of linear wave modes exp(ikx) can be applied. The system (9) shall be solved by a separation approach using the eigenvalues σ q . The Fourier modes are eigenfunctions of the Laplacian, leading to the following equation with the sum over all squared wave numbers q 2 = k 2 i , which is used as an index. This leads to the eigenvalue problem where the matrix Df q is given by For the determination of the eigenvalues σ q of the matrix Df q the roots of the quadratic polynomial must be determined. The eigenvalues σ qi are given by It follows that the mode which belongs to the wave number q is asymptotically stable if and only if Remember, that the condition can also be formulated in terms of determinant and trace of the original ODE system, i.e.
We are interested in conditions for the destabilization of a mode q. Equation (16a) is satisfied because equation (8a) is valid and D 1 , D 2 , q 2 > 0. The only way for destabilization is to violate condition (16b), thus we look for a mode which fulfills The left hand side of (18) defines a quadratic polynomial in q 2 , All modes q with p 2 (q 2 ) < 0 are unstable. The quadratic polynomial admits a minimum at which constitutes the "most unstable" Fourier mode. Inserting the relation (20) into (18) yields the condition or in the reformulated version There is a chance to find an unstable mode q if (21) holds, see figure 1. The conditions (8a) and (21) can be satisfied when a 11 and a 22 have opposite signs.

The trivial equilibrium of the generic cloud model
We now show that the trivial equilibrium of the generic cloud model (1) in case of a stable fix point never leads to Turing instabilities. We start with the generic cloud model (1), leading to the fix point.
Since this implies no cloud, just rain, in the atmospheric layer, this state is called trivial equilibrium. Actually, this fix point is only valid for linear stability analysis for values of the exponents γ ≥ 1 and β c ≥ 1. Otherwise, the initial value problem with q c = 0 is not uniquely solvable, since the right hand side of the ODE is not Lipschitz continuous. For discussions of such cloud models and possible extensions to uniqueness, see the recent study by [Hanke, M, and Porz, N. (2020)]. For linear stability analysis we have to consider the Jacobian Df | (qce,qre) . As discussed by [Rosemeier, J. et al. (2018)] the Jacobian always has the form whereas a 22 = −dζ B d ζ−1 ζ < 0. Actually, for γ > 1 and β c > 1 we obtain a 12 = 0, otherwise a 12 > 0; for details, see calculations in appendix A. Nevertheless, it is clear that the eigenvalues σ i are given by σ 1 = a 11 , σ 2 = a 22 < 0 and thus det(A) = a 11 a 22 = σ 1 σ 2 , and Tr(A) = a 11 + a 22 = σ 1 + σ 2 , respectively. For a stable fix point, both (real) eigenvalues must be negative, leading to the criteria (8); this might be full-filled for the choice of parameters c < a 1 in case of γ = 1; otherwise the fix point cannot be stable (see also appendix A). However, the stable fix point can not lead to Turing instabilities via destabilisation. The criterion for the existence of destabilisation (21) can be reduced to the following form:

unstable modes
Abbildung 1: Quadratic polynomial (19). No unstable modes when the minimum of p 2 is positive (blue line). Unstable modes are possible when the minimum of p 2 is negative (orange line).
Since a 11 = σ 1 < 0 and a 22 = σ 2 < 0, this leads to a contradiction. This proves that the trivial fix point (if it exists) cannot be destabilized by diffusion, and thus it cannot serve for Turing instabilities. From a physical point of view, in this situation the source for cloud droplets as represented by the term cq c is too weak, collision processes (A 1 , A 2 ) reduce the cloud water such strongly that diffusion cannot change the quality of the stable fixpoint (i.e. no cloud with rain).
If the parameters c, a 1 are chosen such that a 11 = σ 1 = c − a 1 > 0, an unstable equilibrium state can be obtained. Physically, this means that condensation and diffusional growth is much stronger than autoconversion, i.e. more cloud mass is generated by condensation than is lost by collision processes. This situation usually occurs in scenarios with a persistent updraught, leading to a steady source of supersaturation and thus cloud droplet formation/growth. In case of an unstable trivial equilibrium, diffusion might lead to stabilisation of the system: the modified Jacobian at the equilibrium state looks like The second eigenvalue σ q2 = a 22 − D 1 q 2 is still negative, since D 1 , q 2 > 0. The first (real) eigenvalue can be negative for Fourier modes q fullfilling the following condition Thus, the absolute value of the diffusion constant D 1 > 0 decides about the stability of the modes.

A case without destabilization
We set γ = 1 and β c = 1 in the system (1) and show that it is not possible to destabilize an asymptotically stable equilibrium of this model by diffusion terms with arbitrary coefficients D 1 , D 2 > 0. The special cloud schemes of the operational numerical weather prediction model COSMO ( [Doms, G. et al. (2011)]) of the German weather service (DWD) and the research model by [Wacker, U. (1992)] admit this special form of the cloud scheme examined in the sequel. Particularly we consider Beside the trivial equilibrium state (see discussion in section 4) a non-trivial equilibrium of the system (28) is given by Note, that for the existence of this (non-negative) equilibrium state two conditions must be fullfilled, i.e.
c > a 1 , and d c − a 1 a 2 ζ βr > B Physically, this means that as before the cloud droplets source cq c is stronger than the sink of autoconversion. Additionally, the rain flux from above B must not be too strong, otherwise no equilibrium state is reached, i.e. the rain will collect almost all cloud droplets. Again we compute the Jacobian at the equilibrium state where a 11 = c − a 1 − a 2 q βr re = 0, (32a) a 21 = a 1 + a 2 q βr r,e = c > 0 (32c) As a 11 = 0, condition (8a) gives a 22 < 0; note, that condition (8b) is fullfilled. On the other hand in this case, condition (21) reduces to This yields a contradiction. Thus, for schemes with γ = β c = 1 pattern formation via Turing instabilities is impossible. Generally, it is of interest if and when the matrix entry a 22 changes its sign, since Turing instabilities can only occur if a 22 > 0. For this purpose we further investigate a 22 in detail For the sign of a 22 the following term must be examined We can now conclude that for β r ≤ ζ we do not predict a pattern forming system. This condition holds for the Wacker and COSMO schemes. The trivial equilibrium state is given by It is unstable for the COSMO and Wacker schemes, and several modes q might be stabilized by the diffusion terms, depending on the diffusion constant D 1 . However it is not unstable for an arbitrary choice for the prefactors. The stable case was already discussed in section 4, it cannot trigger Turing instabilities.

A cloud scheme with pattern formation
In this section we show an example for a cloud scheme of the form (1) which allows pattern formation via Turing instabilities. Again the method described in section (3) is applied. We consider the cloud schemeq The exponents are integers; this permits analytical investigations of the quality of steady states. Note, that the constant rain flux from above (i.e. the term B) is omitted for simplification; we will discuss the inclusion of this term in section 7 for a special set of parameters. Actually, the general qualitative behaviour of the ODE system does not change, if this term is included, but it changes the conditions for the occurrence of Turing instabilities. We first show that system (37) has a stable equilibrium which can be destabilized by diffusion terms, with an appropriate choice of the parameters a 1 , a 2 , c, d, D 1 , D 2 > 0. The corresponding reaction-diffusion system has the formq We start the analysis by computing the non-trivial equilibrium states of (37) which satisfy The mass concentration of cloud droplets q c can be expressed as a function in q r Adding (39a) and (39b) and applying (40) yields which is a cubic equation in q r using the unique real root of eq. (41), we end up with the equilibrium state, which is given by With the equations (42) we obtain the first constraint on the parameters, since q ce and q re should be positive, thus c > a 1 .
As before, condensation is dominant over autoconversion of cloud droplets.
In the next step the Jacobian Df is computed and evaluated at the equilibrium solution (q ce , q re ) a 11 (q ce , q re ) a 12 (q ce , q re ) a 21 (q ce , q re ) a 22 (q ce , q re ) .
For allowing Turing instabilities, the conditions (8) and (21) must be fullfilled. Especially, the equilibrium state must be stable (i.e. condition (8) is valid). For this purpose, we have to find constraints for the parameters a 1 , a 2 , c, d of the model. Furthermore condition (21) can be satisfied when a 11 and a 22 have opposite signs.
The relations (45a), (45b) and (45c) are a consequence of (43). For a positive sign of a 22 , eq. (45d) can be reformulated as the condition This a further constraint on the prefactors. Additionally, the trace of J must be negative when (8a) is supposed to hold, i.e.
a 11 (q ce , q re ) + a 22 (q ce , q re ) = a 1 − c When d is chosen small we obtain for the argument of the square root of eq. (21) a 11 (q ce , q re )a 22 (q ce , q re ) − a 12 (q ce , q re )a 21 (q ce , q re ) = 3d(c − a 1 ) > 0.
Thus with an appropriate choice of d according to (48) we can satisfy the conditions (8), and thus (21) can be fullfilled for a proper choice of the diffusion constants D 1 , D 2 . In summary, we have derived three limiting conditions (43), (46), and (48), which are illustrated in figure 2. For values of c and d in the blueish domain of the parameter space, the equilibrium state is stable and allows Turing instabilities, generally. From eq. (48) as well as from the phase diagram in figure 2 we see that the strength of the sedimentation (parameter d) plays a major role for the existence of Turing instabilities. If sedimentation is too strong as compared to condensation and collision terms, diffusion is not effective enough to distribute the cloud spatially for generating instabilities. Remarks: Abbildung 2: Constraints on the prefactors. We assume that a 1 is given.
1. If we investigate the first equation of the generic ODE System (1a) we realize that there is a conserved quantity. If the condensation term in the original formulation is slightly extended to the final equation i.e. the condensation and the autoconversion term have the same exponential behaviour, then we observe a conserved quantity in the system for the steady state dqc dt = 0 This combination of variables q c , q r is always constant, i.e. it only depends on the chosen parameters for condensation c and collision a 1 , a 2 , respectively.
2. In case of γ = 1, we obtain for the conserved quantity in equation (51) This leads to a strong simplification of the first column in the Jacobi matrix Df , since both entries become constant (a 11 = const., a 12 = const.). This property is true for all values of β r > 0. Thus, the decision about the Turing instability in equation (21) is reduced to the investigation of the sign of a 22 , since det(J) < 0 and a 11 = const. < 0.
3. Including the rain flux from above B into the first equation of the ODE system leads to the modificationq Thus, the determination of the non-trivial fixpoint becomes more complicated. Following the strategy at the beginning of the section, adding the two ODE equations leads to finding real roots of the general cubic polynomial Using Cardano's formulas (see, e.g., [Bosch, S. (2018)] chapter 6), we can determine the roots of the polynomial directly using the following terms: whereas the parameter ∆ decides about the quality of the roots (e.g. one real root and two complex conjugates for ∆ > 0). One real root is given by Since the sign of parameter ∆ decides about the number of real roots, there is a general bifurcation at B 1 , which can be calculated using equation (57): However, the condition (21) for the existence of a Turing instability might be violated at different values of B as can be seen in the numerical simulations in the next section. The equilibrium states can be inserted into the Jacobi matrix for determining the eigenvalues. However, the terms become too complicated, thus it is not possible to determine analytically the condition for Turing instabilities in dependence of all parameters. Using the conserved quantity above, we have only to consider the entries a 12 (B), a 22 (B) for the possibility of Turing instabilities. The sign of a 12 (B) is again the key parameter, and depends on the rain flux B. Actually, for a certain setting of parameters c, a 1 , a 2 , d we will determine the qualitative behaviour and the possibility of Turing instabilities numerically (see next section).
4. The choice of exponents β c > 1, β r > 1 in the accretion term is motivated by existing models as e.g. the operational weather forecasting model IFS ([IFS DOCUMENTATION]), which contains such exponents. However, since the representation of collision processes in bulk models is not well-defined from a basic theory, there is actually no restriction of the choice of parameters. Maybe the parameters also vary for different environmental regimes.

Numerical simulations of cloud pattern
Setup: We carry out 1D and 2D numerical simulations for investigating the special case of a cloud model allowing Turing instabilities, as discussed in section 6. A pseudo-spectral method was applied for the numerical solutions of the system, see appendix B. We choose a domain length of L = 50 for the 1D case and a quadratic domain with length of L = 50 for the 2D case, respectively. In both cases, the domain is cyclic as assumed in the linear stability analysis above. For the simulations, we specify the parameters as follows: a 1 = 1, a 2 = 1, c = 5, d = 0.1, D 1 = 10 3 and D 2 = 10 −1 . In this scenario, the non-trivial fixpoint is asymptotically stable and the parabola p 2 defined by equation (19) is negative for wave numbers q 2 = 4π 2 L 2 n 2 with n ∈ {2, . . . , 7}. Therefore, these modes give rise for linear instability and thus lead to Turing instability. As initial condition we prescribe the equilibrium of the ODE system with spatial, normally distributed perturbations with amplitude of order 0.01. In a first step, the system (37) is simulated, i.e. there is no rain flux from above (B = 0). In a second step we will discuss the impact of the rain flux on the pattern formation in the simulations.
Results of 1D simulations without rain flux (B = 0): First, we investigate the numerical simulations in one spatial dimension. In figure 3 the time evolution of the two variables q c (left panel) and q r (right panel) are shown. The x-axis represents the spatial extension of the 1D domain (with cyclic boundary conditions), The spatial structure is forming out of the noise, i.e. the destabilized modes suddenly grow to larger sizes until they are saturated (and thus stopped) by the nonlinear terms. Their spatial distribution slightly changes during time until at around t ∼ 1200 the situation is consolidated, i.e. the pattern stays quite stationary. In figure 4 the simulations at times t = 20/200/2000 are shown. Here, the evolution can be seen clearly as well as the final "wavy" structure at t = 2000. Note, that the variables q c and q r have contrary behaviour: for high values of q c the rain variable q r is quite small and vice versa. This can be explained by the collision terms, which act as in a generalized predator-prey system. If the predator population (i.e. the rain) is small, the cloud can survive and grow to larger values due to condensation only, since autoconversion is weak. If the rain becomes larger, it reduces the prey (the cloud water) due to collisions.
Using Fourier analysis (not shown) we see that only a part of the Fourier spectrum has reasonable amplitudes, whereas higher modes are of very low amplitudes. However, we do not see the distinct spectrum as predicted by linear stability theory. The reason for this is the nonlinear interaction of the different modes, which smears out the modes and form a more continuous spectrum in contrast to the predicted discrete spectrum. Nevertheless, we see that only a small part of the Fourier spectrum is present in the simulations.
Results of 2D simulations without rain flux (B = 0): In a second simulation we use the 2D setup with white noise as before to investigate pattern formation in a 2D domain.
Qualitatively, we see the same behaviour for the 2D simulations of a quadratic domain of length L = 50 with cyclic boundary conditions. After a short time, the simulation leads to growing unstable modes, which are then saturated by nonlinear terms in the model; these modes form spatial structures, which change only slightly over time until they stay stationary. Thus, pattern formation due to Turing instabilities can be observed as expected from theory. The structures in cloud water q c are less pronounced than in the rain water q r . Nevertheless, the spatial pattern remain stationary, even for longer times.
Results of 1D simulations for including the rain flux (B > 0): In a last series of simulations we investigate the impact of the rain flux B, which was set to zero in section 6 for simplification of the analysis. Using the fixed parameters c, a 1 , a 2 , d we can calculate the roots of the cubic polynomial determining the fix points. One real root can be calculated by q re = u + v + B 3d . The bifurcation value can be calculated as B 1 ∼ 1.10521. Actually, we additionally find that the eigenvalues σ i , i = 1, 2 are always negative for 0 < B < B 1 . Thus, the fixpoint is always asymptotically stable for rain flux from above in the relevant parameter range.
As described in section 6, the entries a 11 < 0, a 21 > 0 of the Jacobi matrix are constant, the entries a 12 (B), a 22 (B) depend on the rain flux. For fullfilling the criterion for Turing instabilities (21), the sign of Abbildung 4: Spatial variation of the variables q c (top row) ans q r (bottom row) for times t = 20 (left), t = 200 (middle), and t = 2000 (right), respectively. Note the different scaling of the y-axis. Actually, at t = 20 there is almost no variation of q c , q r visible, whereas at t = 2000 the change in q c , q r is obvious.
Abbildung 5: Spatial distribution of cloud water q c in 2D for different simulation times (t = 1/10/60/120). Note, that the pattern is already forming at times t ∼ 10. For longer times, the pattern stays stationary although the absolute variation of the cloud water variable is very small over the whole 2D domain Abbildung 6: Spatial distribution of rain water q r in 2D for different simulation times (t = 1/10/60/120). The spatial structure for this variable is more pronounced than for the cloud water, i.e. the spatial variation of q r is quite large. We confirm these findings with a series of numerical simulations using different values 0 < B < 0.17 for the set of parameters as specified at the beginning of the section. As predicted, for values B < 0.137 we find Turing instabilities, whereas for B > 0.137 there are no Turing instabilities. In Figure 7 the simulations at time t = 2000 (i.e. steady state) depending on the parameter B are shown. The absolute values of the pattern in q c and q r slightly vary with changing B; however, the quality of the pattern remains the same until values B ∼ B 2 are reached. Passing this values, a homogeneous state in both variables can be seen and no pattern formation occurs. Note, that the boundary B 2 is not sharp, since the occurrence of the Turing instability depends also on the values of D 1 , D 2 . Due to the large ratio of these coefficient, the transition in the simulations is very close to B 2 .

Summary and conclusion
In this study we investigate a generic cloud model for pure liquid clouds about the possibility of Turing instabilities for forming spatial pattern. This kind of investigation is carried out for the first time for cloud pattern formation. For the theoretical and numerical investigations, the generic cloud model as formulated in the former study [Rosemeier, J. et al. (2018)] is extended by diffusion terms, consisting of Laplacians in spatial directions. The model is analysed using linear stability theory around the steady states of the underlying ODE system. The model is a two-dimensional ODE system; analytical conditions for the existence of Turing instabilities can be determined. Since the model contains quite complex nonlinear terms with several parameters determining the overall quality of the steady states, it is very hard to find general conditions for the existence of Turing instabilities. However, the generic model always admits a trivial fix point in terms of "no clouds, just rain falling through the layer". This fix point could be either stable or unstable, depending on the set of parameters. However, even in the stable case, the steady state cannot be destabilized by diffusion; thus, this state does not admit Turing instabilities. In addition, we can specify a class of models, which do not allow Turing instabilities at all. This class is characterized by a linear autoconversion term A 1 = a 1 q c and also a linear contribution of cloud water q c in the accretion term A 2 . Well-known cloud models as the standard COSMO cloud scheme ( [Doms, G. et al. (2011)]) or the research model by [Wacker, U. (1992)] belong to this class. However, we can also provide an example of a cloud model, which allows Turing instabilities. If the exponents in the accretion parameterisation are chosen to be larger than 1 (in this case β c = β r = 2), the model allows Turing instabilities and thus pattern formation. These theoretical findings can be supported by numerical simulations in one and two spatial dimensions. The inclusion of rain flux from above turns out to be an additional restriction for the instabilities. If the rain flux becomes too large (i.e. if it surpasses a certain threshold B > B 2 ), the criterion for the existence of Turing instabilities is violated, as can be seen also in a series of 1D numerical simulations. This observation leads to the interpretation that collision processes in combination with the sedimentation of cloud particles play the major role for pattern formation; only if these processes can interact in a proper nonlinear way, Turing instabilities are possible. A strong rain flux from above can prevent the formation of cloud pattern; this might be explained by stronger collision terms, which finally almost extinct the cloud droplet population, thus diffusion can not counteract this process.
We can conclude that the generic cloud model admits Turing instabilities in special cases. However, several standard cloud models as used in research and operational weather forecasts do not admit pattern formation due to Turing instabilities. It is still unknown how pattern in clouds form, especially, which processes lead to the emergence of cloud structures. The use of diffusion terms is motivated by parameterisation of subgrid scale processes; maybe this approach is too simple for representing the underlying processes in a meaningful way. On the other hand, it might be that pattern formation in clouds is dominated by Turing instabilities; in this case it might be a major drawback to use cloud models, which do not allow this type of pattern formation. Since pattern formation in clouds is far away from being understood, this has to be taken into account, which might have impact on the choice of cloud models for further investigations. For the investigation of cloud pattern, more theoretical studies are needed for a better understanding of the underlying processes and their interaction, which leads to the emergence of cloud structures.
Remember, that the Jacobian of the trivial state is only defined for values γ ≥ 1, β c ≥ 1; otherwise the initial value problem is not solvable, since the right hand side of the ODE system (1) is not Lipschitz continuous. For the entries in the matrix (61), we have to discriminate between different cases. First, we determine the value of a 11 .
• For γ = 1 we obtain a 11 = c − a 1 . In this case, the trivial fix point can be stable (c < a 1 ) or unstable (c > a 1 ) • For γ > 1 we obtain a 11 = c. In this case, the trivial fix point is always unstable.
Second, the entry a 21 is investigated.
In any case, the entry a 21 does not affect the stability of the trivial fix point.

B Pseudo-spectral method
The pseudo-spectral method is applied to the following type of equationṡ x = L(x) + R(x). (62) As the system 2 is considered, The linear operator L and the reaction term R admit the form L(q c , q r ) = c + D 1 ∇ 2 q c D 2 ∇ 2 q r and R(q c , q r ) = −a 1 q γ c − a 2 q βc c q βr r a 1 q γ c + a 2 q βc c q βr r − dq ζ r + B .
For the spatial discretization a Fourier expansion applies y = ϕ(t) exp(ik n x).
Thus we obtain a system of ordinary differential equationṡ ϕ n (t) = l n ϕ n (t) + R n (t), where l n is the matrix The linear operator R can be expressed with Fourier modes R(q c (t), q r (t)) ≈ N 2 n=− N 2 R n (t) exp(ik n x).
The ETD2 scheme is a two step method. The first step was computed with the exponential integrator scheme of first order (ETD1 scheme) u m+1 = u m exp(kτ ) + N m exp(kτ ) − 1 k .