**1. Introduction**

80 Advanced Fluid Dynamics

Wu, P.-K., Kirkendall, K. A., Fuller, R. P., & Nejad, A. S. (1997), Breakup Processes of

pp. 64 -73.

Liquid Jets in Subsonic Cross-flows. *Journal of Propulsion and Power*, Vol. 13, No. 1,

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).

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***<sup>z</sup>* is the unit vector in the *z* direction. The following dimensionless numbers have been introduced: the Prandtl number *Pr* <sup>=</sup> *<sup>ν</sup>*/*κ*, and the Rayleigh number *<sup>R</sup>* <sup>=</sup> *<sup>g</sup>α*�*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

<sup>83</sup> Influence of Horizontal Temperature Gradients

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* = *a*¯ and

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,

For temperature at the bottom we consider a constant horizontal temperature difference, i.e.

with *<sup>θ</sup>*1(*r*) = <sup>1</sup> <sup>−</sup> *<sup>r</sup>δ*/<sup>Γ</sup> and a second order polynomial which matches the linear profile such

The dimensionless equations and boundary conditions contain five external parameters:

Regarding boundary conditions, several conditions can be considered like allowing flow through the boundaries. At the lateral inner wall *r* = *a*¯ the velocity is zero and an insulating

At *r* = *a*¯ + Γ, the lateral outer wall, a constant radial velocity is assumed and an insulating

On the top surface, the velocity is zero and the temperature is *T* = *T*0, that after rescaling

*ur* = *u<sup>φ</sup>* = *∂ruz* = *∂r*Θ = 0, on *r* = *a*¯ and *r* = *a*¯ + Γ. (4)

*∂zur* = *∂zu<sup>φ</sup>* = *uz* = Θ = 0, on *z* = 1, (5)

*ur* = *u<sup>φ</sup>* = *uz* = 0, on *z* = 0. (6)

Θ = *θ*1(*r*) on *z* = 0. (7)

*ur* = *u<sup>φ</sup>* = *uz* = *∂r*Θ = 0, on *r* = *a*¯. (8)

*ur* = *u<sup>φ</sup>* = *uz* = Θ = 0, on *z* = 1, (10)

*∂rur* = *∂ru<sup>φ</sup>* = *∂ruz* = *∂r*Θ = 0, on *r* = *a*¯ + Γ. (9)

, where *ν*(Θ) = *ν*0*e*−*ηR*Θ.

*r* = *a*¯ + Γ the velocity is zero and an insulating wall is considered,

a linear profile. Here the horizontal temperature gradient appears,

div *<sup>ν</sup>*(Θ) *ν*<sup>0</sup> · 

**2.1 Contained fluid**

and at the bottom

*R*, Γ, *Pr*, *δ*, and *η*.

wall is considered,

wall is considered,

become,

**2.2 Not contained fluid**

<sup>∇</sup>**<sup>u</sup>** + (∇**u**)*<sup>t</sup>*

that *∂rθ*1(*r*) = 0 on *r* = *a*¯ and *r* = *a*¯ + Γ.

on Convective Instabilities with Geophysical Interest

#### **2. 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 �*Tv* = *T*max − *T*0, �*Th* = *T*max − *T*min and *δ* = �*Th*/�*Tv*.

Fig. 1. Physical setup for the cylindrical annulus.

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** = (*ur*, *uφ*, *uz*) 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**� = *d***u**/*κ*, *p*� = *<sup>d</sup>*<sup>2</sup> *<sup>p</sup>*/ (*ρ*0*κν*), <sup>Θ</sup> <sup>=</sup> (*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0) /�*Tv*. Here **<sup>r</sup>** is the position vector, *<sup>κ</sup>* the thermal diffusivity, *<sup>ν</sup>* the kinematic viscosity of the liquid, and *ρ*<sup>0</sup> the mean density at temperature *T*0. The domain is [*a*¯, *a*¯ + Γ] × [0, 1] × [0, 2*π*] where Γ = *l*/*d* is the aspect ratio and *a*¯ = *a*/*d*.

