Influence of Horizontal Temperature Gradients on Convective Instabilities with Geophysical Interest

Since Bénard’s experiments on convection and Rayleigh’s theoretical work in the beginning of the twentieth century (1)-(2), many experimental, theoretical and numerical works related to Rayleigh-Bénard convection have been done (3)-(10) and different problems have been posed depending on what is to be modelled. Classically, heat is applied uniformly from below and the conductive solution becomes unstable for a critical vertical gradient beyond a certain threshold. A setup for natural convection more general than that of uniform heating consists of including a non-zero horizontal temperature gradient which may be either constant or not (11)-(29). In those problems a clear difference is marked by the fact that the fluid is simply contained (11)-(19), where stationary and oscillatory instabilities appear depending on the multiple parameters present in the problem: properties of the fluid, surface tension effects, heat exchange with the atmosphere, aspect ratio, dependence of viscosity with temperature, etc., and the case where the fluid can flow throughout the boundaries (29), where vortical solutions can appear reinforcing the relevance of convective mechanisms for the generation of vertical vortices very similar to those found for some atmospheric phenomena as dust devils or hurricanes (29)-(31). The case where the fluid is simply contained displays stationary and oscillatory instabilities. This problem has been treated from different points of view: experimental (11)-(18) and theoretical, both with semiexact (20)-(21) and numerical solutions (40)-(28). This case contains applications to mantle convection when the viscosity is large (45; 52) or it depends on temperature (19). There are not experiments yet for the case where a flow throughout the boundaries is allowed, only observations of atmospheric phenomena (30; 33; 34; 36; 37), and theoretic numerical results (29; 31). In this work we will review this physical problem, focusing on the latest problems addressed by the authors on this topic, where a non-uniform heating is considered in different geometrical configurations, and we will show the relevant results obtained, some of them in the context of interesting atmospheric and geophysical phenomena (30; 36; 37). 5


Theoretical formulation of the problem
The physical set-up (see figure 1) consists of a horizontal fluid layer in a rectangular domain (19; 45) or a cylindrical annulus (18; 25; 28) between two vertical walls at r = a and r = a + l. The depth of the domain is d (z coordinate). At z = 0 the imposed temperature gradient takes the value T max at a and the value T min at the outer part (a + l). The upper surface is at temperature T = T 0 . We define △T v = T max − T 0 , △T h = T max − T min and δ = △T h /△T v . From now on we will consider an annular domain, therefore cylindrical coordinates will be used in the following. The formulation in a rectangular domain and coordinates would be similar. In the governing equations, u =( u r , u φ , u z ) is the velocity field, T is the temperature, p is the pressure, r is the radial coordinate, and t is the time. They are expressed in dimensionless form after rescaling: r ′ = r/d, t ′ = κt/d 2 , u ′ = du/κ, p ′ = d 2 p/ (ρ 0 κν) , Θ = (T − T 0 ) /△T v . Here r is the position vector, κ the thermal diffusivity, ν the kinematic viscosity of the liquid, and ρ 0 the mean density at temperature T 0 . The domain is [ā,ā + Γ] × [0, 1] × [0, 2π] where Γ = l/d is the aspect ratio andā = a/d. The system evolves according to the mass balance, energy conservation and momentum equations, which in dimensionless form (with primes now omitted) are, ∂ t u + (u ·∇) u = Pr −∇p + ∇ 2 u + RΘe z , where the operators and fields are expressed in cylindrical coordinates and the Oberbeck-Boussinesq approximation has been used (25), i.e. density is constant except in the term of gravity, where a linear dependence on temperature is considered. Here e z is the unit vector in the z direction. The following dimensionless numbers have been introduced: the Prandtl number Pr = ν/κ, and the Rayleigh number R = gα△Td 3 /κν, which represents the effect of buoyancy and in which α is the thermal expansion coefficient and g the gravitational acceleration. In the case of variable viscosity the laplacian operator in Eq. (3) takes the form div ν(Θ)

Contained fluid
Regarding boundary conditions, several conditions can be considered such as that one where flow through the boundaries is not permitted. For instance at the lateral walls r =ā and r =ā + Γ the velocity is zero and an insulating wall is considered, On the top surface, the vertical velocity is zero, the normal derivatives of the rest of components of the velocity are zero and the temperature is T = T 0 , that after rescaling become, and at the bottom For temperature at the bottom we consider a constant horizontal temperature difference, i.e. a linear profile. Here the horizontal temperature gradient appears, with θ 1 (r)=1 − rδ/Γ and a second order polynomial which matches the linear profile such that ∂ r θ 1 (r)=0onr =ā and r =ā + Γ. The dimensionless equations and boundary conditions contain five external parameters: R, Γ, Pr, δ, and η.

