Tornado-type stationary vortex with nonlinear term due to moisture transport

Abstract. Tornado vortex is believed to be essentially nonlinear phenomenon; and the puzzle to choose the nonlinear term(s) responsible for its formation is still unresolved. In the present work we consider the nonlinear term associated with atmosphere humidity, by introducing variable temperature gradient depending on the vertical velocity of the fluid. Such term is able to yield energy to the system and is very suitable for such a problem. Other nonlinear terms are neglected, assuming slow rotation, or in other words a "weak" tornado approximation. We consider one-dimensional radial boundary problem, and use a modificaiton of shooting method to satisfy boundary conditions at large radii. Obtained numerical solutions of the nonlinear differential equation qualitatively agree with the observed atmosphere vortices (tornados, tropical cyclones). The obtained results show general possibility of existence of unstable motion even in convectively stable atmosphere stratification.


Introduction
Formation mechanisms and internal structure of severe atmospheric vortices have been important questions over the past decades.In spite of their importance for both theory and applications, there is no satisfactory theory of tornado structures (Doswell and Burgess, 1993) and recently there is a tendency towards collection of observational data or numerical simulation of the phenomenon (Bluestein et al., 1997;Bluestein and Pazmany, 2000;Bluestein and Weisman, 2000;Lehmiller et al., 2001).
The problem of tornado formation appears to be complicated from both principal and technical points of view.Tornado descends from a long-living rotating thunderstorm cloud.Having appeared from the cloud the funnel narrows in its upper part and dynamically stretches in the bottom part, however the middle part of the tornado column remains cylindrical and almost unchanged during the entire lifetime of the funnel.Such structure is especially pronounced for weak tornadoes (Davies-Jones, 2001), i.e. for tornadoes with high ratio of the height to the diameter.Mathematical model of the middle part of tornado should have a steady-state solution or solution with very low growth rates.
Correspondence to: P. B. Rutkevich (peter home@gmail.com)It is commonly accepted that tornado solution cannot be derived from linear system of equations (continuity, Navier-Stokes, entropy) and therefore it is believed that tornado is described by a set of nonlinear hydrodynamic equations.However, it is still unclear which non-linear processes are responsible for its formation.Indeed, the inertia non-linear terms (v∇)v lead to energy dissipation, therefore cannot be the energy source of the system.Other non-linear terms have the same problem, since they do not introduce energy in the system.Renno and Ingersoll (1996) described convection in the atmosphere as a heat engine between warm surface (due to solar radiation) and cold troposhpere, this assumption implies heat accumulation and release by means of black-body radiation.However, continuous and powerful source of energy observed in tornadoes and cyclones has to come from a source much more efficient than the black body radiation or heat conductivity.One of possibilities is to use latent heat of water vapor as the source, since vapor content in the air is widely spread in the atmosphere, phase transitions occur fast, and the specific condensation energy is high.In this work we imply a mechanism of the latent heat release by varying the vertical temperature profiles in different parts of the tornado structure.Unlike a heat engine, the process is not reversible.Indeed, at different points the tornado has either upward or downward air flows.These flows naturally contain saturated vapor from the surface or dry air from the Published by Copernicus Publications.P. B. Rutkevich and P. P. Rutkevych: Tornado-type stationary vortex with nonlinear term due to moisture transport troposphere, respectively.The upwards fluid flow transports the vapor, which constantly condenses and releases the heat at all levels.Therefore the vertical temperature profiles in these two flows are essentially different.As it will be shown later, such system has a stable configuration with strong vertical and azimuthal velocity components.