The system evolves according to the mass balance, energy conservation and momentum equations, which in dimensionless form (with primes now omitted) are,

$$\nabla \cdot \mathbf{u} = \mathbf{0},\tag{1}$$

$$
\partial\_l \Theta + \mathbf{u} \cdot \nabla \Theta = \nabla^2 \Theta,\tag{2}
$$

$$\left(\partial\_t \mathbf{u} + \left(\mathbf{u} \cdot \nabla\right) \mathbf{u} = Pr \left(-\nabla p + \nabla^2 \mathbf{u} + R \Theta \mathbf{e}\_z\right)\right) \tag{3}$$

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***<sup>z</sup>* is the unit vector in the *z* direction. The following dimensionless numbers have been introduced: the Prandtl number *Pr* <sup>=</sup> *<sup>ν</sup>*/*κ*, and the Rayleigh number *<sup>R</sup>* <sup>=</sup> *<sup>g</sup>α*�*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 *<sup>ν</sup>*(Θ) *ν*<sup>0</sup> · <sup>∇</sup>**<sup>u</sup>** + (∇**u**)*<sup>t</sup>* , where *ν*(Θ) = *ν*0*e*−*ηR*Θ.

#### **2.1 Contained fluid**

2 Will-be-set-by-IN-TECH

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 �*Tv* = *T*max − *T*0, �*Th* = *T*max − *T*min and *δ* = �*Th*/�*Tv*.

ez

g

a


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** = (*ur*, *uφ*, *uz*) is the velocity field, *T* is the temperature, *p* is the pressure, *r* is the radial coordinate, and *t* is the time. They are

*<sup>d</sup>*<sup>2</sup> *<sup>p</sup>*/ (*ρ*0*κν*), <sup>Θ</sup> <sup>=</sup> (*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0) /�*Tv*. Here **<sup>r</sup>** is the position vector, *<sup>κ</sup>* the thermal diffusivity, *<sup>ν</sup>* the kinematic viscosity of the liquid, and *ρ*<sup>0</sup> the mean density at temperature *T*0. The domain is

The system evolves according to the mass balance, energy conservation and momentum

∇ · **u** = 0, (1) *<sup>∂</sup>t*<sup>Θ</sup> <sup>+</sup> **<sup>u</sup>** · ∇<sup>Θ</sup> <sup>=</sup> <sup>∇</sup>2Θ, (2)

−∇*<sup>p</sup>* <sup>+</sup> <sup>∇</sup>2**<sup>u</sup>** <sup>+</sup> *<sup>R</sup>*Θ**e***<sup>z</sup>*

er

l

e

d

� = *κt*/*d*2, **u**� = *d***u**/*κ*, *p*� =

, (3)

**2. Theoretical formulation of the problem**

T0

Tmax

Tmin

expressed in dimensionless form after rescaling: **r**� = **r**/*d*, *t*

[*a*¯, *a*¯ + Γ] × [0, 1] × [0, 2*π*] where Γ = *l*/*d* is the aspect ratio and *a*¯ = *a*/*d*.

equations, which in dimensionless form (with primes now omitted) are,

*∂t***u** + (**u** · ∇) **u** = *Pr*

Tv

Th

Fig. 1. Physical setup for the cylindrical annulus.

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* = *a*¯ and *r* = *a*¯ + Γ the velocity is zero and an insulating wall is considered,

$$
\mu\_r = \mu\_\phi = \partial\_r \mu\_z = \partial\_r \Theta = 0,\text{ on }\ r = \vec{a}\text{ and }\ r = \vec{a} + \Gamma.\tag{4}
$$

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,

$$
\partial\_2 \mathfrak{u}\_{\mathbb{Z}} = \partial\_2 \mathfrak{u}\_{\Phi} = \mathfrak{u}\_{\mathbb{Z}} = \Theta = 0,\text{ on }\ z = 1,\tag{5}
$$

and at the bottom

$$
\mathfrak{u}\_{\mathfrak{r}} = \mathfrak{u}\_{\mathfrak{\phi}} = \mathfrak{u}\_{\mathfrak{z}} = 0,\text{ on }\mathfrak{z} = 0.\tag{6}
$$

For temperature at the bottom we consider a constant horizontal temperature difference, i.e. a linear profile. Here the horizontal temperature gradient appears,

$$
\Theta = \theta^1(r) \quad \text{on} \ z = 0. \tag{7}
$$

with *<sup>θ</sup>*1(*r*) = <sup>1</sup> <sup>−</sup> *<sup>r</sup>δ*/<sup>Γ</sup> and a second order polynomial which matches the linear profile such that *∂rθ*1(*r*) = 0 on *r* = *a*¯ and *r* = *a*¯ + Γ.