Not contained fluid
Regarding boundary conditions, several conditions can be considered like allowing flow through the boundaries. At the lateral inner wall r =ā the velocity is zero and an insulating wall is considered, At r =ā + Γ, the lateral outer wall, a constant radial velocity is assumed and an insulating wall is considered, On the top surface, the velocity is zero and the temperature is T = T 0 , that after rescaling become, and at the bottom u r = ∂ z u φ = u z = 0, on z = 0. (11) For temperature at the bottom we consider a variable horizontal temperature gradient through imposing a Gaussian profile as in Ref. (28), The dimensionless equations and boundary conditions contain five external parameters: R, Γ, Pr, δ, and β.

Metodology: search for solutions and their linear stability
We look for stationary axisymmetric solutions of the problem, then, the equations to be solved are where ∇ * =(∂ r ,0,∂ z ), together with the corresponding boundary conditions. The time independent solution U b (r, z) to the stationary problem obtained from equations (1)-(3) by eliminating the time dependence, is called basic state. It is a non-conductive state (u = 0) as soon as δ = 0. The basic state is considered to be axisymmetric and therefore depends only on r − z coordinates (i.e. all φ derivatives are zero). The velocity field of the basic flow is restricted to u =(u r , u φ = 0, u z ).
A linear stability analysis of the stationary solutions is performed. Fixed (Γ, δ, Pr, β), the solution U(r, z, t)=(u, Θ, p)(r, z, t) of (1)-(3) at given R is expressed as where U b (r, z) is the base flow for the given (R, Γ, δ, Pr, β) andŨ(r, z) refers to the perturbation. We have considered Fourier mode expansions in the angular direction, because along it boundary conditions are periodic. Introducing (16) into the full system (1)-(3) and linearizing the resulting system, the following eigenvalue problem in λ is obtained: where ∇ k =(∂ r , ik, ∂ z ), with the corresponding boundary conditions. The instability is achieved when the real part of the eigenvalue with maximum real part, λ max (R), changes from a negative value to a positive one as R increases, for a specific wave number k. The critical value of R for which λ max (R, k)=0 is denoted by R c and the critical wave number, minimum k for which the bifurcation occurs, by k c .

Numerical methods
The numerical method is described in detail and tested in Refs. (25; 28). The nonlinearities in the basic state problem are solved with a Newton method. Each step of the Newton method and the linear stability analysis have been numerically solved with a Chebyshev collocation method as explained in Refs. (28; 39; 48). The problem is posed in the primitive variables formulation, and the use of the same order approximations for velocity and pressure in the Chebyshev collocation procedure introduces spurious modes for pressure that are solved by adding convenient boundary conditions (43; 44). In the resulting linear problems any unknown field x is expanded in Chebyshev polynomials The corresponding expansions for the four different fields are introduced into the Newton linearized version of equations (13)- (15) and the corresponding boundary conditions and evaluated at the Chebysehv Gauss-Lobatto collocation points (r j , z i ), Some care is necessary in the evaluation rules at the boundaries as explained in Refs. (28; 48). At each iteration of the Newton method a linear system of the form AX = B is derived, where X is a vector containing P = 4 × L × N unknowns and A is a full rank matrix of order P × P. This can be solved with standard routines. The algorithm starts with an approximation to the solution x 0,LN and the iteration procedure is applied until the stop criterion ||x s+1,LN − x s,LN || < 10 −9 is satisfied. The same discretization is used for the eigenvalue problem (17)- (19) with the corresponding boundary conditions. In this way it is transformed into its discrete form by expanding the perturbations in a truncated series of Chebyshev polynomials (20) as performed for the basic state. The evaluation rules are detailed in Ref. (48). Therefore, the eigenvalue problem in its discrete form is, where w is a vector which contains Q unknowns and C and B are Q × Q matrices, with Q = 5 × L × N. QZ or Arnoldi algorithms are used to solve the eigenvalue problem (42). σ are the eigenvalues and w are coefficients in the Chebyshev basis of the corresponding eigenfunctions.
The discrete eigenvalue problem (23) has a finite number of eigenvalues σ i . The stability condition must now be imposed upon σ max where σ max = max(Re(σ i )), bearing in mind In order to test convergence of the method we include, as an example, the calculation of the critical value of the bifurcation parameter, R c , and the critical wave number, k c , for different order expansions in the Chebyshev approximation in the contained fluid case. And we benchmark the method and code to ensure the correctness of the results. Table 1 shows these results for the contained fluid case. When the orders L and N are increased, the critical values tend to a determined value, convergence is very good and for L = 24 and N = 14 the results are sufficiently accurate, in fact they are exact to the thousandth. The values L = 24 and N = 14 can be considered in the computations. In a convergence test comparing the critical R c obtained at different order expansions, the relative difference between the expansions at 26 × 18 and 24 × 16 is found to be less than 10 −4 . There are three significant digits in this calculation. The benchmarking of the method can be done with results in Ref. (48). The critical wave number for Γ = 2.936, η = 0.0862 and δ = 0isk c = 0, so for these values of the parameters the results reported in Ref. (48) are recovered. For this value of the aspect ratio the bifurcation corresponds to a mode 2 in the x direction and the bifurcation takes place at the same value R c = 73.5.  (20) for a 3D fluid with constant viscosity, η = 0, aspect ratio Γ = 2.936 and δ = 0.1.

