Scale selection in low-order spectral models of two-dimensional thermal convection

The non-linear dynamics of two-dimensional thermal convection is investigated with the help of two low-order spectral models. The first model is the three-coefficient Lorenz model, in which the aspect ratio has not been used in scaling the variables. This means that the aspect ratio is a free parameter in addition to the Rayleigh number and the Prandtl number. The stability and the energetics of steady-state convection as a function of the three free parameters are investigated. The second model is a ten-coefficient model which describes the evolution and non-linear interaction of two horizontal scales of motion (with vertical wave number equal to one). It is found that the circulation mostly chooses a steady-state equivalent to the steady-state convective solution of the Lorenz model, i.e., a roll with a definite aspect ratio, even when one or more of the seven other coefficients would grow according to linear theory. It is also found that in most cases, the larger scale is favoured over the smaller scale, as the Rayleigh number is increased, especially at low Prandtl number. This agrees qualitatively with laboratory experiments. The paper is, however, principally meant as a contribution to the physical and mathematical insight into, as well as an illustration of, the non-linear dynamics of two-dimensional thermal convection. DOI: 10.1111/j.1600-0870.1984.tb00262.x


Introduction
It has been noted by among others &gel (1962), that thermal convection at moderate Rayleigh numbers is usually dominated by one harmonic, while the well-known linear theory of Rayleigh (1916) allows for a relatively broad spectrum of exponentially growing disturbances.If all these disturbances would indeed grow, the convection pattern would not be so simple as it apparently is.The spectacular selection of one scale of motion must be a non-linear process.The non-linear terms somehow act to damp all but one of the linearly unstable disturbances.
According to the Rayleigh theory, the most preferred aspect ratio (horizontal wavelength, A, divided by height, h) of convection cells and rolls should be about 3, depending somewhat on boundary conditions (Chandrasekhar,196 I).This is, however, not what is observed in the atmosphere, nor in the laboratory.Convection cells and rolls are observed daily on satellite photographs.They mostly form over the ocean to the west of mid-latitude depressions.In these regions, polar air flows out over relatively warm water.In the first instance, two-dimensional roll cloud patterns are observed.The aspect ratio of these rolls increases gradually as the air travels southward over the ocean.Eventually the rolls become unstable and three-dimensional cells appear.The aspect ratio of these cells has a mean value of 30, while rolls With an aspect ratio of 7 have been observed just before they become unstable (Walter, 1980) (see Fig. I).These values are clearly much larger than predicted by the Rayleigh theory.However, Rayleigh theory is only valid at the onset of convection, when the non-linear advection terms in the governing equations are small.As time increases, the non-linear terms are of increasing importance and they may distribute energy to other scales of motion.The increase in aspect ratio in atmos-Tellus 36A (1984). 5AIR FLOW

-1
Alh-2.0 hlh -6-7 Alh -15 Fig. 1.North-south cross-section view of a typical polar outbreak over a warm ocean showing aspect ratios of rolls and cells for different distances from the ice edge (adapted from Walter, 1980).
pheric convection has, however, been attributed by many workers (e.g.Sheu er 01.. 1980) to anisotropy of eddy diffusivity, large-scale subsidence and latent heat release.Yet, in laboratory experiments of thermal convection, where all of these physical factors are absent, the aspect ratio of rolls increases to more than double the critical aspect ratio as the Rayleigh number is increased (Willis et al., 1972; Koschmeider, 1974).It is therefore very likely that the non-linear terms also play a significant r6le in the scale-enlargement process.Segel (1962) investigated the non-linear interaction of two disturbances in two-dimensional convection theoretically and found that the equilibrium state contained only one of the two linearly unstable disturbances.He called this a "pure state".A "mixed" equilibrium state, containing both disturbances, could only result if one of the disturbances was linearly stable.However, it is not clear from Segel's results which non-linear processes are responsible for this selection.
In this paper, I intend to investigate and illustrate the problems of non-linear scale selection with the help of simplified non-linear models of roll convection.The two basic partial differential equations, describing two-dimensional Rayleigh-Benard convection, will be transformed to wave number space, leaving us with an infinite set of ordinary differential equations describing the time evolution of the amplitude of the stream function and the temperature of an infinite number of wave vectors.From this system, a limited number of equations, describing the evolution of the most energycontaining components is selected.This system of equations is called a low-order model.The steady Tellus 36A (1984).5 state solutions to low-order models are also steady state solutions to higher-order models.The uncertainty lies in the stability of these solutions to infinitesimal perturbations in the neglected components.
Low-order models of the atmosphere were originally conceived to illustrate some of the effects of non-linearity (Lorenz, 1960; Kiillen and Wun Nielsen, 1980).but are now increasingly used to study specific problems such as blocking (Charney and Devore, 1979).The reason that many meteorologists have taken their refuge in low-order models without knowing for sure whether the neglected components are important or not, is that they appear to be successful in, at least qualitatively, explaining and clarifying certain non-linear phenomena in meteorology.By defining a low-order model, we can isolate a certain physical effect from the total system.In this way, we will investigate the interaction of one scale of motion with the temperature stratitication (Section 6) and the interaction of two scales of motion and the temperature stratification (Sections 7-10).Because of their relative simplicity, low-order models give considerably more qualitative insight into the problem than, for example, the perturbation approach (Malkus and Veronis, 1958; Kuo and Platzman, 1961) or the complicated numerical computations based on a Galerkin method reported by Clever and Busse (1974).A very restrictive assumption which is made in the study by Kuo and  Platzman (1961) and many other such studies is that only one fundamental wave mode is linearly unstable, while all other wave modes are excited by the non-linear interactions.A relatively strongly forced system, such as the atmosphere, obviously allows for many linearly unstable wave modes which can grow independently of each other.This introduces the possibility of multiple steady states.The selection of a particular steady state is then dependent on the history and the stability properties of the possible steady states.
The first model we will study is the threecodlcient Lorenz (1963) model.This is the simplest possible non-linear model of convection.It models the interaction of one harmonic with the temperature stratitication.In spite of the fact that an immense body of literature on this model has appeared (e.g.Sparrow, 1982), some interesting properties with physical relevance have not yet been investigated, primarily because Lorenz partially scaled the aspect ratio out of the equations, so that no meaningful results can be deduced about the scale dependence of growth rates and energetics without rescaling the equations.
The second model is a ten-cdticient model with which the interaction of two horizontal scales of motion can be investigated.Agreement is found with &gel' s (1962) results in that the flow tends to choose a state in which all modes not pertaining to a certain, non-linearly selected, scak of motion decay to zero, even if they would grow according to linear theory.The steady-state convective solution is almost always equivalent to the steady-state convective solution of the Lorenz model, i.e., a symmetric roll circulation with a definite aspect ratio.Which mode is selected depends strongly on the external parameters of the problem (the Prandtl number and the Rayleigh number) and the initial conditions.There is a preference for a larger scale at high Rayleigh number and low Prandtl number.