Formulation of the problem
We start from the full system of hydrodynamic equations for compressible fluid, consisting of continuity equation, Navier-Stokes equations, and entropy balance (Landau and Lifshitz, 1987): where thermal diffusivity coefficient χ=k/(c P ρ) is used instead of the thermal conductivity k, and c P is specific heat capacity of air at constant pressure; ν is the kinematic turbulent viscosity, e z Ω is the vector of angular velocity of the rotating mothercloud, and v is the fluid velocity.Although we describe compressible fluid, we have neglected in Eq. ( 2) the term with bulk viscosity for simplicity.We introduce cylindrical coordinates with unit vector e z directed upward.Following the standard linearization procedure (Landau and Lifshitz, 1987), density, pressure and temperature are represented as ρ(r,z)=ρ 0 (z) + ρ 1 (r,z), P(r,z)=P 0 (z) + P 1 (r,z) and T (r,z)=T 0 (z) + T 1 (r,z), respectively.Here the steady-state adiabatic profiles of thermodynamic functions are denoted by index 0, and their small perturbations by index 1.The linearization procedure requires perturbations to be small compared to the steady-state values i.e. "weak" tornado.
The steady-state profiles are determined in an assumption of zero velocities, implying the temperature profiles are close to adiabatic.The corresponding set of equations contain steady-state Navier-Stokes equations ∇P 0 + e z gρ 0 = 0, the perfect gas law P 0 = ρ 0 RT 0 , and the steady-state entropy balance equation: ∇ 2 T 0 = 0. Solution of the latter equation is a linear temperature profile, which we postulate as where T 00 =T 0 (0) is the temperature value near the surface, γ a =−g/c P is a well known neutral adiabatic gradient of dry atmosphere, and γ is the deviation of the chosen steady-state gradient γ from the neutral profile γ a , g is the gravity acceleration.According to our steadystate assumption the deviation is small: γ γ a , and positive value of γ corresponds to convectively stable atmosphere.Substituting temperature profile (Eq.4) into steadystate Navier-Stokes equations one can derive the pressure P 0 (z)=P 00 [T 0 (z)/T 00 ] −g/(Rγ) ≈P 00 [T 0 (z)/T 00 ] κ/(κ−1) , and density ρ 0 (z)=ρ 00 [T 0 (z)/T 00 ] −g/(Rγ)−1 ≈ρ 00 [T 0 (z)/T 00 ] 1/(κ−1) , profiles where κ = c P /c V is the ratio of specific heats.Notice the density derivative where c(z)= √ κRT 0 (z) is the sound velocity.Linearization with the steady state Eqs.( 4)-( 5) leads to the following linearized system of hydrodynamic equations: To derive the last equation we have used the perfect gas entropy as s=c P ln(T/T 00 )−Rln(P/P 00 ), substituted here the pressure and temperature profiles, expressed the temperature through the pressure and density using the perfect gas law, and finally substituted this entropy into Eq.( 3).Now we are constructing a steady-state (∂/∂t=0), isotropic (∂/∂θ=0), and vertically homogeneous (∂/∂z=0) solution of the linearized system (Eqs.6-8).It is convenient to express velocity in terms of potential, toroidal, and poloidal fields (Rutkevich and Rutkevych, 2009): where Φ(t,r), ψ(t,r) and ϕ(t,r) are the potentials of the above mentioned velocity components.In terms of these variables the set of Eqs. ( 6)-( 8) becomes where ∂ ∂r is the radial component of the Laplace operator.The system of ordinary differential Eqs. ( 10)-( 14) describes the radial profiles of the thermodynamic parameters and velocities at any given height z=z 0 , i.e. at given density ρ 0 (z 0 ) and sound velocity c(z 0 ).
As it has been already mentioned in the introduction, only one non-linear term is considered, namely the dependence of where γ 0 >0 is the temperature gradient of stable dry atmosphere, coefficient γ 1 <0 accounts for the temperature profile in moist air, and v T is the velocity normalization coefficient.
For sufficiently high vertical velocity v z the total temperature gradient (Eq.15) becomes negative (i.e.convectively unstable).