The dimensionless equations and boundary conditions contain five external parameters: *R*, Γ, *Pr*, *δ*, and *η*.

#### **2.2 Not contained fluid**

Regarding boundary conditions, several conditions can be considered like allowing flow through the boundaries. At the lateral inner wall *r* = *a*¯ the velocity is zero and an insulating wall is considered,

$$
\mu\_r = \mu\_\phi = \mu\_\mathbb{Z} = \partial\_\mathbb{r} \Theta = 0,\text{ on }\ r = \mathbb{\bar{a}}.\tag{8}
$$

At *r* = *a*¯ + Γ, the lateral outer wall, a constant radial velocity is assumed and an insulating wall is considered,

$$
\partial\_r \mu\_r = \partial\_r \mu\_\phi = \partial\_r \mu\_z = \partial\_r \Theta = 0,\text{ on }r = \overline{a} + \Gamma.\tag{9}
$$

On the top surface, the velocity is zero and the temperature is *T* = *T*0, that after rescaling become,

$$
\mu\_r = \mu\_\phi = \mu\_z = \Theta = 0,\text{ on } z = 1,\tag{10}
$$

number *k*. The critical value of *R* for which *λ*max(*R*, *k*) = 0 is denoted by *Rc* and the critical

<sup>85</sup> Influence of Horizontal Temperature Gradients

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

wave number, minimum *k* for which the bifurcation occurs, by *kc*.

on Convective Instabilities with Geophysical Interest

unknown field **x** is expanded in Chebyshev polynomials

**x***LN* =

evaluated at the Chebysehv Gauss-Lobatto collocation points (*rj*, *zi*),

*rj* <sup>=</sup> cos *<sup>j</sup>* <sup>−</sup> <sup>1</sup>

*zi* <sup>=</sup> cos *<sup>i</sup>* <sup>−</sup> <sup>1</sup>

*L*−1 ∑ *l*=0

*<sup>L</sup>* <sup>−</sup> <sup>1</sup> <sup>−</sup> <sup>1</sup>

*<sup>N</sup>* <sup>−</sup> <sup>1</sup> <sup>−</sup> <sup>1</sup>

*N*−1 ∑ *n*=0 *a***x**

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

> *π*

> > *π*

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 **<sup>x</sup>**0,*LN* and the iteration procedure is applied until the stop criterion ||**x***s*+1,*LN* <sup>−</sup>

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

where *w* is a vector which contains *Q* unknowns and *C* and *B* are *Q* × *Q* matrices, with *Q* =

QZ or Arnoldi algorithms are used to solve the eigenvalue problem (42). *σ* are the eigenvalues

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

and *w* are coefficients in the Chebyshev basis of the corresponding eigenfunctions.

*lnTl*(*r*)*Tn*(*z*). (20)

, *j* = 1, ..., *L*. (21)

, *i* = 1, ..., *N*. (22)

*Cw* = *σBw*, (23)

**3.1 Numerical methods**

**<sup>x</sup>***s*,*LN*|| <sup>&</sup>lt; <sup>10</sup>−<sup>9</sup> is satisfied.