, To basic equations and boundary conditions
We shall consider convection in two dimensions where w is the stream function, g is the acceleration due to gravity, a is the thermal expansion coefficient, v and k are the (constant) eddy diffusion coefficients of momentum and heat, respectively, W = J(0, w) and Y = J(V2w, v), J being the Jacobian operator and W is advection of temperature and Y is advection of the y-component of the vorticity (=VY).Eq (2.3) and (2.4) are based on the Boussinesq approximation (e.g.Spiegel and Veronis, 1960).Their derivation can be found in Saltzman (1962).
We will assume stress-free boundary conditions at L = 0 and L = h, i.e.Lateral boundary conditions are cyclic.

Spectral repmentation
In this section a spectral representation of the eqs.(2.3) and (2.4) will be derived.From this representation, truncated models, describing the evolution in time of only a few scales of motion will be derived.We will first cast the equations in non-dimensional form, by writing the variables of the problem in terms of non-dimensional variables, denoted by primes, as follows: Let us assume that the stream function wand the temperature 8 can be represented as a sum of Fourier components having a fundamental wavelength nLlh = 2nla in the x-direction and 2n in the z-direction as follows: ( H a ) where / is the horizontal wave number, n is the vertical wave number and a is equal to 2hlL. a is therefore inversely proportional to the aspect ratio of the domain, L l h .and 8 are formally defined in the range 0 Q x < 2nla and -n < z Q n.The region of physical interest is 0 Q z < n.We have included a mirror image domain (-n < z < 0 ) because we   (1984). 5 want to describe, among other things, a fundamental flow with half a wave length in the vertical, while retaining orthogonality of the basis functions.The coefficients wr.,(l) and 8,&) are assumed to be real.A factor i has been incorporated in (3.5b) so that ly and 8 of like horizontal wave number will be 90° out of phase, which is required in the linear solution, but not necessarily in the non-linear solution.
A convenient simplification of the notation, introduced by Kuo and Platzman (1961).results if we regard the pair of integers (1.n ) as a vector.The notation 8, where a = (I, n) may then be used in place of Sub, I) = exp I i(1a.x+ nz)}.

The expansion can then be written as
We further define the function where lu means a sum over all integral lattice points in the plane of a = (/,n).The orthogonality of the S,, may be expressed in the form where the asterisk designates a complex conjugate, while 6 is the Kronecker delta.The integration extends over the region 0 < x < Znla, -n < z < + n and do is an area element divided by the total area 4n21a.Substitution of (3.7a.b) into the governing equation (3.2a.b) and projection of the result on a certain wave vector y gives a set of simultaneous ordinary differential equations for the joint evolution of and 8,of the form where k! = a2 f; + n; and Ra = R / ( d ) .By scaling length with h l n and defining a modified Rayleigh number, Ra. as above, we shall see that all factors n disappear out of the spectral equations.
The heat advection spectrum W , and the vorticity advection spectrum V , are found by projecting W and Y on the wave vector y as follows: W,= 1 -is: Wdo, (3.15) From (3.13) it can be seen that L;@ = -L*. (3.16) From this, we can deduce that the interaction of a wave component with itself to itself is impossible.
The boundary condition at L = 0 is fulfilled only if

