Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core

Physics of the Earth and Planetary Interiors, 82(1994) 235—259 Elsevier Science By.

235

Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core Philippe Cardin

“,

Peter Olson

Department of Earth and Planetary Sciences, The Johns Hopkins University, Baltimore, MD 21218, USA (Received 11 August 1992; revision accepted 16 November 1993)

ABSTRACT

Results of a combined laboratory and numerical study of fully developed thermal convection in a rapidly rotating spherical shell are presented. We determine the effects of rotation and a solid inner core on high Reynolds number buoyancy-driven flow in a fluid with the geometry of the Earth’s liquid outer core. Centrifugal acceleration is used as a substitute for the spherical gravity field of the core. Experiments are made with Prandtl number & = 7, Ekman numbers E ~ 2 x 10—6 and Rayleigh numbers RaT ~ 3 x iOu. Numerical calculations of finite-amplitude convection are made using a quasi-geostrophic 5. Near the onset of convection, the flow model with the same Prandtl number, and with Rar ~ 5 x iO~and E ~ 2 x i0 motion consists of periodic columnar cells aligned parallel to the rotation axis and arrayed around the axial cylinder tangent to the inner sphere. The cells spiral and propagate in the direction of rotation, the planform in the equatorial plane resembling a drifting pinwheel. At Rayleigh numbers greater than a few times critical, the periodic pinwheel planform is replaced by chaotic, time-dependent convection consisting of numerous columnar vortices that fill the spherical shell. The vortices are driven by ribbon-shaped plumes concentrated near the equatorial plane. Reynolds stresses derived from the vortices and large-scale temperature gradients produced by the plumes generate a large-scale zonal flow that is retrograde (westward) near the inner boundary tangent cylinder, and prograde (eastward) in an equatorial band near the outer boundary. Our results suggest that convection in the Earth’s core consists of irregularly distributed columnar vortices plus a secondary zonal flow.

1. Introduction

It is generally accepted that convective motions provide the kinetic energyfor magnetic field generation in the Earth’s liquid outer core (see reviews by Stevenson (1981) and Jeanloz (1990)). The structure of convection in the core is governed by the interaction of buoyancy forces (both thermal and chemical), the spherical core—mantle and inner-core boundaries, radial gravity, rotation, and magnetic fields. Fluid mechanical interactions between these effects are nonlinear and highly complex, and, not surprisingly, very little is *

Corresponding author at: Département TAO, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris 05, France.

known about rotating spherical convection in the range of parameters relevant to the core. In this paper, we analyze the structure of high Rayleigh number convection in a rapidly rotating spherical shell with the geometry of the Earth’s core. We focus on convection at higher Rayleigh number and lower Ekman number than considered in previous studies, a regime where the convection is fully developed, highly nonlinear and characterized by a large Reynolds number, but nevertheless dominated by rotational effects. (These dimensionless numbers are defined in Table 2 below, which also lists their estimated values in the core.) We argue that high Reynolds number convection is more relevant than the regime of small Reynolds number convection,

0031-9201/94/$07.00 © 1994 — Elsevier Science B.V. All rights reserved SSDI 0031-9201(93)02919-0

236

P. Cardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

which in the past has served as the basic model of core fluid motions. Rotating spherical convection is only the first step toward a complete model of core flow. There are additional forces affecting core motions that are not included in our model. Most important among these is the Lorentz force, related to interaction of the magnetic field with electrical currents in the fluid. We note that it is possible to achieve a large Lorentz force in laboratory experiments using liquid metals such as gallium, and, with mechanical forcing, laboratory-scale models can reach moderate magnetic Reynolds numbers. However, we do not attempt to do so here. Instead, we focus on the purely mechanically driven convection. We also do not include in this study the effects of the luni-solar precession of the solid mantle which provides an additional source of core motions (Vanyo et al., 1994). In our experiments and numerical calculations, thermal convection is driven by heat flux at the inner boundary of the spherical shell. Heat transport from the solid inner core to the fluid outer core occurs in conjunction with crystallization of the inner core. Thermodynamic calculations of the inner-core solidification process indicate that this heat transport is sufficient to generate thermal convection in the liquid outer core provided the crystallization rate of the inner core exceeds a critical value, of the order of several hundred megagrams per second. Most calculations of heat loss from the core predict crystalization rates higher than this critical value, and so indicate an important role for thermal convection (Verhoogen, 1980; Cardin and Olson, 1992; Buffett et al., 1992). Temperature differences are not the only source of buoyancy in the core. Convection driven by compositional buoyancy, resulting from chemical differentiation at the inner-core boundary, is an expected consequence of inner-core growth, and is often considered the most likely cause of core motions (Loper, 1978; Gubbins et al., 1979; Loper and Roberts, 1981; Fearn, 1989; Braginskii, 1990). It had been suggested that the structure of compositional convection might be very different from the structure of thermal convection in a rotating sphere, implying that study of

thermal convection is irrelevant to the problem of flow in the core. However, our experiments with both types of flow indicate that the structure of thermo-chemical convection is very similar to purely thermal convection over the range of Ekman and Rayleigh numbers considered in this paper (Cardin and Olson, 1992).

2. Previous work on rotating spherical convection The previous work on rotating spherical convection can be classified into (1) theoretical calculations of the onset of thermal convection, (2) laboratory experiments on weakly nonlinear convection, appropriate for Rayleigh numbers very near the critical value, and (3) numerical calculations. Stability analyses of the onset of convection (Roberts, 1968; Busse, 1970a; Durney, 1970; Soward, 1977) have demonstrated that thermal convection in a rapidly rotating sphere first occurs as a circular, periodic array of columns aligned parallelto the rotation axis and extending across the sphere in the axial direction. Depending on the Prandtl number of the fluid, the array consists of propagating elliptical columns, prograde spiral columns or stationary columns attached to the inner sphere (Yano, 1992). Busse (1975) used results from stability theory of fluids with high Prandtl number to argue that the columns would be the propagating type in the Earth’s core, and would form just outside an imaginary axial cylinder enclosing the inner core called the tangent cylinder. Experiments on rotating spherical convection with radial buoyancy forces were made in the microgravity environment of Spacelab by Hart et al. (1986a). Although these experiments reached the fully developed regime, they were restricted to fairly large Ekman number values. The Earthbound technique for laboratory simulation of rotating spherical convection was developed by Busse and Carrigan (1976a,b), and since then has been used in numerous applications (Busse and Hood, 1982; Carrigan and Busse, 1983; Chamberlain and Carrigan, 1986; Cardin and Olson, 1992);

P. Cardin, P. Olson /Physics of theEarth and Planetary Interiors 82 (1994) 235—259

237

(Cordero and Busse, 1992). Rapid rotation imparts a rigidity to the fluid along cylindrical surfaces parallel to the rotation axis, and inhibits the bending of axial fluid columns. Because of this axial rigidity, the dynamically important part of the buoyancy force is the component perpendicular to the rotation axis. Buoyancy forces perpendicular to the rotation axis are achieved experimentally by rotating the fluid sufficiently fast that a large centrifugal acceleration is produced. In the rotating frame the centrifugal acceleration

Previous experimental studies have been restricted to slow rotation rates, or have focused on the weakly nonlinear regime of convection just beyond the critical Rayleigh number, where the motion is periodic or nearly periodic in the equatonal plane, and the Reynolds number is small. However, in the Earth’s core the Reynolds number is large, of the order of 108, implying that the Rayleigh number is far beyond the critical value. It is therefore necessary to determine how the flow structure changes at large Reynolds number

acts as an apparent gravity, directed away from the rotation axis. The basic temperature gradient in the experiment is reversed, relative to the temperature gradient in the core, for the centnifugal buoyancy force in the experiment to be properly directed. Heat loss at the inner sphere in the experiment corresponds to heat production at the inner-core boundary, and heat production at the outer surface in the experiment corresponds to heat loss at the core—mantle boundary. Carrigan and Busse (1983) used this method to verify some of the predictions from marginal stability theory. First, they found that motion first occurs as slowly propagating columnar rolls arrayed periodically around the inner-sphere tangent cylinder, and that the critical Rayleigh number and the wavenumber of the rolls increases with Ekman number, in qualitative agreement with the predictions from the linear theory. They also observed that the rolls tilt in the prograde direction, forming a spiral or pinwheel pattern in the equatorial plane. Carrigan and Busse did not obtain a close quantitative agreement between theory and experiment for the critical Rayleigh number, which they attributed to the presence of a thermal wind in the experiment, driven by the radial gradient in the conductive temperature distribution interacting with the centrifugal acceleration. More recent experiments sought to minimize the effects of this thermal wind. Chamberlain and Carrigan (1986) experimentally simulated internal heat generation in a sphere without an inner core. Cordero and Busse (1992) used an inverted hemisphere geometry with a rotation rate that provided a nearly radial buoyancy force. Both of these studies obtained slowly propagating periodic rolls.