discrete form is,

5 × *L* × *N*.

and at the bottom

$$
\mu\_r = \partial\_{\overline{z}} \mu\_\phi = \mu\_z = 0,\text{ on } z = 0. \tag{11}
$$

For temperature at the bottom we consider a variable horizontal temperature gradient through imposing a Gaussian profile as in Ref. (28),

$$\Theta = 1 - \delta(e^{\left(\frac{1}{\beta}\right)^2} - e^{\left(\frac{1}{\beta} - (\frac{\epsilon - 4}{\Gamma})^2 \frac{1}{\beta}\right)^2})/(e^{\left(\frac{1}{\beta}\right)^2} - 1) \text{ on } z = 0. \tag{12}$$

The dimensionless equations and boundary conditions contain five external parameters: *R*, Γ, *Pr*, *δ*, and *β*.

#### **3. Metodology: search for solutions and their linear stability**

We look for stationary axisymmetric solutions of the problem, then, the equations to be solved are

$$\nabla^\* \cdot \mathbf{u} = \mathbf{0},\tag{13}$$

$$\mathbf{u} \cdot \nabla^\* \Theta = \nabla^{\*2} \Theta,\tag{14}$$

$$\left(\mathbf{u}\cdot\nabla^{\ast}\right)\mathbf{u} = \operatorname{Pr}\left(-\nabla^{\ast}p + \nabla^{\ast2}\mathbf{u} + R\Theta\mathbf{e}\_{z}\right),\tag{15}$$

where ∇<sup>∗</sup> = (*∂r*, 0, *∂z*), together with the corresponding boundary conditions.