Energetics
The total kinetic energy is defined as  Adding together the energy contents of all the wave components and using the symmetry relations (3.17a.b) and (3.18a.b) and assuming no mean horizontal flow.we obtain the total kinetic and potential energies as follows: (4.8) The first term on the right-hand side of (4.9) represents the potential energy of the mean temperature stratification, which is not directly available for transfer into kinetic energy.The Tellus 36A (1984). 5second term on the right-hand side of (4.12) is the available potential energy, AP.
Before concluding this section, we mention an interesting parameter which measures the importance of convection, namely the Nusselt number.Nu.N u is the ratio of the actual steady state heat transport rate to the conductive heat transport rate.It is given by (Kuo and Platzman,196  where the angular brackets denote a horizontal average.Substituting into (4.10) the expansion (3.5b) and using (3.I7b) we obtain, N u = I + I x he,. . .

The stability of the state of rest (4. I I )
It is well-known (see e.g., Kuo and Platzman, 1961) that the system (3.9a.b) becomes linearly unstable for all those modes which have a so-called critical modal Rayleigh number, Racay, which satisfies the criterion (5.1) Ra,. is a minimum when a/,. = ";.Id2 and n;.= 1.
This means that a roll with a horizontal lengthscale which is 2 \/i times greater than h is the first unstable convection roll when Ra is increased from zero.It becomes unstable at Ra = 2714.
We say that all wave modes y, for which Racay < Ra, are sevexcited.These modes can draw energy directly from the mean unstable temperature stratification.Other modes will only be able to receive energy through the non-linear terms.It can be seen from (5.1) that all modes with n,.> I have a critical modal Rayleigh number which is at least one order of magnitude greater than the modes with ni.=I (the so-called "single vertical modes").
The first unstable mode with ny = 2 corresponds to al; = 2 1 4 and Ra,? = 108.Therefore, because they are self-excited, the single vertical modes will, in general, play an active r61e at relatively low Ra, while all the double vertical modes (n, = 2) are passive for Ra < 108, because they can only receive energy from the former modes through the non-linear terms.Higher order vertical modes (n, > 2) are progressively less important, because dissipation increases with nf while forcing does not increase.

The interaction of one scale of motion with the mean temperature stratification: the Lorenz model
In this section, we will investigate the interaction of one scale of motion with the mean temperature stratification.To this end, we truncate the expansions (3.5a.b) by basically considering only the following wave vectors: (I, I) in the W-expansion and (I, I ) and (0,2) in the &expansion.Of course, we have also to consider the mirror images in the Iand n-axis of these wave vectors, but the coefficients corresponding to these wave vectors can be eliminated with the help of symmetry relations (3.17a.b) and (3.18a.b).We will therefore derive ordinary differential equations describing the time-evolution of the following three coefficients: The model which results is equivalent to the Lorenr (1963) model.It describes a completely symmetric roll circulation with an aspect ratio ( L l h ) of 2/a (see Fig. 3).We can adjust the aspect ratio of the roll (by varying a), but not its shape.From eqs. (3.9a, b), we can easily deduce the non-linear equations giving the time evolution of X, Y and 2.
The only non-linear interactions involved are the following: (I, -I ) and (0.2) interact to contribute to the advection term corresponding to wave vector (I, I), and (I, I ) and (-1.I) interact to contribute to (0.2).Other interactions are not permitted due to selection rule (3.15).We ignore all interactions with wave vectors outside the truncation.In this way, the non-linear terms are energy conserving, which is a property of all low-order spectral models in which advection terms are the only nonlinearities (Lorenz, 1982).We then obtain the following model equations: The system is forced through the first term on the right-hand side (r.h.s.) of (6.Ib).This feeds into the equation of motion through the first term on the r.h.s. of (6.la).This in turn leads to stabilization of the mean vertical temperature profile through the non-linear term in (6.1~).as &.* represents the departure of the vertical temperature stratification from the initial linear variation.Through the non-linear term in (6.Ib).this inhibits further growth of the temperature perturbation Y or of available potential energy.In this way, the system is driven to steady-state convection in which a non-linear mean temperature profile is maintained against dissipation.Eqs.(6.la, b, c) differ slightly from the original Lorenz equations because Lorent used a to scale temperature and stream function.He thereby partly scaled the aspect ratio out of the problem.
To find the steady-state solutions of (6.la, b, c), we set all time derivatives equal to zero.We see directly that one possible steady state is (X, Y, 2) = (0, 0, 0).i.e., the state of rest.The stability of rest can be investigated by linearizing the equations (6.la.b, c) around this stead,y state, which can then be written compactly as X = MX, where the dot represents a time derivative and . .where r = (Ra -Ra,)1'2.These two solutions represent steady-roll circulations with opposite flow directions.We will call this steady circulation a "Lorent roll".It is actually equivalent to what Segel (1962) called a pure state.Solutions (6.4a,b) are also exactly equivalent to the second-order steady-state convective solution obtained in the perturbation approach by Kuo and Plattman ( 196 I).
The steady-state kinetic energy is found by substituting the steady state value of X into (4.8): where Ra, is given by (6.3).The steady-state kinetic energy is therefore a function of u and Ra and independent of u.It is a straightforward matter to calculate the position of the maximum of K on the u-axis as a function of Ra.The result is plotted in Fig. 4. Also shown is the position of the maximum of the growth rate 1, as a function of u and Ra for u = I, u = 10 and u = 50 (there is only one maximum in 1, for each value of Ra!).The maximum growth rate curves are symmetric about u = I, that is, the curve for u = a,, is the same as the curve for u = l/uw Even though the maximum growth rate shifts to higher u (lower aspect ratio) as Ra is increased, the maximum value of the steady-state kinetic energy shifts to lower a (higher aspect ratio) as Ra is increased.This can be understood by noting that the non-linear terms in (6.lb.c) contain a factor a. The consequence of this is that a roll with low u does not stabilize the mean temperature profile as quickly as a roll with high u.In other words, the non-linear terms moderate the exponential growth rate less effectively in larger scales of motion.At Ra = 100 (in the atmosphere this would be R = 1 0 0 ~' 5 lo4, which is a typical order of magnitude of Ra in a polar outbreak (Krishnamurti, 1975)), the maximum in K is located at u = 1/3, which corresponds to an aspect ratio of 6.If the circulation chooses its aspect ratio such that the kinetic energy in it is a maximum, this could explain the increase of the aspect ratio of rolls as Ra increases.This decrease in u in conjunction with an increase in Ra is also observed in polar outbreaks of air over a warm ocean (Walter, 1980) (see Fig. I).As the air flows out over the ocean, the depth of the convective layer, h, increases steadily due to penetrative convection.Since Ra is proportional to h4, Ra increases as the air travels farther over the ocean.
Walter reports values of Uh increasing from 2.8 to about 7 as cold air travels southwards over the Bering Sea.Eventually the rolls go over into three dimensional cells whose aspect ratio keep increasing.At this stage, latent heat forcing due to condensation and evaporation undoubtedly plays a significant r6le.To find out whether the flow indeed chooses the larger scale, we have to study the competition between different scales of motion.This is not possible with the Lorent model because it contains only one scale of motion.By substituting the steady state value of Y into (4.9).we find the steady state available potential energy: TeUus 36A (1984).5 If we substitute (6.5) and (6.6) into (6.7).we find that I uRa, e=-.
Therefore, kinetic energy is produced most efficiently from available potential energy if a = l / f i ( L / h = 2 fi) for all Ra.
Following the theory of Malkus (1954).it has often been assumed that the preferred aspect ratio is that which maximizes the heat flux, or Nusselt number.Substituting Z = rY2 into (4.1 I) we get the following expression for Nu: This is exactly equivalent to the expression for Nu obtained from second-order perturbation theory by Malkus and Veronis (1958).From (6.9).we see that Nu is maximized when a = I / v ~ for all Ra.Thus, the heat flux will not be maximized (in two-dimensional convection) if that scale is selected which maximizes the kinetic energy.It will, however, be maximized if that scale is chosen which makes the conversion of potential energy into kinetic energy most efficient.The steady state heat transport is also independent of u.
The above heuristic selection principles cannot give us a definite answer to the question of the preferred scale of motion.They give different answers, although a larger scale than the first self-excited scale (a = ~/fi) is sometimes preferred, while a smaller scale is never preferred.Another approach to this problem is stability analysis.The stability of steady convection to infinitesimal perturbations can be investigated by The eigenvalues of M as a function of u and Ra can be determined numerically by the standard method of Eberlein (see Wilkinson and Reinch, 1971).The result for u = 10 (which is the U-value which was chosen by Lorenz) is presented in Fig. 6.Only regions where all real parts of the eigenvalues are negative (stable steady convection) or regions where at least one real part of an eigenvalue is positive (unstable steady convection) are shown.We find that steady convection is unstable when Ra 2 167 for a roll with a = I/&.Lorenz analytically found a value of p = 24.74(Ra I 4 p ) for the onset of instability of steady convection at a = l / f i and u = 10, which means that our result is in agreement with that of Lorenz.We also find that the minimum value of Ra (=163) for instability of steady convection is located at a = 0.77 ( L / H = 2.6).The curve separating stable steady convection from unstable steady convection is dependent on u, while the curve separating steady convection from rest is not.It is also evident from Fig. 6 that the Lorenz roll with a larger aspect ratio than 2 4 is stable for larger Ra.steady convection is stable for all Ra is verified in Fig. 7.
As was noted before, the curve separating steady from unsteady convection is dependent on u.The value of a where steady convection first becomes unstable as Ra is increased, urn*, is therefore also a function of u and is plotted in Fig. 8, as well as the value of Re, Ramin, at which this happens.We see that urnin shifts to higher wave numbers as u is decreased and it is very tempting to say that amIn eoes to V f i as u goes to infinity.
However, the numerically calculated points are not located exactly on the smooth curve.This could be a numerical inaccuracy, although until now, all Lorenr's analytical results have been verified with relatively very high accuracy.Besides, the values of Ramin do lie on a smooth curve.Nevertheless, we 0 1 0 2 0 3 0 4 0 5 0 U Fig. 8. amin (0) (the value of u at which steady convection first becomes unstable when Ra is increased) and Ramin (x) (Ra-value at which this happens) as a function of u.
can say that for all u < 50, amln > I/&.We can thus conclude that a Loren2 roll with an aspect ratio just greater than 2 f i will probably be favoured over a Lorenz roll with an aspect ratio just smaller than 2 f i as Ra is increased.
For Ra > 167, u = 10, a = l / f i , for example, there are no steady states.All trajectories must however stay in a finite region of (X, Y, 2)-space, because the system is dissipative.Lorenr in- vestigated the behaviour of the system at Ra = Rae, u = 10 and u = I/@ He found that the solution was non-periodic and that it stayed in a complicated region in (X, Y, 2)-space.This region, later termed a "strange attractor", looks roughly like a two-dimensional surface.Any solution wanders through this region and, given sufficient time, will eventually pass arbitrarily close to every point in this region.Dynamical processes which exhibit this property are sometimes called "chaotic".Lorenz concluded that these findings were not very hopeful for medium-to long-range deterministic weather prediction.

Tellus 36A (1984). 5
Almost all literature on the Lorenz model has concentrated on the chaotic properties of its solution.

A low-order model including interaction between different scales of motion
In this section, we will derive a set of equations describing the interaction between two horizontal scales of motion or, more specificdly, the interaction between two single vertical modes with neighbouring horizontal wave numbers, equal to j and j + 1.According to the selection rule (3.15).
these wave modes, however, cannot exchange energy with each other without the mediation of a double vertical mode.The only modes which are excited due to the interaction of (j, I ) and ( j + I, I) are (I, 2) and (2 j + I, 2).We will therefore also include these wave vectors in our model.It should be stressed that we have assumed no inflow at the lateral boundaries, which means that the mean shear coefficients, vO.,,, must be equal to zero.We will also take (0.4) in addition to (0.2) for the temperature, because this gives a more realistic representation of the vertical temperature stratification, especially at high Ra.It also ensures that each wave mode has a direct stabilizing influence on the mean temperature stratification.O0., and OOs3 have been left outside the truncation because they separately imply that the mean temperature perturbation 8 is non-zero.Moreover, these modes can only be excited when we choose j such that we have a single vertical mode and a double vertical mode with identical horizontal wave numbers.This is only the case when j = 1.Note that we have not confined ourselves exclusively to even-parity modes (if (I, + n,) is even, then y has even parity; if (I, + ny) is odd, then y has odd parity).In many non-linear studies (e.g., Kuo and Platzman,196 I ; Veronis, 1966) the ( I, I ) mode (an even parity mode) is considered to be the only self-excited mode, while all other modes are created by non-linear interactions.According to (3.15).the odd-parity modes would then never be excited.Therefore odd-parity modes are ignored in these studies.An effect which is thereby neglected is that there could very well be odd-parity modes which are self-excited.In the following sections, we will Tellus 36A (1984).5 investigate the interaction between two secondorder Kuo and Platzman (1961) solutions of opposite parity.
The coefficients are defined as follows: The system of ordinary differential equations describing the evolution in time of these ten coefficients is: The first four equations describe the velocity field and the last six equations describe the temperature field.A good method to check whether we have made any mistakes in deriving these equations is to verify that the non-linear terms conserve energy as defined in Section 4.
There are three types of non-linearities in this model; (i) advection of vorticity or redistribution of kinetic energy between different scales of motion, described by the non-linear terms in (7.1) to (7.4); (i) advection of temperature, which can be split into two types: (a) redistribution of available potential energy between different scales of motion, described by the first four non-linear terms in eqs.(7.5) and (7.6) and the first two non-linear terms in eqs.(7.7) and (7.8);(b) stabilization of the mean temperature stratification, described by the remaining non-linear terms in (7.5) to (7.8) and all the non-linear terms in (7.9) and (7.10).
This model therefore contains all the different physical mechanisms which are responsible for non-linear behaviour in thermal convection.The Lorenz model only contains the third type of non-linearit y.
Note that a (the aspect ratio of the domain) and j (the wave number of the parameters.If we let a = 1/2 ( j + I, I) or (2, I) wave mode corresponds to the first self-excited mode with Ra, = 6.75.The (j, I)or (I, I) wave mode has double the aspect ratio and becomes linearly unstable at Ra = 11.39.The two double vertical modes (I,2) and ( 2 j + 1,2) or (3.2) have critical modal Rayleigh numbers of 561.5 and 119.65, respectively.Because these wave modes are linearly damped up to such high Rayleigh numbers, they will not play such an active r6le as the two single vertical modes.They actually serve as a (non-linear) syphon of energy between the two competing single vertical modes.According to (6.3). the (2, I) wave mode will have the greatest growth rate for all values of Ra.In Section 6, we saw that as the Rayleigh number is increased, the maximum in the kinetic energy shifts to rolls with greater aspect ratio.We also saw that a roll with large aspect ratio was innerly stable for larger Ra than a roll with small aspect ratio.From this, we suspect that at a certain Ra, the (I, I) wave mode will be favoured over the (2, I) wave mode and will absorb the energy in the (2, I) wave mode through the non-linear terms, because, in contrast to the Lorenz model, this model has non-linear terms which describe energy transfer from one scale to the other.
The If we choose a = I / 3 d and j = 3 and ignore the coefficients D, H and J, the system reduces to the 7-coefficient Saltzman ( 1962) model.Since the system of eqs.(7.1)-(7.10) is highly non-linear, it is not easy and straightforward to determine all steady states and their stability.We can however, by simplifying the equations in a certain way (Section 8) and by numerical experiments (Section 9).gain some insight into the qualitative behaviour of this model and see which steady states are important.After that (in Section 10) we will determine the stability of these solutions as a function of u and Ra for certain choices of a and j, more accurately.

An approacb to non-linear scale selection
We will first investigate the qualitative behaviour of the first type of non-linearity (vorticity advec-TeUus 36A (1984).5 tion), which is dominant when u is smaller than one, using some ideas developed by Haken ( 1978).This will give us some qualitative idea of how non-linear scale selection works.Eqs.(7.1H7.4),describing the velocity field, can be written as follows: where lA, l a , & and l o are so-called relaxation "constants", through which we have parameterized the linear terms, namely, the dissipation and the linear effect of the temperature field on the velocity field.If we do not ignore the non-linear terms in (7.4)-(7.10),we cannot say much about the relaxation "constants", which are then far from constant.If we do ignore these terms, we can say that, as long as Ra is smaller than the critical modal Rayleigh number of either of the double vertical modes and greater than the critical modal Rayleigh number of either of the single vertical modes, 1, and I , are positive and & and l o are negative.If we assume u = 1 / 2 4 and j = 1, then this assertion is valid for 11.39 < Ra < 119.65.If we assume u = l / f i a n d j = 1, it is valid for 13.5 < Ra < 108.5.We can then apply a so-called "adiabatic elimination technique" developed by Haken (1978).This consists of eliminating the fast relaxing variables (in this case ( 2 j + I)' ( k j -&: ) p1 and p2 are certainly positive when k i > kj.This condition is always satisfied for the choices of a and j we will henceforth make, i.e., choices such that the scales we consider have aspect ratios in the neighbourhood of 2 4 .For instance, i f j = 1, then a < 1 for p I and p2 to be positive.Eqs.(8.7) and (8.8) describe the competition of the order parameters A and B among each other.Looking at these equations, we notice the remarkable fact that A, if great enough, can cause B to damp out, while B cunnor do the same thing with A. On the contrary, if B becomes large enough, it can excite A and therefore supress itself.From this, we conclude that it is very likely that there is a Rayleigh number range (high Ra) in which the largest scale ( A ) will take over all the energy from the smaller scale (B).In Haken's words, we say that the mode A then "slaves" the other modes.The fact that one mode can slave all the other modes, even if some of them are linearly unstable, is an intrinsically non-linear phenomenon which accounts for the surprising degree of order in many non-linear systems which are kept far from equilibrium.Eqs.(8.7) and (8.8).
however, only say something about how kinetic energy redistribution between different scales of motion by vorticity advection influences scale selection.We therefore expect that these equations are qualitatively valid only for low 6.An illustration of these equations will be shown in Section 9.
In the context of this section, the papers by Fj0rtoft (1953) and Merrilees and Warn (1975) on the changes in the spectral distribution of kinetic energy for two-dimensional non-divergent flow should be mentioned.Because both enstrophy and kinetic energy are conserved in two-dimensional non-divergent flow, more kinetic energy usually flows from an intermediate scale to a larger scale of motion than from an intermediate scale to a smalkr scak.Therefore the larger scale is favoured over the smaller scale, which is in agreement with our result.However, our result is based on a completely different argument.It is crucial in our deduction that some modes should be linearly damped, i.e. have negative A's.This is not needed in the inviscid theory of Fj0rtoft and its correction by Merrilees and Warn.

Some numerical integrations
We now have a qualitative idea of how nonlinear scale selection could work.In this section we turn to numerical integration of the system (7.1)-(7.10),showing that the larger scale indeed takes over at higher Raykigh numbers and illustrating the qualitative behaviour of the competition eqs.(8.7) and (8.8).
Eqs. (7.1)-(7.10)withj = I and u = 1 / 2 4 are integrated using Young's (1968) method A (doubk forward, centred) with a non-dimensional time step Al = 0.005 and with all coefficients initially equal to zero.We performed an integration of 1200 non-dimensional time-units and let Ra increase linearly in time from 0 to 200.We also imposed a random forcing in the form of finite-amplitude perturbations on the components A, B, C and D. Each time step a perturbation equal to & was added to these four components.The sign of the perturbation was picked at random.This numerical experiment is a crude approximation to what happens in a polar outbreak and also to the laboratory experiments performed by Willis et al. (1972).although the time scales are different.The Prandtl number is fixed at one, which is a realistic value for the atmosphere.In Fig. 9, the amplitudes of the stream function field wave vectors are plotted as a function of the Rayleigh number (or time).We see that there are three discreet transitions.In the beginning, only the smaller scak ( B ) grows to a finite value, while the others stay very close to zero, i.e. are slaved by B. The system in fact settles into a steady state which is equivalent to the Lorenz roll solution (6.5a) with u = I/& (aspect ratio equal to 24).The slaving of A by B must be due to temperature advection.Between Ra = 66 and Ra = 76, the system suddenly evolves to a new solution which is equivalent to the Lorenz roll solution (6.5a) with u = 1 / 2 4 (aspect ratio equal to 4\/2).All stream function amplitudes except A (the larger scale) decay to zero, even B, which according to our linear analysis should have the greatest growth rate.Finally, at Ra 5 120, the critical modal Rayleigh number of the (3, 2) wave mode (coefficients D and H ) is reached.This apparently renders the existing pure large-scale circulation unstable, because B, C and D start to grow slowly, although A stays dominant until the end of the integration at Ra = 200.
Evidently the system in most cases prefers a pure state.In Section 10, we will proceed to determine the exact stability as a function of Ra and u of the two pure states which were selected in the above integration.
Before that, a nice illustration of the qualitative behaviour of the competition eqs.(8.1) and (8.8) will be presented.A numerical integration at low u(=O.I) and constant Ra(= 10) without random forcing was performed.The initial conditions consisted of small perturbations of the stream function field only, given by  dimensional time units.At Ra = 10, the larger scale ( A ) is linearly stable because it has a critical modal Rayleigh number of 11.39.This means that AA is slightly negative.As is positive because B has a critical modal Rayleigh number of 6.75.Therefore A is completely dependent on B for its energy and can indeed grow at the expense of B if B is large enough.However, it quickly depletes its own source of energy and starts to decay.This gives B another chance to grow, which, however, again leads to growth of A and subsequent decay of B.
The nearly periodic solution which ensues is analogous to the famous predator-prey model of Lotka and Volterra (e.g.Haken, 1978).

Stability analysis of the pun states
In this section, we will determine the exact stability as a function of Ra and u of the two Lorenr roll solutions or pure states possible in the 10-coefficient model.This can be done in the familiar way of linearizing the equations around the steady state and looking for Ra-and U-values for which all real parts of the eigenvalues are negative.
The 10 by 10 matrix M corresponding to eqs.where the steady state is given by (6.5a).i.e., Lmatrix J M = with r = (Ra -RaJ1l2 and Ra, = k;/(a* j3, or by (6.5b).but this should, and indeed does not, make any difference in the eigenvalues.
If we are determining the stability of the smaller scale (the (j + 1, I) wave mode), the 6 by 6 matrix has the form, * [:2] -,,i The 3 by 3 matrix in M is given by (6.10).The matrix in the middk of M correspond to the stability of the particular Lorenz roll to infinitesimal perturbations in other scales of motion.We will call these "outer" perturbations.
The last eigenvalue is equal to -16, which expresses the fact that the vertical temperature stratification coefficient 5(6)o4) acts in stabilizing way on the flow.We therefore only have to determine the eigenvalues of the 6 by 6 matrix.
If we are determining the stability of the larger scale (the (j, I) wave mode), the 6 by 6 matrix has the form: where the steady state is given by (6.5a).i.e., with r = (Ra -RaJ1'* and Ra = k",/(u*(j + The eigenvalues of MA and Me are determined numerically by the same method as used in Section 6.Again, Rayleigh number ranges are determined for which the solution is stable (all real parts of the eigenvalues are negative) or unstable (at least one eigenvalue has a real part which is positive).The  It turns out that there is a large region in parameter space in which at least two completely different stable flow patterns are possible.It depends on the initial condiions to which of the two solutions the system will evolve.
Another interesting result, which is in agreement with the qualitative result of Section 8, is that at low u and sufficiently high Ra, only the larger scale is stable.At u = I, the smaller scale becomes unstable to infinitesimal outer perturbations at Ra = 64.This is in agreement with the first numerical integration discussed in Section 9 (see Fig. 9).Also the point at which the larger scale becomes unstable at Ra = I 19 is correctly reproduced in this integration.This last bifurcation point coincides with the linear marginal instability of the coefficients corresponding to the wave vector Rayleigh numbers of 6.75 and 13.5, respectively.According to (6.2) (linear theory), the (2, I) wave mode becomes the fastest growing mode when Ra z 54.I, but the steady state kinetic energy of the large-scale Lorent roll is always greater than that of the small-scale Lorenr roll.We again determine the signs of the real parts of the eigenvalues of the matrices MA and Mg, but this time with u = l / f iand j = 1.The result is shown in Fig. 13.We again find a region in parameter space in which both the small-scale and the large-scale Lorenr rolls are stable.The large-scale Loren2 roll is the only stable solution for relatively low u and relatively high Ra.This is again in qualitative agreement with the competition eqs.

Summary and conclusions
We have studied the non-linear behaviour of two-dimensional thermal convection as a function of Rayleigh number and Prandtl number with the help of low-order spectral models.Although these models are a crude approximation to the full non-linear equations (for instance, in the atmosphere there is a continuous spectrum of wave numbers instead of a discrete spectrum), they do permit a fairly rigorous investigation of the different types of non-linearities which are responsible for various observational facts such as scale selection and enlargement.The solutions of these low-order models display a wide variety of possible types of behaviour.We have ignored the possibility of chaotic solutions and concentrated on the "orderly'' solutions and the mechanisms which can lead to transition from one such solution to the other.We have divided the non-linearities into three types, namely, redistribution of kinetic energy between different scales of motion, redistribution of available potential energy between different scales of motion and stabilization of the mean temperature stratification.The model described in Section 7 contains all these non-linear mechanisms.The Lorenz model (Section 6) contains only the last type, because it describes only one scale of motion.In spite of this, we can still infer from the Lorenz model that a mode with a larger aspect ratio than that of the first self-excited mode will probably be preferred over a smaller scale as Ra is increased.This, because the maximum steady-state kinetic energy of convection shifts to higher aspect ratios as Ra increases, in spite of the fact that, according to linear theory, the maximum growthrate shifts to lower aspect ratios as Ra increases.The reason for this is that the mean temperature stratification is stabilized less quickly in larger scales of motion.
Baroclinicly unstable waves show the same behaviour (Hart,198 I).According to linear theory, the growth-rate of unstable modes is greatest in large zonal wave numbers (tn 2 IS), but the nonlinear regime is characterized by tn < 10.In this case, the basic zonal shear apparently plays the same r6le as the temperature stratiilcation in the case of convection.The non-linear interaction of the growing wave with the mean zonal shear produces a correction to the mean zonal shear, thereby reducing the available potential energy Tellus 36A (1984).5 until there is a non-linear equilibration.This equilibration goes least quickly in larger scales of motion.
The question remains as to how realistic the Lorenz model is quantitatively.This can best be Seen by comparing the Nusselt number as given by (6.9) with observations.The problem is that there are hardly any observations of Nu He concluded that his calculations were reasonably accurate because they agreed to within 1% of the next higher approximation and because they agreed "almost exactly" with the experimental results.However, these experimental results say nothing about the absolute value of Nu and can also be very well fitted to eqs.(6.9) for Ra < 100.
The question of the heat transport in a layer with stress-free boundary conditions is thus still open.In terms of heat transport, the Lorenz model is not so unrealistic for Ra < 100 ( R < 10') as is commonly assumed.
The second model permits two horizontal scales of motion (with vertical wave number equal to one).
With this model, we have tried to find an answer to the following question: how docs it come about that convection at moderate Rayldgh numbers is usually dominated by one harmonic, whik many modes are linearly unstabk with growth rates which do not dfler much?The mode which wins can somehow take most advantage of the nonlinearities, leading to the complete decay of the other modes.Somehow, most researchers have taken this fact for granted.Intuitively we anticipate that the system should strive to the simplest possible steady-state solution, because this is probably the solution in which kinetic energy or heat transport or efficiency are maximized, or available potential energy is minimized.But how this comes about is not trivial.In the ten-coefficient model, we indeed find that for Ra c 100, the selected steady-state solution is almost always equivalent to the steady-state solution of the simplest possible non-linear model of convection: the Lorenz model.Segel (1962) callad this a "pure state".
We have determined the stability of two possible Lorenr roll solutions to perturbations in other modes for 0.1 < u < 50 and 0 < Ra < 200.We find that the larger scale is favoured over the smaller scale as Ra is increased, especially for u < 1. Vorticity advcction, especially, destabilizes the smaller scale of motion.The effect of redistribution of availabk potential energy on scale selection is more difficult to investigate, because it involves the non-linear interaction between two scalar fields, resulting in a greater number of non-linear terms in the spectral equations.We also find that there is a rather large region in parameter space (especially when u > I) in which both pure Lorenz roll solutions with different aspect ratios are stable to infinitesimal perturbations in the other modes.It then depends on the initial conditions or history as to which solution is chosen.
In other words, hysteresis is possible in the model.The roll aspect ratio as a function of the Rayleigh number for u = 6.7 and u = 0.71 as observed in the laboratory by Willis et al. (1972) is plotted in Fig. 14.It can be seen that when Ra was steadily increased in time (as in the numerical integration shown in Fig. 9 ) at u = 6.7, the aspect ratio (at the same Ra) differed from the case when Ra was decreased in time.This did not happen when u = 0.7 I. Thus, hysteresis is observed at u = 6.7 but not at u = 0.7 1.It can also be seen that the increase in aspect ratio, as Ra increases, goes faster at u = 0.71 than at u = 6.7.Our model results are compatible with these qualitative observations.
We have not investigated the stabTty of convection rolls to perturbations in the third dimension.This implies that the third dimension is not always needed for spontaneous changes to occur in the diameter of rolls.Therefore, the qualitative observation by Willis er al. (1972), who noted that ". . .mechanisms which can change the diameter of rolls imply the existence of structure which is three-dimensional", does not always have to hold.

Aclmowledgements
The many helpful comments and suggestions given by Erland Kiillen, Leo Maas, Hans Oerlemans, Cor Schuurmans and Huib de Swart are gratefully acknowledged.

(
x, L), between two rigid horizontal boundaries.At the lower boundary ( L = 0), the temperature To is maintained and at the upper boundary (z = h), the lower temperature To -AT is maintained.It is convenient to divide the temperature field into two parts in the following way:where The parameter r is called the lapse rate; cp is the specific heat of dry air at constant pressure; 6is the deviation from the linear temperature profile that would prevail in the absence of convection.Two equations describing two-dimensional dry convection are the vorticity equation and the thermodynamic equation, respectively: scalar form this means that I,= I. + Is, ny= na+ n, q.
that the horizontal dependence of w and 8 is represented by a sine and a cosine series, relations (3.I7a, b) and (3.18a, b) to eliminate the negative wave number amplitudes out of the truncated versions of eqs.(3.9a.b).

Fig. 2 .
Fig. 2.Qualitative picture of the dependence of the dissipation and the kinetic energy generation on the horizontal wave number.Only in a finite band of wave numbers around the critical wavenumber I, does the kinetic energy generation exceed the dissipation.

Fig. 3 .
Fig. 3. Stream-line panern of the circulation described by the Lorenr model.

A
This eigenvalue can be interpreted as the (exponential) growth rate of the intensity of convection when no non-linearities are taken into account.It is a complicated function of u, Ra and u.The growthrate is positive when The equations (7.2) have two other steady-state solutions

Fig. 4 .
Fig. 4.Position of the maximum growth rate (for dflwent a-values) and the maximum steady-state kinetic energy on the a-axis as a funaion of Ra.

Fig. 6 .
Fig. 6.The stability of steady convection as a function of u and Ra at u = 10.
10-coefficient model reduces to the Lorenz model if we ignore the coefficients B, C, D, F, G, H and J or A, C, D, E, G, H and J.It thus contains two possible Lorenz roll solutions simultaneously.
C and D) by assuming that they follow the self-excited modes (A and B) exactly (adiabatically).This means that C and D are always in equilibrium with A and B. In TeUus 36A (1984).5that case, we can assume that the time derivatives of C and D are zero.This leads to D are not self-excited, they are completely dependent on A and B for their energy.Therefore, they can never take the initiative, and thus must "obey" the orders of the so-called "order parameters" A and B. Neverthekss, C and D play a crucial r6e in that they serve as the non-linear coupling between the order parameters.In other words, all energy flowing from A to B, or vice versa, has to flow through C and D.Elimination of C and D in (8.1) and (8.2) with the help of (8.5) and (8.6) leaves us with, dB ---( l a -P Z A')

Fig. 9 .Fig. 10 .
Fig. 9.The coefficients A (solid line), B (broken line, long dashes), C (broken line, medium size dashes) and D (broken line, short dashes) as a function of Ra (or time) in a numerical integration of the loooeRicicnt model (with perturbations (=lo-') imposed on the four coefficients shown, with u = V 2 4 .j = I, u = I), in which Ra was increased linearly in time from 0 to 200 over I200 time units.TeUus 36A (1984).5 The time step was & = 0.01.The solution is shown in Fig. 10.It is almost periodic with a relatively long period of about 200 non-Tellus 36A (1984). 5 linearized around either one of the two steady states has the following structure: Lorenz roll, the stagiity of which we are determining.The eigenvalues of the 6 by 6 Fig. 11.The stability of the two possible Lorenz roll solutions to the 10-coefficient model as a function of u and Ra when u = 1/2\/2 and j = 1.Broken (solid) lines separate innerly (outerly) stable from innerly (outerly) unstable regions for the two Lorenz rolls.Thick lines correspond to the large scale Lorenz roll (A: aspect ratio = 4fi) and thin lines correspond to the small scale Loren2 roll (B: aspect ratio = 2fi).The letters A and B denote the stability of that particular mode in that particular region of parameter space.In the hatched regions there is no stable Lorenz roll.

(
2j + 1, 2). or (3, 2).which have a critical modal Rayleigh number of 119.65.Fig.I I shows the relative stability of two Lorenr rolls, one of which has double the aspect ratio of the other.With the scale parameters u andj we can, however, not only vary the absolute aspect ratio but also the relative aspect ratio of the two possible Loren2 roll solutions.The question is, whether the stability diagram changes qualitatively if we vary the relative scales.If we assume u = I/(( j + ~)fi), then the smaller scale ( B ) always has an aspect ration of 2 f i .The ratio, s of the aspect ratio of the larger scale to that of the smaller scale is given by j + I s=-.j We w i l l vary s (or j ) and see if the lines of marginal outer stability (the solid lines in Fig.I I ) change place significantly in parameter space.We repeated the stability analysis shown in Fig.I 1 (corresponding to s = 2 ) for s = 1.5, s = 1.25 and s = 1.1.The result is shown in Fig. 12.It turns out that the picture does not change much qualitatively as s decreases from 2 to 1.1.The lines of marginal outer stability of the larger scale (labeled 2 and 3 in Fig. 1 I ) hardly change location.A dip appears at u z 6 in the line of marginal outer stability of the smaller scale (labelled I in Fig.I I) if s = 1.25 and at u z 7 if s = 1.1.Nevertheless, as is the case with s = 2 and s = 1.5.as u increases, both lines go to a certain asymptotic value of Ra which lies close to the critical modal Rayleigh number of the ( 2 j + 1.2) wave mode.Before concluding this section, we will determine the stability diagram of the two Lorenr rolls in the 10-coefficient model when u = I/d and j = 1.In that case, the ( j , I) or (],I) wave mode (A) corresponds to a roll with an aspect ratio of 2 4 (the first self-excited mode) and the ( j , I) or (2, I) wave mode (B) corresponds to a roll with an aspect ratio of fi.(A) and (B) have critical modal

Fig. 12 .Fig. 13 .
Fig. 12.The lines of marginal outer stability of the two Lorem roll solutions of the 10-coelficient model (already drawn in Fig.I I for s = 2) for s = 1.1, s = 1.25, s = 1.5 and s = 2, where the smalkr scale has M aspect ratio equal to 2 4 .

Fig. 14 .
Fig. 14.The (smoothed) diameter of laboratory convection rolls (as a fraction of the critical diameter, A,, at Ra,) as a function of Ra, for u = 0.71 (---) and u = 6.7 (4.steadily increasing: +, Ra steadily decreasing), as measured by Willis et al. (1972).The boundaries were not stress-free, so A, = 2.016 h and Rae = 17081~' = 17.6 (Chandrasekhar, 1961). 36A The stability of the two possible Lorenz roll solutions to the 10-coefficient model as a function of u and Ra when u = 1/2\/2 and j = 1.Broken (solid) lines separate innerly (outerly) stable from innerly (outerly) unstable regions for the two Lorenz rolls.Thick lines correspond to the large scale Lorenz roll (A: aspect ratio = 4fi) and thin lines correspond to the small scale Loren2 roll (B: aspect ratio = 2fi).The letters A and B denote the stability of that particular mode in that particular region of parameter space.In the hatched regions there is no stable Lorenz roll.
for stress-free boundary conditions.Malkus (198 I) has estimated theoretically that Nu for stress-free boundary conditions is about 13% greater than Nu for rigid boundary conditions.At Ra = 100 (R = 9700).Nu has been measured to be about 2.2 for air (u = 0.71) and rigid boundary conditions (Chandrasekhar, I96 1 ) .Therefore, Nu would be approximately equal to 2.5 for stress-free boundary conditions at the same parameter values.According to the Lorenz model, the maximum value Nu can attain at Ra = 100 is Nu = 2.86.But if the aspect ratio of the rolls has increased to 4 4 , for example, then, according to (6.9).Nu = 2.77.On the other hand, Veronis (1966) has calculated values for Nu of about 4.8 for Ra = 100.His calculations were based on a spectral expansion in which all odd-parity wave modes were excluded.