and the convection becomes fully developed, but still dominated by the Coriolis acceleration. In particular, it is important to establish whether the periodicity or the columnar structure is replaced first as the Rayleigh number is increased. Numerical calculations of three-dimensional thermal convection at Ekman numbers of i03 or greater indicate that the columnar structure is either absent altogether or disappears along with the periodicity as the Rayleigh number increased beyond the critical value (Glatzmaier, 1985a,b; Hart et a!., 1986b; Sun et al., 1993). A somewhat different picture emerges from calculations made at lower Ekman numbers. Using a truncated Galerkin method, Zhang and Busse (1987) and Zhang (1991, 1992) have obtained periodic columnar convection patterns for Rayleigh numbers up to approximately twice the critical value, for Ekman numbers as low as i05. Using a fully three-dimensional model, Glatzmaier and Olson (1993) found chaotic but essentially columnar convection at an Ekman number of iO~ and a Rayleigh number of iO~,which is approximately 50 times critical. As a test of the validity of the experimental technique, Glatzmaier and Olson also compared the structure of convection driven by radial buoyancy with convection driven by centrifugal buoyancy, at the same Rayleigh number. Although the two flows are not identical, they have many characteristics in common, particularly outside the tangent cylinder of the inner core. Both exhibit columnar flow outside of the tangent cylinder, driven by thermal plumes lying near the equatorial plane. Glatzmaier and Olson found that superposition of numerous plumes produces an equatorial torus of anomalous ternperature fluid. They also obtained for each case

238

P. Cardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

generally retrograde zonal flow, consisting of a thermal wind (with axial shear) driven by the equatorial torus and a columnar geostrophic wind driven by Reynolds stresses derived from the smaller-scale convection. Although three-dimensional numerical models offer very powerful tools for investigation of notating spherical convection, they are currently limited to rather modest Rayleigh and Ekman number combinations, and it is unlikely that this situation will change very much in the near fitture. The structural complexity of convection at low Ekman number and highly supercnitical Rayleigh numbers in a rotating fluid sphere cannot be adequately resolved with the computer power currently available. However, as we show here, laboratory experiments can easily reach more extreme combinations of Rayleigh number and Ekman number than are at present accessible to calculation, and thereby more closely represent the flow regime in the core. 3. Columnar convection in the core Is it reasonable to expect that convection in the core is columnar? Columnar flow results from a geostrophic force balance, in which the Ekman number, Rossby number and Elsasser number are small compared with one. In addition, axial columns (columns parallel to the rotation axis) cannot form if the fluid is thermally or chemically stratified. Finally, the ratio of buoyancy to Coniolis forces in the fluid must be small, otherwise the columns will be destroyed by the shearing motion in the large-scale thermal wind. As indicated in Table 2 below, the Rossby number is small and the Ekman number is extremely small in the core, so inertia is secondary to the Coniolis acceleration, and viscous forces are probably negligible except in thin boundary layers. The Elsasser number, the ratio of Lorentz force to Coniolis acceleration, depends on the strength of the magnetic field. If the strength of the toroidal field inside the core is comparable with the strength of the poloidal field, then the Elsasser number is small, less than 0.1. Alternatively, a strong toroidal field, say 0.01 T or larger,

implies an Elsasser number greater than one. Because the toroidal field strength is essentially unconstrained, we can only estimate a lower bound for the Elsasser number of about 0.01. Columnar convection is favored when the Elsasser number is less than one and is therefore expected in the core provided the total field is not much stronger than the poloidal component. The situation regarding stable stratification is less clear. There has been speculation about the existence of stable chemical stratification, and the uncertainty in seismic models of the core permits a weak stratification, but there is no strong evidence in favor of it. It is easier to judge whether thermal effects are capable of destroying columnar flow, because we can estimate the ratio of the thermal buoyancy force in the fluid to the Coniolis acceleration. The magnitude of thermal buoyancy forces is agT’, where a is the thermal expansion coefficient, g is gravity and T’ is the magnitude of lateral temperature fluctuations in the fluid. V can be estimated indirectly, using the average convective heat flux q. Temperature fluctuations are related to the convective heat flux according to q pC~uT’, where p. C,, and u are fluid density, specific heat and velocity, respectively. In order of magnitude, the ratio of buoyancy force to Coriolis acceleration for thermal convection in the core is then agq C flu2

(1) p where fi is the angular velocity of rotation. If we use q = 0.01 Wm2 for the convective flux denived for inner-core growth by Cardin and Olson (1992), a typical fluid velocity of u = 5X iO~m 5 (Bloxham et al., 1989), assume a iO~K1, and use p = io~kg m3, g 10 m ~2, c~, = 620 J kg’ K’, and fi = 7.29 x iO~s1, then A 10-2 in the core. A roughly similar value is obtained for compositionally driven core convection. The expression for compositional buoyancy flux given by Cardin and Olson (1992) results in A 0.05. We conclude that it is likely that buoyancy forces are small compared with the Coriolis acceleration in the outer-core fluid. Except for the highly uncertain magnitude of the Lonentz A

P. Cardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

force, and the speculative possibility of chemical stratification, the conditions necessary for columnan convection in the core seem to be satisfied. There is evidence in support of columnar convection in the core from the structure of the geomagnetic field on the core—mantle boundary. Maps of the magnetic field on the core—mantle boundary derivedby Gubbins and Bloxham (1987) show that most of the dipole field is due to flux concentrated in a couple of high-latitude patches. Flux patches in the northern hemisphere have antisymmetnic counterparts (with reversed flux) at the same longitude and the same latitude in the southern hemisphere. The patches of concentrated field are located where the inner-core tangent cylinder intersects the cone—mantle boundary. It has been suggested that the flux patches are caused by convection rolls positioned around the inner-cone tangent cylinder. Models of flow at the top of the core have been derived by inverting geomagnetic secular variation data using the frozen flux assumption. Some of these models include structures which are suggestive of columnar convection (see reviews by Bloxham et al. (1989) and Bloxham and Roberts (1991)). For example, Hulot et al. (1990)

239

have shown that the assumption of tangentially geostrophic flow, together with the frozen-flux approximation, yields a pattern of circulation in the outermost part of the core which is nearly mirror-symmetric with respect to the equatorial plane. Hulot et al. pointed out that equatorial mirror-symmetry in the circulation is consistent with convection columns aligned parallel to the rotation axis. It is worth emphasizing that geomagnetically derived core flow patterns are not periodic in longitude. 4. Experimental methods Figure 1 is a schematic illustration of the experimental apparatus. A plastic outer sphere of 30 cm diameter and a copper inner sphere of 10 cm diameter are mounted concentrically on a stainless steel shaft of 2.2 cm diameter. A flywheel is also mounted on the shaft, to damp high-frequency fluctuations in the notation rate. The shaft is bored along its axis, providing a conduit for the cooling fluid used to control the temperature of the copper sphere. Steady, ventical counter-clockwise rotation of the shaft is

IUNf

TOP ~

SIDE VIEW

cold bath

gD

Fig. 1. Schematic diagrams of the apparatus in plan and section views, showing the experimental geometry, temperature boundary conditions, and viewing orientations.

244)

P. GardEn, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

maintained by a d.c. motor with a specially designed feedback control system. The feedback system operates by comparing the frequency of an analog encoder attached to the motor axle with the frequency of an independent harmonic oscillator, and generating an error signal propontional to the difference between the two frequencies. The current into the armature of the motor is passed through a transistor and chopped by an amount proportional to the error signal. This control device, when used with an adjustable mechanical brake, allows us to maintain the notation rate constant to one pant in iO~during an expeniment. Thermal convection is produced by imposing a fixed temperature difference between the inner and outer spheres. The temperature at the inner boundary is controlled by circulating cooling fluid at uniform temperature along the axis of the shaft and through the copper sphere. Rotaryjoints are used to connect the spinning shaft to the cooling fluid reservoir. Temperature variations on the inner-sphere boundary are less than 0.05 K. The outer sphere is maintained at approximately room temperature, using the air currents genenated by the spinning apparatus. Water is used as the working fluid in the spherical shell. The pentinent thermodynamic and transport properties of water, together with the physical parameters of the experimental device, are listed in Table 1.

To begin each experiment, the apparatus is notated in the clockwise direction at a steady angularvelocity for at least 30 mm, for the working fluid to reach solid body rotation. Experiments were made with rotation rates ranging from 250 to 400 rev min At 400 rev min 1, for example, the centrifugal acceleration is 25 g at the outer equator, and 8 g at the inner equator. After solid body notation is reached, the cinculating bath is switched on and the inner-sphere temperature is maintained at a preselected value, either above or below the initial temperature of the working fluid. In addition, the room temperatune is carefully monitored so that the temperatune difference between the inner and outer boundary is known. Thermocouples are mounted in the notating shaft at entrance and exit points of the inner sphere, providing a measurement of the temperature rise in the cooling fluid as it passes through the inner sphere. The temperature rise is converted to average heat flux, using the heat capacity and flow rate of the cooling fluid. This measurement allows us to determine, for each experiment, the Rayleigh number based on the time-averaged heat flux at the inner sphere. We use the Rayleigh number based on heat flux for extrapolation of our results to the core. Two methods of flow visualization are used. In the first method, a high-intensity, high-speed flash is collimated into a vertical sheet of parallel light ~.

TABLE 1 Physical parameters Parameter R

0 R1 D

gD

a v IC

p k

Q

Definition

Experiment

Core

Outer radius Inner radius Shell depth Rotation velocity Temperature difference Gravity at outer radius Gravity at inner radius Gravity at R = D Thermalheat Specific expansivity Kinematicviscosity Thermal diffusivity Meandensity Thermal conductivity Convective heat flux