The time independent solution *Ub*(*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** = (*ur*, *u<sup>φ</sup>* = 0, *uz*).

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

$$\mathcal{U}(r,z,t) = \mathcal{U}^b(r,z) + \tilde{\mathcal{U}}(r,z)e^{i\mathbf{k}\phi + \lambda t},\tag{16}$$

where *Ub*(*r*, *z*) is the base flow for the given (*R*, Γ, *δ*, *Pr*, *β*) and *U*˜ (*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:

$$
\nabla^k \cdot \vec{\mathcal{U}} = 0,\tag{17}
$$

$$
\lambda \not\Theta + \tilde{\mathcal{U}} \cdot \nabla^k \Theta^b + \mathcal{U}^b \cdot \nabla^k \tilde{\Theta} = (\nabla^k)^2 \tilde{\Theta},\tag{18}
$$

$$
\lambda \mathcal{U} + \left( \mathcal{U} \cdot \nabla^k \right) \mathcal{U}^b + \left( \mathcal{U}^b \cdot \nabla^k \right) \mathcal{U} = \text{Pr} \left( -\nabla^k \vec{p} + (\nabla^k)^2 \mathcal{U} + R \vec{\Theta} \mathbf{e}\_z \right), \tag{19}
$$

where <sup>∇</sup>*<sup>k</sup>* = (*∂r*, *ik*, *<sup>∂</sup>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 *Rc* and the critical wave number, minimum *k* for which the bifurcation occurs, by *kc*.

#### **3.1 Numerical methods**

4 Will-be-set-by-IN-TECH

For temperature at the bottom we consider a variable horizontal temperature gradient through

The dimensionless equations and boundary conditions contain five external parameters:

We look for stationary axisymmetric solutions of the problem, then, the equations to be solved

The time independent solution *Ub*(*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

A linear stability analysis of the stationary solutions is performed. Fixed (Γ, *δ*, *Pr*, *β*), the

where *Ub*(*r*, *z*) is the base flow for the given (*R*, Γ, *δ*, *Pr*, *β*) and *U*˜ (*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

<sup>∇</sup>*<sup>k</sup>* · *<sup>U</sup>*˜ <sup>=</sup> 0, (17)

*<sup>λ</sup>*Θ˜ <sup>+</sup> *<sup>U</sup>*˜ · ∇*k*Θ*<sup>b</sup>* <sup>+</sup> *<sup>U</sup><sup>b</sup>* · ∇*k*Θ˜ = (∇*k*)2Θ˜ , (18)

*U*˜ = *Pr*

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

*U*(*r*, *z*, *t*) = *Ub*(*r*, *z*) + *U*˜ (*r*, *z*)*e*

linearizing the resulting system, the following eigenvalue problem in *λ* is obtained:

where ∇<sup>∗</sup> = (*∂r*, 0, *∂z*), together with the corresponding boundary conditions.

solution *U*(*r*, *z*, *t*)=(**u**, Θ, *p*)(*r*, *z*, *t*) of (1)-(3) at given *R* is expressed as

where <sup>∇</sup>*<sup>k</sup>* = (*∂r*, *ik*, *<sup>∂</sup>z*), with the corresponding boundary conditions.

*ur* = *∂zu<sup>φ</sup>* = *uz* = 0, on *z* = 0. (11)

∇<sup>∗</sup> · **u** = 0, (13)

**<sup>u</sup>** · ∇∗<sup>Θ</sup> <sup>=</sup> <sup>∇</sup>∗2Θ, (14)

*ikφ*+*λt*

−∇*<sup>k</sup> <sup>p</sup>*˜ + (∇*k*)2*U*˜ <sup>+</sup> *<sup>R</sup>*Θ˜ **<sup>e</sup>***<sup>z</sup>*

−∇<sup>∗</sup> *<sup>p</sup>* <sup>+</sup> <sup>∇</sup>∗2**<sup>u</sup>** <sup>+</sup> *<sup>R</sup>*Θ**e***<sup>z</sup>*

− 1) on *z* = 0. (12)

, (15)

, (16)

, (19)

and at the bottom

*R*, Γ, *Pr*, *δ*, and *β*.

are

imposing a Gaussian profile as in Ref. (28),

Θ = 1 − *δ*(*e*

( 1 *β* ) 2 − *e* ( 1 *<sup>β</sup>* <sup>−</sup>( *<sup>r</sup>*−*a*¯ <sup>Γ</sup> )<sup>2</sup> <sup>1</sup> *β* )2 )/(*e* ( 1 *β* )2

**3. Metodology: search for solutions and their linear stability**

(**u** · ∇∗) **u** = *Pr*

basic flow is restricted to **u** = (*ur*, *u<sup>φ</sup>* = 0, *uz*).

*λU*˜ + *<sup>U</sup>*˜ · ∇*<sup>k</sup> U<sup>b</sup>* + *<sup>U</sup><sup>b</sup>* · ∇*<sup>k</sup>*  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

$$\mathbf{x}^{LN} = \sum\_{l=0}^{L-1} \sum\_{n=0}^{N-1} a\_{ln}^{\mathbf{x}} T\_l(r) T\_n(z). \tag{20}$$

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 (*rj*, *zi*),

$$r\_{\dot{j}} = \cos\left(\left(\frac{j-1}{L-1} - 1\right)\pi\right), \quad \dot{j} = 1, \ldots, L. \tag{21}$$

$$z\_i = \cos\left(\left(\frac{i-1}{N-1} - 1\right)\pi\right), \quad i = 1, \ldots, N. \tag{22}$$

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 **<sup>x</sup>**0,*LN* and the iteration procedure is applied until the stop criterion ||**x***s*+1,*LN* <sup>−</sup> **<sup>x</sup>***s*,*LN*|| <sup>&</sup>lt; <sup>10</sup>−<sup>9</sup> 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,

$$
\mathbb{C}w = \sigma Bw,\tag{23}
$$

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

**r**

 **r**

�0.06 �0.04 �0.02 0 0.02 0.04 0.06

�1 �0.5 <sup>0</sup> 0.5 <sup>1</sup> �1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

�1 �0.8 �0.6 �0.4 �0.2 <sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> �1

 **z**

<sup>87</sup> Influence of Horizontal Temperature Gradients

�0.8 �0.6 �0.4 �0.2 0 0.2 0.4 0.6 0.8 1

Fig. 2. Basic state for Γ = 2.936, *η* = 0.0862, *δ* = 0.1 and *R* = 72.65. On the left velocity field

*δ* = 0.1. 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 *kc* = 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 *kc* = 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

 **z**

�0.8 �0.6 �0.4 �0.2 0 0.2 0.4 0.6 0.8 1

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

Fig. 3. Growing mode or eigenfunction at the instability threshold for Γ = 2.936, *δ* = 0, *η* = 0.0862 and *Rc* = 72.504. On the left velocity field **u**. On the right Isotherms of

�1 �0.5 <sup>0</sup> 0.5 <sup>1</sup> �1

**r**

�0.8 �0.6 �0.4 �0.2 0 0.2 0.4 0.6 0.8 1

**z**

�1 �0.8 �0.6 �0.4 �0.2 <sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> �1

**r**

�0.8 �0.6 �0.4 �0.2 0 0.2 0.4 0.6 0.8 1

to totally 3D structures.

temperature Θ.

angular velocity.

**4.2 Not contained fluid**

**u**. On the right Isotherms of temperature Θ.

on Convective Instabilities with Geophysical Interest

**z**

that if *σ*max < 0 the stationary state is stable while if *σ*max > 0 the stationary state becomes unstable. The control or bifurcation parameter is the Rayleigh number *R*. For fixed values of the parameters, in those cases Γ, *Pr*, *δ*, *β* or *η*, the critical values are the minimum value of *Rc* for which there exists a value of *k*, *kc*, such that *σ*(*Rc*, *kc*) = 0.

In order to test convergence of the method we include, as an example, the calculation of the critical value of the bifurcation parameter, *Rc*, and the critical wave number, *kc*, 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 *Rc* obtained at different order expansions, the relative difference between the expansions at <sup>26</sup> <sup>×</sup> 18 and 24 <sup>×</sup> 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 *δ* = 0 is *kc* = 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 *Rc* = 73.5.


Table 1. (*kc*, *Rc*) for different order expansions in *L* and *N* in the Chebyshev expansion (20) for a 3D fluid with constant viscosity, *η* = 0, aspect ratio Γ = 2.936 and *δ* = 0.1.

#### **4. Numerical solutions with geophysical interest**

#### **4.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 at infinite *Pr* number, aspect ratio Γ = 2.936, *η* = 0.0862, *R* = 72.650 and **z**

6 Will-be-set-by-IN-TECH

that if *σ*max < 0 the stationary state is stable while if *σ*max > 0 the stationary state becomes unstable. The control or bifurcation parameter is the Rayleigh number *R*. For fixed values of the parameters, in those cases Γ, *Pr*, *δ*, *β* or *η*, the critical values are the minimum value of *Rc*

In order to test convergence of the method we include, as an example, the calculation of the critical value of the bifurcation parameter, *Rc*, and the critical wave number, *kc*, 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 *Rc* obtained at different order expansions, the relative difference between the expansions at <sup>26</sup> <sup>×</sup> 18 and 24 <sup>×</sup> 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 *δ* = 0 is *kc* = 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

*N* = 12 *N* = 14 *N* = 16 *N* = 18

*L* = 14 (2.5, 1203.91) (2.5, 1210.00) (2.5, 1212.73) (2.5, 1208.70) *L* = 16 (2.5, 1220.00) (2.5, 1214.00) (2.5, 1214.92) (2.5, 1214.94) *L* = 18 (2.5, 1220.10) (2.5, 1225.00) (2.5, 1224.92) (2.5, 1224.07) *L* = 20 (2.5, 1220.10) (2.5, 1224.15) (2.5, 1224.90) (2.5, 1224.90) *L* = 22 (2.5, 1220.20) (2.5, 1224.92) (2.5, 1224.92) (2.5, 1224.92) *L* = 24 (2.5, 1220.20) (2.5, 1225.00) (2.5, 1225.00) (2.5, 1224.92) *L* = 26 (2.5, 1220.20) (2.5, 1225.00) (2.5, 1224.92) (2.5, 1224.92)

Table 1. (*kc*, *Rc*) for different order expansions in *L* and *N* in the Chebyshev expansion (20)

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 at infinite *Pr* number, aspect ratio Γ = 2.936, *η* = 0.0862, *R* = 72.650 and

for a 3D fluid with constant viscosity, *η* = 0, aspect ratio Γ = 2.936 and *δ* = 0.1.

**4. Numerical solutions with geophysical interest**

for which there exists a value of *k*, *kc*, such that *σ*(*Rc*, *kc*) = 0.

the same value *Rc* = 73.5.

**4.1 Contained fluid**

Fig. 2. Basic state for Γ = 2.936, *η* = 0.0862, *δ* = 0.1 and *R* = 72.65. On the left velocity field **u**. On the right Isotherms of temperature Θ.

*δ* = 0.1. 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 *kc* = 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 *kc* = 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.

Fig. 3. Growing mode or eigenfunction at the instability threshold for Γ = 2.936, *δ* = 0, *η* = 0.0862 and *Rc* = 72.504. On the left velocity field **u**. On the right Isotherms of temperature Θ.

#### **4.2 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.

on Convective Instabilities with Geophysical Interest


<sup>89</sup> Influence of Horizontal Temperature Gradients

 -

-

t = 0

a) b)

Fig. 5. Clockwise vortex at Γ = 0.5, *δ* = 10 and *R* = 4367. a) Isotherms of Θ; b) meridional velocity (*ur*, *uz*); c) contour plot of the pressure *p* ; d) contour plot of the radial velocity component *ur*; e) contour plot of the vertical velocity component *uz*; f) contour plot of the azimuthal velocity component *uφ*. The contours correspond to equally spaced values within


