Advances in Geosciences Model of oscillatory instability in vertically-homogeneous atmosphere

Existence and repeatability of tornadoes could be straightforwardly explained if there existed instability, responsible for their formation. However, it is well known that convection is the only instability in initially stable air, and the usual convective instability is not applicable for these phenomena. In the present paper we describe an instability in the atmosphere, which can be responsible for intense vortices. This instability appears in a fluid with Coriolis force and dissipation and has oscillatory behaviour, where the amplitude growth is accompanied by oscillations with frequency comparable to the growth rate of the instability. In the paper, both analytical analysis of the linear phase of the instability and nonlinear simulation of the developed stage of the air motion are addressed. This work was supported by the RFBR grant no. 09-05-00374-a.


Introduction
Atmospheric phenomena like tornado and tropical cyclone have been intensively studied, however, their basics is not clearly understood.There exists a great amount of recorded data on these phenomena; however theory of these phenomena is not yet developed.From observations one knows, that tornado is a rotating funnel-shaped cloud which stretches down from the storm-cloud base.Tornadoes can either reach the earth surface which stops further propagation, or can have a fixed shape in the air without visible limitations.Despite of their moderate dimensions they are considered as the most dangerous atmospheric phenomena.It is known from observations that tornado appears in initially motionless air.Convection is the only known instability in initially motionless air; however in the simplest case the convective instability it is not applicable for description of tornadoes.
Correspondence to: P. B. Rutkevich (peter@d902.iki.rssi.ru) The quest for theoretical modelling of tornadoes has a long history with remarkable number of publications (see, for example, reviews by Doswel and Burgess (1993) and by Bengtsson and Lighthill (1982)).A noteworthy work of Zavolzhenskii (2002) describes sinking process of a rotating vertical column into underlying heavy fluid under action of tangent stress at the upper surface.The results show, that the equilibrium state of rotating layer of viscose fluid losses its stability under temperature inversion (i.e.stable) conditions.Formation of the tornado column is a result of complex interactions between different components of air velocity and thermodynamic parameters, and can be a response of stably stratified air to vertical or horizontal air motion, which may have other nature as well.For example, regular motion in thunder cloud, detected by Doppler radars (Rasmussen et al., 1994;Wurman, 2002), are commonly accepted to appear due to the stream flow in the upper layers of the atmosphere.Morgulis and Yudovich (2001) described ideal liquid flow through a two-dimensional area, assuming the normal velocity and curl of the velocity at the boundaries are not equal to zero.They have shown that there is no stationary or oscillatory solution in such flow.Instead, all solutions of the problem infinitely grow.In other words, in such system the inflowing fluid flux generates unlimited rotation of the fluid.
In 2002 authors also investigated a problem of a stationary rotating structure in a linear approximation.The development of the structure was assumed to be associated with air motion in the cloud layer above the system and was introduced into the problem as an external parameter responsible for boundary conditions between the cloud and the stable air below the cloud.Mathematically the problem was based on hydrodynamic equations for stably stratified air with inhomogeneous boundary conditions at its upper edge.The rotation velocity (i.e. the amplitude of the corresponding Bessel mode) of the system appears to be inversely proportional to the viscosity, and since the air viscosity is small the toroidal component of the velocity field becomes predominant.P. B. Rutkevich and P. P. Rutkevych: Model of oscillatory instability in vertically-homogeneous atmosphere In this view, investigation of the full system of hydrodynamic equations is necessary.In the present paper we describe unstable solution of linearized hydrodynamics equations, which can be responsible for tornado formation.We analyze several particularly interesting cases, related to real atmosphere conditions.Numerical solution of nonlinear hydrodynamic equations shows that the amplitudes of the parameters become saturated by Bessel modes of higher order.