0.15 m 0.05 m 0.1 m 25—40 rad s 0—10 K 94—240 m ~2 3 1—80 m ~_2 63—1604mK1 ~_2 2X10 4.l8XlO3Jkg’ K1 1x106 m2 s 1.4 x i0~ m2 s~ lXlO3kgm3 5.9x10’ Js1 m1 K1 —(5—50) W

3.48 x 106 1.22x 106 2.26X 106 7.29x iO~ ? 10.7 4.4 7~1X105 620 ~4X106 4 x 10_6 1X104 30 (0.2—2)x 1012

P. Cardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

241

TABLE 2 Dimensionless parameters Parameter

Definition

Inner/outer radius ratio

R

Prandtl number, F~Ekman number, Rayleigh number,ERaT Rayleigh number, RaQ Reynolds number, Re Rossbynumber,ReE Hartmann number Elsasser number

1 /R0 v/K 2 v/flD ag 3/Kv ag~QD/kKv 0t~TD UD/v U/liD B2D2o-/pv o~B2/flp

of 0.5 cm thickness. The sheet is projected through the sphere along one of the two ray paths labeled A and B in Fig. 1. The distribution of shear planes in the beam is revealed by the variations in intensity of light reflected from rheoscopic flakes seeded in the fluid. This is a well-established tecluiique for convection experiments, and is particularly suitable for visualizing the flow structure in planes parallel to the rotation axis of the sphere. In discussing our results, we refer to this as side view visualization. The second method is used to visualize the planform in the equatorial plane. Fluorescein dye is transferred to the rotating sphere through a rotary joint and injected into the convecting fluid through tubes placed near the inner sphere, as shown in the section view in Fig. 1. The dye is illuminated from the side, using a strobe light source filtered to reduce

Experiments

Outer core

0.33 7io~—io~ ~ 3 x i0~ ~ 7x 10° 102_103 “iO~ 0 0

0.34 0.1—10 10_15_10_13 ? 10~_1029 107_108 ~iO~ 1012 1013 0.1—10

are the Pnandtl number, the ratio of fluid kinematic viscosity to thermal diffusivity: F~=

(2)

—

K

the Ekman number, the ratio of the viscous force to Coniolis acceleration: p

E=

(3)

and the Rayleigh number, a dimensionless measure of the thermal buoyancy: ag ATD3 Ra~= D (4) Ky

Here, I.’, K, and a are the kinematic viscosity, thermal diffusivity and thermal expansion coefficient, respectively, and D = R 0 R1 is the shell thickness (the difference between the outer radius R0 and inner radius R1), ~T = T~, ~‘1 is the difference between the outer boundary 2D andisinner the boundary temperatures, gD = L1D from the centrifugal acceleration atand a distance notation axis. It should be noted that R = D is the midpoint of the spherical shell in the equatorial —

visible wavelengths and transmit uv wavelengths, Filtering is necessary to reduce reflections of visible light from the spherical surfaces, which tend to isobscure the fluorescent dye pattern then photographed fromimage. above,The in what we refer to as the top view, The three important experimental parameters

—

TABLE 3 Summary of experimental results Figure

E

RaT

RaQ

Planform

2(a) 2(b) 3

2.5x106 2.5 x 10_6 2.5 x 10—6 3

5 6x107 6x i0 4x iO° 4 x 10°

No columns Pinwheel

4

~106(0.1Ra,~,, Ra~,, 1) 7 X iO~(50 Ra~, 1(decaying) 8 (50 Ra~,,1) 7 x i0 1)

x 10~

Chaotic vortices + zonal flow Chaotic vortices + zonal flow

242

P. Gardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

plane. Using total heat flux Q in place of ~T, an alternative definition of the Rayleigh number is RaQ=

ag~QD kKv

=NuRaT

(5)

where Nu = Q/kDi~T is the Nusselt number of the convection and k is thermal conductivity of the fluid. Table 2 gives the range of dimensionless parameters we have investigated experimentally, plus estimates of these parameters in the Earth’s core.

—

5. Experimental results Table 3 lists the parameter values of experiments illustrated in Figs. 2—4, along with a qualitative description of the resulting convection patterns 5 1 Convection near the critical Rayleigh number

Our main interest in this paper is the structure of fully developed convection, at highly super critical Rayleigh number Here we briefly de scribe the pattern motion near the onset of convection, to demonstrate compatibility with previously published experiments. Figure 2(b) shows the structure of convection at E = 2.5 x 10—6 and RaT = 1.2 x iO~which is in the neighborhood of the critical Rayleigh number. We have used the following procedure to obtain this pattern. Thermal convection was initiated by maintaining the temperature of the inner sphere several degrees below room temperature for a period of several minutes. Once the convection developed, the inner-sphere temperature was increased slowly, gradually approaching the room temperature. The amplitude of the convection decayed as the driving temperature gradient diminished. The image in Fig. 2(b) represents the flow structure in the final stage of decay, when the temperature difference ~T is nearly zero. The flow shown in Fig. 2(b) consists of a nearly periodic array of axial, columnar vortices arrayed in a circle around the tangent cylinder surrounding the inner sphere. The absence of structure

Fig. 2. Experimental results on the structure of thermal convection near the critical Rayleigh number, at E = 2.5 x 10_6. (a) Side view of subcritical flow (no columns), using beam B shown in Fig. 1; (b) side view of pinwheel convection planform near the critical Rayleigh number using beam B.

elsewhere in the image indicates that convection is limited to the neighborhood of the tangent cylinder. It should be noted that the image is not completely symmetrical with respect to the rotation axis. The apparent width of the columns increases from left to right across the inner sphere. We infer from this variation that each column is elongate in the radial direction and also tilted in the direction of notation, implying that the necklace of columns forms a prograde pinwheel planfonm when viewed from above. This is essentially the pattern described previously by

P. Gardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

Carnigan and Busse (1983) and Chamberlain and Cannigan (1986). We also observe that the pattern is not stationary. Photographs of the same sector of the sphere taken several tens of seconds apart indicate a very slow prognade drift on propagation of the pinwheel pattern. The direction and speed of drift are qualitatively consistent with the propagation rate measured experimentally by Condeno and Busse (1992) in a hemispherical shell with a proportionally larger inner-core size than ours. On the basis the number of columns appearing in this image, we estimate there is a total of about 40 columns arrayed around the tangent cylinder, corresponding to an azimuthal wavenumben of 20. The marginal stability calculations presented in the next section predict an azimuthal wavenumben of about 20 for the preferred planform at this Ekman number. The number of columns varies inversely with notation rate. As the rotation rate is increased, and the Ekman number becomes smaller, we observe that the azimuthal wavelength of the tilted cells becomes shorter. However, the reduction in cell size is a rather weak function of rotation rate, and oven the restricted range of rotation speeds we have investigated the variation in column diameter is small. As the previous experimental work by Busse, Carnigan and coworkers has already determined the critical Rayleigh number for the onset of convection, we have not made a full exploration of that condition in our apparatus. However, to check our interpretation that the first appearance of columns corresponds to the onset of convection, we have determined the existence of a minimum Rayleigh number for column formation. Figure 2(a) shows a side view of the fluid at a subcritical Rayleigh number, approximately 106. Columnar flow is cleanly absent at this Rayleigh number. 5.2. Observations at supercritical Rayleigh numbers

At larger Rayleigh numbers, the convection is fully developed and the columnar motion fills the sphere. In this regime, the spatially periodic planforms found near the critical Rayleigh number

243

evolve into a complex pattern consisting of irregulanly distributed vortices driven by thermal plumes. The pattern of motion is strongly time dependent in this regime. A lange-scale, zonal motion is also present as a secondary flow. These structural changes mark the transition from low Reynolds number convection dominated by viscous forces, to high Reynolds number convection in which inertia of the fluid is important. The flow structure in the fully developed columnar regime is shown in Figs. 3(a)—3(d), from an experiment with E = 2.5 x 10—6 and RaT = 7 X 108 50 Ra~~1, where Ra~~1 denotes the critical Rayleigh number. Figures 3(a) and 3(b) are side views, using beam orientations A and B, respectively, as shown in Fig. 1. These images demonstrate that columnar flow occupies the entine spherical shell. There are several columnar vortices in each radial section, rather than the single circle of vortices observed near the critical Rayleigh number. Columns are present inside the tangent cylinder, although the flow associated with them is less vigorous than the flow outside the tangent cylinder. Whereas there exists a charactenistic dimension to the columns, it is clear the columns have a variety of cross-section dimensions, and their spatial distribution is not penodic. Also, we observe that the column distnibution in these images changes within a period of a few seconds, indicating that the flow structure is variable on that time-scale. The planform of motion in the equatorial plane in the fully developed columnar regime is shown in Figs. 3(c) and 3(d). These figures show dye patterns as seen from above, in top view. Figure 3(c) was taken about 2 mm after a band of dye was released from a probe located near the rotation shaft and just above the inner sphere. Figure 3(d) is the same view of the same experiment, taken about 10 mm later than 3(c). The dye is transported from the probe toward the equator in the boundary layer surrounding the inner sphere. The cold and dense fluid in the boundary layer separates from the inner sphere as it nears the equator, forming ribbon-shaped plumes centered in the equatorial plane. The structural evolution of these plumes is revealed by the dye streaks in Fig. 3(c). As the fluid in each plume moves radi-