their ranges of [-9:1] for <sup>Θ</sup>, [-0.02:20.4] ·103 for *<sup>p</sup>*, [-9.3:4.3] for *ur*, [-0.6:14.5] for *uz* and

Fig. 6. Track of a particle in the fluid for the stable clockwise vortex. The values of the parameters are Γ = 0.5, *δ* = 10 and *R* = 4367. a) Starting point at (*r* = 0.085, *z* = 0.05, *φ* = 0);


 -

[-26.1:0] for *uφ*. The pressure *p* is determined up to a constant.

t = 0

b) starting point at (*r* = 0.31, *z* = 0.05, *φ* = 0).

 - - - - - - -


Fig. 4. Growing mode or eigenfunction at the instability threshold for Γ = 2.936, *δ* = 0.1, *η* = 0.0862 and *Rc* = 72.650. On the left velocity field **u**. On the right Isotherms of temperature Θ.

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 *kc* = 0 and *u*˜*<sup>φ</sup>* � 0 is a 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,

$$\frac{dr}{dt} = \mathfrak{u}\_r(r, z),\tag{24}$$

$$\frac{d\phi}{dt} = \mu\_{\phi}(r, z),\tag{25}$$

$$\frac{dz}{dt} = \mathfrak{u}\_z(r, z). \tag{26}$$

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 *uz* 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 *ur* becomes positive and *uz* 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<sup>φ</sup>* is stronger and the spiral up motion of the particle starts as soon as the particle begins to move.