Oscillatory instability
Typically the term "instability" describes the fact, that a solution of linear equation has an exponentially growing solution.Usually instability appears through interaction of two fields that appear in a feedback loop, increasing each other, like convection instability.See, for example, Landau and Lifshitz (1987).In the simplest case, there exists a field y 1 (t) increasing another field y 2 (t), while the second field increases the first one.Mathematically this can be described by the following system: here we assume the coefficients a and b are positive, and primes in the left-hand side of the equations stand for the time derivatives.From Eq. ( 1) one straightforwardly gets the second-order differential equation: Solution of this equation has one increasing term with growth rate √ ab.In case of convective instability in incompressible fluid (Landau and Lifshitz, 1987) the two fields involved in the feedback loop are the temperature and the poloidal velocity field (more details on definition of poloidal and toroidal velocity components will be presented in Sect.3).The instability growth rate depends on vertical temperature profile gradient, which can be negative or positive, corresponding to convection or internal wave, respectively.
Similarly, in case of three fields in a feedback loop there also can appear a situation when the existence of one field increases the second field, existence of second field increases the third, and the existence of the third field increases the first one.Mathematically this model can be described by three first-order equations, with more intricate feedback loop: where we assume a, b and c to be positive.System (4) can be reduced to the following third order equation: and the solution describes an instability with the growth rate equal to 3 √ abc.As it has been mentioned above, convection involves interaction of the temperature and of the poloidal fields in the feedback loop.The only third field that can be introduced into the system is the toroidal field (since the potential velocity field and the pressure field describe sound waves and almost never are considered separately).Therefore, in order to find a new instability we have to add the toroidal field into the system.However, for real atmosphere the toroidal field enters the system together with Coriolis parameter and the equation of the type (5) takes the form: The three independent solutions of this equation are: where i stands for the imaginary unit: i 2 =−1.One can see, that two unstable solutions (8) and ( 9) have oscillating behaviour.A similar conclusion can be drawn for higher order system of equations.