244

P. Cardm, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

~C~Cs

OCs

a~ a,~ Cs

8: 2 oCs ~~ao

_______

£~‘ Cs CS

.

~ n.E aC) ~Cs

V

n..

C)

~

~

.~

.~

~

2 .s ~ ~ C)

C)

2 .E ~ ~)

.n

.a .E “C

P. Gardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

ally outward from the inner sphere, it tends to spinal in the retrograde sense. This occurs because the vortex on the retrograde side of the plume is enhanced by the Coriolis acceleration, whereas that on the prograde side is diminished by the Coniolis acceleration. The spiraling action

245

limits the radial distance each cold plume travels, and eventually all the cold plumes originating from the inner sphere are turned into retrograde vortices. These cold vortices are arrayed in an irregular band aborit an imaginary axial cylinder with a diameter a few centimeters larger than the

Fig. 4. Experimental results on the structure of thermal convection for Pr—•• 7, RaT = SORa~~ and E = 3~10~.Top view showing convection pattern at different times in the equatorial plane, using fluorescein dye released above inner core. Photographs were taken at the following times (min:s) after the start of dye release: (a) 5:59; (b) 10:50; (c) 11:12; (d) 13:26; (e) 34:00; (0 36:00.

246

P. GardEn, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

inner sphere, as can be seen in Fig. 3(c). After a few tens of seconds, a secondary set of cold plumes develops from the first band of vortices, The secondary plumes transport cold fluid funther outward, until they too are turned into netnograde vortices by the Coriolis acceleration. The result of this turning action is a chaotic pattern of columnar convection, consisting of several uneven vortex bands, as shown in Fig. 3(d). The bands are irregular, as is the distribution of vortices within each band. Consequently, the flow structure is highly time dependent. The ensemble of ribbon-shaped plumes defines a lange-scale torus of cold fluid centered on the equatorial plane. Figures 4(a)—4(f) are from an experiment with E = 3 X 10—6 and Ra~ 50 Racrit. The images were made by a steady release of dye from the probe above the inner sphere, begun about 5 mm after the convection was initiated. This sequence illustrates the radial transport processes, including plume formation and the development of the vortex band structure (particularly in Figs. 4(c), 4(d) and 4(f)). The single dye filament visible in Fig. 4(a) shows first-, second-, and third-generation cold plumes and the vortices derived from them. We note that only retrograde, cold vortices are marked by the dye in Fig. 4. There are also prograde-notating, warm vortices which originate as radially inward moving warm plumes. The positions of warm plumes and warm vortices are indicated by holes in the dye pattern in Figs. 4(a)—4(f). Using successive photographs of the dye pattern taken a few seconds apart, we have estimated the maximum radial velocity of the cold plumes. The velocity of a typical cold plume is approximately 0.5 cm s~ at this Rayleigh number. However, plume formation is highly intermittent. The average radial transport velocity, as indicated by the rate of advance of the dye front, is much less, perhaps 0.05 cm s~.The swirling velocity within the vortices is probably the most energetic component of the motion, with sustamed velocities typically 0.5—1 cm s~. In addition to columnar plumes and vortices, a lange-scale zonal flow is present when the convection is fully developed. The general pattern of zonal motion can be inferred from the distribution of dye in Figs. 4(a)—4(f). Dye is injected

above the inner sphere, near the left edge of the images. As dye flows toward the equator in the cold thermal boundary layer surrounding the inner sphere, it is deflected in the retrograde dinection. The point at which dye reaches the equator and separates from the inner sphere is typically 20—30° retrograde of the injection point. The retrograde drift of the dye continues after it becomes incorporated into the cold plumes and vortices. From the orientation of the dye streak in Fig. 4(b), we infer that the angular speed of retrograde zonal flow is greatest near the inner sphere, and decreases with radius. Comparison of photo pairs taken sequentially in time indicates that the zonal flow angularvelocity is typically 10° min’, which corresponds to a linear velocity of about 0.02 cm s~ near the inner sphere. This is smaller than the swirling velocity in the vortices, which we have estimated to be about 0.5 cm s~ at this Rayleigh number. The speed of the retrograde zonal velocity decreases with radius, and the drift direction becomes prograde near the outer-sphere equator. Evidence for this is seen in Fig. 4(e), where many of the cold plumes near the outer equator are strongly sheared and tilted in the prograde direction. In addition, we have observed prograde drift of dye injected through a probe located at the outer-sphere equator. The zone of prograde motion appears to be limited to a narrow band, approximately 30°wide and 1 cm deep, and symmetric about the equator. In every experiment we observe more retrograde motion than prograde motion. The prograde motion, when present, is confined to the region near the outer-sphere equator, whereas retrograde motion is present throughout the nest of the fluid. Although the zonal flow is ubiquitous when the convection is fully developed, it is clearly a secondary motion compared with the columnar convection. The total kinetic energy in the vortices exceeds the kinetic energy of the zonal flow in all of our experiments, by a factor of at least 10 and typically 100. There are two driving forces for the zonal flow, one kinetic and one thermal. The kinetic driving force is Reynolds stresses derived from the prograde tilt of the convection columns. Zonal

P. GardEn, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

flow driven by this mechanism is cylindrical, and reflects the columnar structure of the convection. As shown previously by Busse (1970b) and Busse and Hood (1982), and also by our numerical calculations in the next section, this component of the flow is retrograde near the inner tangent cylinder and prognade near the outer equator, similar to the pattern we observe. The thermal driving force is the baroclinicity of the thermal torus, and results from superposition of numenous cold plumes in the equatorial region. The torus generates in a retrograde thermal wind with shear in the z-dinection. Previous studies came to different conclusions regarding the relative importance of these two forces. Canrigan and Busse (1983) inferred that the thermal wind component drives the zonal flow, whereas Glatzmaier and Olson (1993) found in their numerical calculations that the two forces are comparable. To establish which of the two driving forces is dominant in our experiments, we made some observations to determine the relative magnitude of the columnar transport vs. the axial shear in the zonal flow. We injected a steady stream of dye onto the inner sphere and observed the motion of the dye in the side view. The combination of column stretching and Ekman pumping gradually draws the dye into filaments parallel to the rotation axis. We measured the bending of these filaments as they were transported by the zonal flow. In every case, the bending is observed to be very small, indicating that the columnar (z-independent) component of the zonal flow is larger than the thermal wind (z-dependent) pant. We also made measurements of the time-average heat loss from the inner sphere Q, as a function of the temperature difference between the inner and outer sphere, ~T. These are obtamed by maintaining L~Tconstant for a lengthy period of time, typically 4 or 5 h. It is difficult to maintain a uniform temperature on the outer boundary in our apparatus because the outer sphere is made of Plexiglas, which is a poor heat conductor. Because the temperature varies on the outer sphere, there is some uncertainty in the actual value of i~Tacross the fluid. We estimate that the absolute value of RaT is uncertain in our experiments by as much as ±10%,although the

a8

50

ROT 5 xl ~8

0

°

30 20

6 xl

4x

~-

~

1 xi 8x1

40

~

247

10

~.o

a

9

~..•

2 xl 0

•

10

10

~ (°c) .

.

Fig. 5. Experimental relationship between outer—inner boundary temperature difference /~Tand total inner-sphere heat flux Q for Pr =7 and E =3 x 106. Opposing axes show the equivalent RaT—RaQ relationship. The dashed line is RaQ = 6RaT.

relative uncertainty between the various expeniments is less. Most of the uncertainty can be eliminated by adopting the Rayleigh number based on total heat flow, RaQ, defined in Eq. (5). Figure 5 shows the correlation between i~Tand Q, and the corresponding relationship between RaT and RaQ, for a number of experiments. The linear relationship Ra 6 Ra (6 1’

is a good approximation to the data, indicating an approximately constant Nusselt number Nu 6. In the Introduction, we argued that thermal convection in the outer cone is likely to be columnan. Part of that argument was based on the parameter A, the ratio of thermal buoyancy forces to the Coriolis acceleration, being much less than one. We now demonstrate that the equivalent parameter in these experiments is also small compared with one, and indeed is comparable with the value of A in the cone. As g = fl2R and q = Q/4irR~in our experiments, (1) can be written aflQ A’ = 2 (7) 8rrpR~C~u

248

P. Gardin, P Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

From the swirling velocity in the cold vortices, u 0.5 cm s~.Using the other parameters listed in Table 1, we find A’ ~ 10—2, As anticipated, thermal buoyancy in the convecting fluid is small compared with the Coniolis acceleration. Three

Here u is fluid velocity, P is reduced pressure, T is temperature, I is a unit axial vector (parallel to the rotation axis fi), and g is the gravity vector. These equations have been made dimensionless using the spherical annulus gap D = R0 R1 as

consequences follow from the fact that A’ is small.atFirst, the numbers convection columnar, even Rayleigh fan remains beyond the critical

the length the vorticity diffusion t~= 2/v as thescale, time-scale, and PrAT as thetime temperD ature scale, where i~T= T 1 is the temperature difference across the annulus. The parameters Pr, Ra~and E are defmed by (2)—(4). In this model, we include only a centrifugal gravity, corresponding to a constant body force directed perpendiculan to the rotation axis, so that

value. Second, it insures that the equatorial component of gravity is the dynamically important component, which permits us to model core convection using the centrifugal acceleration. Finally, because the value of the parameter A’ in the experiments is comparable with its estimated value in the core, we conclude that our experiments correctly model the balance between Conohs and buoyancy forces in the core.

6. Numerical model of columnar convection

—

—

g F (11) where P is the unit radius vector in the (r, 0, z) cylindrical coordinate system. The vorticity equation, obtained by taking the curl of (8), has the following r, 0, and z components: —

8u

The columnar structure observed in the expeniments suggests that rapidly notating spherical convection can be closely approximated by flow with uniform vonticity in the z-dinection. Models

+

-~

(u

.

V)Wr

(a

—

V)Ur

-

—

aw 0

of this sort have already been used for convection in a notating cylindrical annulus (On and Busse, 1987; Brummell and Hart, 1993). Here, we derive the equations of motion for a spherical shell with centrifugal buoyancy, subject to the condition that the ratio of buoyancy to Coniolis force is small. We show that the equations of motion reduce to columnar flow when both E and A are small and when the mass flux across the equator is zero. In a system of coordinates notating with the spherical boundaries, equations for conservation of momentum, mass and heat in the fluid are, using the Boussinesq approximation

2u

—E’VP— Ra~gT+V V u= 0

(8)

=

.

8T -~--+(uV)T=Pr~V

2T

8u9 + (is

—

‘

V) to0

(9) (10)

(to

—

.

V) u9

13

a +

(u V)w ‘

at =

2E

toO +

a —~-~

—

R

~2

—

—

z

Z

~

+

—i

to

Ra aT —i —

r

az

a —

( )

2E~—i

az (14)

80

We expand the fluid velocity and vorticity in powers of the parameter A defined in (1): u=u0+Au1+... w=w0+Aw1+...

—+(u’V)u+2E~IXu

2C0r

2E-~--~ = V (12)

(15)

(16) We substitute (15) and (16) into (12)—(14) and collect terms with like powers of A and E. The lowest-order equations are simply au

-~--—=0 0

(17)

P. Cardin, P Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259

The restriction of zero transport across the equaton then implies u~0= 0, and consequently (17)

where

yields

r~=

u0=u~(r,0)

(18)

w~=Iw(r, 0)

(19)

where the subscript e denotes the equatorial yelocity, the projection of the velocity vector in planes parallel to the equator. We defme columnar flow to be fluid motion satisfying the conditions (18) and (19). At the next order in A and E, Eqs. (12)(14) are (20) —

2A

Equation (22) is the predictive equation for the axial (columnar) vorticity and (21) is the thermal wind equation for the deviations from columnar flow. The thermal wind u 91 is higher order in A than the columnar velocity Ue. Following the method used by On and Busse (1987), we average Eqs. (10) and (22) along columns in the axial direction. This results in an column-average temperature,

I

T dz (23) 2LJL where 2L(r) is the height of the tangent cylinder at radius r. As to and Ue do not depend on z, the ~

terms in (10) and (22) involving only these vanables are not affected by averaging. However, terms involving u~ are affected, as this component of velocity depends on z. Averaging the third term in (22) yields

~I

L -

u~1(z=L)

=

(26)

=L)

~~UrO(Z

The column-averaged vonticity and heat equations then become

ato 2~(r) T•+(Ue’Ve)to_EL~•~~Ur=

~2w V

~

+

RaT r aT 80

au~1

A

EL~u z1]~L 2 ~~777U,0

—

(27)

a? (22)

A

is the slope of the outer spherical boundary. We have used the condition that the normal component of the velocity vanishes at r = r0, leading to the following relationships on the outer boundary:

and

2A au~

8w

=

(25)

dr

(21)

8z

~+(ue.Ve)w_~~!=Ve2w+RaT~•j

T

dL —

—ERaT 8T

—

az

249

where 82

V2e

=

8r2

—

+

1 & r8r

— —

1

82

(29)

+ —i r

In deriving (27) and (28) we have ignored axial diffusion of heat and Ekman boundary layer transport at the outer sphere. We have also assumed that there is zero heat transfer across the equator. It is convenient to introduce a streamfunction description of the columnar velocity field. As Ve’ue=O,wehave Ue = Ve x (I~r) (30) so that = V~

~

(31)

—

Substituting (30) and (31) into (27) and (28), we obtain (1

8w -~-

(24)

(28)

~_+(ue~Ve)?=Pr_iVe2?

+

~-~:-

2’q(r) 1

OçIi 1 8w \

8,~v8w -

Ra8T

=V~w+—~.~-

-

açli

EL(r) r80

(32)

250

P GardEn, P Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259

and

a?

7. Linear stability analysis +

(1 a~~a?a~fr1a? ~

—

=

Pr

1VC2?

(33)

The height of the tangent cylinder at the outer sphere is given by 2L

=

2(r~

—

(34)

r2)i~’2

dL

di

rotating spherical shell, using linearized versions of Eqs. (31)—(38). The steady-state conduction solution to (33), satisfying the boundary conditions (38) is 1 ln(r/r 1)

?~= Pr ln( r0/r1)

and the slope of the outer boundary is

=

vection outside the tangent cylinder in a rapidly We have computed the onset of columnar con-

(39)

r =

—

(r~

—

)

(35)

~2 1/2

We assume infinitesimal perturbations, with the

form We adopt the following boundary conditions on temperature and fluid velocity. At the tangent cylinder of the inner sphere, we assume zero shear stress and constant average temperature,?. At the outer-sphere equator we assume zero slip and constant average temperature. Both boundaries are assumed impermeable. Thus

~( r~)= ‘I’( r1)

=

0

(36) 2 a~

w(r0)= —-~—~-;w(r1)= ~

(37)

(v,, 7’)

[‘P(r), 0(r)] exp(iao+o’t)

=

where a is angular wavenumben, and (41) includes the growth nate °~rand the frequency o~, of the disturbance. Substituting (39)-.(41) into (31)—(33), and neglecting all the terms involving products of disturbance quantities yields the following linearized perturbation equations: 2+1 2 (2a 2

and

r

(38) Equations (31)—(33), along with the geometrical coefficients (34) and (35), and the boundary conditions (36)—(38) are the equations of motion for columnar convection in a spherical shell. These equations apply outside the tangent cylinder. ~ side the tangent cylinder, where fluid columns touch the inner sphere, the definitions of tangent cylinder height 2L must be modified, and the different values of the boundary slope ~ at the top and bottom of each column must be included in the model equations. As the volume of fluid inside the tangent cylinder is small compared with the total volume of the spherical shell, because convection is more difficult to excite near the poles of a notating sphere, and because the centrifuge model is a poor approximation to a self-gravitating sphere in this region, we do not compute the flow inside of the tangent cylinder.

(40)

r

‘2a2+1

u\

I a4

r/ o~a2

H,

4a2

—

2w

)

+

r~

+ r2

+

E(r~ r2) —

íaRaT 1

=

(42)

r

and ®~+ 1

—~

—

/~

r

+

—no a2\

=

__________ r2(ln r ia 0

—

ln r~)

r

(43) Primes in these equations denote derivatives with respect to radius. The boundary conditions neduce to ~I’(r~)= ‘P(r0)

=

0(r1)

=

0(r0)

=

0

(44)

P Gardin, P. Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259

and

i””’~

1’ 1 4 /

a’!’

a2~I, 18’!’ (45) 2 r 8r r=r =0 8r r=r0 8r Equations (42)—(45) were solved numerically with specified values of E and Pr to obtain the —

251

=~———

critical Rayleigh number ~ the critical wavenumber ~ and the critical frequency ~ The numerical technique we use has been described fully by Candin and Nataf (1991). Briefly, the Rayleigh number and the wave frequency at the onset of convection for each wavenumben a were obtained at specified Prandtl and Ekman numbers by setting the disturbance growth rate 0~r to zero, and solving the system of Eqs. (42) and (43) with boundary conditions (44) and (45) as a sixth-order, complex eigenvalue problem. The critical Rayleigh number ~ corresponds to the minimum RaT value, ~ is the corresponding critical frequency, and ~ is the critical angular wavenumber. For Ekman numbers greater than iO—~,we used a uniform 100-point grid extending oven the radial interval (r-,I r At smaller Ekman numbers, where the flow is confined to the neighborhood of the tangent cylinden, we found it necessary to restrict the domain of calculation to a smaller interval, (r ).

1,

rm~),

where r
a

0crit crit 3

5 6 10 13 21 26 35

4.2 13 19 56 88 254 400 1152

___________________________________________

-

N,

/

experimerit5

106

~ 1~

-

no

.~

I

columns

1 o~

1 0~ 1 0_2 .

.

‘‘

1

0~ Ekman

Number

1 0_6

.

Fig. 6. Regime diagram for thermal convection in a rapidly rotating spherical shell for a Pr =7 fluid. The solid curve is the critical Rayleigh number for convection, determined using the columnar flow model. Range of experimental results is given by the rectangular shaded region. Parameter combinations of several calculations are shown with asterisks. Dashed curve is the stability curve from Busse’s (1970a) narrow cylinder model, modified according to the analysis in Section 6. The regime of columnar convection lies above the stability curve and includes the region where the parameter A from Eq. (1) is small. The regime where A is large and convection is not columnar is indicated schematically as three-dimenstona.

ous curve in Fig. 6, and the corresponding variations in the (nondimensional) critical wavenumber and critical frequency with Ekman number are shown in Fig. 7. The small Ekman number the following approximate results have been fitted torelationships: power laws, yielding ~

4E615

(46)

~

0.615E°3

(47)

—i 3

—0.13E / (48) Also, we find that the radial scale of the flow 8 .

.

the distance from the inner-sphere tangent cylin-

252

P Cardin, P Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259 ii’

are too weak to overcome the resistance to column stretching, and convection is suppressed. We can use these concepts to relate the results of our stability calculations to Busse’s (1970a) analytical theory, derived for a thin cylindrical shell with uniformly sloping upper and lower boundaries. We identify the thickness of Busse’s shell with 6, the nondimensional radial lengthscale in our solution. Then the parameters for the thin cylindrical shell become local parameters in the spherical geometry. Denoting local panametens with subscripts 6, Busse’s (1970a) expression

‘i’

1 o~

102

f req u enc~,_/

-

• ~

for the critical Rayleigh number can be written

4—~”~venumber

-

10

~

-

10

_______________________________________________________ 1 ~_2 1 o~Ekman 1 o~Number 1 O~ 1 0_6 Fig. 7. Variation of dimensionless azimuthal wavenumber, propagation frequency, and cell radial length at the critical Rayleigh number as a function of Ekman number for a Pr =7 fluid calculated using the columnar flow model.

Ra8~~~ = ~ ~ (1 Pr) (50) The relationships between spherical and local Rayleigh and Ekman numbers are 3RaT Ra8 = 6 and E~= 6

(51)

(52)

2~

We substitute (51) and (52), along with the relationship (49), into (50), obtaining den to the first zero of the eigenfunction, depends 2i/2

on the Ekman number approximately as 6

3Ei/5

Ra~~6 = (49)

The Coniolis acceleration acts to stabilize the fluid, reducing the scale of motion and delaying the onset of convection. With decreasing Ekman number, the axial rigidity of fluid columns increases, and radial motion of the nearly rigid fluid columns is inhibited by the spherical outer boundary. At small Ekman numbers, fluid columns are restricted to small radial displacements, and this effectively limits the radial extent of the convecting region. Instead of filling the entire spherical shell, convection near the critical Rayleigh number is confined, at small Ekman numbers, to a thin cylindrical shell around the inner core, where the stretching of columns is least. Outside this narrow shell, buoyancy forces

~

4/3~

Pr ~4/3 7/i5 + Pr) E~

(53)

For the parameters in our experiment, (53) is approximately ~

8E

i7/i5

(54)

Formula (54) is shown by the dashed line in the regime diagram (Fig. 6). It is clearly an excellent fit to the marginal stability curve calculated for the sphere, especially for small Ekman numbers. We find that Busse’s (1970a) results for the thin cylinder model agree with our results for the onset of convection in a spherical shell with Pr = 7 fluid, provided the thickness of the cylinder in his model is interpreted as the radial length-scale of the flow, and not the thickness of the entire spherical shell.

P. Cardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

8. Finite-amplitude calculations The general characteristics of the experimental flows near the onset of convection—columnar motion limited to the tangent cylinder surrounding the inner core—are similar to the linear stability calculations. However, to make a meaningful comparison with the experiments on convection at large Reynolds number, it is necessary to solve the full nonlinear Eqs. (31)—(38) at Rayleigh numbers fan beyond the critical value, The parameter range of the experiments, mdicated by the shaded region in Fig. 5, includes Ekman numbers of 6 X 10-6 and smaller. We have found that it is very difficult to make neliable numerical calculations in this range, because the fine scales of motion require more spatial and temporal resolution than we can achieve with the computational resources available to us. Accordingly, we have concentrated on understanding the dynamics of high Rayleigh number thermal convection in a rapidly rotating sphere in the more accessible parameter range iO~~ E ~ 2 x iO~. We have obtained time-dependent numerical solutions to the column-averaged equations of motion (31)—(38) using finite difference methods on grids with uniform spacing (sr, ~0). We use second-order central difference formulae for spatial derivatives, except for advective terms, which are calculated using up-wind differences. The solution is advanced in time using an explicit Runge—Kutta time step, also correct to second order. At each time step, the vonticity and ternpenatune fields are advanced using (32) and (33),

and then the streamfunction is updated, by solving (31) using successive over-relaxation. We found it necessary to restrict the time step to about 0.1t~,where t~is the Counant time based on convective velocity. Very short time steps were required to suppress spurious numerically excited wave motion in the solution, particularly at the lowest Ekman numbers. Table 5 lists the parameters of calculations illustrated with figures. These calculations were continued until they reach either steady state or, in cases where the flow is intrinsically time dependent, a statistical thermal equilibrium. The only exception is the highest Rayleigh number case at E = 2 X iO~,which may not have reached thermal equilibrium. Parameter combinations for several of these calculations are plotted on the regime diagram in Fig. 6. 8.1. Solutions for slightly supercritical Rayleigh number

Beyond the critical Rayleigh number the pattern of convection begins to deviate from charactenistics of the flow at marginal stability. Figure 8 shows contours of streamfunction, temperature and vorticity at RaT= 1.1Ra~~~ for E = iO~,E = iO~,and E = 2 x iO~from the finite difference calculations. The proportional relationship between cell size and Ekman number is evident. At these Rayleigh numbers the flow remains penodic in longitude and confined to the region near the tangent cylinder. It should be noted that the streamfunction consists of cells that are tilted in

TABLE 5 Calculation parameters Figure 8(a) 8(b) 8(c)

9(a) 9(b)

9(c) 10(a) 10(b) 10(c) 11

E

xRa~~~

Grid

Planform

1x10’ 1x104 2x105 lx iO~ ix iO~’ 2x iO~ lx iO~

1.1 1.1 1.1 3 3 3 10 10 10 200

40x150 4tJxl5O 50x300 40 x 150 40x 150 50x300 40x 150 40x150 50x300 40x 150

Pinwheel Pinwheel Pinwheel Aperiodic pinwheel Aperiodic pinwheel Aperiodic pinwheel Chaotic Chaotic Chaotic Highly chaotic

ixio4 2x105 lx i04

253

Zonal/total energy 4 5x10 1x103 1x103 0.05 0.09 0.24±0.08 0.35±0.05 0.6±0.1

254

P Cardin, P Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259

tent with the calculations. Under these condi-

7

(~) (:I~II)~

Hons,(32~educesto

(55)

The streamfunction and temperature perturbation for this flow can be written

where a(r) is the phase function measuring the tilt of the cells in the equatorial plane. Substituting(56) and (57) into (55) yields the tilt relation-

o

9’ w Fig. 8. Contours of streamfunction, temperature and axial vorticity for thermal convection in a rapidly rotating spherical shell viewed from above for Pr =7 and RaT 1.lRaerit at Ekman numbers E= ~ iO~ and 2x105, calculated by finite difference using the columnar flow model. Rotation of the sphere is counter-clockwise in this view. Numerical grid

tan(a)

ELa3 2~r3

=

(58)

The trajectory of the cells is defined by the constant phase curve 3

a

-i

—

2

2

Ea (r

04 r ) (59) 2r According to (59), the angle of a cell 0 is negative, so the cells tilt in the prognade direction. In addition, Eq. (59) indicates that the tilt angle increases with radius, which means that the trajectory of a fluid parcel is a prograde spiral. —

in the pinwheel planform regime. parameters are listed in Table 5. All of these calculations are

the prograde direction, forming a spiral or pinwheel pattern. The temperature perturbations, in contrast, do not show the same degree of tilt. The tilted, spiral pattern of motion in the pinwheel planform is a consequence of the Coniohis acceleration acting on fluid columns. The Coniolis acceleration modifies the vorticity of each column as it moves toward or away from the notation axis, and the fluid adjusts by altering its trajectory in the (r,0)-plane, forming spinalshaped cells, We can demonstrate that the Coniolis acceleration produces the spinal effect by determining analytically the planform of steady-state periodic cells for which viscous, Coniolis and buoyancy forces are in balance. To simplify the analysis, we adopt the following approximations. First, we let 8/8r 1/r(a/aO), which implicitly assumes the cells are narrow. We also allow the flow to be tilted, but not the driving temperature perturbations. Both of these approximations are consis.~

0

=

—

—

=

—

a tan

8.2. Transition to chaotic planforms

As the Rayleigh number is increased further beyond the critical value, the periodicity in the pinwheel planform is replaced by a chaotic pattern consisting of irregularly distributed vortices plus a lange-scale zonal flow. The transition from the pinwheel to chaotic flow is illustrated in Fig. 9—13. Figures 9 and 10 show contours of streamfunction, temperature and axial vorticity at three Ekman numbers, E = i0~, iO~and 2 x iO~, for Rayleigh numbers equal to 3Ra~~~ and 10Ra~~, respectively. Figure 11 shows contours of streamfunction, temperature and axial vorticity at various times for E = i0~at RaT = 200Ra~~~. In Fig. 12 we show a snapshot of the radial profile of the zonal flow for this case. Figure 13

255

P. Gardin, P Olson /Physics ofthe Earth and Planetary Interiors 82 (1994) 235—259

shows time series of the total and zonal kinetic energy. At RaT = 3Racnt, the deviation from a penodic planform remains small for the entire calculation, but at RaT = 10Ra~~ the flow becomes fully chaotic after a dimensionless time t = 0.06, corresponding to only a few eddy circulation times. The average prognade tilt of the convective periodicity number eddies persists is gone. found atInnearer RaT addition, = 10Ra~~~, to the there critical but is a the zonal Rayleigh spatial flow at 10Ra~~~, which is retrograde near the innercore tangent cylinder and prograde near the equator. When RaT>> ~ the zonal flow increases in strength until its kinetic energy becomes comparable with the kinetic energy in the convective eddy motion. The zonal flow draws the convective plumes into a thin tendril in the prognade direction, as illustrated in Fig. 11 for the case Ra = 200Racrjt and E = iO~.The vorticity associated

,

c

• ~ •

Li ii

w

~ ~ 9’

Fig. 10. Contours of streamfunction, temperature and axial vorticity for thermal convection in a rapidly rotating spherical shell viewed from above for Pr = 7 and RaT = iORaerit, at Ekman numbers E = i0~, i04 and 2 x i0~, calculated by finite difference using the columnar flow model. Rotation of the sphere is counter-clockwise in this view. Numerical grid parameters are listed in Table 5. All of these calculations are in the chaotic regime.

~0

‘ci

—II~ Li

9,

Fig. 9. Contours of streamfunction, temperature and axial vorticity for thermal convection in a rapidly rotating spherical shell viewed from above for Pr 7 and RaT = 3Raerit, at Ekman numbers E = i0’, i04 and 2x i0~, calculated by finite difference using the columnar flow model. Rotation of the sphere is counter-clockwise in this view. Numerical grid parameters are listed in Table 5. All of these calculations are in the pinwheel planform regime, but deviations from cornplete periodicity are evident,

with each plume is localized in discrete centers along the trajectory of the plume. This is consistent with the behavior of the vortices observed in the experiments at highly supercnitical Rayleigh number (see Figs. 3 and 4). The magnitude of the zonal flow and its variation with radius are shown in Fig. 12 at a particular time in the RaT = 200Raerjt and E = iO~calculation. There is continual exchange of energy between the zonal flow and the convection columns. The zonal flow extracts energy from the convection and also and tendsmixing to suppress the into convection by shearing the plumes the fluid. Figure 13 shows time series of kinetic energy for the 200Ra~~~, E = iO~case. The bold continuous curve is the total (nondimensional) kinetic energy, the thin continuous curve is the kinetic energy in the zonal flow, and the dashed curve is the percentage of the total energy in the zonal flow, as a function of time. Small arrows indicate the times displayed in Fig. 11. Kinetic energy

256

P Gardin, P Olson/Physics of theEarth and Planetary Interiors 82 (1994) 235—259

~

tO.I~

~•

Di

•

5.0x104 1 .0 xl 0

~

0

t0.1315

Fig. 11. Contours of streamfunction, temperature and axial vorticity for thermal convection in a rapidly rotating spherical shell viewed from above for Pr 7, Rar =200 Raerjt, and E = i0~, at various times. Rotation of the sphere is counterclockwise in this view. Numerical grid parameters are listed in Table 5. This calculation is in the chaotic regime.

peaks correspond to episodes of rapid plume growth. The zonal energy tracks the total energy, with a small phase lag, indicating that kinetic

7

200

80 100 60

C

•

~

~

5 ~ 2.5 xl ~ 2.0x10 1.5x10~ a, >-

~

•

• 20 40

~

C)

thermal 0.110 in ain 0.120 rapidly rotating 0.130energy 0.140 spherical shell in at Fig. 13. convection The variation thetime kinetic with time Ra = 200 ~ Pr = 7, and E = l0~. Bold line: total energy; thin line: zonal energy; dashed line: percentage of total energy in the zonal part. Arrows mark the times of images in Fig. 11.

energy is transferred from the plumes to the zonal flow. It has been shown previously that this columnar zonal flow is maintained by Reynolds stresses, which are present in the fluid as a result of the tilt of the convective eddies (Busse and Hood, 1982; Zhang, 1991). Adapting the method of Busse and Hood to our problem, we can account for the generally westward nature of the zonal stress in the zonal flow balances the Reynolds stress produced by a periodic array of steady-state flow using a simplified model in which viscous tilted cells. To simplify the analysis, we adopt a thin= gap geometry, in which r r, = x and rO y. The zonal average ofwe theset0-component of —

_____________________________________

0 >

82U ar (60) 8x28x the momentum equation (Eq. (8)) then becomes where U = Ku 9> is the zonal (0-averaged) az-

‘5 0

N ‘5

y//

0 —200

______________________________ 0.6

0.8

1.0 radius

1.2

1.4

Fig. 12. The variation of zonal (azimuthally averaged) velocity with radius in the equatorial plane of the rotating spherical shell in chaotic convection at Ra = 200 Racrit, Pr = 7, and E = iO~. The scale factor for dimensionless velocity is the Reynolds number, UD/v.

imuthal velocity, and r denotes the *y-component (i.e. the r0-component) of the Reynolds stress, ~

T=

(61)

The streamfunction for an array of tilted cells can be written as i/i

=

iIi0(x) sin(ay + Ia)

(62)

P Gardin, P

Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

Negative values of the product ak correspond to prograde tilt. The Reynolds stress produced by (62) is 2 =

—

(63)

ak~Iio

2 We let ~1’o( x)

=

‘P0sm( m~rx)

(64) Substituting (63) and (64) into (60) yields, after integration U=

ak ‘~‘o~1 2 ~

—

x

1

,

+ .~—sm(2mirx)

(65)

We assign a negative value to ak, cornespondmg to prognade tilted cells. Then, according to (65) the zonal flow is negative (retrograde) and strongest relatively near the more innerprognade sphere (small x),outer and becomes near the sphere. 9. Comparison of laboratory experiments and numerical calculations

257

convection (the supercnitical regime). Periodic, cellular convection is the transition between these two broad regimes, and occurs oven a relatively narrow range of Rayleigh numbers. We can derive general scaling laws for length and velocity scales in the high Reynolds number regime, to nationalize the similarities and differences between the laboratory experiments and the columnar flow calculations, and also to extrapolate to conditions in the outer cone. We let 6 represent the length-scale for the diameter of an eddy in the equatorial plane. The equatorial velocity and axial vorticity are then related by Ue o6. A balance between advection of vorticity, vortex stretching and vortex generation by buoyancy in Eq. (32) requires that the following approximate relationships must hold: 2 t~e Ra (66) ‘~ 6 EL T 6 As in Section 4, we can express the (dimensionless) temperature fluctuations in the fluid in terms of the heat transport, using Nu —

~

T’

(67)

=

u~Pr2

In both the laboratory experiments and the numerical calculations, we find that the pinwheel pattern is stable to only 10Ra~,~, approximately. Why is the regime of periodic, cellular convection so limited? The combined effects of rotation and sphenicity provide the answer. Rotation in the spherical geometry greatly increases the critical Rayleigh number, so that a large amount of po-

Solving for the eddy length and velocity scales yields

tential energy is stored in the fluid. At the onset of convection, the steeply sloping spherical boundaries insure that numerous small-scale vontices form. The spatial organization of a pattern consisting of many vortices is easily disrupted, especially at high Rayleigh number. Because it is impossible to preserve periodicity of so many vortices when the buoyancy force is finite, the flow evolves into a chaotic pattern at a slightly supencritical Rayleigh number. In summary, two dynamic regimes cover large areas of the parameter space for thermal convection in a rapidly rotating sphere: no convection (the subcritical regime) and chaotic, columnar

Ue

3/5

6

— —

~-

~ k Pr2 )

i/5

‘~

C~

~68

d an EL

-i--—

i/5

RaQ ~2/5 ~.-~--~-)