Convection without rotation
In non-rotational (Ω→0) and incompressible (c→∞) fluid the system Eqs.( 10)-( 14) reduces to one equation for the poloidal velocity potential ϕ: Other velocity components, pressure and density can be derived from the solution of Eq. ( 16).Now we substitute the temperature gradient (Eq.15) in Eq. ( 16), express v z in terms of ϕ by means of Eq. ( 9), and obtain the following nonlinear equation or in dimensionless variables where r = λx and ϕ = C f , and Note, that according to this definition the coefficient C is negative.The unknown function f (x) has to be a regular function of its argument, therefore solution of Eq. ( 18) can be expanded in a Taylor series near x=0: where only two coefficients are independent.Indeed, one can easily check that coefficients a 4 , a 6 , etc, can be recursively expressed in terms of parameters A and a 2 using Eq. ( 18).Construction of the solution (Eq.20) for large x is computationally intensive due to fast divergence of the series, therefore we use 4-th order Runge-Kutta method to build a numerical solution of Cauchy problem for the 4-th order ordinary Eq. ( 18), using the values of function f (x) and its derivatives at point x=0.1, as calculated from the Taylor expansion (Eq.20).It appears impossible to apply Runge-Kutta method from x=0 due to singularity of the Eq. ( 18).For arbitrary values of parameters A and a 2 solution (Eq.20) diverges even for moderate x.However, the localized solution of our interest here should have a non-zero value near the axis, gradually vanish at medium distance from the axis, and uncontrollably grow far away from the axis due to numerical errors.We denote x m (A,a 2 ) the maximum value of dimensionless radius in the far-away range x > 10, where the solution is small enough Here b is a fixed number in the range 0 < b < 1, and the "faraway" radius x = 10 corresponds to r = 700 m, which we suppose sufficiently large.Let us denote A 0 and a 2,0 the values of parameters A and a 2 , corresponding to the solution (Eq.20) vanishing at x → ∞.It is clear that function x m (A,a 2 ) defined above should go to +∞, as A→A 0 , a 2 →a 2,0 .
Equation (18) has to be solved numerically, therefore the exact values (A 0 ,a 2,0 ) are hardly obtainable, however, one can approximate the solution by choosing parameters (A,a 2 ), which give large enough value for the function x m (A,a 2 ).In this case the best solution would be the solution with the biggest value of x m .In our numerical scheme the coefficients A and a 2 were varied simultaneously (see Fig. 1) to obtain maximum value of x m (A,a 2 ).The sharp peak on the plot corresponds to the best approximation of solution f (x,A 0 ,a 2,0 ) of the nonlinear equation ( 18).The position of the peak is at A≈ − 3.3641719 and a 2 ≈ − 0.265424, giving the maximum dimensionless radius x m ≈25.The corresponding solution of dimensionless poloidal potential f (x) is shown in Fig. 2. Corresponding velocity profiles (v r ,v θ ,v z ) are plotted in Fig. 3 for the following values of parameters: It should be noted here, that strictly speaking the radial and azimuthal velocity components (v r and v θ ) are not solutions of the full system of hydrodynamic Eqs. ( 10)-( 14), since in non-rotational case the system significantly simplifies.However, one can formally derive the profiles of all   velocity components from the known profile of f (x), using definition (Eq.9) and system (Eqs.10-14) as: Azimuthal component of the velocity is proportional to the rotation speed of the mothercloud.The profile of v θ in Fig. 3 has been plotted for Ω = 0.1 s −1 .For the chosen values of parameters the radial and azimuthal components are significantly smaller than the radial component.This fact is in accordance with the initial assumption of non-rotational fluid (Ω→0).Noteworthy, that the radial velocity of the fluid v r is positive near the axis, indicating continuous expansion of the gas during its upwards movement.Another noteworthy conclusion is that the dimensional parameter λ decreases with decreasing temperature (Eq. 19).This means that for other constant conditions the radius of the tornado column decreases with hight.However, this effect is of the order of 1-2% for temperature variation of ±5% (i.e.tornado height below 2 km).

Convection with rotation
Taking into account compressibility of the gas (i.e.c ∞) and its rotation (Ω 0), one can reduce the linearized system (10)-( 14) to a single high-order ordinary differential equation for the poloidal field ϕ.Again we assume vertically homogeneous (∂/∂z=0), isotropic (∂/∂θ=0), and stationary (∂/∂t=0) system.Using the same dimensionless units (Eq.19) the equation reads as P.B. Rutkevich and P.P. Rutkevych: Tornado-type stationary vortex with nonlinear term due to moisture transport 5    where dimensionless parameter Q = (κ − 1)(2Ωgλ 3 ) 2 ν −2 c −4 characterizes the Coriolis force and compressibility.Clearly, in the case Q=0 one obtains Eq. ( 18).Due to 6-th order of Eq. ( 22), the expansion coefficients a 6 , a 8 etc. can be derived from given A, a 2 , a 4 , therefore the solutions now form a three-parameters set f (x;A,a 2 ,a 4 ), making the numerical search procedure significantly more complicated.To solve Eq. ( 22) we apply the procedure described in Sect.3, including Taylor series expansion in the vicinity of x=0 and Runge-Kutta numerical solution from x=0.1.Due to higher noise in the function, the values of the three independent parameters have been varied to achieve the smallest possible value of integral x 2 x 1 f (x) 2 dx, instead of the maximum x m described earlier.This choice of the convergency condition reflects the fact, that the absolute value of the function f should be the smallest in the medium range of radii.Figure 4 shows the best numerical solution, which we have found in such a way, using x 1 =8, x 2 =14, and given Q=0.1.The solution corresponds to the following set of parameters: A≈11.29547, a 2 ≈ − 5.454639, a 4 ≈ − 2.179057.One can see, that the found solution diverges at x m ≈13.
The radial and vertical velocity components can be obtained from the potential profile f (x) using Eq. ( 21), while for the azimuthal profile one has to solve differential equation Velocity profiles are plotted in Fig. 5.    of the fluid.In our simulation the azimuthal velocity changes the sign and infinitely increases, however we suppose this is an artifact of the numerical method.The expected convergency region of the azimuthal velocity is in the range x ≥ x 2 , where the profile f (x) used in Eq. ( 23) has large errors.The actual profile of azimuthal velocity, which can be regarded as fine-tuning of the solution, is a subject of further study.The coefficient λ used for dimensionless radii (Eq.19) equals to 70 m for the chosen values of parameters.Therefore one can conclude, that the main upward air stream in Fig. 3 has radius slightly more than 100 m.Although this value is relatively large for a tornado, the model does not include the rotation, which is likely to sharpen the structure.Indeed, solution with rotation in Fig. 5 has the radius of the upwards flow less than 70 m, which is very close to a typical tornado size (Davies-Jones, 2001).
The absolute value of velocity of 3 ms −1 in non-rotational case is reasonable, since effectively this is the velocity of a free localized convection in atmosphere.Account of rotation in the system leads to dramatic growth of velocity, especially its vertical component.This suggests that a model of tornado must include other non-linear dissipating terms to suppress the fluid flow.