Analytical solution
One can expect that the complete system of atmosphere hydrodynamic equations of the fifth order (with respect to time) also may contain an effect shown in the previous section.The complete set of equations include Navier-Stokes equation, continuity and equation for entropy.In the initial state the thermodynamic parameters T 0 (z), P 0 (z) and ρ 0 (z) are assumed to follow adiabatic (neutral) stratification.Using representations for pressure P =P 0 +P 1 and density ρ=ρ 0 +ρ 1 , where P 1 P 0 and ρ 1 ρ 0 , and v for velocity, and eliminating the temperature T , one can obtain a linearized hydrodynamic system for v, P 1 , and ρ 1 : Here t is the time, g is the gravitational acceleration, e z is a vertically directed unit vector along the z axis, 2 is the Coriolis parameter, ν and χ are the gas viscosity and thermal diffusivity, respectively, c P and c V are the heat capacities of the air at constant pressure and constant volume, c is the sound velocity.We can use P =ρRT in the ideal gas approximation, and consider linear initial temperature profile dT 0 /dz=T 0 ( +γ ).Here =−g/(c P T 0 ) is the neutral gradient in adiabatic atmosphere.As it was already mentioned above, positive values of γ >0 correspond to internal waves, and negative values γ <0 correspond to convection.Using axial symmetry of the problem we introduce cylindrical coordinates with the z axis directed along the gravitational force.It is convenient (Lupyan et al., 1992) to use the gas velocity representation as a sum of potential, toroidal and poloidal velocity fields: here the potentials of the mentioned fields (t, r), ψ(t, r) and ϕ(t, r) are scalar functions of time and coordinates, and r=re r +ze z .From Eq. ( 13) one can show, that in axiallysymmetrical vertically-homogeneous problem the potential field corresponds to the radial velocity, the toroidal field ψ stands for the azimuthal velocity, and the poloidal field φ for the vertical velocity component.Applying (∇•), (e z • ∇×), and (e z • ∇×∇×) operators to Eq. ( 10) one can obtain equations for velocity potentials.In these variables the system (10)-( 12) becomes: where =∇ 2 is Laplace operator in cylindrical coordinates, and ⊥ is the radial part of this operator.c 2 =c P /c V •RT 0 is the sound velocity in the air.
After axially-symmetrical Hankel transforms with respect to horizontal coordinate F (r, t)=F (t) F r (r)J 0 (kr)rdr and Fourier transform for time F (t) exp(xt)dt, one obtains dispersion relation between horizontal wave number k and growth rate x: where the dispersion relation coefficients f 0 (k)• • •f 5 (k) as functions on the wave number k are presented below: Here we have introduced the ratio of specific heats κ=c P /c V , and assumed ν=χ.Equation ( 19) has been solved numerically for typical Earth's atmosphere parameters.Figures 1 and 2 show real parts of all five solutions of the dispersion relation ( 19) for different values of the temperature stratification gradient γ .One can see that for unstable negative stratification (γ <0) in Fig. 1a there is only one positive root (convection), other negative curves correspond to the second "convective" root, sound waves (two curves) and potential vorticity (fifth curve).For neutral (γ =0) stratification Fig. 1b two convective solution transforms into two unstable (oscillating) solutions.At small positive stratification in Fig. 2a the (oscillatory) instability still exists though the growth rate and wave number diminish with increasing γ , and after certain positive value of the temperature gradient all the solutions become stable (internal wave), as one can see in Fig. 2b.Sound waves and potential vorticity are almost unchanged for different γ .
It is noteworthy, that the behaviours of the dispersion curves at small wave numbers for the case of neutral and negative stratifications are different.One can see in Fig. 1a and 2a that these dispersion curves have different derivatives at k→0.Below we performed analysis of these solutions in the vicinity of k=0.
For wave number equal to zero the Eq. ( 19) essentially simplifies and becomes:   and its roots are: Detailed analysis of Eq. ( 19) with neutral stratification γ =0, shows, that in vicinity of k=0 this equation can be rewritten as One can conclude, that the expansion of the growth rate x as function of the wave number should be performed in powers of k 2/3 (see Arnold et al., 1988, for more details).The following approximate values for the growth rate can be obtained:  where we have introduced parameter G=G(c, g, κ, ): Note that G does not depend on viscosity.The wave number k max corresponding to the maximum growth rate of the instability accurate within two terms of expansion ( 21) is and the growth rate x max within the same accuracy becomes: It should be noted that the Coriolis parameter for tornado vortex does not mean the Earth rotation, but the rotation of the mother cloud, which can significantly exceed the natural Earth's rotation.The effect of the rotation (i.e. the Coriolis parameter 2 ) of the mother cloud on the growth rate is shown in Fig. 3.One can see, that the instability appears only   in the presence of rotation, and for >0.02 s −1 the growth rate is almost constant.The chosen range of ν between 10 and 100 m 2 /s is the typical value of turbulent viscosity in the atmosphere.
At large values of the Coriolis parameter the expression ( 21) can be further simplified leading to maximum wave number Using numerical values c= √ κRT =335.4 m/s for the sound velocity at T =280 K, ν=χ=100 m 2 /s for viscosity and thermal diffusivity, κ=c P /c V =1.4 for specific heats ratio, and g=9.8 m/s 2 for the gravitation acceleration, which are almost invariable in low atmospheric layers, one can estimate the maximum wave number k max (ν)=1/(23 √ ν) m −1 , that corresponds to the radius R max = 90 √ ν m.Under these assumptions the maximum growth rate of the instability takes the form: Substituting the values for the sound velocity and the gravitation acceleration yields x max ≈2.916×10 −3 s −1 .Though according to ( 24) and ( 27) the growth rate does not depend on the air viscosity, this fact appears as a result of the used two-terms approximation ( 21) and ( 25), respectively.Figure 4a compares the approximate theoretical curve ( 21) and the exact solution obtained numerically from ( 19).Another noteworthy result of the growth rate expressions ( 21   and ( 25) is that the viscosity ν appears in the numerator of the fractions in the RHS of these equations, we will discussed this later.
The used approximation of neutral air stratification makes the overall dispersion relation rather simple, though one has to use inconvenient power expansion k 2/3 series.Analytical analysis of stable air stratification γ >0 is more relevant and makes possible the expansion in terms of more natural k 2 : however, the expressions for coefficients A 1 and A 2 are far more complicated and are not presented here.The comparison of the approximation (28) and the exact numerical solution for the stable temperature profile is shown in Fig. 4b.

Effect of nonlinear terms
Linear analysis presented in the previous section properly described the beginning stage of the instability, while in its developed stage the nonlinear terms of Navier-Stokes equations must be taken into account.The obtained growth rates and typical wavelengths provide an estimation for the typical size of the structure, since the main impact into the entire system should come from the mode with maximum growth rate.At the developed stage, account of nonlinear terms is especially important, since the interaction between the Bessel modes will describe the fine structure of the motion.The nonlinear terms appear due to (v∇)v in (10), and ∇(ρv) in ( 11).Nonlinear term in equation ( 18) have been neglected, due to small deviation of temperature profile from adiabatic profile γ .In the system ( 14)-( 17) additional terms: ..
are to be added in the left hand sides of the corresponding equations.Using cylindrical symmetry of the problem we assume the solution as a series of orthogonal Bessel functions.Each mode can be described in a linear approximation, and depending on the wavelength will either decay or grow, according to (19), and this has been confirmed by numerical simulation described below.
The system of Eqs. ( 14)-( 18) has been solved numerically in a cylinder of radius R, assuming F (r, t)= i Fi (t)J 0 (k i r), where F stands for each of the fields , ψ, ϕ, P 1 and ρ 1 , J 0 is Bessel function of zeroth order, k i is the wavelength found from the boundary condition J 1 (k i R)=0, and subscript i corresponds to the i-th mode.The amplitudes F (t) have been calculated using Runge-Kutta 4-th order method with time step 10 −2 s.The results have shown that the growth rate of each single Bessel mode is in agreement with the analytical analysis (19).
It is noteworthy that there is a threshold value of the cylinder radius R th , corresponding to negligibly-small growth rate value.A solution in cylinder with radius exceeding the threshold value R>R th will exponentially grow, while any smaller radius R<R th will lead to a vanishing solution.The value of R th corresponds to the threshold k th , which is the solution of equation f 0 (k)=0 from (24).
However, one can see in Fig. 5 that accounting only one single Bessel modes with nonlinear terms does not provide stabilization of the solution.After certain time (approx.T =27 h) the amplitude growth is faster than exponential (Fig. 5 shows only the poloidal field).Therefore one can conclude, that further expansion of the modes is required.Indeed, four Bessel modes can successfully describe the entire system near the threshold radius.One can notice saturation of the amplitude growth caused by non-linear interaction with four modes taken into account in Fig. 6a.In the picture, black line corresponding to the first mode in a system of four Bessel modes is compared to the solution of linear equation for the first Bessel mode (red line).At the initial stage both solutions are almost identical, however, interaction of the first (growing) mode with other (decreasing) modes results in a balance and thus forms an equilibrium amplitude.One can see the oscillatory behaviour of the solution even in its developed stage.Figure 6b shows time dependence of the second Bessel mode.It is expectable, that the energy transfer between the modes should be directed from increasing modes (i.e.modes with small wave number k) to decreasing modes (with bigger k).In equilibrium the energy transfer rate from increasing to decreasing modes is compensated by the decay rate of the latters.Mathematically, an unlimited number of available higher-order (decreasing) modes suggests that they should be able to suppress any driving mechanism and thus form a stable configuration.

Conclusions
The presented model of growing oscillations in atmosphere shows a general possibility of existence of unstable motion in stable air stratification.Equation (19) gives the relation between the growth rate and the wave number.The fastest growing mode corresponds to the maximum of the function x(k), thus one can estimate the typical size of the formed structure.The main parameter responsible for the instability is the initial temperature profile γ .
The analyzed analytical solutions in two different approximations (γ =0 and γ >0) show the following two main features of the motion.First, the maximum growth rate is almost  independent of the viscosity, while the corresponding wave number is proportional to 1/ √ ν (further analysis shows that the maximum growth rate x max slightly increases with the viscosity, which possibly can be explained as a response to significant change of k max ).And second, the growth rate increases with the Coriolis parameter value.Obtained expression for wave number k max in neutral stratification corresponding to maximum growth rate suggests the expectable value of the unstable area R max =µ i /k max , where µ i is the i-th root of Bessel function J 1 .
In case of convectively stable temperature profile (γ >0) the situation becomes more complicated, and can be summarize as follows.The growth rate of the instability decreases with increasing absolute value of γ , and entirely disappears after certain threshold value.In the Earth's atmosphere the stratification is often near the values sufficient for the oscillatory instability, therefore this instability might be responsible for the beginning of the severe weather phenomena in real atmosphere.
We can associate the results of the oscillatory instability with the thin and tall tornado vortex, though it should not be forgotten that in the present investigation the vertical depen-dence is omitted.A real tornado in its upper part is usually funnel-shaped, and at its lower part the vortex touches the ground, therefore the vertical dependence cannot be omitted.
Since all hydrodynamic equations together with all the parameters proved to be relevant, the problem appeared to be rather complicated even in the linear stage.Numerical simulation of the nonlinear approximation of a single Bessel mode have shown that the growth rate in not limited by nonlinear terms.However, accounting several modes (we used up to four modes, three of them with negative growth rate) saturates the oscillations and shows a reasonable mechanism for the amplitude limitation, through energy transfer to decreasing modes.Similar mechanism of energy transfer is usually used in turbulence theory (Monin and Yaglom, 1967).Accounting sufficient number of decreasing modes one can expect to saturate the amplitude of the stable solution.However, accounting a big number of modes in numerical simulation requires big computational resources.

Fig. 1 .
Fig. 1.Five roots (only real parts are shown) of dispersion relation(19), corresponding to the full system of axially symmetric vertically homogeneous hydrodynamic equations (10)-(12) for negative and zero values of temperature gradient γ: a) γ = −1 × 10 −5 m −1 has one positive (red) root, b) γ = 0 has two positive (red and blue) roots.Other (negative) roots correspond to the second convective root, sound wave and potential vorticity.