**z**

8 Will-be-set-by-IN-TECH

**z**


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 *kc* = 0 and *u*˜*<sup>φ</sup>* � 0 is a non-zero azimuthal velocity component. The fluid inside the annulus begins to move in the

The linear stability analysis of the vortical structures shows that there is a wide range of

The track of a particle in the vortex can be obtained by integrating the evolution of the element

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 *uz* 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 *ur* becomes positive and *uz* 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<sup>φ</sup>* is stronger and the spiral up motion of

*dr*

*dφ*

*dz*

Fig. 4. Growing mode or eigenfunction at the instability threshold for Γ = 2.936, *δ* = 0.1, *η* = 0.0862 and *Rc* = 72.650. On the left velocity field **u**. On the right Isotherms of

 - - - - - - - -

parameters for which this state is stable.

of fluid which follows the velocity field,


temperature Θ.

**z**

azimuthal direction, rotating around the inner cylinder.

the particle starts as soon as the particle begins to move.

**r**

 - - - - - - - -

**r**

*dt* <sup>=</sup> *ur*(*r*, *<sup>z</sup>*), (24)

*dt* <sup>=</sup> *<sup>u</sup>φ*(*r*, *<sup>z</sup>*), (25)

*dt* <sup>=</sup> *uz*(*r*, *<sup>z</sup>*). (26)

 - - - - - - -