1

(69) ‘11 where c1 and c2 are constants. In deriving (68) and (69) we have used the relationship Nu = RaQ/RaT. From the experimental results we find c~ 5 and c2 1. According to the above scaling laws, both the fluid velocity and the column diameter increase with Rayleigh number, and decrease with Ekman number. These are the qualitative trends found in both the experiments and in the calculations at high Reynolds number. Formulae (68) and (69) also provide reasonably good quantitative fits to the numerical calculations with large Reynolds c2

258

P Gardin, P Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

number. For example, the mean fluid velocity in calculation at E = iO~and Ra = 10Ra~~shown 1 in Fig. 10 is in the range 20—25, and the eddy diameter, on the basis of the vorticity pattern, is about 1/5. Assuming that the experimentally denived relation RaQ 6RaT is valid, and setting L/2fl = 1, formulae (68) and (69) give 6 0.24 and Ue 25. These relationships also predict that the column diameter and convective velocity are small near the outer boundary, where the slope parameter t~ becomes large, which is consistent with our experimental results.

buoyancy forces are small compared with the

Coriolis acceleration in the fluid, which, we have argued, is also likely to be the case in the outer core. We note that the symmetry properties of chaotic columnar convection are consistent with the equatorial mirror-symmetry inferred for flow in the outermost core by Hulot et ah. (1990), using geomagnetic secular variation data. Chaotic columnar convection is also consistent with the non-periodic character of magnetic flux bundles on the core—mantle boundary inferred by Gubbins and Bloxham (1987). It is interesting to extrapolate our results to 10. Implications for the outer core cone conditions, to estimate the column diameter and fluid velocity expected for thermal convecWe have used laboratory and numerical expention in the core if there were no magnetic field. iments to determine the structure of chaotic conUsing the parameters in Tables 1 and 2, (68) and vection in a rapidly rotating spherical shell. We (69) predict u~ 108_109 for the Reynolds numfind that thermal convection is columnar, for ben and 6 10— 2~These are equivalent to fluid Rayleigh numbers up to 200 times the critical velocities and column diameters of about 0.1 cm value. When the Rayleigh number exceeds Ra~~1 s’ and 20 1cm, respectively. by a factor of 10, approximately, the periodic The fluid velocity predicted by this model is pinwheel pattern found at lower Reynolds numcomparable with the fluid velocities derived from bens is replaced by a highly time-dependent, geomagnetic secular variation (Bloxham et al., chaotic pattern, consisting of an irregular distri1989). However, the column diameter is much too bution of columnar vortices driven by plumes in small to correspond to any fluid motion detected the equatorial region, plus a zonal flow. The with secular variation. Indeed, the most signifizonal flow is dominantly retrograde but becomes cant difference between our results and structure prograde near the the outer-sphere equator, and of cone flow inferred from the geomagnetic secuis driven primarily by Reynolds stresses derived lan variation is that the eddies are much larger from the prognade tilt of the columnar convecdiameter in the cone than purely hydnodynamical tion. convection predicts. This difference in scale may Our results strongly indicate that core convecbe due to the action of Lonentz forces in the core, tion is turbulent. We have found that the regime which we have not considered in this paper. of laminar, periodic convection in a notating sphere is restricted to the immediate neighborhood of the critical Rayleigh number, where the Reynolds number of the motion is small. The fact Acknowledgments that the Reynolds number for the core is large indicates that laminar convection is highly imThis research has been supported by the Naprobable in the outer core. In contrast to the very tional Science Foundation through Grant restricted laminar regime, our experiments and EAR8916152. We thank Monitz Heimpel and calculations indicate that the chaotic columnar James Kinkhin for their assistance in this project, regime is found oven a very large range of Harvey Palmer for designing the speed control Rayleigh number values, in situations where the device used in the experiments, and F. Busse, C. Reynolds number of the motion is large, as it is in Cannigan and two anonymous referees for their the cone. Columnar convection occurs when helpful comments.