Fig. 5 .
Fig. 5. Numerical simulation of nonlinear Navier-Stokes equations, one Bessel mode.One can see, that the solution tends to increase infinitely.

Fig. 5 .
Fig. 5. Numerical simulation of nonlinear Navier-Stokes equations, one Bessel mode.One can see, that the solution tends to increase infinitely.

Fig. 5 .
Fig. 5. Numerical simulation of nonlinear Navier-Stokes equations, one Bessel mode.One can see, that the solution tends to increase infinitely.

Fig. 5 .
Fig. 5. Numerical simulation of nonlinear Navier-Stokes equations, one Bessel mode.One can see, that the solution tends to increase infinitely.

Fig. 6 .
Fig. 6.Time dependence of (a) first and (b) second Bessel modes of poloidal field ϕ in a four-modes approximation.In figure (a) linear solution ϕ1,lin(t) (red) is compared to the same mode ϕ1,nonlin (black) reduced due to interaction with higher modes.

Fig. 6 .
Fig. 6.Time dependence of (a) first and (b) second Bessel modes of poloidal field ϕ in a four-modes approximation.In figure (a) linear solution ϕ 1,lin (t) (red) is compared to the same mode ϕ 1,nonlin (black) reduced due to interaction with higher modes.
Rutkevich and P. P. Rutkevych: Model of oscillatory instability in vertically-homogeneous atmosphere mind a big amount of nonlinear terms, one should consider numerical simulation as the best opportunity for further analysis.