Fig. 5. Clockwise vortex at Γ = 0.5, *δ* = 10 and *R* = 4367. a) Isotherms of Θ; b) meridional velocity (*ur*, *uz*); c) contour plot of the pressure *p* ; d) contour plot of the radial velocity component *ur*; e) contour plot of the vertical velocity component *uz*; f) contour plot of the azimuthal velocity component *uφ*. The contours correspond to equally spaced values within their ranges of [-9:1] for <sup>Θ</sup>, [-0.02:20.4] ·103 for *<sup>p</sup>*, [-9.3:4.3] for *ur*, [-0.6:14.5] for *uz* and [-26.1:0] for *uφ*. The pressure *p* is determined up to a constant.

Fig. 6. Track of a particle in the fluid for the stable clockwise vortex. The values of the parameters are Γ = 0.5, *δ* = 10 and *R* = 4367. a) Starting point at (*r* = 0.085, *z* = 0.05, *φ* = 0); b) starting point at (*r* = 0.31, *z* = 0.05, *φ* = 0).

We have distinguished two cases, a first one where the fluid is simply contained in a domain,

<sup>91</sup> Influence of Horizontal Temperature Gradients

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

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

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,

[2] Lord Rayleigh. On convective currents in a horizontal layer of fluid when the

[3] S. Chandrasekhar. *Hydrodynamic and Hydromagnetic Stability*. Dover Publications, New

[4] E. Bodenschatz, J. R. de Bruyn, G. Ahlers, and D. S. Cannell. Transitions between

[5] S.W Morris, E. Bodenschatz, D. S. Cannell and G. Ahlers. Spiral defect chaos in large aspect ratio Rayleigh-Bénard convection. *Phys. Rev. Lett.* 71, 2026-2029, 1993. [6] E.L. Koschmider. Bénard Cells and Taylor Vortices. Cambridge University Press, 1993. [7] M. Assenheimer, V. Steinberg, "Transition between spiral and target states in

[8] E. Bodenschatz, W. Pesch and G. Ahlers. Recent developments in Rayleigh-Bénard

[9] B. B. Plapp, D. A Egolf, E. Bodenschatz and W. Pesch. Dynamics and selection of giant

spirals in Rayleigh-Bénard convection. *Phys. Rev. Lett.* 81, 5334-5337, 1998.

temperature is on the under side. Phil. Mag. 32, pp. 529-46 (1916).

patterns in thermal convection. *Phys. Rev. Lett.* 67, 3078-3081, 1991.

Rayleigh-Bénard Convection". *Nature* 367 (6461) 345 (1994).

convection. *Annual Review of Fluid Mechanics* 32, 709-778, 2000.

and a second one where the fluid can flow throughout the boundaries.

of waves and the totally three dimensional patterns.

on Convective Instabilities with Geophysical Interest

[1] H. Bénard *Rev. Gén. Sci. Pures Appl.* 11, 1261 (1900).

phenomena as an instability.

**7. Acknowledgements**

which include RDEF funds.

York, 1981.

**8. References**