P. Gardin, P. Olson /Physics of the Earth and Planetary Interiors 82 (1994) 235—259

References Bloxham, J. and Roberts, P.H., 1991. The geomagnetic main field and the geodynamo. Rev. Geophys. Space Phys. Suppl., pp 428—432. Bloxharn, J., Gubbins, D. and Jackson, A., 1989. Geomagnetic secular variation. Philos. Trans. R. Soc. London, 329; 415—502. Braginskii, S., 1990. Towards a realistic theory of the geodynamo. Geophys. Astrophys. Fluid Dyn., 60: 89—134. Brummell, N.H. and Hart, i.E., 1993. High Rayleigh number /3-convection. Geophys. Astrophys. Fluid Dyn., 68; 85—132. Buffett, B.A., Huppert, H.E., Lister, J.R. and Woods, A.W., 1992. Analytical model for solidification of the Earth’s core. Nature, 356; 329—331. Busse, F.H., 1970a. Thermal instabilities in rapidly rotating systems. J. Fluid Mech., 44: 441—460. Busse, F.H., 1970b. Differential rotation in stella convection zones. Astrophys. J., 159; 629—639. Busse, F.H., 1975. A model of the geodynamo. Geophys. J. R. Astron. Soc., 42; 437—459. Busse, F.H., and Carrigan, CR. 1976a. Laboratory simulation of thermal convection in rotating planets and stars. Science, 191; 81—83. Busse, F.H. and Carrigan, C.R. 1976b. Convection, induced by centrifugal buoyancy. J. Fluid Mech., 62: 579—592. Busse, F.H. and Hood, L.L, 1982. Differential rotation driven by a convection in a rapidly rotating annulus. Geophys. Astrophys. Fluid Dyn., 21; 59—74. Cardin, P. and Nataf, H.-C., 1991. Thermal coupling in layered convection; marginal stability analysis. J. Phys., Paris, II, 1; 599—622. Cardin, P. and Olson, P., 1992. An experimental approach to thermochemical convection in the Earth’s core. Geophys. Res. Lett. 19: 1995—1998. Carrigan, C.R. and Busse, F.H., 1983. An experimental and theoretical investigation of the onset of convection in rotating spherical shells, J. Fluid Mech., 126; 287—305. Chamberlain, J.A. and Carrigan, C.R., 1986. An experimental investigation of convection in a rotating sphere subject to time varying thermal boundary conditions, Geophys. Astrophys. Fluid Dyn., 41: 17—41. Cordero, S. and Busse, F.H., 1992. Experiments on convection in rotating spherical shells: transition to a quasiperiodic state, Geophys. Res. Lett., 19: 733—736. Durney, B., 1970. Nonaxisymmetric convection in a rotating spherical shell. Astrophys. J., 161: 1115—1127. Fearn, D.R., 1989. Compositional convection and the Earth’s core, In: F.J. Lowes, D.W. Collinson, J.H. Parry, S.K. Runcorn, D.C. Tozer and A. Soward, Geomagnetism and Paleornagnetism. Kluwer Academic, Dordrecht, 335—346. Glatzmaier, G.A., 1985a. Numerical simulations of stellar convective dynamos, II: Field propagation in the convection zone. Astrophys. J., 291: 300—307. Glatzmaier, G.A., 1985b. Numerical simulations of stellar convective dynamos, III: At the base of the convection zone. Geophys. Astrophys. Fluid Dyn., 31: 137—150. Glatzmaier, G.A. and Olson, P., 1993. Highly supercritical