Conclusions
The proposed one-dimensional vertically-homogeneous model cannot be applied to the tornado funnel near the mothercloud, or to the turbulent area near the ground surface.The model is applicable at the middle level of a weak tornado, where the shape of the column changes slowly in space and time.The model shows that the system of linearized hydrodynamic equations can have a stable solution with only one non-linear term.Since the full system of equations has several important non-linear terms (including viscosity and centrifugal force) this paper cannot be regarded as the final solution of the hydrodynamic system or as a complete theory of tornado.However, the results show, that a steady-state configuration is possible even in convectively-stable atmosphere.
Proposed non-linear effect related to asymmetric vapor transport in a "weak" tornado-like structure introduces the source of energy, which can be the true mechanism in natural intensive vortices.Obtained structures have all features of real tornadoes/cyclones, however estimations show that the size of the solution is closer to the size of tornado.We expect, that accounting other non-linear terms (which are usually dissipating) in the full system of equations will lead to stabilization of the solution and especially far from the axis.
Edited by: D. Giaiotti Reviewed by: two anonymous referees

Figure 1 .Figure 2 .
Figure 1.Profile of the x m versus the parameters a 2 and A. The peak position indicates the best approximation of the solution of the nonlinear equation (18) in non-rotational fluid.

Figure 1 .
Figure 1.Profile of the x m versus the parameters a 2 and A. The peak position indicates the best approximation of the solution of the nonlinear Eq. (18) in non-rotational fluid.

Figure 2 .Figure 4 .
Figure 2. The best approximation of dimensionless poloidal velocity profile f (x).The values of the parameters A and a 2 in the Taylor expansion (20) correspond to the position of the peak in Fig. 1.

Figure 2 .
Figure 2. The best approximation of dimensionless poloidal velocity profile f (x).The values of the parameters A and a 2 in the Taylor expansion (Eq.20) correspond to the position of the peak in Fig. 1.

Figure 1 .
Figure 1.Profile of the x m versus the parameters a 2 and A. The peak position indicates the best approximation of the solution of the nonlinear equation (18) in non-rotational fluid.

Figure 2 .Figure 3 .
Figure 2. The best approximation of dimensionless poloidal velocity profile f (x).The values of the parameters A and a 2 in the Taylor expansion (20) correspond to the position of the peak in Fig. 1.

Figure 4 .
Figure 4. Profile of dimentionless potential field f (x) with rotation.The values of parameters are the same as the non-rotational case, except the Coriolis parameter Ω=0.23s −1 .

Figure 3 .
Figure 3. Profiles of velocity components in a non-rotational approximation.Values of parameters are listed at the end of Sect. 3. Notice a factor in front of the radial velocity component v r .

Figure 4 .
Figure 4. Profile of dimentionless potential field f (x) with rotation.The values of parameters are the same as the non-rotational case, except the Coriolis parameter Ω=0.23s −1 .

Figure 4 .Figure 5 .
Figure 4. Profile of dimensionless potential field f (x) with rotation.The values of parameters are the same as the non-rotational case, except the Coriolis parameter Ω=0.23 s −1 .

Figure 5 .
Figure 5. Profiles of velocity components with rotation, parameter values are the same as in Fig. 4.
Rutkevich and P. P. Rutkevych: Tornado-type stationary vortex with nonlinear term due to moisture transport 79 the temperature profile on the vertical velocity of the fluid.In this work we use a simple approximation Adv.Sci.Res., 4, 77-82, 2010www.adv-sci-res.net/4/77/2010/P. B.
The air motion has several important features, namely updrift in the central area, relatively small radial velocity, and azimuthal velocity has maximum at approximately x≈8.The azimuthal velocity in axially-symmetrical system is the parameter with the slowest radial convergency, because it is limited only by the viscosity Adv.Sci.Res., 4, 77-82, 2010 www.adv-sci-res.net/4/77/2010/P. B. Rutkevich and P. P. Rutkevych: Tornado-type stationary vortex with nonlinear term due to moisture transport 81