Contained fluid
In references of small cells the case of large viscosity (or Prandlt number) could be considered as an approximation to mantle convection. The largest value of Prandlt number considered in the experiments is Pr = 60 in Ref. (52), in this case boundary layer waves are observed.
Numerical results with infinite Pr number are reported in (45). In this case only stationary patterns of rolls perpendicular to the temperature gradient are observed. Also it is of interest the case of variable viscosity dependent on temperature, this case is plenty of references, but all of them consider homogeneous heating without horizontal temperature gradients (35; 36; 47). The only reference in which those gradients are taken into account in a variable viscosity case is (19). Some numerical solutions obtained in the case considered in Ref. (19) are presented in figure 2 Figure 2 shows that the structure of the velocity field is more localized close to the zone where the temperature is higher, i.e, at r = −1. The presence of the horizontal gradient generates convective basic states, that were conductive without the horizontal gradients. In Ref. (19) it is shown the fluid motion is produced in the region where viscosity is lower.
Regarding the instabilities, in the case of large Γ the influence of the horizontal temperature gradient is considerable, the problem is nearly two-dimensional (2D) in the uniform heating case, but it is three-dimensional (3D) with the horizontal temperature gradient. Figure 3 shows the growing mode or eigenfunction in the case Γ = 2.936, δ = 0 and η = 0.0862, the critical wave number in this case is k c = 0, so a 2D structure appears after the bifurcation. Figure  4 shows the growing mode or eigenfunction in the same case as before, but with horizontal gradient δ = 0.1, the critical wave number in this case is k c = 1.7, so a 3D structure appears after the bifurcation. Also we can observe from figures 3 and 4 that the bifurcating pattern is more structured in the r − z plane for δ = 0 and becomes more structured in the y − z plane for δ = 0. Hence, a horizontal temperature gradient gives rise to thermal plumes which bifurcate to totally 3D structures.

Not contained fluid
This case is plenty of references of direct simulations solving numerically the partial differential equations (31; 34). But under the instability or bifurcation perspective the case in which the fluid can flow through the boundaries is only treated in reference (29). In that paper we show that a vortical structure appears after a stationary bifurcation of a state without angular velocity.  A numerical solution obtained in the problem considered in Ref. (29) is presented in figure 5. Figure 5 shows the profiles of temperature, pressure and velocity components corresponding to the clockwise vortex for Pr = 0.7, Γ = 0.5 and δ = 10 at R = 4367. This vortex appears after a bifurcation of a basic state with zero azimuthal velocity (see Ref. (29)). The main feature of the new steady flow emerging from the convective instability with k c = 0 andũ φ = 0i sa non-zero azimuthal velocity component. The fluid inside the annulus begins to move in the azimuthal direction, rotating around the inner cylinder.
The linear stability analysis of the vortical structures shows that there is a wide range of parameters for which this state is stable. The track of a particle in the vortex can be obtained by integrating the evolution of the element of fluid which follows the velocity field, In our simulations we observe a spiral upward motion of the particle around the inner cylinder, which implies a transport of mass in the azimuthal direction. Starting from below, the particle goes up, moves towards the inner cylinder and rotates around it. The combination of these movements gives the spiral trajectory shown in figure 6, where the trajectory of a particle in the fluid is presented for Γ = 0.5, δ = 10 and R = 4367 at two different initial conditions. Starting from a point close to the bottom plate but near the inner cylinder, where the effect of u z is stronger than the effect of u φ , the particle goes up very fast without much turning around the inner cylinder. This can be appreciated in figure 6 a) where the starting point considered is (r = 0.085, z = 0.05, φ = 0) in [0.06, 0.56] × [0, 1] × [0, 2π]. When the particle reaches the upper part of the structure it describes wider circles around the inner cylinder as u r becomes positive and u z is very small at those levels (see figure 5). Figure  6 b) shows the effect of localizing the starting point further from the inner cylinder, e.g. at (r = 0.31, z = 0.05, φ = 0). In this case, the effect of u φ is stronger and the spiral up motion of the particle starts as soon as the particle begins to move.

