This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright
Transcript
This article appeared in a journal published by Elsevier. The attachedcopy is furnished to the author for internal non-commercial researchand education use, including for instruction at the authors institution
and sharing with colleagues.
Other uses, including reproduction and distribution, or selling orlicensing copies, or posting to personal, institutional or third party
websites are prohibited.
In most cases authors are permitted to post their version of thearticle (e.g. in Word or Tex form) to their personal website orinstitutional repository. Authors requiring further information
regarding Elsevier’s archiving and manuscript policies areencouraged to visit:
Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Contents lists available at ScienceDirect
Physics of the Earth and Planetary Interiors
journa l homepage: www.e lsev ier .com/ locate /pepi
On the rise of strongly tilted mantle plume tails
C.A. Mériauxa,∗, J.A. Mansourb, L.N. Moresib, R.C. Kerrc, D.A. Mayd
a Instituto Dom Luiz, Departamento de Engenharia Geografica, Geofisica e Energia, Faculdade de Ciencas, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugalb School of Mathematical Sciences, Monash University, Victoria 3800, Australiac Research School of Earth Sciences, Australian National University, Canberra, 0200 ACT, Australiad Geophysical Fluid Dynamics, Department of Earth Sciences, ETH Zürich, Zürich, Switzerland
a r t i c l e i n f o
Article history:Received 6 May 2010Received in revised form 13 October 2010Accepted 15 October 2010
The rise of an initially horizontal, buoyant cylinder of fluid through a denser fluid at low Reynolds num-ber is used to look at the ascent of strongly tilted mantle plumes through the mantle. Such ascentsare characterized by (1) the growth of instabilities and (2) the development of a thermal wake down-stream. Three-dimensional numerical experiments were carried out to examine these features. An hybridparticle-in-cell finite element method was used to look at the rise of non-diffusing cylinders and, a stan-dard finite element method was used to look at the diffusing case. First the experiments show that thetimescale of the fastest growing instability vary with the Rayleigh number and the viscosity ratio. Inparticular the growth rate decreases as the Rayleigh number decreases, in agreement with our analysisof the laboratory experiments of Kerr et al. (2008). Second the experiments show that the length of thethermal wake increases with the Rayleigh number but the change in viscosity has almost no influenceon the wake length. Applied to strongly tilted mantle plumes we conclude that such plumes cannot beunstable given the plume timescales. We also discuss the application of this conclusion to weakly tiltedplumes. Besides, this study allows to predict that mantle plumes are unlikely to have developed a signif-icant thermal wake by the time they reach the surface. Finally, the resolution that is required to allow forthe growth of mantle plume tails by combined diffusion and thermal entrainment is shown to representa challenge for the large scale mantle convection simulations.
Thermal plumes are an intrinsic feature of convection. Theirlong-lived shape can loosely be described by a cylindrical conduit.Provided that they are not subjected to a shearing current, plumeswill rise opposite to gravity, whereas, in the presence of a velocityshear, they will tilt.
The idea that plumes are deflected from the vertical has espe-cially been portrayed in the case of the Earth’s mantle convectionbecause of the large scale plate-driven shear flow (Steinberger andOConnell, 1998). In particular, laboratory experiments (Skilbeckand Whitehead, 1978; Whitehead, 1982; Olson and Singer, 1985;Kerr and Lister, 1988) showed that gravitational instabilities ofcompositional plume conduits could occur if the conduits wereeither horizontal or tilted by a background shear by more than60◦ from the vertical. In the laboratory experiments by Kerr andLister (1988), the diffusion was compositional and not thermal,which resulted in slight diffusion only. The nature of these grav-
itational instabilities were analyzed by Lister and Kerr (1989) ina linear stability analysis of a buoyant cylinder of non-diffusingfluids. In particular, the authors found that a cylinder of radiusa is gravitationally unstable, with the most unstable wavelength� given by � ≈ �∞a, with �∞ = 6.28 when the viscosity contrastM = 1 and �∞ = 8.09 when M> ∼10. The authors further pro-posed that the timescale of the fastest growing mode � was givenby � ≈ c∞(�0/��ga), where �0 is the outer viscosity, �� is the den-sity difference between the outer fluid and the fluid within thecylindrical region, g is the gravitational acceleration and c∞ is aconstant, which equals 9.35 whenM = 1 and converges to a limitvalue of 2.92 as M> ∼10 (see Fig. 9 of Lister and Kerr, 1989).Recently, the effect of thermal diffusion on the gravitational sta-bility of a rising horizontal cylindrical region of buoyant viscousfluid was examined by Kerr et al. (2008). At large viscosity ratios,the authors found that the instability is unaffected by diffusionwhen the thermal Rayleigh number Ra = ��ga3/��0, defined by thetemperature-driven density contrast �� between the warm cylin-drical region and cool ambient fluid, the thermal diffusivity � andthe viscosity �0 of the ambient fluid, is greater than about 300.When Ra is less than 300, diffusion significantly increases the timefor instability, as the rising fluid region needs to grow substantially
64 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
by both diffusion and entrainment before it becomes unstable. Forinstance, when Ra was less than about 140 and within an avail-able rise height of about 40 times the cylinder radius, the risingregion of fluid was unable to grow sufficiently and the instabilitywas not seen. In particular a dimensionless time for gravitationalinstability at which the amplitude of the gravitational instabilitieshad grown to be equal to the cylinder diameter was measured.This characteristic time for instability was shown to increase atlow Rayleigh numbers (Fig. 4 of Kerr et al., 2008) and suggesteda trend for a decrease in growth rates at low Rayleigh numbers.Besides gravitational instabilities, Kerr et al. (2008) showed that athermal wake left behind the buoyant cylinder develops during therise. This downstream loss of buoyancy during a viscous rise of athermally buoyant body was first predicted by Whittaker and Lister(2008).
Applied to the Earth’s mantle, the study of Kerr et al. (2008)presented two predictions for strongly mantle plume tails. First,strongly tilted plume tails were to be stable thanks to diffusionin both the upper and lower mantle. This conclusion was essen-tially drawn from the low range of Rayleigh numbers Ra thatwas associated with each plume (4.8 ≤ Ra ≤ 104). Importantly thisrange relied on estimates of the plume radius a that was derivedfrom a pipe flow model in which the plume buoyancy flux wasgiven by models of swell uplift of the lithosphere (Davies, 1988;Sleep, 1990). Second Kerr et al. (2008) predicted the loss of someof the plume buoyancy flux from the tail in a thermal wakeduring the rise. This prediction implicitly meant that estimatesof plume buoyancy fluxes based on mantle swells at the sur-face were likely to represent underestimated values of the plumeflux.
The objective of this paper was to further investigate the rise ofstrongly tilted thermal plumes with an emphasis on the timescalesfor the growth of the instability and the development of the ther-mal wake. First we looked at the growth of the instabilities througha numerical study with experiments considering the case of the riseof an initially horizontal cylinder in a three-dimensional domain. Akey to this work was to work in a 3D domain because the features ofthe Rayleigh–Taylor instability to be investigated were intrinsicallythree-dimensional. We thus studied the growth rates using twotypes of numerical experiments depending on whether we con-sidered non-diffusing or diffusing cylinders (i.e. infinite and finiteRa, respectively). The rise of non-diffusing cylinders was simulatedusing a hybrid particle-in-cell finite element method while simula-tions of the diffusing cylinders were simulated via a standard finiteelement method. We also analyzed the growth rates of the insta-bilities for 17 of the laboratory experiments by Kerr et al. (2008) forcomparison. Second we looked at the development of the thermalwake during the rise. Implications for both the stability and loss ofbuoyancy of strongly tilted mantle plumes were further analyzedin view of the mantle timescales. In particular we will show that thetimescales provide strong evidence of (1) the stability of stronglytilted mantle plumes and (2) the unfeasibility of mantle plumesto have significantly lost their buoyancy by the time they reach thesurface. We also discuss the application of this conclusion to weaklytilted plumes.
In Section 2, we describe the numerical problem. Our numericalexperiments at infinite Rayleigh number are presented in Section3, while we detail the experiments at finite Rayleigh numbers inSection 4. We present the characterization of the growth rate for thegravitational instabilities that were obtained from the laboratoryexperiments of Kerr et al. (2008) in Section 5. Further discussion isgiven in Section 6. Section 7 discusses the implications for mantleplumes. Section 8 draw the conclusions. A dimensional analysis inthe finite domain is presented in Appendix A. Appendix B detailsthe methods of analysis of the growth rate. The thermal boundarylayer is finally discussed in Appendix C.
Fig. 1. Schematic of the numerical model.
2. Numerical problem description
We consider a cylinder of fluid of radius a, that is initially laid outhorizontally in a fluid delimited by a box (Fig. 1). The fluid withinthe cylindrical domain has an initial temperature T = 1, while theouter fluid is at temperature T = 0. We assume that the density �within the box follows
� = �0(1 − ˛T), (1)
where �0 is the initial density of the fluid outside the cylinder and˛ is the coefficient of thermal expansion. Driven by the densitydifference �� = �0 − �, which initially is ��0 = ˛�0, the cylinderrises with a vertical velocity that is proportional to
U = g��a2
�0. (2)
We assume that the viscosities of the outer fluid �0 and the cylin-der inner fluid �1 are either constant or following an Arrheniusrheology (Reese et al., 1999):
�(T) = �∗0 exp
(E
T + 1
). (3)
The parameters E and �∗0 were chosen so that �(0) = 1 and �(1) =
1/M. We note that the total buoyancy per unit length of the cylin-der is B = ��ga2.
The system is fully described by Stokes equations for incom-pressible creeping flow and the thermal energy equation. Whenthe fluid properties are non-dimensionalized with respect to thoseof the outer fluid �0 and �0 and �, the lengths with respect to theradius of the cylinder and time with respect to �0/��0ga, thoseequations take the non-dimensionalized form:
−∇(�D) + ∇p = Tn, (4)
∇ · u = 0, (5)
where Dij =(∂ ui/∂ xj + ∂ uj/∂ xi) is the dimensionless strain rate ten-sor with dimensionless velocity u, p is the dimensionless dynamicpressure, and n = g/g is a unit vector with g the gravitational accel-eration vector, and
∂T
∂t+ u · ∇T = 1
Ra∇2T. (6)
The system is found to be described by two dimensionless num-bers that are the Rayleigh number
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 65
and the viscosity ratio
M = �0
�1. (8)
In the limit of infinite Rayleigh number, thermal diffusion can beneglected (see Eq. (6)) and, the rise is equivalent to the rise of thecompositionally buoyant cylinder fluid. When thermal diffusion ispresent, the radius a of the buoyant cylinder fluid increases duringthe rise while �� decreases. We note that a is yet used in referenceto the initial radius of the cylinder throughout this paper.
3. Numerical experiments at infinite Rayleigh number
3.1. Model setup and numerical method
The models were set in a box of dimensions Lx, Ly, and Lz inthe x, y and z directions respectively. Model cylinders of radius awere initially laid out across the box in the x direction, centeredin the box in the z direction and at a given height, which wastypically toward the bottom of the box, in the y direction. Cylin-ders that are either horizontal or slightly tilted from the horizontalalong the x direction were considered. Free-slip velocity bound-ary conditions were applied on all walls within the model domain.The calculations were performed using a hybrid particle-in-cellfinite element method (Moresi et al., 2003, 2007) together with amultigrid solver and implemented in the StGermain/Underworldframework (http://www.underworldproject.org/). The visualiza-tion used the gLucifer’s architecture (Duboz et al., 2005). Thecomputational domain had Nx, Ny, and Nz elements in the x, y andz directions respectively. The values of Nx, Ny, and Nz were chosento ensure that the geometry of the cylinder was well resolved inall directions (Li/Ni ≤ a in each i direction). The initial distributionof particles was such that each element contained approximately30 particles. This number, together with a random (non-uniform)particle coordinate distribution helped to ensure that voids due tostretching of the fluid would not occur during the calculations. Boththe density and viscosity ratio between the cylinder (inner fluid)and the outer fluid were initially prescribed. All calculations wererun for times that are long compared to the timescale for growth� given by Lister and Kerr (1989). The calculations were typicallyperformed on 32 processors with a memory usage of 50 Gb usingthe SGI Atlix of the NCI facility (http://nf.nci.org.au/). Table 1 detailsthe experiments that were carried out.
3.2. Wavelengths of instabilities
Gravitational instabilities developed in all experiments exceptexperiments 10I and 16I (see Table 1). Fig. 2 shows a typical devel-opment of the instabilities during the buoyant cylinder rise forexperiment 14IPT. The instability becomes evident after it has risenover a distance which is approximately twice the cylinder diame-ter. The wavelength of the instabilities is found to depend on theradius of the cylinder. For instance, four plumes developed in exper-iment 1I (Fig. 3a) while they reduced to two in experiment 6IPT(Fig. 3b), where the radius of the cylinder was twice the radius usedin experiment 1I.
When the viscosity ratio M = 1, we find that the mean plumespacing �num = �num/a is � = 10.21 ± 1.4. This mean plume spac-ing appears slightly larger than the value of 6.28 predicted by thelinear stability of Lister and Kerr (1989). A similar larger meandimensionless plume spacing of 9 ± 1 at Ra = ∞ and M = 1 wasfound in the study of Kerr et al. (2008) and invoked to result from thesuppression of weak plumes overwhelmed by their rapidly grow-ing neighbours. When the viscosity ratioM ≥ 10, the mean plumespacing is �num = 10.60 ± 1.9 whenM = 50. This value agrees rea-sonably well with the value of 8.09 predicted by the linear stability
Fig. 2. Snapshots of experiment 14IPT at dimensionless times t = t/(�0/��0ga) =0, 40, 80 and 120 from a front and a perspective view. The passive tracers are notshown in these frames.
of Lister and Kerr (1989) and the result of Kerr et al. (2008), whoreported a mean value of �/a of 12.7. Overall our wavelengths areindependent of the viscosity ratio.
3.3. Boundaries and development of instability
As shown in Appendix A from dimensional analysis, the distanceof rise necessary to the development of the instability depends onthe rise velocity and the growth rate. These two quantities are func-tionals that especially depend on the dimensions of the box relativeto that of the cylinder. So one expects that the dimensions of thebox may influence whether or not instabilities will occur withinthe available rise height. Fig. 4 illustrates the effect of the verticalboundary parallel to the cylinder axis between experiment 16I, 17Iand 18I. Reducing the domain of the box in the direction z per-
66 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 3. Three rising cylinders at infinite Rayleigh number. A perspective and sideviews are shown. Frames a)–c) show experiment 1I, experiment 6IPT and experi-ment 8IPT, respectively. Frames a) and b) show the influence on the wavelengths ofthe radius of the cylinder. Frames a) and c) show that the number of plumes doesnot depend on the viscosity ratioM.
pendicular to the cylinder axis allows the instabilities to grow ina shorter rising distance. In the first instance, this results from theeffect of the boundary on the rise velocity. Fig. 5 shows the heightof the cylinder as a function of time for experiment 16I, 17I an 18Ibefore the development of the instabilities. The velocity profiles arelinear and a linear regression on each curve indicates that the rise
velocity decreased by a factor 24 between experiment 16I and 17Iand, by a factor 49 between experiment 16I and 18I. Complemen-tary experiments with non-diffusing buoyant disks showed that therise velocity of the disk converges to a solution which is no longerinfluenced by the lateral boundaries when the distance of the diskcentre to the lateral boundary d satisfies d/a ≥ 40. From the com-parison of experiments 17I and 18I, it was difficult to further inferif there was also a change in growth rate as the vertical bound-aries of the box in the direction perpendicular to the cylinder axiswas made smaller. It was indeed not obvious to assess when theinstabilities were equally developed in experiment 18I comparedto experiment 17I due to asymmetry of the bump in experiment 17Iand the shift of the unstable wavelength in the box in experiment18I. A proper analysis of the growth rate with the influence of theboundary is to be found in the next section.
3.4. Timescale for growth
The analysis of the timescale for growth was made using seven-teen passive tracers that were initially equally distributed in spaceevery 0.25 non-dimensional unit length along the axis of the cylin-der. As the position of each tracer was recorded at each time step,we used the profiles of their positions to determine the growthrate. We detail in Appendix B the three methods which we used forthe analysis. The first method relies on a discrete Fourier transformof the tracer profiles and the growth of the dimensionless ampli-tude A = A/a of the dominant growing mode. The second method isbased on a measure of the dimensionless growing amplitude of thedisturbance A = A/a and its fit with an exponential curve. The thirdmethod estimates the slope of the natural logarithm of a correctedamplitude versus time. The amplitude is measured as in the secondapproach, which also provides the factor of correction.
Table 2 gives dimensionless timescales for the growth of fewexperiments at infinite Rayleigh number. The three approachesused to characterize the growth differ by less than 8% of error.First the estimates whenM = 1 have a mean value of 9.88, whichis consistent with the predicted value by Lister and Kerr (1989) of�G = 9.35. Second, the dimensionless timescale for growth is foundto decrease with an increase in the viscosity ratioM as predicted byEq. (1). Third, there is a consistent increase in growth rate �−1
G as Lz/ais smaller. We note that this increase in �−1
G appears to depend onthe viscosity ratioM. The increase in growth rate as Lz/a is reducedis found to be ≈17%, when M = 1, but it appears to be less than≈10%, whenM = 50, which is within the accuracy of our estimates.
4. Numerical experiments at finite Rayleigh number
4.1. Model setup and numerical method
The models were again set in a box of dimensions Lx, Ly, andLz in respectively the x, y and z directions. Model cylinders wereinitially laid out horizontally across the entire length of the box inthe x direction, centered in the box in the z direction and at a givenheight, which was typically toward the bottom of the box, in the ydirection. Free-slip velocity boundary conditions were applied onall walls.
These calculations were performed using a finite elementmultigrid solver and implemented in the StGermain/Underworldframework with the advection-diffusion solver based on Brooksand Hughes (1982). We did make use of a dual mesh approach toease the calculations. The method consisted of using two grids Nx,Ny, and Nz in respectively the x, y and z directions for the discretiza-tion of the differential momentum and energy equations. Given thatthe temperature field typically varies rapidly across the thermalboundary layer, a dense mesh is needed for this field to be accu-
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 67
Table 1Variables of the experiments at infinite Rayleigh number. The dimensionless wavelength �num = �num/a was estimated from the average of the wavelengths �num when morethan a wavelength could be seen. Otherwise �num = �num/a was based on the measured single wavelength. We note that there was no instability that could be observedwithin the given heights in experiments 10I and 16I. For experiments 11I, 13I, 19IPT and 20IPT only a bound for �num is given because if one instability was seen in thedomain, the length of the domain was too short to observe a complete wavelength. Experiments that were performed with passive tracers are indicated by the extra label PT.
rately represented. However, the momentum is fairly constant atthis scale which allow for a less dense sampling. Since the solutionof the momentum equations claims the majority of the memoryusage, we also reduced the number of unknowns in solving for thevelocity using a coarse mesh while the energy equation was solvedat the finer scale that the thermal boundary layer required usingvelocities values that were interpolated at the temperature meshpoints.
The root mean square (RMS) of the velocity Vrms was checkedto be identical between our single mesh and two mesh test case.Note that the size of the coarse grid was chosen to be twice thesize of the finer grid which was constrained by the length scale ofdiffusion (see Appendix C for further discussion). The visualizationagain used the gLucifer’s architecture.
For all calculations, the initial temperature of the thermallybuoyant fluid region was prescribed as a Gaussian tube with a
Fig. 4. Effects of the boundary on the growth rate. a) Experiment 16I, b) experiment 17I and c) experiment 18I. Experiments 16I, 17I and 18I have similar parameters andonly differ by the length Lz (see Table 1). The times shown in a)–c) are dimensionless times t = t/(�0/��0ga).
Author's personal copy
68 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 5. Velocity profiles of the vertical position of the cylinder as a function of timefor experiments 16I (dotted line), 17I (dashed line) and 18I (solid line). Note theincrease in velocity for experiment 18I, which corresponds to the rise of the unstableplume.
sinusoidal disturbance which follows:
T = a1 exp
(−r2yz
2c21
), (9)
ryz =√
(y − b1)2 + z2, (10)
b1 = a2 sin(b2x − c2), (11)
where a1, b1, c1, a2, b2 and c2 are constants to be chosen. Forexample, T can be smoothed without a sinusoid in choosing a2 = 0.c1 represents the initial smoothing length in temperature. Eq. (9)implies that the effective radius of the Gaussian cylindrical regionis then given by aeq =
√2c1. We discuss in Appendix C that such an
initial condition is equivalent to impose an already developed ther-mal boundary layer compared to a step in temperature. In doing so,the numerical grid needs only to resolve the thermal length c1.
During the calculations, a number of variables were recorded ateach time step. These variables include the total buoyancy B in thebox defined by
B = ˛�g
∫TdV, (12)
the vertical position of the centre of buoyancy yb by calculating
yb(t) = ˛�g
B
∫yTdV, (13)
Table 2Dimensionless timescale for growth �G measured at infinite Rayleigh number. �DFT
Gis the value derived from the Discrete Fourier Transform of the tracer profiles. �EFA
Gis the value determined from an Exponential Fit of the Amplitude versus time and�LLA
Gis the value derived from the Linear regression of the natural Logarithm of the
corrected Amplitude (see Eq. (21)) at times where the amplitudes A = (A/a) grewup to 1. For these runs, the correction ˛ was small with an average value of about0.0013.
Fig. 6. Growth of an initial disturbance from experiment 1FPT. Frame a) is the Gaus-sian tube initial condition. Frame b) is the Gaussian tube after a dimensionless timet = t/(�0/��0ga) = 46. The successive black dots are the passive tracers, which areused for the analysis of the growth. The single black dot in the right panel frame bindicates the position of the centre of buoyancy yb defined equation (13). The colorbar represents the non-dimensional temperature. (For interpretation of the refer-ences to color in this figure legend, the reader is referred to the web version of thearticle.)
where y is the vertical coordinate, the maximum temperatureTmax(t) in the box and the temperature of the centre of buoyancyTyb
(t) during the rise (Whittaker and Lister, 2008).
4.2. Growth of a disturbance
We captured the growth of an initial disturbance. At a Rayleighnumber of 56 in a box of dimensions 4 × 16 × 4, a Gaussian dis-turbance of length scale c1 = 0.2 (i.e. aeq = 0.28) with a1 = 1 anda sinusoid characterized by a2 = 0.05, b2 = 1.94 and c2 = /2 grew(experiment 1FPT, Table 3; Fig. 6). A priori the choice of wavelengthof the disturbance 2/b2 equal to 3.24 was expected to be within therange of wavelengths allowing growth based on the scaling foundin the experiments at infinite Rayleigh number. Using an aver-
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 69
Table 3Experiments at finite Rayleigh number. Passive tracers were used. The dimensionless timescale for growth �G with given dimensionless numbers. �DFT
Gis the value derived
from the Discrete Fourier Transform of the tracer profiles. �EFAG
is the value determined from an Exponential Fit of the Amplitude versus time and �LLAG
is the value derived fromthe Linear regression of the natural Logarithm of the corrected Amplitude (see Eq. (21)) at times where the amplitudes A = (A/a) grew up to 1. For these runs, the corrections˛ were about 10 times larger than those of the non-diffusing runs.
Experiment a Ra M �0/��ga c1 b2 Lx × Ly × Lz �DFTG
age of �num = �num/a = 10.21 (see Table 1), � = 10.21 × 0.28 ∼ 2.86.The calculation was done using the dual mesh approach withthe velocity and temperature grids of size 64 × 256 × 64 and128 × 512 × 128, and 96 × 384 × 96 and 192 × 768 × 192 on respec-tively 32 and 64 processors and 25 and 115 gigabytes of memoryusage. The sizes of the temperature grids implied that Li/Ni was
Fig. 7. a) Plot of the dimensionless amplitude A = A/a of the dominant Fourier modeas a function of dimensionless time from experiment 5FPT (Table 4). The dashed lineis the exponential fit, which has been determined from dimensionless amplitudesless than 1. The fit indicates a dimensionless timescale for growth of 41.70. b) Plotof the natural logarithm of dimensionless amplitude A + ¯ , where ¯ is the dimen-sionless factor of correction defined in Appendix B, as a function of dimensionlesstime from experiment 5FPT (Table 4). The value of ¯ determined via the exponentialfit on estimated amplitudes is 0.029082. The dashed line is to indicate the slope of0.0291 corresponding to the dimensionless growth rate, which is equivalent to adimensionless timescale for growth of 34.38. (For interpretation of the referencesto color in this figure legend, the reader is referred to the web version of the article.)
respectively 0.0312 and 0.0156 which was about 1/6 and 1/12 timesthe smoothing length c1. The two runs, which used a total walltimeof 24 h and 170 h, gave similar results of Vrms, B, yb, Tmax and Tyb
.As in some of the runs at infinite Rayleigh numbers, seventeen
passive tracers that were initially equally distributed every 0.25along the axis of the tube (see Fig. 6) were used to analyze thegrowth rate. The profiles of the tracer positions were here used tomeasure the amplitude of the initial disturbance as a function oftime.
Fig. 8. Plot of the dimensionless growth rate �DFTG
as a function of the wavenumberb2 of the initial sinusoid. a) WhenM = 1 and b) whenM = 50. (For interpretation ofthe references to color in this figure legend, the reader is referred to the web versionof the article.)
Author's personal copy
70 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 9. Rise and mergence of a disturbance. Experiment 15FPT at dimensionless times t = t/(�0/��0ga) = 0, 45.7, 91.3 and 136.9. The single black dot in the right panel foreach time indicates the position of the centre of buoyancy yb defined equation (13). The color bar represents the non-dimensional temperature. (For interpretation of thereferences to color in this figure legend, the reader is referred to the web version of the article.)
Fig. 7a shows the dimensionless amplitude A = A/a of the dom-inant Fourier mode as a function of dimensionless time t =t/(�0/��0ga) with its exponential fit, which has been deter-mined from dimensionless amplitudes less than 1. The fit indicatesa dimensionless timescale for growth of �DFT
G = 41.70. Fig. 7bshows the natural logarithm of the dimensionless correctedamplitude A + ¯ , where ¯ is the factor of correction defined inAppendix B, versus dimensionless time t = t/(�0/��0ga). Thegrowth rate of the logarithm of the dimensionless correctedamplitude tends to be linear, which reflects an exponential
growth of the instability. We estimate the exponential dimen-sionless growth rate while the amplitude grows from 0 to1. This corresponds to a dimensionless timescale for growth�LLA
G = 34.38.
4.3. Dependence of the timescale for growth with the wavelengthof the initial sinusoid
The dependence of the timescale for growth as a function ofthe wavenumber b2 of the initial sinusoid was examined through
Author's personal copy
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 71
Fig. 10. Rise of a disturbance at Ra = 600 and M = 50. a) Experiment 25FPT. b) Experiment 19FPT. Left panel is the temperature field. Right panel is the viscosity field. Thesingle black dot in the left panel for each time indicates the position of the centre of buoyancy yb defined equation (13). The wavenumber b2 allows growth in a) while it isdamped in b). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
experiments 2FPT–14FPT when M = 1, and experiments 19FPTto 29FPT when M = 50. Fig. 8a and b shows the dimensionlessgrowth rate �−1
G that was measured from the discrete Fourier trans-form approach as b2 varies. In both cases, a range of wavenumberspermits growth around a wavenumber, which leads to maximumgrowth. This wavenumber of maximum growth depends on theviscosity ratio M. In Fig. 8a, the growth rate peaks around b2 = 1.When b2 > 2 or < 1, the growth rate drops. In Fig. 8b, the growthrate peaks around b2 = 0.4. When b2 > 1 or < 0.5, the growth ratedrops. In experiments 15FPT and 20FPT, we did not estimate thegrowth rate because we observed a mergence rather than a growthof the initial disturbances. This is illustrated in Fig. 9 for experiment15FPT. Fig. 10 shows the growth and damping of the initial distur-bance when b2 = 0.45 (Fig. 9a) and b2 = 7.78 (Fig. 9b), respectively.For experiments 3FPT, 4FPT and 21 FPT, the Fourier decompositionof the disturbance was not reconstructing the profiles very well sowe did not measure the growth rate. We note that these cases arein the limit of the b2 region of favourable growth.
4.4. Dependence of the timescale for growth on dimensionlessnumbers
A few experiments (see Table 3) allow to further look at theinfluence of our dimensionless numbers on the dimensionlesstimescales for growth in addition to the benchmark cases at infiniteRayleigh number discussed earlier (see Table 2). Table 3 reportson the dimensionless timescales for growth �G that were foundaccording to the method of analysis. First, we find that the dimen-
sionless timescale for growth increases as the Rayleigh number Radecreases when M = 1. Second, the dimensionless timescale forgrowth decreases as the viscosity ratioM decreases. We find thatthe timescale growth rate of the diffusing cylinders when M = 1is about three times larger than the one which is obtained whenM ≥ 10. This actually extends the conclusions of Lister and Kerr(1989) for non-diffusing cylinders to diffusing cylinders.
4.5. Downstream thermal wake
A thermal wake developed during the rise of the diffusing cylin-ders. Fig. 11 shows the length of the wake at similar dimensionlesstimescales for four runs. Frames a) and b) of Fig. 11 show a cross sec-tion of the cylinders for experiments 5FPT and 16FPT at t = 54.68.Frames c) and d) of Fig. 11 show a cross section of the cylindersfor experiments 18FPT and 24FPT at t = 17.14. Importantly, whileexperiments 5FPT and 16FPT differed by their Rayleigh number,experiments 18FPT and 24FPT differed by their viscosity contrast.The length of the wake was measured as the distance between thezero of the centreline temperature in the tail and the position ofmaximum temperature on the axis. From experiments 5FPT and16FPT, we first found an increase by a factor 1.10 of the length ofthe wake with the increase of Rayleigh number by a factor 2.82.This increase is consistent with a recent scaling proposed for thelength of the wake lw , in which is lw ∝ (Ra ln Ra)1/2 (Whittakerand Lister, 2008). Using this scaling, a factor 1.88 of increase wasexpected for the increase of the wake length between Ra = 56 and
Author's personal copy
72 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 11. Cross sections of the cylinders showing the length of the thermal wake(vertical dashed line). Frames a) and b) are experiments 5FPT and 16FPT at t = 54.68.Frames c) and d) are experiments 18FPT and 24FPT at t = 17.14. The length of thewake was measured as the distance between the zero of the centreline temperaturein the tail and the position of maximum temperature on the axis. The single black dotin each panel indicates the position of the centre of buoyancy yb defined equation(13). The color bar represents the non-dimensional temperature. (For interpretationof the references to color in this figure legend, the reader is referred to the webversion of the article.)
Ra = 158. From experiments 18FPT and 24FPT, we then found thatthe change in viscosity has almost no influence on the wake lengthwhich differed by only 1% between M = 1 and M = 50. This out-come had only been suggested in the study of Whittaker and Lister(2008).
Some buoyancy appears to be lost by the cylinder. To analyzethe distribution of the buoyancy, we considered times when thecentre of buoyancy is a good proxy delimiting the buoyant cylinderand the trailing halo. The buoyancy in the wake was thus identifiedas the region of the box, whose vertical positions y were less thanthat of the centre of buoyancy yb and, where temperatures wereless than the temperature of the centre of buoyancy Tyb, as shownin Fig. 12. We thus found that 36% of the total buoyancy is located inthe wake region in experiment 1FPT of Ra = 56 andM = 1 at t = 46(Fig. 6). At Ra = 600 andM = 50 (experiment 25FPT, Fig. 10a), 36%and 37% of the total buoyancy were found in the wake region atrespectively at t = 9.02 and t = 16.05.
Fig. 12. Temperature distribution for experiment 1FPT at t = 46. Black dots are thevertical positions y of a sampled grid cell as a function of temperature in the cell. Thetwo red lines are the ordinates of the centre of buoyancy yb and its temperature Tyb
.The wake region is defined by y < yb and T < Tyb
. (For interpretation of the referencesto color in this figure legend, the reader is referred to the web version of the article.)
5. Growth rates from laboratory experiments
In Kerr et al. (2008), the time for gravitational instability wasquantified as a function of Rayleigh number (see Fig. 4 of Kerr et al.,2008). It was measured, from the front view, as the time ti at whichthe amplitude of the gravitational instabilities had grown to beequal to the cylinder diameter. Here, we looked at the growth rateof the instabilities as we analyzed the initial exponential growth ofthe amplitudes when amplitudes are smaller or equal to the radiusof the cylinder. Our methodology followed several steps. First weextracted frames every second from our front view videos. Second,each image of the front views was scaled to have the same verti-cal and horizontal scales. A minor correction of the horizontal (lessthan 1◦) was applied to the images in some of our experiments.We then used the software LoggerPro to digitize the profiles of thecylinders along which we added tracers as shown in Fig. 13. All the
Fig. 13. Eleven digitized profiles obtained for experiment 5 of Kerr et al. (2008, seeTable 4) at a Rayleigh number of 343. The profiles are shown from the bottom tothe top at times t = 50, 53, 56, 59, 62, 65, 68, 71, 74, 77 and 80 s. The horizontal axisis the position x along the tank and the vertical axis is the vertical position y of thecylinder at time t.
Author's personal copy
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 73
Fig. 14. Growth analysis of experiment of Kerr et al. (2008, see Table 4) at a Rayleighnumber of 343. a) Plot of the dimensionless amplitude A = A/a, defined as the stan-dard deviation of the y position of the profiles between x = 12 and 25 cm seen inFig. 10, as a function of dimensionless time. The dashed line is the exponential fit,which is determined from dimensionless amplitudes less than 1. The fit indicatesa dimensionless timescale for growth of 7.94. b) Plot of the natural logarithm ofthe dimensionless corrected amplitude A + ¯ as a function of dimensionless timet. ¯ is the dimensionless factor of correction defined in Appendix B. The standarddeviation of the y position of the profiles between x = 12 and 25 cm seen in Fig. 10was again used as a proxy for the amplitude. The last data points are not used forthe regression as the growth is not likely to represent the initial exponential growthbut rather the linear growth which is expected as amplitudes increase. The dashedline is to indicate the slope of 0.12 corresponding to the dimensionless growth rate,which is equivalent to a dimensionless time scale for growth of 8.35. (For interpre-tation of the references to color in this figure legend, the reader is referred to theweb version of the article.)
profiles were ultimately analyzed to estimate the amplitude of thegrowing disturbances while the dimensionless amplitude A = A/awas less than 1.
The growth rate was estimated in looking for the exponen-tial fit of the dimensionless amplitude A = A/a < 1 as a function ofdimensionless time t = t/(�0/��0ga). A = A/a < 1 was defined as thestandard deviation of the y position of the profiles for one plume.For instance the plume selected for Fig. 14a was the one betweenx = 12 and 25 cm as seen in Fig. 10. In this example, the fit indicatesa dimensionless timescale for growth �G of 7.94. We note that theorigin of time was chosen as the time at which the cylinder wasemplaced at x = 25 cm in the tank. In Fig. 14b, we also show the linearregression on the natural logarithm of the dimensionless correctedamplitude A + ¯ versus the dimensionless time t = t/(�0/��0ga).The constant ¯ is the factor of correction defined in Appendix B.Errors on the growth rate reported Table 4 were estimated bysampling different plumes along the cylinder. Fig. 15 shows thedimensionless characteristic time for growth (i.e. the inverse of thegrowth rate) �G that was estimated for 17 of the experiments byKerr et al. (2008). The results of Table 4 primarily show that thegrowth rates �−1
G increase with the Rayleigh number.
6. Discussion
There is qualitative agreement between the growth rate derivedfrom the numerical experiments and the laboratory study of Kerret al. (2008). For instance, the numerical values of �EFA
G appears to be5% larger than an interpolated value �EFA
G at Ra = 600 from the lab-oratory experiments at Ra = 620. The difference is within the errorbar of ≈15% for the estimate of the experimental value.
Table 4Dimensionless characteristic times for growth �G of experiments by Kerr et al. (2008)based on the value �EFA
Gdetermined from an Exponential Fit of the Amplitude versus
time. The error was estimated by sampling a different plume.
To some degree, the errors in our estimates hide the differences,which could have been expected between the numerical and lab-oratory experiments. For instance, a difference could be attributedto the aspect ratio. The aspect ratio Lz/a of the laboratory tankwas Lz/a = 31.58–66.67, which is 3–6 times larger than ours. Soour current estimates of the characteristic time for growth shouldbe smaller but the difference should be marginal because we onlyfound 10% increase in growth rate at most in the infinite Rayleighnumber experiments as the aspect ratio Lz/a was made 2 timessmaller whenM = 50.
Another difference between our numerical diffusing experi-ments and the laboratory experiments is the initial tilt. The initialtilt was small but varied from 0.2 to 4.1◦ in Kerr et al. (2008), whilein the present study, there was no initial tilt. Both compositionalexperiments on the rise of horizontal cylinders (Kerr and Lister,1988) and tilted ones (Whitehead, 1982) led to chains of unstableplumes which does not suggest a dramatic change of behaviour dueto tilt. Besides, where the cylinder is horizontal, it rises verticallythrough the host fluid as the buoyancy force is simply perpendic-
0 200 400 600 800 1000 12000
2
4
6
8
10
12
14
16
18
20
τ G
Ra
Limit at Ra =∞
Fig. 15. Plot of the characteristic dimensionless time for growth �EFAG
as a functionof Rayleigh number Ra. Note the increase of the dimensionless characteristic timefor growth (i.e. the decrease of the growth rate) as Rayleigh numbers are less than≈200. The black dots are from our analysis of the experiments of Kerr et al. (2008).The data at Ra = 600 represented by a triangle symbol is a data obtained with thenumerical experiments. The dashed line represents the limit found for �DFT
Gat Ra = ∞
andM = 50 in experiment 5IPT.
Author's personal copy
74 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 16. Stand alone cylinders at t = 21.89. a) The cylinder was initially horizontal. b) The cylinder was initially tilted 2.5◦ to the horizontal.
ular to the axis of the cylinder. This vertical rise of the buoyantcylinder through the host fluid leads to Rayleigh–Taylor instability.As the cylinder is inclined to the horizontal increases, the buoy-ancy force splits into two components: a diminished componentof the buoyancy force perpendicular to the axis of the cylinderand a new component of the buoyancy force causing the fluidto flow rather along the cylinder axis. So one would expect thatthe horizontal rise is the most favourable case to the formationof creeping plumes by gravitational instability. Furthermore theflow along the cylinder can be described by the pipe flow model.Such a model predicts that a volumetric flux of viscous fluid Qwould flow through a more viscous fluid in a vertical conduit ofradius of a = (8iQ/��g)1/4,where i is the viscosity of the fluidinside the conduit. When inclined of angle � to the horizontal, thesame given flux Q would flow along a conduit of increased radiusa = (8iQ/��g sin �)1/4 as the component of the gravity force isg sin �. Consequently, we could speculate that, when the cylinder istilted from the horizontal, it is still unstable because of the compo-nent of the buoyancy force perpendicular to the axis of the cylinder.However the buoyancy force along the cylinder axis might allowthat small disturbances migrate toward each other due to the flowalong the axis and feed fewer major plumes.
When the cylinder is tilted, the flow along the axis causesan inflow/outflow at the lateral walls, which corresponds to aboundary condition that is no longer compatible with the free-slipboundary conditions used in this study. So to get an insight intothe behaviour of tilted cylinders, we performed one simulation atRa = ∞ andM = 50 of a long stand-alone cylinder, which was tilted� = 2.5◦ to the horizontal. For comparison, we also simulated theequivalent flat one. The cylinders of radius a = 0.2 and length L = 44were laid in a box with Lx = 48. The cylinders were made long to limitthe effects of the free ends of the cylinder. Fig. 16 shows the ini-tial growth of the flat and tilted stand-alone cylinders at t = 21.89.Apart from the effects of the free ends of the cylinder, the behaviouris very similar: both shows the development of a series of creep-ing plumes. However, the experiments suggest that the plumes aremore frequent and evolved in the flat case than in the tilted case.It is interesting to note that the instabilities of both free ends growfaster in the tilted case.
Beyond the growth rate, this study shows a consistent 35% lossof buoyancy in the wake at times when the wake is well devel-oped. To some degree this estimate results from our initial thermal
condition, which imposes an already developed thermal boundarylayer. When no thermal boundary layer is initially prescribed (e.g.condition 2 of Appendix C), the loss is not as large. An additionalexperiment at Ra = 158, a = 1 andM = 1 with a step in temperaturefrom 1 to 0 as initial condition gave 15% of loss at t = 14 and 29%at t = 29. By comparison, in experiment 16FPT of Rayleigh numberRa = 158, 35% of the total buoyancy was found in the wake regionat t = 26. Even though our estimates are apparently larger thanthe 7% loss of buoyancy estimated t = 60 for the rise of a cylin-der of Rayleigh number Ra = 80 in an infinite medium (Kerr et al.,2008). The difference is likely due to the very different numeri-cal schemes used. In the early study, the diffusion of buoyancy isrepresented by a Brownian motion, while in this study it is solvedusing an explicit advection-diffusion solver (Brooks and Hughes,1982). In this regard, we performed 2D simulations to check thatthe diffusion was sufficiently resolved to give reliable estimates ofthe buoyancy loss. These 2D test runs were built to represent thevertical cross section of the cylinder with the initial temperaturestep but at successive finer resolutions. Table 5 details the buoy-ancy estimates in the 3D experiment and, in the 2D experiments.The estimates of the loss of buoyancy for the 3D and 2D do not differby more than 5%.
In the future, it would be useful to further improve on theperformance of the code, as the calculations are currently ratherexpensive. In particular, one improvement of the Underworldframework, which would help is the ability to use a structured gridwith non-uniform element sizes in the simplest approach by hav-ing the mesh at the resolution of the thermal boundary layer byhaving the outer region of coarser resolution. With a lesser cost,experiments could be considered more easily.
Table 5Measure of buoyancy at dimensionless time t = 14 for a cylinder at Ra = 158, a = 1andM = 1 with a step in temperature as initial condition and, the 2D vertical crosssection of the cylinder at successive finer resolutions. The cylinder was in a box oflength 4. Hence the factor 4 in buoyancy between the 3D and 2D case. The resolutionis Li/Ni . The total buoyancy in the domain is B while Bw is the buoyancy in the wakedefined by T < 0.45 and y < 3.43.
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 75
Table 6Mantle plume timescales: (1) time estimates for instability growth �G = c�0/��ga,where the factor c is estimated from Fig. 15 and Table 4, and (2) rise time throughthe mantle tr . ��g is set to a typical value of 300 kgm−2 s−2. The thermal diffusivity� is set to a value of 7.5 × 10−7 m2 s−1 in the upper mantle and 1.5 × 10−6 m2 s−1 inthe lower mantle (Steinberger and Antretter, 2006). The rise time tr is calculatedusing the Stokes velocity vs = ks��ga2/�0 with ks = 0.13(�0/�p)0.62, where �p is theplume viscosity (Kerr and Lister, 2008) and, a length scale H of 600 and 2100 km forthe upper and lower mantle respectively. We used �0/�p = 10 in the upper mantleand �0/�p = 100 in the lower mantle (Steinberger and Antretter, 2006).
One conclusion of Kerr et al. (2008) is that strongly tilted plumeswere stable because of thermal diffusion. This conclusion wasessentially drawn from estimates of the Rayleigh numbers for man-tle plumes, which were less than 140. As Ra ∝ a3, values for Racritically depend on the values used for the plume radii a, whichcan be simply derived from the pipe flow model using the buoyancyflux Bf attributed to plumes. In doing so, Ra ∝ B3/4
f. This approach
provides estimates of plume radii at the source (i.e. the core-mantleboundary) before a thermal halo caused by diffusion and thermalentrainment develops around the conduit. So estimates of plumeradii (≤52 km) are about twice as less as those from tomographicmodels (≈100 km) in the upper mantle (Montelli et al., 2004;Wolfe et al., 2009). Whereas a direct measure of plume radii at thecore-mantle boundary is not accessible, buoyancy fluxes of mantleplumes are known up to some degree. While Kerr et al. (2008) used avalue of 6200 kg s−1 for Hawaii (Davies, 1988), a value of 8700 kg s−1
(Sleep, 1990) has also been proposed. This latter value would beequivalent to a factor 1.28 increase in terms of plume Rayleighnumber. More, plume heat fluxes for individual plumes at mid-mantle depths were found to be are greater than those predictedfrom the buoyancy flux at the surface within a factor 2–6 (Noletet al., 2006). Besides some of the mantle plume buoyancy is lost dur-ing ascent to form a downstream thermal wake (Kerr et al., 2008).If so, buoyancy fluxes would be underestimated when measuredfrom surface topographic swells. In other words plume Rayleighnumbers could be larger. The overall incertitude in the buoyancyflux could result in unstable mantle plumes. From this study wecan actually disregard this possibility and draw the conclusionthat strongly tilted mantle plumes are to be stable by examiningtheir characteristic times. The growth timescales for plumes in theupper or lower mantle and with radius 50 and 100 km togetherwith the plume rise timescales through the upper and lower man-tle are given in Table 6. When compared, the growth timescalesare larger than the rise timescales, which clearly eliminates thepossibility of unstable plumes within the mantle. One exceptionis the case of plumes that would have a diameter of 100 km inthe upper mantle. Such plumes would likely develop instabilitiesin a rise height H of about 4 times the plume radius (Kerr et al.,2008), whereas the available rise height in the upper mantle wouldbe 6 times the plume radius. The buoyancy flux that would feedsuch a plume would be 9 × 103 kg s−1. This is larger than current
estimates of buoyancy flux for Hawaii, which has the largest buoy-ancy by in large for all plumes. A buoyancy flux of 9 × 103 kg s−1
would only be realistic if it represented the Hawaiian plume withapparent surface buoyancy of 6300 kg s−1 but, which would havelost 30% of its source buoyancy. Such scenario seems very unlikelygiven the absence of geophysical observations that could suggestthe presence of secondary instabilities under Hawaii. In particu-lar, unlike the Louisville hotspot, the age progression along theHawaiian volcanic chain is well described by the Pacific platemotion.
7.2. Influence of the tilt
While long-lived mantle plume tails are sheared by surround-ing mantle flow, the degree of tilt of the plume tails varies both inspace and time. For instance, Table 2 of Steinberger (2000) showthat tilts mantle plume conduit are often substantial with maxi-mum tilts ranging from 3 to 94◦. In this context, horizontal plumeconduits are thus strictly speaking a limit case. We have seen thathorizontal plume conduits are not unstable although they are themost in favour of instability compared to tilted conduits whosecomponent of buoyancy force responsible for instability is dimin-ished. So weakly tilted plume conduits should be even less likely tobecome unstable provided they do not encounter a discontinuityin the mantle.
7.3. Thermal wake
The development of a thermal wake and the loss of buoy-ancy from mantle plume tails essentially depend on the rise timeavailable. For instance, we showed that a strongly tilted plume ofRa = 158 could have lost 30% of its buoyancy after a dimensionlessrise time of t = 30. By comparison, characteristic dimensionlessrise times for the whole mantle would range between 20 and 40(Table 6). This suggests that mantle plumes could lose some buoy-ancy during the rise implying that estimates of plume radius andRayleigh number would be underestimated. A 30% loss of buoyancywould correspond to a plume radius and plume Rayleigh numberunderestimated by about 8% and 22% respectively. However sucha shift in Rayleigh number would not be sufficient to allow thegrowth of instabilities.
7.4. Large-scale mantle convection simulations
Mantles plume tails have a characteristic thermal boundarylayer that scales as ı ∝ aRa−1//2, where Ra is the plume Rayleighnumber. As shown in Appendix C, this characteristic length scale isactually quite constraining in terms of resolution for any numer-ical simulations. In particular the criterion is more constrainingthan the condition which is usually used in the large-scale con-vection models. In such models the finest resolution is chosenaccording to the thickness of the boundary layer which scales as� ∝ HRa−1/3
c , where H is the height of the box and Rac is definedby Rac = ˛g��TH3/��0. Note that we use here for simplicity thepower 1/3 but this scaling was shown to be more subtle by Zhong(2005). At Rac = 106 within a box of H = 1, � would thus need to beof order 0.01. The thermal plumes developed by such a convect-ing system would be characterized by a plume Rayleigh number ofRa = (a/H)3Rac ≈ 296, assuming a ratio of about 1/15 between theradius of a plume and the depth of the mantle. The plume wouldthus be resolved when the resolution is ı = 1/15 × 296−0.5 ≈ 0.004that is 2.5 times smaller than �. This severe requirement wouldget worse for even larger Rayleigh number convection models dueto the power law in −1/2 scaling the plume thermal boundarylayer. Therefore, if they generate plumes, typical convection models
Author's personal copy
76 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
have yet to resolve the growth by combined diffusion and thermalentrainment of mantle plumes.
8. Conclusions
In this paper, we have used laboratory and numerical exper-iments to investigate the rise and instabilities developed by aninitially horizontal, buoyant cylinder of fluid rising through a denserfluid at low Reynolds number. We have then examined the impli-cations for mantle plumes.
The numerical study shows that diffusing cylinders have a finiterange of wavenumbers for which there is growth of a thermaldisturbance, with a maximum that demonstrates the presence ofa most unstable wavenumber. Based on both the laboratory andnumerical experiments, we show that the growth rate increases asthe Rayleigh number increases. That growth rate increases withthe viscosity ratio both in the case of non-diffusing and diffus-ing cylinders. Numerical experiments for non-diffusing cylindersreveal that a greater proximity of the cylinder to the vertical bound-aries parallel to his axis results in a decrease of the rise velocity buta slight increase of the growth rate. The rise of the diffusing cylin-ders is associated with the development of a thermal wake whoselength increases with the Rayleigh number but does not dependon the viscosity ratio. Applied to the timescales of the mantle, thisstudy shows that mantle plumes cannot be unstable. The study alsosuggests that mantle plumes would not lose significant buoyancyduring the rise. Finally, the resolution that is required to allow forthe growth of mantle plume tails by combined diffusion and ther-mal entrainment is finally shown to represent a challenge for thelarge scale mantle convection simulations.
Acknowledgements
CM especially thanks editor Mark Jellinek for handling thispaper and one anonymous reviewer for his very sound com-ments. The calculations presented here were conducted on theSGI Altix AC and SGI XE clusters at the NCI National Facility(http://nf.nci.org.au/). We would like to thank David Singleton andMargaret Khan from the NCI National Facility for their support. CMwas supported by the Merit Allocation Scheme, which the NationalComputational Infrastructure at the Australian National Univer-sity has established to support computationally intensive projects(http://www.nci.org.au/). CM acknowledges financial support fromLM’s ARC grant DP0449979.
Appendix A. Dimensional analysis in a finite domain
Consistent with dimensional analysis, in our finite domain, therise velocity u is given by a functional relation of the form
u = UF(
Ra,M,Li
a
)= g��0a2
�0F
(Ra,M,
Li
a
), (14)
where Li refers to the domain length in each spatial direction i andF is an unknown function. Similarly the growth rate � = 1/� follows
� = ��0ga
�0G(
Ra,M,Li
a
), (15)
where G is another unknown function. So the distance of rise hnecessary to the development of instabilities satisfies
h = u� = aF(
Ra,M,Li
a
)G−1
(Ra,M,
Li
a
)= aH
(Ra,M,
Li
a
).(16)
Appendix B. Methods of analysis of the growth rate
The instability of the cylinder interface classifies as aRayleigh–Taylor instability under constant gravitational accelera-
tion (Taylor, 1950; Chandrasekhar, 1961; Kull, 1991). The unstablemodes are usually assumed to arise from an initial surface displace-ment 0, while the fluid is initially at rest, or from an initial surfacevelocity ∂t 0. The radius of the cylinder can then be consideredas the combination of its mean value plus the perturbation. It iswritten
r = a + = a + Af (x, y, z), (17)
where A is a function of time and represents the amplitude of theperturbation and, f(x, y, z) is the spatial disturbance (e.g. the Gaus-sian tube of the diffusive numerical experiments). The fluid velocityat the interface in the vertical direction is then
v = a + Af (x, y, z). (18)
One may show then that the amplitude of the perturbation satisfiessolutions whose dependence on time t is exponential.
In the numerical experiments, we analyze the growth from thepassive tracers and, although a disturbance of the cylinder inter-face may initially be applied, the initial disturbance of the tracers isactually zero. In addition, v(0) is not zero. We can therefore assumethe vertical velocity of the disturbance decomposed into its Fouriercomponents k is given by
A(t) =∑
k
Ak(0) exp(
t
�G(k)
), (19)
and the vertical displacement of the disturbance
A(t) =∑
k
Ak(0)�G(k)(
exp(
t
�G(k)
)− 1
), (20)
or, for a single mode k,
ln(
Ak(t) + Ak(0)�G(k))
= t
�G(k)+ ln
(Ak(0)�G(k)
). (21)
In this paper, we found that instability actually resulted from adominant growing mode k. For simplicity, �G throughout this paperis thus implicitly the characteristic time for growth of the dominantmode.
The analysis of the timescale for growth of the numerical exper-iments is made using passive tracers that are initially equallydistributed along the axis of the cylinder. As the position of eachtracer is recorded at each time step, the successive profiles of theirvertical positions represent the function A(t) (see Fig. 17). Two inde-pendent approaches can then be adopted to estimate the growthrate �−1
G in the range where amplitudes are smaller than the radiusof the cylinder a. In approach 1, a Discrete Fourier Transform (DFT)of the tracer profiles is performed which provides the differentmodes of displacement and their associated amplitudes as a func-tion of time. A fit of the dominant amplitude versus time witha function of the form ˛(exp (�t) − 1) then gives the growth rate� . In approach 2, one assumes that it is the most unstable wave-length that grows and dominates over the other wavelengths. Theamplitude of growth A(t) can then be estimated as being half thedifference between the maximum and the minimum of the verticalpositions of the tracers in each profile. A subsequent fit of A(t) witha function of the form ˛(exp (�t) − 1) gives the growth rate � andthe constant ˛. A third approach can also be applied, but this onedepends on approach 2. In approach 3, the amplitude of growth isestimated as in approach 2 but one plots the natural logarithm ofa corrected amplitude (see Eq. (21)) versus time using the valueof ˛ = A(0)�G determined in approach 2. The growth rate of thelogarithm of the amplitude tends to be linear, which reflects theexponential growth of the instability. A fit of that slope then givesthe growth rate.
In this paper, the analysis of the growth rate is made in dimen-sionless form with length scale a and timescale �0/��0ga. Fig. 18
Author's personal copy
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 77
Fig. 17. Three profiles of the 17 tracers in experiment 14IPT. Tracers are shown bycrosses at their horizontal and vertical positions x and y in the box. From bottom totop, the three profiles that are shown corresponds to the profiles at dimensionlessmeasured amplitudes A = A/a = 0.4, 1 and 2 respectively.
is an example of the DFT procedure applied to experiment 14IPT.Fig. 18a shows the tracer profiles at three times and those pro-files, which we reconstructed from the Fourier decomposition. Thepower spectrum of the tracer profiles, which is illustrated Fig. 18b,shows that growth occurs for a dominant wavenumber. Fig. 18cshows the dimensionless amplitude A = A/a of this dominant modeversus dimensionless time t = t/(�0/��0ga) and the exponentialfit found for 0 ≤ A < 1. We estimate a dimensionless growth rate �G
of 0.097.Fig. 19 shows the results of approach 2 and 3 based on the ampli-
tude proxy. Fig. 19a shows the dimensionless amplitude versusdimensionless time and its associated exponential fit for 0 ≤ A < 1.We estimate a dimensionless growth rate of 0.095 in this case and
Fig. 18. Discrete Fourier Transform analysis of the growth rate applied to experi-ment 14IPT. a) The tracer profiles at three times (continuous line) and those profiles,which we reconstructed from the Fourier decomposition (crossed line). b) The powerspectrum of the tracer profiles, which shows that growth occurs for a dominantwavenumber. c) The dimensionless amplitude of the dominant mode versus dimen-sionless time (continuous line) and the exponential fit found for 0 ≤ A < 1 (dottedline). We estimate a dimensionless growth rate of 0.097. (For interpretation of thereferences to color in this figure legend, the reader is referred to the web version ofthe article.)
Fig. 19. a) Plot of the dimensionless amplitude A = A/a versus dimensionless timefrom experiment 14IPT. The dashed line is to indicate its associated exponential fitfor 0 ≤ A < 1. We estimate a dimensionless growth rate of 0.095 and ¯ = 0.001171.b) Plot of the natural logarithm of dimensionless amplitude A + ¯ as a function ofdimensionless time from experiment 14IPT (Table 1). The dashed line is to indicatethe slope of 0.091 corresponding to the dimensionless growth rate. (For interpreta-tion of the references to color in this figure legend, the reader is referred to the webversion of the article.)
find that the dimensionless value of ˛, ¯ = ˙A(0)�G , is 0.001171. Wenote that the coefficients of fits were systematically determinedwith 95% confidence bounds. Fig. 19b shows the natural logarithmof the amplitude plus ¯ = 0.001171 versus time and the linearregression found for A < 1, which gives a dimensionless growth rateof 0.091.
The DFT based method is by far the most complete approach toestimate the growth rate. It is rather a straight forward analysis inthe context of the numerical experiments in which tracers are ini-tially equally distributed along the unstable interface but it is not sofrom the tracer profiles extracted from the laboratory experiments.In this later case, an additional step would require the interpolationof the profiles before they can be sampled at equally spaced inter-vals. In this paper, only approach 2 was applied to the laboratorydata of Kerr et al. (2008).
In the context of the laboratory experiment, it could seem betterto suppose that the fluid is initially at rest, in which case the solutionfor A is of the form:
A(t) ∝(
cosh(
t
�G
)− 1
). (22)
We actually found that, when the amplitude is represented by Eq.(22), the growth rates derived for the laboratory experiments arenot very different from those derived using Eq. (20). Fits can befound either using Eqs. (20) or (22) and, respectively, the two esti-mates of �G appear to differ by a minimum of 2% to a maximumof 26%. Only for the Ra = 170 and Ra = 158 experiments, �G derivedfrom using Eq. (22) is found to be about twice the value given inTable 4. Whatever the case may be, the two sets of estimates showa very similar curve of �G versus Ra. In particular the growth ratedecreases as the Rayleigh number decreases.
Appendix C. Thermal boundary layer
When the initial temperature of the thermally buoyant fluidregion is prescribed as a Gaussian tube, the initial disturbance intemperature decreases smoothly from the axis of the cylinder inevery radial directions within a distance which is given by the
Author's personal copy
78 C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79
Fig. 20. Two temperature profiles for each experiment TBL1–TBL4 as referred in Table 7. The temperature profiles are shown within the first time steps as Tmax ∼ 1, and afterTmax has dropped to 0.5. The distances ı at which the extrapolation of the linear profile (dashed lines) of the temperature profile equals the central mean of temperature isindicated in each frame. The crosses indicate the temperature recorded at the nodes of the grid.
effective radius aeq =√
2c1. This condition (condition 1) effectivelymimics a thermal boundary layer (TBL).
Another approach could be applied to set the initial condition ofthe thermally buoyant cylindrical region. A temperature T = 1 canbe simply prescribed in a cylindrical fluid region of radius a (seeFig. 1) while the outer region is at T = 0. When using this later initialcondition (condition 2), the TBL is a step in temperature from 1 to0. The values of Nx, Ny, and Nz have then to be chosen to ensurethat the thermal halo, which will develop around the cylinder dueto diffusion, is resolved. This requirement turns out to be very con-straining. For large Rayleigh numbers, the thickness ı of the TBL fora thermally buoyant drop in an unbounded region is characterizedby
ı = k1
(�a
U
)1/2= k1aRa−1/2, (23)
where k1 is a constant to be determined (Clift et al., 1978; Griffiths,1986; Griffiths and Campbell, 1991). As ı is characterized by Eq.(23), one must satisfy in each i direction
Li
Ni= c3ı = c3k1aRa−1/2, (24)
where the constants k1 and c3 need to be determined.To estimate k1c3 we performed four experiments of increasing
grid size labelled TBL1–TBL4 at a given Rayleigh (see Table 7). Fig. 20shows 2 temperature profiles for each experiment TBL1–TBL4.The temperature profiles are shown within the first time steps asTmax ∼ 1, and after Tmax has dropped to 0.5. The thickness of the
Table 7Numerical experiments for the scaling of the thermal boundary layer.
Experiment a Ra aRa−1/2 Lx × Ly × Lz Nx × Ny × Nz Li/Ni
TBL was then defined as the distance ı at which the extrapolationof the linear profile of the temperature profile equals the centralmean of temperature. Fig. 20 shows that (1) the thermal bound-ary grows with time, which implies that it is better resolved withtime (2) only for experiment TBL4, the TBL is initially resolved bytwo and a half cells. Therefore we can infer that k1c3 ∼ 1 and oneneeds to satisfy Li/Ni ≈ aRa−1/2 in each i direction when the initialstep condition is used. However this condition is a very difficultone to satisfy in regard to the numerical resources. For instancea thermally buoyant cylinder of Rayleigh number 200 and radius0.2 in a region of dimensions Lx × Ly × Lz = 4 × 8 × 2 would requireNx × Ny × Nz = 284 × 568 × 142 grid cells.
Incidentally, using two runs 1 and 2 of similar ��, g and Lx
but of initial condition 1 and 2 respectively, one can calculatea2 = a1(B2/B1)1/2, where B1 and B2 are the buoyancy of run 1 and2. As an example, we measured B2 = 1.63865 with a Gaussian tubewithout sinusoid but c1 = 0.2, and B1 was 3.28547 for a cylinderfluid of radius a1 = 0.4. This gives a1(B2/B1)1/2 = 0.2825 which showsconsistency with the scaling aeq =
√2c1 = 0.283.
References
Brooks, A.N., Hughes, T.J.R., 1982. Streamline upwind /Petrov-Galerkin formulationsfor convection dominated flows with particular emphasis on the incom-pressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. 32,199.
Chandrasekhar, S., 1961. Hydrodynamics and Hydromagnetic Stability. Oxford atThe Clarendon Press, p. 429.
Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops and Particles. Academic Press,New York, p. 246.
Davies, G.F., 1988. Ocean bathymetry and mantle convection: 1. Large-scale flowand hotspots. J. Geophys. Res. 93, 10,467.
Duboz, C., Turnbull, R., Stegman, D., Moresi, L., Lo, A., Quenette, S., 2005. gLucifer:Next Generation Visualization Framework, APAC05.
Griffiths, R.W., 1986. Thermals in extremely viscous fluids, including the effects oftemperature-dependent viscosity. J. Fluid Mech. 166, 115.
Griffiths, R.W., Campbell, I.H., 1991. On the dynamics of long-lived plume conduitsin the convecting mantle. Earth Planet. Sci. Lett. 103, 214.
Kerr, R.C., Lister, J.R., 1988. Island arc and mid-ocean ridge volcanism, modelled bydiapirism from linear source regions. Earth Planet. Sci. Lett. 88, 143.
Author's personal copy
C.A. Mériaux et al. / Physics of the Earth and Planetary Interiors 184 (2011) 63–79 79
Kerr, R.C., Mériaux, C., Lister, J.R., 2008. Effect of thermal diffusion on the sta-bility of strongly tilted mantle plume tails. J. Geophys. Res. 113, B09401,doi:10.1029/2007JB005510.
Kerr, R.C., Lister, J.R., 2008. Rise and deflection of mantle plume tails. Geochem.Geophys. Geosyst. 9, Q10004, doi:10.1029/2008GC002124.
Lister, J.R., Kerr, R.C., 1989. The effect of geometry on the gravitational instability ofa buoyant region of viscous fluid. J. Fluid Mech. 202, 577.
Kull, H.J., 1991. Theory of the Rayleigh–Taylor instability. Phys. Rep. 206, 197.Montelli, R., Nolet, G., Dahlen, F.A., Masters, G., Engdahl, E.R., Hung, S., 2004. Finite-
frequency tomography reveals a variety of plumes in the mantle. Science 303,338.
Moresi, L.N., Dufour, F., Mühlhaus, H.B., 2003. A Lagrangian integration point finiteelement method for large deformation modeling of viscoelastic geomaterials. J.Comput. Phys. 184, 476.
Moresi, L.N., Quenette, S., Lemiale, V., Mériaux, C., Appelbe, B., Mühlhaus, H.B., 2007.Computational approaches to studying non-linear dynamics of the crust andmantle. Phys. Earth Planet. Int. 163, 69.
Nolet, G., Karato, S., Montelli, R., 2006. Plumes fluxes from seismic tomography. EarthPlanet. Sci. Lett. 248, 685.
tion and magmatic resurfacing on Venus. Icarus 139, 67.Skilbeck, J.N., Whitehead, J.A., 1978. Formation of discrete islands in linear island
chains. Nature 272, 499.
Sleep, N.H., 1990. Hotspots and mantle plumes: some phenomenology. J. Geophys.Res. 95, 6715.
Steinberger, B., OConnell, R.J., 1998. Advection of plumes in mantle flow: implica-tions for hotspot motion, mantle viscosity and plume distribution. Geophys. J.Int. 132, 412.
Steinberger, B., 2000. Plumes in a convecting mantle: models and observations forindividual hotspots. J. Geophys. Res. 105, 11, 127.
Steinberger, B., Antretter, M., 2006. Conduit diameter and buoyant rising speed ofmantle plumes: implications for the motion of hot spots and shape of plumeconduits. Geochem. Geophys. Geosyst. 7, Q11018, doi:10.1029/2006GC001409.
Taylor, G., 1950. The instability of liquid surfaces when accelerated in a direc-tion perpendicular to their planes. Proc. Roy. Soc. Lond. Math. Phys. Sci. 201,192.
Whitehead, J.A., 1982. Instabilities of fluid conduits in a flowing Earth—areplates lubricated by the asthenosphere? Geophys. J. Roy. Astron. Soc. 70,415.
Whittaker, R.J., Lister, J.R., 2008. The self-similar rise of a buoyant thermal in veryviscous flow. J. Fluid Mech. 606, 295–324.