259

thermal convection in a rotating spherical shell: centrifugal vs. radial gravity. Geophys. Astrophys. Fluid Dyn., 70: 113—136. Gubbins, D. and Bloxham, J., 1987. Morphology of the geomagnetic field and implications for the geodynamo. Nature, 325: 509—511. Gubbins, D., Masters, T.G. and Jacobs, J.A., 1979. Thermal evolution of the Earth’s core. Geophys. JR. Astron. Soc., 59: 57—99. Hart, i.E., Toomre, J., Deane, A.E., Hurlburt, N.E., Glatzmaier, G.A., Fichtl, G.H., Leslie, F., Fowlis, W.W., and Gilman, P.A., 1986a. Laboratory experiments on planetary and stellar convection performed on Spacelab 3. Science, 234: 61—64. Hart, i.E., Glatzmaier, G.A. and Toomre, J., 1986b. Spacelaboratory and numerical simulations of thermal convection in a rotating hemispherical shell with radial gravity, J. Fluid Mech., 173: 519—544. Hulot, G., Ic Mouel, J.-L. and Jault, D., 1990. The flow at the core—mantle boundary: symmetry properties, J. Geomagn. Geoelectr., 42: 857—874. Jeanloz, R., 1990. The nature of the Earth’s core. Annu. Rev. Earth Planet. Sci., 18: 357—386. Loper, D.E., 1978. Some thermal consequences of a gravitationally powered dynamo. J. Geophys. Res., 83: 5961—5970. Loper, D.E. and Roberts, P.H., 1981. A study of conditions at the inner core boundary of the Earth, Phys. Earth Planet. Inter: 24: 302—307. Or, A.D. and Busse, F.H., 1987. Convection in a rotating cylindrical annulus, Part 2: Transition to asynunetric and vacillating flow, J. Fluid Mech., 174: 313—326. Roberts, P.H., 1968. On the thermal instability of a self-gravitating fluid sphere containing heat sources. Philos. Trans. R. Soc. London, Ser. A, 263: 93—117. Soward, A.M., 1977. On the finite amplitude thermal instability of a rapidly rotating fluid sphere. Geophys. Astrophys. Fluid Dyn., 9: 19—74. Stevenson, Di., 1981. Models of the Earth’s core. Science, 214: 611—619. Sun, Z.-P., Schubert, G. and Glatzmaier, G.A., 1993. Transition to chaotic thermal convection in a rapidly rotating spherical shell. Geophys. Astrophys. Fluid Dyn., 69: 95— 131. Vanyo, J., Wilde, P., Cardin, P. and Olson, P., 1994. Experiments on processing flows in the Earth’s liquid core. Science, submitted. Verhoogen, J., 1980. Energetics of the Earth. National Academy of Sciences, Washington, D.C., 135 pp. Yano, J.L., 1992. Asymptotic theory of thermal convection in rapidly rotating systems. J. Fluid Mech., 243: 103—131. Zhang, K., 1991. Convection in a rapidly rotating spherical shell at infinite Prandtl number: steadily drifting rolls. Phys. Earth Planet Inter., 68: 156—169. Zhang, K., 1992. Spiralling columnar convection in rapidly rotating spherical shells. J. Fluid Mech., 236: 535—556. Zhang, K. and Busse, F.H., 1987. On the onset of convection in rotating spherical shells. Geophys. Astrophys Fluid Dyn., 39: 119—147.