Contained fluid
Regarding the numerical solutions found in this case, the horizontal temperature gradient generates convective states and tends to concentrate motion near the warmer wall. This fact coincides with experiments in Refs. (38; 40; 41) and is consistent with previous numerical results reported in (27; 45). The temperature dependent viscosity localizes motion near the region of lower viscosity, i.e., the bottom plate. This also coincides with experiments in Refs. (36; 46; 49) and numerical results in Ref. (48). It is remarkable that the horizontal temperature gradient favours a threedimensional structure after the bifurcation, while the pattern continues being axisymmetric after the bifurcation in the only vertical gradient case.

Not contained fluid
As detailed in Ref. (48) we found qualitative similarities between the vortical structures computed numerically and some meteorological phenomena such as dust devils and cyclones. One of the main characteristics of dust devils is a low-pressure region in the center of the dust devil which coincides with the dust devil's warm core (33). This is also observed in our numerical vortices (see figure 5 c). Regarding temperature, in dust devils, the maximum temperature deviation from the environment temperature (i.e. the temperature furthest from the dust devil center) occurs at the lowest levels. This feature is observed in the temperature profile of our vortices. The experimental measures provided in Ref. (33) show that there is radial inflow at the lower levels of the dust devil and radial outflow in the upper levels. It is also shown that the vertical velocity reaches highest values and then falls off rapidly as the radius is increased. These features are appreciated in the profile of u r and u z shown in figures 5 d) and e). The trajectory of particles around the inner cylinder described in this section appears to be very similar to the trajectory of particles of air (or dust) in a dust devil, characterized by a spiral up motion (33). Other more complex meteorological phenomena such as cyclones also present these structural characteristics. It is known that the center (eye) of a cyclone is the area of lowest atmospheric pressure in the region, which corresponds to a warm core in some kind of cyclones (e.g. tropical and mesoscale) (31; 34). This coincides with that observed in figures 5 a) and 5 b). Regarding the motion in cyclones, it is observed the inward flow next to the surface, strong upward motion in the eyewall and outflow in a layer near the top of the storm (31; 34). This characteristic is described in the combined effect of the radial and vertical velocity components observed in our vortices as pointed out above (see figures figures 5 c, d) and e)). In cyclones, a counter-clockwise motion (clockwise in the southern hemisphere) is observed around the center of the storm, stronger just above the surface in a ring around the center and sligther as we go up from the surface (31; 34). That coincides with the effect of the azimuthal velocity component observed in the vortices we have computed numerically responsible for the movement of the particles around the inner cylinder.

Conclusions
In this work we have reviewed the influence of horizontal temperature gradients on convective instabilities, focusing on results with geophysical interest.
We have distinguished two cases, a first one where the fluid is simply contained in a domain, and a second one where the fluid can flow throughout the boundaries. In the first case three subcases can be grouped. The case corresponding to small cells and small Pr number displays stationary and oscillatory instabilities depending on the multiple parameters present in the problem: properties of the fluid, surface tension effects, heat exchange with the atmosphere, aspect ratio, dependence of viscosity with temperature, etc. This problem has been treated from experimental, theoretical and numerical points of view.
The cases corresponding to small cells and large or infinite Pr number are closer to mantle convection. Boundary layer waves are observed in experiments and 3D stationary patterns of rolls perpendicular to the temperature gradient appear numerically. Finally for the case of infinite Pr number with temperature dependent viscosity, the closest to mantle convection, 3D stationary patterns concentrated in the region of lower viscosity and waves for larger values of the R number appear. Summarizing, horizontal temperature gradients favour the presence of waves and the totally three dimensional patterns. The problem where the fluid can flow throughout the boundaries has been treated usually as direct numerical simulations. For the first time it has been studied under the perspective of instabilities or bifurcations in Ref. (29). In this reference vortical solutions, very similar to those found for some atmospheric phenomena such as dust devils or hurricanes, appear after a stationary bifurcation. This is a powerfull and simple explanation of those atmospheric phenomena as an instability.

Acknowledgements
This work was partially supported by Research Grant MICINN (the Government of Spain) MTM2009-13084 and CCYT (Junta de Comunidades de Castilla-La Mancha) PAI08-0269-1261, which include RDEF funds.