26. A. Limare, B. Kenda, E. Kaminski, E. Surducan, V. Surducan, C. Neamtu, Transient convection experiments in internally-heated systems, submitted to MethodsX, 2019.

Brevets

1. A. Limare, E Surducan, V Surducan, C Neamtu and E di Giuseppe, Microwave heating device, RO131921 A2, EP2947961 A1, WO2015177244 A1.

2. V Surducan, E Surducan, A. Limare, Stabilizing and controlling block meant to supply magnetron filament current, RO130089 A2.

Communications dédiées au grand public

1. Modélisation expérimentale de rivières. Géochronique 124, 23-39, F. Métivier, E. La-jeunesse, O. Devauchelle, A. Limare, L. Guérit, M. Houssais, G. Seizilles, P. Meunier, L. Malverti, H. Bouquerel, 2012.

2. Participation au film documentaire de Eric Pagès pour PLANETE +BAC : Science de la vie et de la Terre : La terre en mouvement 2013.

3. Formation en Techniques de visualization des écoulements de 4 artistes, étudiants en Master "Art et Science" de LadHyX, Ecole Polytechnique France.

4. Interview à LiveScience, December 18, 2015 : Microwavable Mantle : Physicists Nuke Mock Earth Layer, for Science.

5. Conference ‘Fête de la Science at IPGP’, November 2019, Tokyo.
## Liste des publications avant IPGP (1992-2006)

1. A. Frandas, M. Chirtoc, R. M. Candea, Photopyroelectric (PPE) Phase Spectroscopy for Sensitive Optical Absorption Coefficient and Thermal Diffusivity Measurement, Spectroscopy Letters, 25(8), 1155-1163, 1992

2. D.Dadarlat, R.M.Candea, A. Frandas, M.Chirtoc, D.Bicanic, PPE detection of 3D-2D Antiferromagnetic Phase Transition of Single Crystal La2CuO4, Photoacoustic and Photothermal Phenomena II, Ed. D.Bicanic, Springer Verlag 1992

3. D.Dadarlat, A. Frandas, I.Bratu, Photopyroelectric Spectroscopy of H2O-D2O Mixtures, Infrared Physics, 33(6), 575-579, 1992

4. D.Dadarlat, M.Chirtoc, R.M.Candea, A. Frandas, D.Bicanic, The Anomalous Behaviour of the PPE Signal in the Critical Region of a PhaseTransition, a New Support for the General Theoretical Model, Photoacoustic and Photothermal Phenomena II, Ed. D.Bicanic, Springer Verlag 1992

5. R.M.Candea, A.Lupu, M.Chirtoc, D.Dadarlat, A. Frandas, Photopyroelectric Spectroscopy of the Hydrogen Uranyl Phosphate, Spectroscopy Letters - 26(5), 923-934, 1993

6. D.Dadarlat, A. Frandas, Inverse Photopyroelectric Detection of Phase Transitions, Applied Physics A 56, 235-238, 1993

7. D.Dadarlat, D.Bicanic, H.Visser, A. Frandas, Combined Standard (SPPE)-Inverse (IPPE) Photopyroelectric Configuration for Measurement of Thermal Parameters of Liquid Fatty Acids, Revue de Physique IV, C7-483, 1994

8. D.Dadarlat, D.Bicanic, H.Visser, F.Mercuri, A. Frandas, Phototpyroelectric Method for Determination of Thermophysical Parameters and Detection of Phase Transitions in Fatty Acids and Triglycerides. Part I: Principles, Theory and Instrumentation Concepts, Journal of the Amer. Oil Chemists' Soc. 72(3), 273, 1995

9. D.Dadarlat, D.Bicanic, H.Visser, F.Mercuri, A. Frandas, Phototpyroelectric Method for Determination of Thermophysical Parameters and Detection of Phase Transitions in Fatty Acids and Triglycerides. Part II: The Temperature Dependence of Thermophysical Parameters, Journal of the Amer. Oil Chemists' Soc. 72(3), 281, 1995

10. A. Frandas, H.Jalink, R.Turcu, M.Brie, The Impulse Photopyroelectric Method for Thermal Characterization of Electrically Conducting Polymers, Applied Physics A 60, 455, 1995

11. D.Dadarlat, A. Frandas, M.Marinelli, F.Mercuri, D.Bicanic, A New Photopyroelectric Scheme Suitable for Phase Transitions Investigations: The Front Configuration with Semitransparent Sensor, Applied Physics A 61, 183, 1995

12. R.Turcu, M.Brie, A. Frandas, S.Pruneanu, Optical Studies on Free-standing Polypyrrole Films by the Photopyroelectric Method, Applied Physics B 62, 499-502, 1996

13. H. Jalink, A. Frandas, R. van der Schoor, D. Bicanic, New photothermal cell equipped with Peltier elements for phase transition studies, Review of Scientific Instruments 67(11), 3990, 1996

14. A. Frandas, H.Jalink, R. van der Schoor, Low frequency photoacoustics for monitoring the photobaric component in vivo of green leaves, Photosynthesis Research 52: 65-67, 1997

15. M. Chirtoc, I. Chirtoc, A. Frandas, D. Dadarlat, D. Bicanic, Rayleigh-Bénard-Marangoni Instabilities During Evaporation of Aqueous Alcohol Solutions, Detected by Pyroelectric Sensors. Instrumentation Science and Technology, 26, (1), 133 (1998)

16. A. Frandas, H. Jalink, R. van der Schoor, Monitoring the photobaric component of green leaves by low frequency photoacoustics, Photoacoustic and Photothermal Phenomena, AIP Conference Proceedings 463, F. Scudieri and M. Bertolotti (eds.), 594-596, 1998
17. A. Frandas, M.Chirtoc, D.Paris and M.Egee, Influence of humidity on thermophysical properties of starch sheets, Photoacoustic and Photothermal Phenomena, AIP Conference Proceedings 463, F. Scudieri and M. Bertolotti (eds.), 640-642, 1998

18. A. Frandas, V. Surducan, G. Nagy, D. Bicanic, The colligative properties of fruit juices by photopyroelectric calorimetry, Photoacoustic and Photothermal Phenomena, AIP Conference Proceedings 463, F. Scudieri and M. Bertolotti (eds.), 667-669, 1998

19. D. Paris, C. Bissieux, J.S.Antoniow, M.Chirtoc, A. Frandas, M.Egee, PPE study of water migration in starch sheets under moist atmosphere, Photoacoustic and Photothermal Phenomena, AIP Conference Proceedings 463, F. Scudieri and M. Bertolotti (eds.), 679-681, 1998

20. H. Jalink, A. Frandas, R. van der Schoor, R.J. Bino, Chlorophyll fluorescence of the testa of Brassica Oleracea as an indicator of seed maturity and seed quality. Scientia Agricola, 55, 88-93, 1998.

21. H. Jalink, R. van der Schoor, A. Frandas, J.G. van Pijlen, R.Bino, Chlorophyll fluorescence of Brassica Oleracea seeds as a new non-destructive marker for seed maturity and seed performance. Seed Science Research, 8, 437-443, 1998

22. A. Frandas, D. Dadarlat, M. Chirtoc, H. Jalink, D. Bicanic, D. Paris, J.S. Antoniow, M. Egee, and C. Ungureanu, Photothermal method using a pyroelectric sensor for thermophysical characterization of agricultural and biological samples, Proc. SPIE, Vol. 3405, 1039 (1998); DOI:10.1117/12.312710

23. D.Dadarlat, G.Nagy, U.Stanescu, C.Bele, M.Chirtoc, A. Frandas, I.Bratu, V.Surducan, E.Surducan, D.Bicanic, Combined photothermal and classical methods for the study of early spoilage of some edible oils, Advances in Food Sciences, 21 (5/6), 177-183, 1999

24. A. Frandas, D. Paris, C. Bissieux, M. Chirtoc, J.S. Antoniow, M. Egée, Classical and photopyroelectric study on optical and thermophysical properties of starch sheets function of water content and temperature, Applied Physics B71, 69-76, 2000

25. A. Frandas, T. Duvaut, D. Paris, Inverse analysis of water profile in starch by non contact photopyroelectric method, Applied Physics B71, 77-84, 2000

26. T. Duvaut, A. Limare, J.M Bachmann, Photothermal device for water dynamics measurement and thermophysical caracterisation. Application on textile samples. Rev. Sci.Instrum. 73 (3), 1299-1303, 2002

27. A. Limare, T. Duvaut, J.F. Henry, C. Bissieux, Non-contact photopyroelectric method applied to thermal and optical characterization of textiles. Four-flux modeling of a scattering sample. Int. J. of Thermal Science, 42 (11), 951-961, 2003

28. A. Limare, T. Duvaut, J.M Bachmann, Photopyroelectric measurement of the thermal properties of knitted textile samples. Influence of composition, structural parameters and water content. Int. J. of Thermal Science, 42 (11), 963-972, 2003

29. J. Bouajila, A. Limare, C. Joly, P. Dole, Lignin Plasticization to Improve Binderless Fiberboard Mechanical Properties. Polymer Engineering and Science 45(6) : 809-816., 2005

30. J. Bouajila, P. Dole, C. Joly, A. Limare, Some laws of a lignin plasticization, Journal of Applied Polymer Science, 102 (2) : 1445-1451, 2006
## Sélection d'articles

1. Davaille A., A. Limare, F. Touitou, I. Kumagai, J. Vatteville, Anatomy of a laminar starting thermal plume at high Prandtl number, Experiments in Fluids, DOI:10.1007/s00348-010-0924-y, 2010.

2. Androvandi S., A. Davaille, A. Limare, A. Foucquier, C. Marais, At least three scales of convection in a mantle with strongly temperature-dependent viscosity, Physics of the Earth and Planetary Interiors 188, 132-141, 2011.

3. A. Limare, M. Tal, M. D. Reitz, E. Lajeunesse, F. Métivier, Optical method for measuring bed topography and flow depth in an experimental flume, Solid Earth, 2, 143-154, doi:10.5194/se-2-143-2011, 2011.

4. A. Limare, K. Vilella, E. Di Giuseppe, C.G. Farnetani, E. Kaminski, E. Surducan, V. Surducan, C. Neamtu, L. Fourel, C. Jaupart, Microwave-heating laboratory experiments for planetary mantle convection, Journal of Fluid Mechanics, 777, 50-67, 2015.

5. K. Vilella, A. Limare, C. Jaupart, C.G. Farnetani, L. Fourel and E. Kaminski, Fundamentals of laminar free convection in internally heated fluids at values of the Rayleigh Roberts number up to $10^9$, Journal of Fluid Mechanics vol. 846, pp. 966-998, 2018.

6. A.Limare, C. Jaupart, E. Kaminski, L. Fourel, C. G. Farnetani, Convection in an internally-heated stratified heterogeneous reservoir, Journal of Fluid Mechanics, 870, 67-105, 2019.
# Anatomy of a laminar starting thermal plume at high Prandtl number

Anne Davaille · Angela Limare · Floriane Touitou · Ichiro Kumagai · Judith Vatteville

Received: 18 January 2010 / Revised: 20 May 2010 / Accepted: 22 June 2010
© Springer-Verlag 2010

**Abstract** We present an experimental study of the dynamics of a plume generated from a small heat source in a high Prandtl number fluid with a strongly temperature-dependent viscosity. The velocity field was determined with particle image velocimetry, while the temperature field was measured using differential interferometry and thermochromic liquid crystals. The combination of these different techniques run simultaneously allows us to identify the different stages of plume development, and to compare the positions of key-features of the velocity field (centers of rotation, maximum vorticity locations, stagnation points) respective to the plume thermal anomaly, for

Prandtl numbers greater than $10^3$. We further show that the thermal structure of the plume stem is well predicted by the constant viscosity model of Batchelor (Q J R Met Soc 80: 339–358, 1954) for viscosity ratios up to 50.

## 1 Introduction

Plumes are localized regions of buoyant fluid released continuously through the effect of gravity. They are, therefore, one of the most common feature of thermal convection, whereby heat transfer occurs by mass transport. Plumes are encountered in a number of industrial and natural situations, such as smoke release from chimneys, waste disposal in oceans, cloud formation processes, and convection in chemical reactors, magma chambers or in solid-rock planetary mantles. Plumes can either be turbulent, whereby their growth is dominated by turbulent entrainment of surrounding material (Morton et al. 1956) or laminar, in which case their growth is due to thermal diffusion, continuous feeding from the source, and laminar entrainment of surrounding material at the rear of a leading vortical head (Griffiths and Campbell 1990; Kumagai 2002).

Turner (1962) was the first to study what he called the starting plume, comprising a steady plume conduit capped by a large buoyant head. He further suggested that the buoyancy of the head increases since the head is fed by the stem as a result of the slower upwards motion of the head compared with the steady stem below it. Since then, much effort has been devoted to understand plume dynamics and to provide scalings for plume ascent velocity (Whitehead and Luther 1975; Shlien 1976; Olson and Singer 1985; Chay and Shlien 1986; Griffiths and Campbell 1990; Moses and al. 1993; Couliette and Loper 1995; van Keken

**Electronic supplementary material** The online version of this article (doi:10.1007/s00348-010-0924-y) contains supplementary material, which is available to authorized users.

A. Davaille (✉) · F. Touitou
Laboratoire FAST, CNRS/UPMC/U-PSud, Bat. 502,
Rue du Belvedere, Campus Universitaire,
91405 Orsay cedex, France
e-mail: davaille@fast.u-psud.fr

A. Limare · I. Kumagai · J. Vatteville
Dynamique des Fluides Géologiques,
IPGP, Univ. Paris 7/CNRS, 4 Place Jussieu,
75 252 Paris cedex 05, France

**Present Address:**
I. Kumagai
Division of Energy and Environmental System,
School of Engineering, Hokkaido University,
Kita-13, Nishi-8, Kita-ku, Sapporo 060-8628, Japan

**Present Address:**
J. Vatteville
Department of Mechanical Engineering, University of Thessaly,
Leoforos Athinon, Pedion Areos 38834, Volos, Greece
1997; Kaminski and Jaupart 2003; Rogers and Morris 2009), for steady-state plume stem structure (Batchelor 1954; Fuji 1963; Shlien and Boxman 1979; Tanny and Shlien 1985; Worster 1986; Olson et al. 1993; Vasquez et al. 1996; Laudenbach and Christensen 2001; Whittaker and Lister 2006a, b), or for plume growth by entrainment of ambient fluid (Griffiths and Campbell 1990; Moses and al. 1993; Couliette and Loper 1995; Kumagai 2002). However, different studies have often proposed different scalings. One reason could be the differences in experimental set-ups, such as the fluid used, the source supplying the buoyancy (e.g. injection of buoyant material through a hose or local heating) or the size of the experimental tank. Another reason is that the definition of the plume often depends on the visualization technique used. Distribution of dye (Whitehead and Luther 1975; Olson and Singer 1985; Griffiths and Campbell 1990; Kumagai 2002; Rogers and Morris 2009), small tracer bubbles (Shlien and Brosh 1979; Tanny and Shlien 1985; Couliette and Loper 1995), shadowgraph (Shlien 1976; Moses et al. 1993), Mach Zender (Boxman and Shlien 1978; Shlien and Boxman 1979, 1981; Chay and Shlien 1986) or differential interferometry (Kaminski and Jaupart 2003), and isotherms with thermochromic liquid crystals (Kumagai et al. 2007, 2008) led to different scaling laws and to diverging explanations about entrainment in the head and tail. Hence, it is essential to investigate how the technical differences and underlying assumptions influence our understanding of plume dynamics.

We focus here on the development of a laminar thermal plume generated from a local plane heat source in a viscous fluid, the viscosity of which depends on temperature. The plume dynamics are characterized using a combination of three different visualization techniques to determine quantitatively the temperature field (with Differential Interferometry DI and Thermochromic Liquid Crystals TLCs) and the velocity field (with Particule Image Velocimetry PIV). Section 2 describes the experimental setup and visualization techniques. The plume ascent is studied in Sect. 3, while Sect. 4 deals with the plume internal structure. The differences in plume shape, structure and uplift determined by each visualization technique are investigated. We also evaluate the influence of the temperature-dependent viscosity on the plume structure.

## 2 Laboratory experiments

### 2.1 Experimental set-up

Figure 1 shows the experimental set-up. The heat source is a Peltier element covered by a 4-mm-thick, 25-mm-

<figure><img src="image_40.png" /><figcaption>Fig. 1 Top view of the experimental set-up</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. 2 Temporal evolution of the electric power supplied to the Peltier device (top) and of its surface temperature (bottom). RUN 2</figcaption></figure>

diameter Cu disk. For each experiment, the Peltier heater is switched on at time $t = 0$ and a constant voltage is applied to it throughout the experiment. Both the ambient temperature $T_0$ and the temperature of the Peltier heater surface $T_{scr}$, measured by thermocouples, and the electric current voltage and intensity, are monitored through time. As shown on Fig. 2, this type of heat source delivers a quasi-constant power $P_{el}$ over the development of the plume. The latter is followed simultaneously through DI, TLCs and PIV. Electric power was systematically varied between 0.6 and 5.8 W. However, since a Peltier heater is a semiconductor, the thermal power $Q$ delivered to the fluid does not scale linearly with the electric power (see the manufacturer, http://www.melcor.com/cpseries.html, for details). Hence, for each experiment, the real thermal power was calculated from the temperature and velocity fields measured in the plume stem (see Sect. 2.5; Table 2).

Two cubic tanks of different sizes were used: a glass tank (60 × 60 × 60 cm³) and a plexiglas tank (15 × 15 × 20 cm³). They were filled with glucose syrup
(GS6075 from Syral). The physical properties of the fluid were measured in the laboratory, except for the specific heat which value is taken from the manufacturer (Table 1). The fluid rheology was characterized with a rheometer RS600 from ThermoRheo: it shows a Newtonian behavior but with a viscosity which depends strongly on temperature $T$ (Appendix 1).

For temperature-dependent viscosity fluids, plume dynamics are described essentially by three dimensionless numbers (i.e. Griffiths 1986). The Rayleigh number compares the driving thermal buoyancy forces to the resisting effects of thermal diffusion and viscous dissipation:

$$Ra = \frac{\alpha \rho g \Delta T_0 a^3}{\kappa v_0} \qquad (1)$$

where $a$ is the heater diameter, $\Delta T_0 = T_{scr} - T_0$ the temperature jump applied by the heater, $g$ the gravitational acceleration, $\alpha$ the thermal expansivity, $\kappa$ the thermal diffusivity and $v_0$ the kinematic viscosity ($v = \eta/\rho$) at the initial temperature of the fluid. The Prandtl number $Pr = v_0/\kappa$ compares thermal and momentum diffusions.

Finally, for temperature-dependent viscosity fluids, another important parameter is the viscosity ratio between the hotter parts of the fluid and the colder parts of the tank $\gamma = v_0/v_{hot}$. Studies in the Rayleigh-Benard configuration have shown that if this viscosity ratio was much greater than 10, motions tend to occur in a sublayer, under a cold conductive lid across which most of the viscosity variations occur (e.g. Morris 1984; Davaille and Jaupart 1993). A large viscosity ratio could therefore potentially change plume morphology and dynamics (Whitehead and Luther 1975; Olson and Singer 1985; Olson et al. 1993).

<table><caption>**Table 1** Physical properties of the glucose syrup at room temperature (20°C) properties of the fluids at 20°C</caption><thead><tr><th>Fluid</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Viscosity η, Pa s</td><td>11.95</td><td>7.35</td></tr><tr><td>Thermal diffusivity κ, m²/s</td><td>1.1 × 10<sup>-7</sup></td><td>1.1 × 10<sup>-7</sup></td></tr><tr><td>Specific heat C<sub>p</sub>, J/kg/K</td><td>2,280</td><td>2,280</td></tr><tr><td>Density ρ, kg/m³</td><td>1,399</td><td>1,389</td></tr><tr><td>Thermal expansion coefficient α, K<sup>-1</sup></td><td>4.15 × 10<sup>-4</sup></td><td>4.22 × 10<sup>-4</sup></td></tr><tr><td>Refractive index temperature variation dn/dT, K<sup>-1</sup></td><td>-2.05 × 10<sup>-4</sup></td><td>-2.05 × 10<sup>-4</sup></td></tr></tbody></table>

<table><caption>**Table 2** Characteristics of the experiments</caption><thead><tr><th>Run</th><th>Tank</th><th>T<sub>0</sub> (°C)</th><th>T<sub>scr</sub> (°C)</th><th>P<sub>el</sub> (W)</th><th>ΔT<sub>0</sub> (°C)</th><th>γ</th><th>Ra</th><th>Pr</th><th>Q(W)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>21.6</td><td>45.5</td><td>1.46</td><td>23.9</td><td>10.8</td><td>2,011</td><td>64,600</td><td>0.72 ± 0.2</td></tr><tr><td>2</td><td>1</td><td>23.0</td><td>43.8</td><td>1.17</td><td>20.8</td><td>8.0</td><td>1,751</td><td>54,750</td><td>0.65 ± 0.2</td></tr><tr><td>3</td><td>2</td><td>21.06</td><td>75.45</td><td>5.80</td><td>54.4</td><td>47.2</td><td>6,543</td><td>44,440</td><td>1.60 ± 0.25</td></tr><tr><td>4</td><td>2</td><td>21.6</td><td>64.3</td><td>3.1</td><td>42.7</td><td>20.6</td><td>5,337</td><td>42,800</td><td>1.10 ± 0.2</td></tr></tbody></table>

The values of these parameters for the experiments discussed in this paper are given in Table 2. The high $Pr$ values ensure the absence of inertial effects in our tank, as for example in the Earth's mantle. The $Ra$ values ensure the formation of plumes (e.g. Griffiths and Campbell 1990). We see also that the viscosity ratio is between 8 and 50. One of the questions to be answered will therefore be whether or not such viscosity ratios are sufficient to have a strong effect on plume dynamics.

## 2.2 Temperature field: differential interferometry

One of the main interferometric tools used in convection studies has been the Mach Zender interferometer, which produces fringe shifts proportional to local temperature (e.g. Gebhart et al. 1970; Boxmann and Shlien 1978). Here, we prefer to use differential interferometry (also called shearing interferometry) because of the simpler setup involved, which enables us to combine it with other techniques. DI measures derivatives of light phase distortions (e.g. Sernas and Fletcher 1970), using the interferences between a beam passing through an object and the same beam laterally displaced. The simplest way of producing these two beams is to use a flat optical plate tilted at a certain angle (Kaminski and Jaupart 2003). So, we used a HeNe ($\lambda$ = 632 nm)-expanded laser beam as a light source and placed behind the tank a tilting flat plate (thickness $t_p$ = 4.16 mm and refractive index $n_p$ = 1.51872) at an incidence angle $i$ = 10.5° (Fig. 1). In the absence of the object, the displacement $\delta$ introduced by the tilting plate is (Kelley and Hargreaves 1970):

$$\delta = 2t_p / \sqrt{\sin i \cdot \cos i \cdot n_p^2 - \sin^2 i} \qquad (2)$$

On the interferogram, the orientation of fringes is perpendicular to the tilt plane and their frequency is related to $\delta$. Therefore the sensitivity of the method can be varied. The interferograms were recorded every 10 s with a digital camera (Canon 8.2 Mpixels). The vertical fringes away from the plume correspond to light propagation through uniform background (Fig. 3a). In the plume region, those fringes are distorted by horizontal temperature gradients. The dimensions of investigated objects are only limited by the diameter of the exit lens of the beam expander (80 mm in our case). For strongly refracting axisymmetric objects,
**Fig. 3** Plume at *t* = 475 s for RUN #1. **a Left:** differential interferometry image. The red line shows the contour of the non-zero thermal anomaly (DI on Fig. 3b, d, e). **Right:** the bright lines are the isotherms. **b** Phase *h* on a horizontal cross-section at 60 mm height. **c** Phase difference *d* between the distorted interferogram *h* and the initial one *h*₀. The black Sh line shows the location of the maximum phase displacement derivative, which would correspond to a shadowgraph contour for a purely cylindrical geometry. **d** refractive index distribution, and **e** temperature field as a function of radius

<figure><img src="image_42.png" /></figure>

it has been shown (Vest 1975) that accurate results can be obtained only if the interferogram is formed with appropriate imaging. The convergent lens introduced into the pathway forms the image of the center plane of the object on the screen (Fig. 1). The positioning of the lens ensures a 1:1 object to image ratio.

### 2.2.1 Evaluation of the interferogram

From the intensity profile of the interferogram along one straight horizontal line (Fig. 3a), the phase function *h* is calculated. The difference in phase between each consecutive fringe is $2\pi$ radians (360°). For the initial interferogram, the phase function $h_0$ is a straight line as a function of the pixel number (Fig. 3b). When an object is present, the phase function will be distorted (Fig. 3b), and the difference between the two lines is $d(y)$. The plume axis is determined as the point where the two phase functions are intersecting and corresponds to phase difference $d(y) = 0$ (Fig. 3b, c). Then, for an optically thin, axially symmetric object such as our laminar thermal plume, the refractive index distribution introduced by the object can be calculated numerically by simple integration (Pretzler et al. 1993) (Fig. 3d):

$$ \frac{2\pi}{\lambda}(n(r) - n_0) = -\frac{1}{\pi\delta}\int_{r}^{\infty} \frac{d(y)dy}{\sqrt{y^2 - r^2}} \quad (3) $$

Then, knowing the temperature dependence of the refractive index, the temperature field can be calculated (Fig. 3e). This procedure was applied to all the images. We chose the incidence angle and the characteristics of the flat plate in order to image accurately most of the plume temperature anomaly. For distances farther than 3 mm from the heat source, the uncertainty is estimated at 0.2°C

or 2% of the total temperature difference. Closer to the heat source, the strong temperature gradients render the fringes too thin to be accurately described.

DI has the advantage to give the complete temperature field upon inversion of the interferograms. Its disadvantages are that: (a) the dimension of the investigated object depends on the expansion beam diameter which is limited by the optics (here only 80 mm diameter), (b) since the test region is imaged, the position of the thermal anomaly has to be well determined in advance and (c) it cannot be used quantitatively for geometries which are neither 2D nor axisymmetric.

In the present study, we used DI as a reference to validate our TLC method (see next section) and to determine how the partial image of a plume given by the TLCs compares with the total temperature field given by DI. Disadvantage (a) is the most limiting problem as we cannot image the plume on its whole ascent. We have therefore either focussed on the source region of the plume, keeping the whole set-up fixed (RUN 1, movie 1), or on the plume cap (RUN 2, movie 2). In the latter, we overcome (a) and follow the plume cap by rising both the laser and cameras at the same velocity as the plume cap (movie 2).

### 2.3 Selected isotherms: thermochromic liquid crystals

The use of thermochromic liquid crystals (TLC) allows temperature mapping on a 2D-plane in the fluid flow without perturbing it (e.g. Dabiri 2009 for a review, Rhee et al. 1984; Dabiri and Gharib 1991; Willert and Gharib 1991). When illuminated by white light, TLC color changes with increasing temperature from invisible to red at low temperatures, passes through green and blue to violet and turns invisible again at high temperatures. The total temperature range $\Delta T_c$ accessible with one TLC is between 1
and 40°C. With a high precision color CCD camera, a very
stable white light and a precise calibration of the TLC color
against true temperature, it is then possible to determine the
temperature field with a 5% precision over ΔT<sub>c</sub> (Fujisawa
and Adrian 1999). Since: (a) thermal convection in viscous
fluids at high Rayleigh number involves typical tempera-
ture heterogeneity of 20–60°C (e.g. Davaille and Limare
2007), and (b) we want a temperature precision better than
0.2°C, we choose to seed the experimental cell with several
types of TLCs, each reflecting light at a different temper-
ature, and to illuminate a cross-section of the tank with a
monochromatic laser light sheet. Therefore, each type of
TLC is responsible for one bright line, which represents an
isotherm (Figs. 5, 7a). This method was previously used to
study two-layer convection (Le Bars and Davaille 2002),

<figure><img src="image_43.png" /><figcaption>**Fig. 4** Thermal plume visualized by TLCs (RUN 4). The *bright lines* are the isotherms. **a** Vertical cross-section along the plume axis. The snapshot is taken with the camera at 90° from the laser sheet. The *red arrows* indicate the position of the horizontal cross-section of (**b**) and (**c**). **b** Sketch showing the laser sheets positions. **c** Snapshot taken from the top at a 20° from the vertical. It allows seeing the horizontal cross-section of the plume stem as well as the vertical cross section above it</figcaption></figure>

<caption>**Fig. 5** Comparison between the TCL images and the temperature field calculated from DI for glucose syrup at different moments of the plume evolution. The points are obtained by intersecting the DI temperature field to the discrete values of the isotherms. *Red line* indicates the DI plume boundary. RUN 1</caption>

plumes out of a viscous thermal boundary layer (Davaille
and Vatteville 2005; Davaille et al. 2008) and thermo-
chemical plumes out of a localized heat source (Kumagai
et al. 2007, 2008). It can be combined with Laser Induced
Fluorescence and Particle Image Velocimetry to obtain
simultaneously the composition and velocity field (Kuma-
gai et al. 2007).

Here, the sugar syrup was seeded with five types of
liquid crystals (20–30 µm aqueous slurries from Japan
Capsular Products, Inc. for the large tank; 40 µm aqueous
slurries from Hallcrest for the small tank). A vertical cross-
section of the experimental tank centered on the plume axis
is illuminated by a 532-nm YAG (Coherent VERDI 2W)
laser sheet 1 mm thick (Fig. 4a). The vertical cross-section
forms a small angle (~15°) with the section imaged by DI
(Fig. 1), so that the results from TLC and DI will be
comparable only if the plume is axisymmetric. This prop-
erty is indeed verified on the horizontal cross-sections of
the plume (Fig. 4b, c). A black and white CCD camera
(LaVision 1,280 × 1,020 pixels) recorded the scattered
light from the liquid crystals (Fig. 1). Five isotherms were
imaged, with a precision for their peak intensity value
of ±0.1°C (see Appendix 2 for details).

2.3.1 Comparison with DI

Figure 5 shows the TLC images superposed with the cor-
responding data for several heights measured by DI. The
agreement between the two methods is good. Even in the
hottest part of the plume, the positions of the TLCs iso-
therms can be superimposed on the same isotherms deter-
mined by DI within 0.1 mm. This is coherent with the fact
that optical distortion should be everywhere negligible in
our plume parameter range (see Appendix 3).

The TLC method has several advantages: (a) direct visualization of some discrete isotherms, (b) a simple setup with a visualization cross-section which can encompass the

<figure><img src="image_44.png" /></figure>
whole experimental tank (a simple cylindrical lens can produce almost any dimensions for the laser sheet), (c) the laser scanning system allows imaging of fully 3D thermal anomalies since no axisymmetric hypothesis is needed (Le Bars and Davaille 2002; Davaille and Vatteville 2005), (d) a simple calibration procedure which can allow isotherms determination with a 0.1°C precision for aqueous solutions (Davaille et al. 2008). The biggest disadvantage of the TLC method remains that only a few isotherms are measured.

## 2.4 Velocity field

The first velocity field measurements in thermals and laminar plumes were done using the lengths of streaks in time-exposure photographs of dust in the fluid (Shlien and Brosh 1979; Tanny and Shlien 1985). Here, the fluid was seeded with tiny bubbles (0.1 mm in diameter) for the large tank, and with hollow glass spheres 10 µm in diameter for the small tank. The viscosity of the fluids is large enough and the particles small enough so that the Stokes velocity of the particles is negligible compared to the velocity of the fluid motions. Therefore, the particles act as passive tracers of the fluid motions. Images were taken every 0.5 s along with the TLCs, with a CCD camera (1,024 × 1,391 pixels). The velocity field was then calculated using two methods: Particle Image Velocimetry (PIV) and individual particle tracking. First, PIV was performed through cross-correlation between images, using the PIV package Davis from LaVision. The size of the interrogation window was 64 pixels with an overlap of 87% or 32 pixels with an overlap of 50%, leading to a spatial resolution of 3.5 mm. Due to the small velocity of our laminar plumes, the time lapse between two correlated images was taken between 1 and 2.5 s. Typical velocities along the plume axis reach 1 mm/s. Due to spurious light scattering on the bubbles, the laser sheet was enlarged in the big tank and the velocity field accuracy was only 30% (RUN 1 and 2). However, the use of smaller glass spheres in the small tank allowed to reach a 5% precision (RUN 3).

However, the PIV technique gives a spatial average of the velocity field and therefore, we could miss the details of the velocity structure around the plume axis where the velocity gradients are larger. In order to check that issue, the velocity was also measured by following individual particles. This method is much more time-consuming than the PIV. Figure 6 shows a comparison between the two techniques for the experiment where the velocities and velocity gradients were the largest (RUN 3), i.e. the case where the discrepancies should be the largest. We see that both techniques capture the essential features of the plume velocity field, but that the largest PIV interrogation window underestimates the velocities around the plume axis in the area close to the heat source.

<figure><img src="image_45.png" /><figcaption>Fig. 6 Comparison between the different determinations of the velocity a along the axis, b at height z = 25 mm, c at z = 70 mm. In green, the PIV was calculated in square interrogation windows of 64 × 64 pixels. In black, the interrogation window was 32 × 32 pixels. The red line shows the values obtained by following individual particles. RUN 3</figcaption></figure>

## 2.5 Estimate of the thermal power injected in the plume in steady state

From the temperature and velocity fields in the plume stem, we can compute the real heat power Q delivered by the Peltier heat source with:

$$Q(z) = 2\pi\rho C_p \int_{r=0}^{\infty} rV_z(r, z)(T(r, z) - T_0)dr \quad (4)$$

Within the uncertainty of the velocity field, we found a constant heat power along the stem. However, this thermal power Q is not directly proportional to the electric power P, but their ratio Q/P decreases as P increases (Table 2). This is due to the non-linear answer of the Peltier heater (see the manufacturer, http://www.melcor.com/cpseries.html, for details).

## 3 Plume ascent

Movies 1–2 and Figs. 7, 8 and 9 show the development of the thermal plume through its temperature and velocity fields. Once the heater is turned on, a thermal boundary layer (TBL) grows over it by heat conduction (Fig. 5a). When the Stokes velocity of this hot pocket becomes greater than its rate of growth by heat diffusion, a thermal starting plume is generated (Whitehead and Luther 1975; Moses et al. 1993; Kaminski and Jaupart 2003). The hot plume has the classical mushroom structure (Figs. 5b–d, 7a), with a nearly spherical head on top of a thin stem. Previous experimental studies (Whitehead and Luther
**Fig. 7** Development of the plume through time. RUN 3. **a** Snapshots at several times. The *bright lines* are the TLCs isotherms. **b** Velocity modulus. **c** Velocity spatial derivative $dV_r/dr(r, z)$. The *insert* shows the stagnation point above the plume head, where $dV_r/dr(r, z)$ should be maximum. **d** Grid deformation calculated from the velocity fields. The *yellow crosses* show the positions of the stagnation point determined in (c).

<figure><img src="image_46.png" /></figure>

1975; Olson and Singer 1985) have shown that this morphology always obtained when the plume material was less viscous than the ambient fluid. This is indeed the case here where the hot material has a viscosity much lower than the ambient fluid. During the plume rise, the velocity maximum along the plume axis is first located in the plume head but later descends into the stem where it remains steady (Figs. 7b, 9a).

To measure plume ascent, we need to consider the plume as an object invading the ambient fluid. Then, there

are several ways to identify the plume by its contours and/or averaged properties.

## 3.1 Plume as a thermal anomaly

Since a plume develops because of its buoyant thermal anomaly, the most natural way to describe a plume is by following the thermal anomaly (DI, TLCs) or by following the associated temperature gradients (as a shadowgraph would do).
<figure><img src="image_47.png" /><figcaption>**Fig. 8** Plume uplift through time. **a** Image intensity along the plume vertical axis; the *bright lines* are the isotherms (labelled in *white*). **b** Superposition of the plume images taken every 5 s between *t* = 100 and *t* = 500 s. The internal structure of the plume conduit, with its 5 isotherms, appears to be steady. The *outer bright line* underlines the</figcaption></figure>

steady growth of the plume head. The positions of the centers of rotation (corresponding to the position of the maximum values of the stream function) are represented in *blue* (early times) to *red* (late times). RUN 3

<figure><img src="image_48.png" /><figcaption>**Fig. 9** **a** Velocity along the plume axis as a function of time. **b** $dV/dr$ along the plume axis as a function of time. **c** Mean horizontal vorticity modulus as a function of depth and time. All three variables have been normalized by their respective maximum values. The *black line*</figcaption></figure>

represents the uplift of the stagnation point A as measured from (b). The *green line* has been obtained by subtracting the plume half radius to the stagnation point height, and we see on (c) that it follows quite well the uplift of the vorticity maximum. RUN 3

The maximum extension of the plume can be determined from the complete temperature field measured by DI. The thermal plume can be defined as the zone where the thermal anomaly introduces significant phase distortion; so the plume boundary is defined as the point where the slope of the phase function changes abruptly (red markers in Fig. 3b, c, e; red line in Figs. 3, 5): in our case, this corresponds (arbitrarily) to a phase shift of $\pi/2$. Because of the non-linear character of the equations, it does not correspond to a given isotherm or to a given temperature anomaly fixed for a given fluid or even along the duration of an experiment. For our experiments, it decreases from 0.8 to 0.2°C relative to the tank initial temperature as time increases.

Fixed temperature anomaly contours can also be obtained from DI or directly from TLCs. Here, our coldest isotherm is +2.4°C higher than the initial temperature of the tank. Its contour will therefore underestimate the plume extension compared to DI.

Temperature gradient maximum can be obtained from the temperature field measured by DI, TLCs, or directly visualized by shadowgraph.

## 3.2 Plume visualized by the velocity field and/or a passive tracer

The distribution of the vertical component of the velocity field is given in Fig. 7b. At first, the velocity maximum is
localized in the plume head; then it moves into the stem. A similar behavior was observed in plumes in nearly constant viscosity silicone oils (Vatteville et al. 2009). Following, the velocity maximum is therefore not a good way to follow the plume uplift. However, previous laboratory experiments reported that a dye carried by the buoyant axial conduit reaches a stagnation point at the top of the head on the vertical axis (point A on Fig. 7c), from where it spreads laterally into thin lamellae (e.g. Griffiths and Campbell 1990; Kumagai 2002). The existence of this stagnation point in the referential of the growing head is also inferred by theoretical studies of isolated thermals (e.g. Griffiths 1986; Whittaker and Lister 2008). Therefore, an interesting quantity to follow is $\partial V_r/\partial r$, which should present a local maximum at point A. This is indeed the case (Figs. 7c, 9b). Moreover, for a stagnation point in an axisymmetric system, conservation of mass implies that

$$ \partial V_z/\partial z = -2\partial V_r/\partial r. \qquad (5) $$

For the plume portrayed in Fig. 7, the ratio between $\partial V_z/\partial z$ and $\partial V_r/\partial r$ is equal to $1.94 \pm 0.32$, which, within error, agrees with the theoretical value of 2.0.

We also used numerical particles as passive markers (to mimic dye) and advected them using the velocity fields calculated by PIV from the laboratory experiment. Initially, they are regularly dispersed in a 2D numerical grid. At each time step, they are advected using the velocity field spatially and temporally interpolated at the position of the particles using a fourth-order Runge-Kutta method (van Keken et al. 1997; F. Touitou, G. Brandeis, A. Davaille, Influence of a strongly temperature-dependent viscosity on mantle plume dynamics and sampling, in prep.). The top of the plume head is characterized by strong horizontal stretching and vertical compression of the initial square grid (Fig. 7d). The top of the plume head followed by a passive dye is therefore well described by the point where $\partial V_r/\partial r$ is maximum on the vertical axis. This gives a useful mean to follow the plume head uplift when only velocity fields are measured.

### 3.3 Plume ascent and growth seen by the different techniques

Figure 10a presents a comparison between the different observations recorded along the plume axis for a typical experiment. In particular, it shows the temporal evolution of the top edge of the plume cap evaluated by the four different methods previously described: the top boundary of the minimum thermal anomaly measured by DI (hereafter DI-top), the maximum temperature gradient calculated from the temperature structure (hereafter shadow-top), the height of each TLC isotherm (hereafter TLC-top), and the location of the stagnation point A (hereafter ptA).

As expected, DI-top, shadow-top and ptA all describe the same temporal evolution. Shadow-top and ptA locations are

<figure><img src="image_49.png" /><figcaption>Fig. 10 a Plume uplift along vertical axis measured with DI maximal extension (dashed line), TLC isotherms (solid line), stagnation point “A” (red disks), and maximum temperature gradient (black squares). b Plume growth measured by DI anomaly maximum radial extent (dashed line) and the coldest TLC isotherm radial extent (solid line). RUN 2</figcaption></figure>

undistinguishable within the 1-mm precision of the experiments. Away from the top boundary, the plume appears to rise at a constant velocity, in agreement with previous studies (e.g. Chay and Shlien 1986; Tanny and Shlien 1985; Moses et al. 1993; Kaminski and Jaupart 2003; Rogers and Morris 2009). Following Batchelor (1954) and Worster (1986), Kaminski and Jaupart (2003) further found that at high Prandtl number in a large box and a nearly constant-viscosity fluid (silicone oil), the “shadow-top” velocity was scaling as:

$$ W_{KJ} = (0.57 \pm 0.02) \frac{\ln \epsilon^{-2}}{2\pi} V_{\text{Batch}}, \qquad (6) $$

where $V_{\text{Batch}} = \sqrt{\frac{\alpha g Q}{\rho c_p v}}$ (Batchelor 1954), and $\epsilon$ is solution of $\epsilon^4 \ln \epsilon^{-2} = 1/Pr$ (Worster 1986). In our experiments, the viscosity depends strongly on temperature and the question is, therefore, which value of the viscosity is relevant? Now, the ascent velocity of the plume head should be controlled by the viscosity of the most viscous fluid, i.e. the ambient fluid, so that the viscosity in Eq. 6 should be $v_0$. Moreover, the Stokes velocity of a less viscous sphere rising in a more viscous fluid depends also on the viscosity ratio (e.g. Whitehead and Luther 1975; Griffiths 1986), so that the velocity of our sugar plumes should finally scales as:
$$W(\gamma) = \frac{5\gamma + 1}{4\gamma + \frac{3}{2}}W_{KJ}, \qquad (7)$$

For RUN 2 (Fig. 10, γ = 8, large box), Eq. 7 predicts $W = 0.682V_{Batch}$, in very good agreement with the measured value ($W = 0.675V_{Batch}$). However, for experiments with larger viscosity ratios and run in a smaller box (e.g. RUN 3, Figs. 7, 8, 9), the measured velocity value ($W = 1.05V_{Batch}$) is systematically greater than the predicted one ($W = 0.705V_{Batch}$). Since confinement is expected to decrease the plume head velocity (e.g. Happel and Brenner 1973), the opposite trend observed here is due to temperature-dependent viscosity effects. They will be studied systematically elsewhere (Touitou et al., in prep.).

The position along the axis of each isotherm follows the same constant uplift until a certain time. As the plume head rises, it cools, and therefore the isotherms are disappearing one by one from the head (Figs. 5, 7, 10). They will eventually be located only in the plume stem. In RUN 1 and 2, only the coldest isotherm follows the three other indicators throughout the plume life (Fig. 10a).

Figure 10b follows the plume head growth with the maximum lateral extensions of the minimum temperature anomaly seen by DI, and of the coldest TLC. The two measurements give different values, owing to the diffusive thermal halo around the plume head. However, up to ~600 s, DI and the coldest isotherm give the same growth rate. For $t > 600$ s, the TLC-derived radius becomes constant and will eventually decrease as the head cools under the action of thermal diffusion.

# 4 Plume internal structure

## 4.1 Plume stem

Figures 10a and 8 show that the plume stem thermal structure remains constant as soon as it is built, therefore confirming an earlier observation by Shlien and Boxman (1981). The stem quickly achieves a steady-state structure which is independent of the dynamics of the ascending head above. Since we have the complete temperature field, we can test the scaling laws which have been proposed for a constant viscosity plume stem (Batchelor 1954; Shlien and Boxman 1979), and see whether they still apply for the temperature-dependent viscosity stem for our range of viscosity ratios (8–50).

### 4.1.1 Thermal structure

Batchelor (1954) considered a model in which a plume rises from a point source of heat of power $Q$ in an infinite fluid with constant properties and high Prandtl number.

Simple scaling arguments suggest that the temperature anomaly $\Delta T(z)$ along the axis ($r = 0$), and the stem radius $R(z)$, scale with the height $z$ above the source as:

$$\Delta T(z) = T(r = 0, z) - T_0 \sim \frac{Q}{\kappa \rho C_p z}, \qquad (8)$$

$$R(z) = \left( \frac{\rho C_p v \kappa^2}{g \alpha Q} \right)^{1/4} z^{1/2} \qquad (9)$$

Thermal diffusion causes the temperature anomaly to decrease with height as $z^{-1}$ and the stem radius to increase as $z^{1/2}$. Experiments in water ($Pr \sim 7$; $\gamma \le 2$; Shlien and Boxman 1979; Moses et al. 1993) and silicone oils ($20 \le Pr \le 10^5$; $\gamma \le 2$; Moses et al. 1993; Kaminski and Jaupart 2003) verify the scalings (Eqs. 8–9) if corrections are made for the finite size of the heat source: this is done by replacing in Eq. 8 and 9 $z$ by $(z - z_0)$ where $z_0$ is the origin of a real point source equivalent to the finite size source.

Figure 11 shows the temperature along the plume axis measured by DI as the plume rises and in steady state when the isotherms have reached their final height (Fig. 8a). From the heat source to right below the plume head, the two types of measurement agree, and Eq. 8 explains the data well, including the temperature of the Peltier heater. Figure 12 further shows that the steady-state structure of the plume stem for all the experiments can be predicted by:

<figure><img src="image_50.png" /><figcaption>Fig. 11 Temperature anomaly along axis for RUN 2. Red stars stand for the position of the TLCs in steady state, and blue stars for the maximum height reached by each TLC. Dashed thin lines represent the temperature profiles measured by DI at different times. The red solid line represents a fit by Batchelor's model (Eqs. 8–10) along the whole plume with $z_0 = -20.8$ mm and $C = 0.22$. The dashed thick blue line follows the decrease in the temperature maximum located in the plume head and can also be predicted with Batchelor's model with $C_{head} = 0.40$ and $z_{head} = 34.6$ mm</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>Fig. 12 a Temperature anomaly along the plume stem axis in steady state for all the experiments. For each experiment, the temperature anomaly of each TLC isotherm has been normalized by the total temperature difference in the experiment. The solid line represents the fit by Eq. 10. Purple diamonds stand for RUN 1, blue crosses for RUN</figcaption></figure>

2, green disks for RUN 3, and red circles for RUN 4. b Velocity anomaly along the plume stem normalized by $V_{\text{Batch}}(z)$ (see Eq. 12) for RUN 2 (blue) and RUN 3 (green). The ratio appears constant far from the bottom and top boundaries

$$ \Delta T(z) = C \frac{Q}{\kappa \rho C_p (z - z_0)} = \frac{-z_0 \Delta T_0}{(z - z_0)} \quad (10) $$

with $z_0 = -21.5 \pm 1.0$ mm. A priori, $z_0$ could depend on the Rayleigh number and on the geometry of the source (e.g. Whittaker and Lister 2006b). However, in the experiments presented here, $z_0$ does not seem to depend on the heat source power, nor on the viscosity ratio.

The experimental constant $C = 0.23 \pm 0.02$ has a bigger uncertainty, due in particular to the large error bars on some of the velocities, and therefore on the calculation of the thermal power $Q$. However, a closer look shows also that Eq. 10 systematically underpredicts (by 5–10%) the axis temperature of the experiment with the biggest viscosity ratio. This might hint that the viscosity ratio begins to play a role when $\gamma \ge 40$.

According to Batchelor’s model, the radial temperature structure of the stem should be given by:

$$ \theta(r, z) = T(r, z) - T_0 = \Delta T(z) f \left( \frac{r}{R(z - z_0)} \right) \quad (11) $$

where $f$ should describe the same function for all the experiments. Figure 13 shows that the positions of the isotherms indeed describe a single curve, while Fig. 14 shows that all the experiments are within 10% of the same curve. So, within the parameter range of the present study, Batchelor’s model (derived for constant viscosity fluids) described the data with a maximum error of 10%.

### 4.1.2 Velocity structure

According to Batchelor’s model for a constant viscosity fluid, the velocity along the plume stem axis should be

<figure><img src="image_52.png" /><figcaption>Fig. 13 Function $f$ from Eq. 11 as a function of the normalized radius $r/(z - z_0)^{1/2}$. For the isotherms, Eq. 10 implies that $f \sim \Delta T_{\text{iso}}(z - z_0)/\Delta T_{0z}$, with $\Delta T_{\text{iso}} = T_{\text{iso}} - T_0$. Triangles stand for isotherm 24.35°C, diamonds for isotherm 27.0°C, crosses for isotherm 31.1°C, circles for isotherm 35.0°C and stars for isotherm 39.5°C. RUN 3</figcaption></figure>

constant, equal to $V_{\text{Batch}}$. However, Fig. 6a shows that the velocity is decreasing with increasing height. This is true for all our experiments. Because the temperature along the stem decreases with increasing height, the viscosity increases within the plume conduit. So that velocity should be reduced accordingly. Figure 12b indeed shows that within our parameter range, the material velocity along the plume axis away from the bottom and top boundaries scales as:
<figure><img src="image_53.png" /><figcaption>**Fig. 14** Function f from Eq. 11 as a function of the normalized radius r/R(z). The power and physical properties of the fluid enter in the calculation of R(z) according to Eq. 8. Triangles stand for isotherm 24.35°C, diamonds for isotherm 27.0°C, crosses for isotherm 31.1°C, circles for isotherm 35.0°C and stars for isotherm 39.5°C. RUN 3 is in green and RUN 4 in red. The black dashed and solid lines represent the radial structure obtained from DI measurements for RUN 1 and RUN 2, respectively. Only one side of the plume has been represented in the latter cases for the sake of clarity</figcaption></figure>

$$W(z) = 0.78(\pm 0.1) \sqrt{\frac{\alpha g Q}{\rho C_p \nu(z)}} \quad (12)$$

The velocity along the axis seems, therefore, proportional to the local Batchelor's velocity scale $V_{\text{Batch}}(z)$ where the viscosity value is taken at $T(z)$. We were not able to find a simple scaling to describe the velocity radial structure. This is due to the effects of confinement and of temperature-dependent viscosity. We shall discuss the effect of confinement in Sect. 4.4.

## 4.2 Plume head

The plume head contains a local maximum in temperature which is located on the plume axis (Figs. 7a, 11). This temperature local maximum along the axis was first reported by Chay and Shlien (1986), but its origin had remained elusive. We see here that for each experiment, we can trace it all the way to the heat source, and from the very early stages of the plume development. It is, therefore, inherited from the initial construction of the head.

This local temperature maximum decreases with time, as the head cools by heat diffusion. As already noticed by Moses et al. (1993), its decrease can also be fitted by Batchelor's model, albeit with different values for C and $z_0$ (Fig. 11).

<figure><img src="image_54.png" /><figcaption>**Fig. 15** Particle velocity $V_z(A)$ at the stagnation point A (black), compared to the uplift velocity $V_A = dz_A/dt$ of the stagnation point (red), and to $V_A$ minus the plume head growth rate (green on Fig. 9c) as measured from Fig. 10b). RUN 3</figcaption></figure>

## 4.3 “Duration” of an isotherm

According to Fig. 8, we see that each isotherm becomes extinct as the rising plume is cooling. The isotherms visualization is, therefore, useful to characterize the plume dynamics only for a certain time. From the two previous sections and Eqs. 7–10, we can determine for how long a given isotherm will be useful to trace the plume head and growth, or to trace the plume conduit. From Eq. 10, the maximum height an isotherm will reach is given by:

$$z_{\text{isoth}} = -z_0 - \frac{z_0 \Delta T_0}{\Delta T_{\text{isoth}}} = -z_0 + C \frac{Q}{\kappa \rho C_p \Delta T_{\text{isoth}}} \quad (13)$$

with $C = 0.24$ and $z_0 = -21.5$ mm for the stem, and $C = 0.40$ and $z_0 = 34.6$ mm for the head.

## 4.4 Anatomy of a thermal starting plume at high Prandtl number

We have seen in Sect. 3.3 that the ascent of the top of the plume could be followed by the uplift of the stagnation point A. It is tempting to assimilate it to the “plume uplift”, implying that at least the bulk of the plume head is rising at this velocity. However, we can notice that the vertical velocity $V_z(A)$ of the particles at point A (calculated from PIV) is always smaller (by 5–10% in our experiments) than the uplift velocity of A, $V_A = dz_A/dt$ (Fig. 15). This is due to the fact that the plume head is also growing by thermal diffusion: as shown in Fig. 15, $V_A$ is in fact the sum of the particles velocity $V_z(A)$ and the diffusive growth rate of the head. So is there a strategic plume feature which would rise at the “true” velocity $V_z(A)$?

In their pioneering study, Shlien and co-authors (Shlien and Brosh 1979; Tanny and Shlien 1985) had defined the plume radius by the distance of the center of rotation $C_R$ (Fig. 16a) from the plume axis, and the plume uplift velocity by the uplift velocity of $C_R$. The centers of rotation
<figure><img src="image_55.png" /><figcaption>**Fig. 16** a Streamlines in the laboratory reference frame on top of the plume image. b Streamlines in the “plume head” reference frame on top of the plume image and c on top of the vorticity field. The yellow</figcaption></figure>

crosses represent the stagnation point A, the stars the rotation centers $C_R$ in the laboratory reference frame, and the yellow open circles point to the vorticity extrema. RUN 3, $t = 225$ s

are the points where the flow velocity is zero in the laboratory reference frame, which corresponds also to extrema of the stream function. Their fluid was water ($Pr \sim 7$) and the circulation was due to inertial effects. $C_R$ were very close to the plume thermal anomaly since the velocity boundary layer thickness was comparable to the thermal boundary layer's, i.e. the plume conduit width.

On the other hand, in an unbounded fluid at infinite Prandtl number, the velocity is the integral of a Stokelet distribution (e.g. Whittaker and Lister 2006a) and should be everywhere upwards. So, the existence of $C_R$ in a laboratory experiment with a large $Pr$ should therefore be entirely the effect of confinement in the experimental tank.

We calculated the stream function $\Psi$ from our velocity fields using (Tanny and Shlien 1985):

$$ \Psi(r, z) = \int_{0}^{r} r \cdot V_z(r, z) \cdot dr. \qquad (14) $$

As expected for high Prandtl number fluids, the velocity boundary layer is much thicker than the plume thermal anomaly, and therefore $C_R$ are well outside the plume thermal anomaly (Fig. 16a). Moreover their positions are indeed very sensitive to the presence of side walls: on Fig. 8b, $C_R$ are quickly confined within a constant distance of the plume axis, due to the small size of the box. By the same token, $C_R$ uplift stops, due to the upper surface effects, much sooner than the plume thermal anomaly (Fig. 8b). Therefore in a high Prandtl number fluid, a thermal plume cannot be adequately characterized by the rotation centers.

We also calculated the vorticity fields, which do not depend on the reference frame in which they are calculated (Figs. 16c, 9c). As already pointed out by Tanny and Shlien (1985), vorticity is diffuse and there is no localized vorticity core. However, the vorticity field presents two extremum in the plume head (“vorticity centers”). Their

ascent velocity is found to be well predicted by $V_z(A)$ (Fig. 9c). If we now compute the stream function in the reference frame of the plume head, i.e. with the velocity field $V - V_z(A)$ (Fig. 16b), we see that the plume thermal anomaly is contained into a closed streamline, and that the two stream function extrema are very close to the vorticity centers. The plume radius determined on the thermal anomaly is now very close to the plume radius defined as the maximum radial extension of the closed streamline.

Figure 16 shows the relative positions of the key-features of the plume thermal anomaly and dynamics. Within our parameter range, the diffusive growth of the plume cannot be neglected compared to advection. It is in particular required to compute the “true” starting plume ascent velocity.

## 5 Conclusions

The simultaneous measurements of the temperature and velocity fields have allowed us to characterize the anatomy of laminar starting thermal plumes in fluids with high Prandtl numbers and strongly temperature-dependent viscosity.

The plume thermal anomaly can be efficiently traced by means of thermochromic liquid crystals “isotherms”. However, a given isotherm will remain useful only for a certain time, depending on the physical properties of the fluid and the heating power of the source. We show that the thermal structure of the plume stem is well predicted by the constant viscosity model of Batchelor (1954) for viscosity ratios up to 50.

It is also possible to follow the ascent and growth of the plume through its velocity field. Key measurements are then the velocity of the stagnation point located on the top of the plume head $V_A$, and the velocity of the particles at this point $V_z(A)$. The former will give the same information
as following dye (which underlines the maximum deformation) or following shadowgraph (which underlines the maximum temperature gradient), while $V_z(A)$ gives the same value as the uplift of the maximum vorticity centers. The latter constitutes the “true” plume head velocity. The difference $V_A - V_z(A)$ represents plume head growth by thermal diffusion. We further show that the stem velocity structure is quite sensitive to the temperature dependence of the viscosity and to confinement.

More work is now needed to study in details the influence of the heat source geometry and side walls on the scalings, and to extend them to large viscosity ratios.

**Acknowledgments** This work has benefited from discussions with Neil Ribe, Eric Mittelsteadt, Peter van Keken, and Béatrice Guerrier. It was funded by program DyETI of INSU/CNRS, the French ANR “BEGDY” and the collaboration between IPGP in Paris and ERI in Tokyo. The manuscript has been improved, thanks to the constructive comments of two anonymous reviewers.

## Appendix 1: Viscosity laws

The fluid viscosity was measured between 15 and 65°C using a rheometer RS600 from ThermoRheo with a cylindrical geometry. We obtained for fluid 1 (RUN 1 and 2):

$$ \eta(T) = \exp(8.283 \times 10^{-4}T^2 - 1.551 \times 10^{-1}T + 5.251) \quad (15) $$

and for fluid 2 (RUN 3 and 4):

$$ \eta(T) = \exp(-7.086 \times 10^{-2}T + 3.408) \quad (16) $$

where T is the temperature in °C.

## Appendix 2: TLC calibration

The TLCs were calibrated by imposing a stable linear temperature gradient to a 10 cm layer of the experimental fluid: the light intensity across a vertical cross-section of the tank is then plotted as a function of the vertical profile of temperature (Fig. 17). Each type of TLC results in an intensity peak for a given temperature. Each bright stripe has a finite thickness because the polymeric capsules enclosing the TLCs introduce some scatter around the temperature at which they respond. The peak maximum gives us the isotherm value; and the width of the peak at half maximum intensity (or bandwidth) gives the local temperature gradient (Table 3). Peaks are well defined for aqueous slurries in sugar syrups, where the precision of the isotherms values reaches ±0.1°C, and the bandwidth is about 0.2°C. The use of the same TLC technique in silicone oils of similar viscosity show much poorer results with bandwidth reaching 1.5°C (Limare et al. 2008; Vatteville

<figure><img src="image_56.png" /><figcaption>**Fig. 17** Light intensity as a function of temperature for the KMN liquid crystal type. Each TLC type results in a peak centered on a given temperature</figcaption></figure>

<table><caption>**Table 3** Characteristics of the liquid crystals (KMN** are from Japan Capsular Products, Inc.; and BM/** are from Hallcrest, Inc.) used in the experiments with their calibration values</caption><thead><tr><th>Product name</th><th>Peak value (°C)</th><th>Bandwidth (°C)</th></tr></thead><tbody><tr><td>KMN37-39</td><td>37.5</td><td>0.20</td></tr><tr><td>KMN34-36</td><td>34.5</td><td>0.20</td></tr><tr><td>KMN31-33</td><td>31.4</td><td>0.14</td></tr><tr><td>KMN28-30</td><td>28.4</td><td>0.17</td></tr><tr><td>KMN25-27</td><td>25.3</td><td>0.23</td></tr><tr><td>BM/24C2W/S40</td><td>24.35</td><td>0.30</td></tr><tr><td>BM/27C2W/S40</td><td>27.0</td><td>0.35</td></tr><tr><td>BM/31C2W/S40</td><td>31.10</td><td>0.25</td></tr><tr><td>BM/35C2W/S40</td><td>35.0</td><td>0.25</td></tr><tr><td>BM/40C2W/S40</td><td>39.5</td><td>0.25</td></tr></tbody></table>

et al. 2009). This may be due to the response of the TLCs coating, which is chemically different in the two cases.

## Appendix 3: Estimate of the optical distortion

Since a thermal anomaly acts as a divergent lens, we estimated its effect on the TLC isotherms positions. For this purpose we have used the formula given by Laudenbach and Christensen (2001) for the deflection angle of a laser beam through an axisymmetric thermal plume for a given distribution of the refractive index:

$$ \alpha(p) = 2p \int_{p_0}^{\infty} \frac{1}{n(p^2 - p_0^2)} \frac{dn}{dp} dp \quad (17) $$

where p is defined as $p = rn(r)$ and $p_0 = rn_0$; and $n(r)$ is the refractive index distribution and $n_0$ the refractive index of the medium at ambient temperature.

The distribution of the refractive index of the plume is known from the DI results. Since the isotherms are
<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /><figcaption>**Fig. 18** Optical distortion induced by the thermal anomaly on the TLC images: a horizontal temperature profiles, b deflection angle, and c shift on the images, as a function of the radial coordinate</figcaption></figure>

visualized in a plane passing through the center of the plume, we therefore considered that the deflection angle of the scattered light from the isotherms was equal to half of the total deflection angle of a virtual beam passing through the entire plume. In our case, this angle remains smaller than $10^{-3}$ radian. This implies that the associated shift in location of a given isotherm would remain smaller than $5 \times 10^{-2}$ mm (Fig. 18).

## References

Batchelor GK (1954) Heat convection and buoyancy effects in fluids. Q J R Met Soc 80: 339–358

Boxman RL, Shlien DJ (1978) Interferometric measurement technique for the temperature field of axisymmetric buoyant phenomena. Appl Opt 17: 2788–2793

Chay A, Shlien DJ (1986) Scalar field measurements of a laminar starting plume cap using digital processing of interferograms. Phys Fluids 29: 2358–2366

Couliette DL, Loper DE (1995) Experimental, numerical and analytical models of mantle starting plumes. Phys Earth Planet Inter 92: 143–167

Dabiri D, Gharib M (1991) Digital particle image thermometry: the method and implementation. Exp Fluids 11: 77–86

Dabiri D (2009) Digital liquid crystal particle thermometry/velocity (DLCPT/V)—a review. Exp Fluids 46(2): 191–241

Davaille A, Jaupart C (1993) Transient high-Rayleigh number thermal convection with large viscosity variations. J Fluid Mech 253: 141–166

Davaille A, Vatteville J (2005) On the transient nature of mantle plumes. Geophys Res Lett 32: L14309 doi:10.1029/2005GL023029

Davaille A, Limare A (2007) Laboratory studies on mantle convection. In: Bercovici D, Schubert G (eds) Treatise of geophysics. Elsevier, Amsterdam, pp 89–165

Davaille A, Androvandi S, Vatteville J, Limare A, Vidal V, Lebars M (2008) Thermal boundary layer instabilities in viscous fluids. In: ISFV13—13th International symposium on flow visualization, July 1–4, 2008, Nice, France, p 12, available at http://www.ipgp.fr/~limare/317.pdf

Fuji T (1963) Theory of the steady laminar natural convection above a horizontal line heat source and a point heat source. Int J Heat Mass Transf 6: 597–606

Fujisawa N, Adrian R (1999) Three-dimensional temperature measurement in turbulent thermal convection by extended range scanning liquid crystal thermometry. J Vis 1: 355–364

Gebhart B, Pera L, Schorr WA (1970) Steady laminar natural convection plumes above a horizontal line heat source. Int J Heat Mass Transf 13: 161–171

Griffiths RW (1986) Thermals in extremely viscous fluids, including the effects of temperature dependent viscosity. J Fluid Mech 166: 115–138

Griffiths RW, Campbell IH (1990) Stirring and structure in mantle starting plumes. Earth Planet Sci Lett 99: 66–78

Happel J, Brenner H (1973) Low Reynolds number hydrodynamics. Noordhoff, Leyden

Kaminski E, Jaupart C (2003) Laminar starting plumes in high-Prandtl-number fluids. J Fluid Mech 478: 287–298

Kumagai I (2002) On the anatomy of mantle plumes: effect of the viscosity ratio on entrainment and stirring. Earth Planet Sci Lett 198: 211–224

Kumagai I, Davaille A, Kurita K (2007) On the fate of thermal plumes at density interface, earth planet. Sci Lett 254: 180–193

Kumagai I, Davaille A, Kurita K, Stutzmann E (2008) Mantle plumes: thin, fat, successful, or failing? Constraints to explain hot spot volcanism through time and space, Geophys. Res Lett 35: L16301 http://www dx.doi.org/10.1029/2008GL035079

Kelley JG, Hargreaves RA (1970) A rugged inexpensive shearing interferometer Appl Opt 9: 948–952

Laudenbach N, Christensen UR (2001) An optical method for measuring temperature in laboratory models of mantle plumes. Geophys J Int 145: 528–534

Le Bars M, Davaille A (2002) Stability of thermal convection in two superimposed miscible viscous fluids. J Fluid Mech 471: 339–363

Limare A, Kumagai I, Vatteville J, Davaille A (2008) Thermal plumes visualisation: differential interferometry versus thermochromic liquid crystals, ISFV13—13th international symposium on flow visualization, July 1–4, 2008, Nice, France, p 12, available at http://www.ipgp.fr/~limare/318.pdf
Morris S, Canright DR (1984) A boundary-layer analysis of Benard convection with strongly temperature-dependent viscosity. *Phys Earth Planet Inter* 36: 355-373

Morton BR, Taylor GI, Turner JS (1956) Turbulent gravitational convection from maintained and instantaneous sources. *Proc R Soc Lond A* 234:1-23

Moses E, Zocchi G, Libchaber A (1993) An experimental study of laminar plumes. *J Fluid Mech* 251: 581-601

Olson P, Singer H (1985) Creeping plumes. *J Fluid Mech* 158: 511-531

Olson P, Schubert G, Anderson C (1993) Structure of axisymmetric plumes. *J Geophys Res* 98: 6829-6844

Pretzler G, Jäger H, Neger T (1993) High-accuracy differential interferometry for the investigation of phase objects. *Meas Sci Technol* 4: 649-658

Rhee H, Koseff J, Street R (1984) Flow visualization of a recirculating flow by rheoscopic and liquid crystal techniques. *Exp Fluids* 2: 57-64

Rogers MC, Morris SW (2009) Natural versus forced convection in laminar starting plumes. *Phys Fluids* 21: 083601

Sernas V, Fletcher LS (1970) A schlieren interferometer method for heat transfer studies. *J Heat Transf* 92: 202-204

Shlien DJ (1976) Some laminar thermal and plume experiments. *Phys Fluids* 19: 1089-1098

Shlien DJ, Boxman RL (1979) Temperature field measurement of an axisymmetric laminar plume. *Phys Fluids* 22(4): 631-634

Shlien DJ, Brosh A (1979) Velocity field measurements of a laminar thermal. *Phys Fluids* 22: 1044-1053

Shlien DJ, Boxman RL (1981) Laminar starting plume temperature field measurement. *Intl J Heat Mass Transf* 24: 919-930

Tanny J, Shlien DJ (1985) Velocity field measurements of a laminar starting plume. *Phys Fluids* 28: 1027-1032

Turner JS (1962) The starting plume in neutral surroundings. *J Fluid Mech* 13:356-368

van Keken PE (1997) Evolution of starting mantle plumes: a comparison between laboratory and numerical models. *Earth Planet Sci Lett* 148: 1-14

van Keken P, King S, Schmeling H, Christensen U, Neumeister D, Doin M-P (1997) A comparison of methods for the modeling of thermochemical convection. *J Geophys Res* 102(B10): 22477-22495

Vasquez PA, Perez AT, Castellanos A (1996) Thermal and electro-hydrodynamics plumes. A comparative study. *Phys Fluids* 8: 2091-2096

Vatteville J, van Keken P, Limare A, Davaille A (2009) Starting laminar plumes: comparison of laboratory and numerical modeling. *Geochem Geophys Geosyst* 10: Q12013 doi:10.1029/2009GC002739

Vest CM (1975) Interferometry of strongly refracting axisymmetric phase objects. *Appl Opt* 14: 1601-1606

Whitehead JA, Luther DS (1975) Dynamics of laboratory diapir and plume models. *J Geophys Res* 80: 705-717

Whittaker RJ, Lister JR (2006a) Steady axisymmetric creeping plumes above a planar boundary. Part I: a point source. *J Fluid Mech* 567: 361-378

Whittaker RJ, Lister JR (2006b) Steady axisymmetric creeping plumes above a planar boundary. Part II: a distributed source. *J Fluid Mech* 567: 379-397

Whittaker RJ, Lister JR (2008) The self-similar rise of a buoyant thermal in very viscous flow. *J Fluid Mech* 606: 295-324

Willert C, Gharib M (1991) Digital particle image thermometry. *Exp Fluids* 10: 181-193

Worster MG (1986) The axisymmetric laminar plumes: asymptotic solution for large Prandtl number. *Stud Appl Math* 75: 139-152
<figure><img src="image_60.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Physics of the Earth and Planetary Interiors

journal homepage: www.elsevier.com/locate/pepi

<figure><img src="image_61.png" /></figure>

# At least three scales of convection in a mantle with strongly temperature-dependent viscosity

Sophie Androvandi<sup>a</sup>, Anne Davaille<sup>b,*</sup>, Angela Limare<sup>a</sup>, Aurélie Foucquier<sup>a,c</sup>, Catherine Marais<sup>a,d</sup>

## ARTICLE INFO

**Article history:**

Received 31 January 2011

Received in revised form 8 July 2011

Accepted 13 July 2011

Available online 23 July 2011

Edited by Keke Zhang

**Keywords:**

Mantle convection

Plumes

Subduction

Asthenosphere

Laboratory experiments

## ABSTRACT

We studied the convective patterns developing at high Rayleigh numbers and intermediate viscosity ratios in a fluid with a strongly temperature-dependent viscosity. Within this sluggish lid regime, three different scales of convection develop. The largest convective scale is cellular, with cold downwelling sheets of viscous fluid encasing hotter, less viscous, parts of the tank. Within each of those cells develop several (typically 3–7) hot 3D upwelling plumes. Upon impinging under the cold thermal boundary layer, each plume in turn generates locally a small ring of cold material which does not reach the bottom of the tank.

Applying those results to the Earth's mantle, we suggest that the large-scale features of mantle convection and the co-existence of several scales of convection, i.e. slabs and plumes, are produced by thermal convection in a mantle with high Rayleigh number and a strongly temperature-dependent viscosity material. This generates a convective pattern in two large-scale cells : the Pacific and the Indo-Atlantic "boxes". Our experiments further suggest that what has been named the two hot superplumes, i.e. the two seismically slow regions encased within the subduction rings, are in fact each constituted of several hot instabilities. Moreover, plumes impacts under the lithosphere should be surrounded by cold rings of small extent. The asthenosphere appears therefore as the graveyard of both small-scale convection and hot plumes generated in the system.

© 2011 Elsevier B.V. All rights reserved.

## 1. Introduction

Owing to the presence of continents and to the chemical heterogeneity and complex rheology of mantle material, several scales of convection co-exist in the Earth's mantle, producing cold linear slabs, mid-ocean ridges, 3D superswells and hot spots. Tomographic images (Fig. 1) show that the Earth's mantle and the hot spots population is divided in two large boxes by subduction (Masters et al., 1982; Molnar and Stock, 1987; Weinstein and Olson, 1989; Grand et al., 1997). In terms of convection, mantle cold downwellings are therefore sheet-like and are delimiting two large-scale cells, in which several hot instabilities are generated to produce hot spots and superswells (e.g. Davaille et al., 2002, 2005; Jellinek et al., 2003; Gonnermann et al., 2004; Schubert et al., 2004; McNamara and Zhong, 2005). However, even though these features have individually been generated and studied in numerical and laboratory experiments, the exact conditions for

their co-existence in a self-consistent convective model have remained elusive.

Numerical studies (e.g. Bercovici, 1996, 1998, 2003; Moresi and Solomatov, 1995; Trompert and Hansen, 1998; Gurnis et al., 1998, 2000; Tackley, 2000a,b,c; Ricard et al., 2001; Stein et al., 2004) have shown that the presence of rigid plates and the strong localization of deformation typical of plate tectonics on Earth are obtained only in fluids with complex rheology (yield-stress, strongly temperature-dependent, and self-lubricating). Nevertheless, if we consider the large-scale feature, that is the upper cold more viscous thermal boundary layer of the mantle (i.e. the lithosphere) participating in the convective overturn of the whole mantle, plate tectonics convection appears like an example of the mobile-lid regime featured in convection in a strongly temperature-dependent viscosity fluid at intermediate viscosity ratios (e.g. Weinstein and Christensen, 1991; Solomatov, 1995; Moresi and Solomatov, 1995; Ratcliff et al., 1997; Kameyama and Ogawa, 2000; Stemmer et al., 2006). But the question is then whether this rheology is sufficient to produce several scales of convection, and under which conditions.

We have therefore run laboratory experiments to explore the influence of a strongly temperature-dependent rheology on the

<footnote><sup>a</sup> Dynamique des Fluides Géologiques, IPGP/Univ. Paris 7/CNRS, Paris, France</footnote>

<footnote><sup>b</sup> Laboratoire FAST, CNRS/UPMC/U-PSud, Orsay, France</footnote>

<footnote><sup>c</sup> LPMC, Ecole Polytechnique/CNRS, Palaiseau, France</footnote>

<footnote><sup>d</sup> PMMH, ESPCI/CNRS/UPMC Paris, France</footnote>

<footnote>* Corresponding author. Tel.: +33 0 1 69 15 80 56; fax: +33 0 1 69 15 80 60.
E-mail address: davaille@fast.u-psud.fr (A. Davaille).</footnote>
<figure><img src="image_62.png" /><figcaption>**Fig. 1.** Variations of S-waves seismic velocities at 2700 km depth, extracted from tomographic model TXBM (Grand, 2002). The circles show the hot spots on the Earth surface. The subduction ring around the Pacific seems to divide the mantle in two “boxes”: Indo-Atlantic and Pacific. However, note that the slow seismic velocity anomaly (in orange-red) under the Pacific is not uniform but presents at least 5 local maxima. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

pattern of thermal convection when the intensity of convection is mantle-like (i.e. high Rayleigh number, $10^6–10^8$). New visualization techniques (Davaille and Limare, 2007) allowed us to measure the detailed structure of the temperature and velocity fields. This work complements previous studies run either in 2D geometry (Christensen, 1984a,b, 1985; Hansen and Yuen, 1993; Moresi and Solomatov, 1995; Kameyama and Ogawa, 2000), or in 3D geometry at lower Rayleigh numbers (Olson et al., 1988; Richter et al., 1983; White, 1988; Ogawa et al., 1991; Christensen and Harder, 1991; Weinstein and Christensen, 1991; Giannandrea and Christensen, 1993; Tackley et al., 1994; Tackley, 1996; Ratcliff et al., 1997; Yoshida and Kageyama, 2006; Stemmer et al., 2006), or with no strongly temperature-dependent viscosity (Davies, 2005; Davies and Davies, 2009), or with less detailed information on the temperature and velocity fields (Weeraratne and Manga, 1998; Manga and Weeraratne, 1999; Manga et al., 2001; Schaeffer and Manga, 2001). We were able to expand the regime diagram of the different convective regimes, which then allows us to explain in a same framework a number of terrestrial observations. The laboratory experiments are described in Section 2, while the implications for the style of mantle convection, and the nature of the asthenosphere are presented in Sections 3 and 4.

## 2. Laboratory experiments

### 2.1. Experimental set-up

Layers of sugar syrup were cooled from above and heated from below at constant temperatures in square-based (30 cm × 30 cm) plexiglas tanks of aspect ratios 2 and 5. The top and bottom boundaries are copper exchangers within which thermostated liquids (cold ethanol and hot water) were circulated. Each experiment was run during several days in order to reach quasi-steady state. To simultaneously measure in situ the temperature and velocity fields, the fluid was seeded with glass spheres (10-µm diameter) and four types of thermochromic liquid crystals (TLC). Each type of TLC brightens at a different temperature. So, we visualize four isotherms which appear as bright lines (39.50 ± 0.25, 31.10 ± 0.25, 24.35 ± 0.30, 10.70 ± 0.60 °C, labelled, respectively A, B, C and D in the figures) when vertical cross-sections of the tank are illuminated

with a 532-µm laser sheet (Davaille and Limare, 2007). Images were taken every half-second by a CCD camera (1280 × 1020 pixels). Thanks to the glass spheres, the velocity field was calculated by cross-correlation between images, using the PIV package Davis from LaVision. We can also scan the tank with the laser sheet in order to determine the 3D structure of the instabilities. In addition, a vertical probe (labelled “T” on the figures) containing 14 thermocouples was used to measure a vertical temperature profile at a given location. Thermocouples were also placed in the lower and the upper plates. We verified that temperature was homogeneous in each copper plate within 0.05 °C. All temperatures were read every 30 s through a scanning voltmeter connected to a computer. The thermocouples allow us to take very long time series (e.g. Fig. 7).

Twenty-six experiments were performed, using three fluids (differing by the amount of solid matter dissolved in water). Their physical properties are given in Annex 1 and Tables 2 and 3. All fluids are in the limit $v/κ ≫ 1$ relevant to the Earth, where $κ$ is the thermal diffusivity and $v$ is the kinematic viscosity (i.e. momentum diffusivity). The intensity of convection is related to the global Rayleigh number, which compares the driving thermal buoyancy forces to the resisting effects of thermal diffusion and viscous dissipation:

$$Ra(v) = \frac{\alpha \rho g \Delta T H^3}{\kappa v} \qquad (1)$$

where $H$ is the layer depth, $ΔT$ is the temperature difference applied across it, $α$ is the thermal expansivity, $ρ$ is the density, $g$ is the gravity acceleration, and $v$, the kinematic viscosity. Since the viscosity of the fluids depends strongly on temperature, the temperature at which the viscosity is taken in Eq. (1) has to be defined; several choices are possible: $Ra_{hot} = Ra(v_{hot})$ with the viscosity taken at the temperature of the bottom hot boundary, $Ra_{cold} = Ra(v_{cold})$ with the viscosity taken at the temperature of the top cold boundary, $Ra_{1/2} = Ra(v_{1/2})$ with the viscosity taken at the averaged temperature between the two boundaries, and $Ra_m = Ra(v_m)$ with the viscosity taken at the mean temperature of the convecting bulk interior of the tank.

In the limit $v/κ ≫ 1$, the system dynamics depend on only two parameters, one of the Rayleigh numbers, -here we shall take $Ra_{hot}^-$, and the viscosity ratio across the tank defined as the ratio of the viscosity of the cold (top) fluid to the hot (bottom) one $γ = v_{cold}/v_{hot}$. The
parameter range of our experiments covers: $10^5 < Ra_{hot} < 10^8$ and $7 < \gamma < 4200$ (Table 1).

## 2.2. Regime diagram

Fig. 2 summarizes the different regimes, based on the laboratory and 3D numerical experiments reported in the literature (for a recent review, see Davaille and Limare, 2007), and on our new data. Their systematics, especially in the mantle perspective (i.e. high $Ra_{hot}$) are best explained in terms of $Ra_{cold}$ and $Ra_{hot} = \gamma Ra_{cold}$, the Rayleigh numbers based on the local properties of the two thermal boundary layers.

### 2.2.1. Nearly constant viscosity

The convective patterns for constant viscosity are well known for both rigid and stress-free boundary conditions (e.g. Krishnamurti, 1970): 2D rolls are stable up to $Ra = 20000$. Then steady cellular patterns develop to become time-dependent above $Ra = 10^5$. Between $Ra = 10^6$ and $10^7$, the cellular character of the pattern is lost in the bulk of the fluid, and hot and cold plumes develop chaotically both in space and time from the TBLS (e.g. Guillou and Jaupart, 1995; Weeraratne and Manga, 1998; Schmalzl et al., 2002; Davies, 2005; Davies and Davies, 2009; Androvandi, 2009). However, close to the boundaries, plumes appear to originate from triple junctions of ridges imbedded into the TBLS (Travis et al., 1990; Houseman, 1990; Parmentier and Sotin, 2000; Davies and Davies, 2009).

### 2.2.2. High viscosity ratios $\gamma > 10^4$: stagnant lid regime

For $\gamma > 10^3$ (Nataf and Richter, 1982; Richter et al., 1983) to $2 \times 10^4$ (Giannandrea and Christensen, 1993; Solomatov, 1995; Stemmer et al., 2006), the “stagnant lid” regime prevails whereby convection develops below the coldest and most viscous part of the upper TBL, which remains immobile (Davaille and Jaupart, 1993). This is consistent with the fact that $Ra_{cold}$ is subcritical (or very close to it): therefore instabilities at the coldest temperature are not able to develop in a finite time. The convective pattern below this “stagnant lid” depends on the local Rayleigh number in the

sublayer and is similar to nearly constant viscosity patterns (Trompert and Hansen, 1998).

### 2.2.3. Intermediate viscosity ratios and low $Ra$

At low $Ra$, 3D-flow patterns such as squares and hexagons are stabilized by the temperature-dependence of viscosity. They exhibit plume-like ascending flow, while the descending flow becomes more similar to a sheet as the viscosity contrast increases. This is observed for rigid boundaries (White, 1988), as well as for free boundaries (Christensen and Harder, 1991; Ogawa et al., 1991). For $Ra_{hot} = 1.5 \times 10^5$ and $\gamma = 50$, Christensen and Harder (1991) found an irregular and time-dependent spoke pattern, with a network of sinking cold sheets at shallow depth breaking up into single falling jets at greater depth. The rising flow also forms a network of ridges near the bottom, which at greater height breaks into plumes. Our experiments with rigid boundaries show that this regime extends to at least $Ra_{hot} \sim 5 \times 10^6$ and $\gamma \sim 4000$ (Fig. 2). A closer look on the convective pattern shows the asymmetry between the hot and cold TBLS: the cold downwellings impose an overall cellular character with cell aspect ratio ranging between 0.7 and 1.5 (Fig. 3b; see also Giannandrea and Christensen, 1993). Moreover, the cold currents present no or only weak time-dependence (Fig. 3a and c), contrary to the hot currents which move around and can even disappear for a while (Fig. 3a and c). However, each cold cell contains only 1 hot plume at a time. This is consistent with the fact that both $Ra_{hot}$ and $Ra_{cold}$ correspond to the cellular or spokes regimes (black triangles on Fig. 2).

For a stress-free upper boundary, data was obtained only for $Ra < 3 \times 10^5$ (Weinstein and Christensen, 1991; Giannandrea and Christensen, 1993): for $\gamma < 5000$, the upper surface material remains mobile and large aspect ratio cells (1.5–3.5) were observed with slab-like stable cold downwellings and narrower hot plumes whose position changes with time. The existence of this large-scale, sluggish lid (or “transitional”) mode is in agreement with earlier predictions (Stengel et al., 1982; Jaupart and Parsons, 1985). Numerical simulations in 2D also reported this regime with large aspect ratio cells, expanding to higher $Ra_{hot}$ (Moresi and Solomatov, 1995; Kameyama and Ogawa, 2000).

<table><caption>**Table 1**
Experimental parameters in steady-state.</caption><thead><tr><th>Exp.</th><th>T<sub>cold</sub></th><th>T<sub>hot</sub></th><th>T<sub>m</sub></th><th>γ</th><th>Ra<sub>1/2</sub></th><th>Ra<sub>m</sub></th><th>Ra<sub>hot</sub></th><th>Ra<sub>cold</sub></th><th>fluid</th></tr></thead><tbody><tr><td>LRa1</td><td>6.3</td><td>51.4</td><td>37.1</td><td>103</td><td>3.0 × 10<sup>6</sup></td><td>6.7 × 10<sup>6</sup></td><td>2.1 × 10<sup>7</sup></td><td>2.1 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa2</td><td>11.3</td><td>51.6</td><td>38.2</td><td>54</td><td>3.5 × 10<sup>6</sup></td><td>6.6 × 10<sup>6</sup></td><td>1.9 × 10<sup>7</sup></td><td>3.6 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa3</td><td>16.0</td><td>51.9</td><td>39.0</td><td>31</td><td>4.0 × 10<sup>6</sup></td><td>6.3 × 10<sup>6</sup></td><td>1.8 × 10<sup>7</sup></td><td>5.7 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa4</td><td>20.7</td><td>52.0</td><td>40.4</td><td>18</td><td>4.3 × 10<sup>6</sup></td><td>6.2 × 10<sup>6</sup></td><td>1.5 × 10<sup>7</sup></td><td>8.6 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa5</td><td>-3.9</td><td>51.0</td><td>35.0</td><td>425</td><td>2.1 × 10<sup>6</sup></td><td>6.7 × 10<sup>6</sup></td><td>2.5 × 10<sup>7</sup></td><td>5.9 × 10<sup>4</sup></td><td>1</td></tr><tr><td>LRa6</td><td>10.0</td><td>50.8</td><td>38.0</td><td>61</td><td>3.2 × 10<sup>6</sup></td><td>6.6 × 10<sup>6</sup></td><td>1.8 × 10<sup>7</sup></td><td>3.0 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa7</td><td>-18.1</td><td>50.1</td><td>34.1</td><td>3484</td><td>1.1 × 10<sup>6</sup></td><td>7.6 × 10<sup>6</sup></td><td>2.9 × 10<sup>7</sup></td><td>7.8 × 10<sup>3</sup></td><td>1</td></tr><tr><td>LRa8</td><td>-6.0</td><td>50.0</td><td>34.9</td><td>541</td><td>1.8 × 10<sup>6</sup></td><td>6.8 × 10<sup>6</sup></td><td>2.4 × 10<sup>7</sup></td><td>4.4 × 10<sup>4</sup></td><td>1</td></tr><tr><td>LRa9</td><td>16.2</td><td>50.7</td><td>38.7</td><td>29</td><td>3.6 × 10<sup>6</sup></td><td>5.9 × 10<sup>6</sup></td><td>1.6 × 10<sup>7</sup></td><td>5.6 × 10<sup>5</sup></td><td>1</td></tr><tr><td>LRa10</td><td>-13.2</td><td>50.0</td><td>33.9</td><td>1673</td><td>1.3 × 10<sup>6</sup></td><td>7.0 × 10<sup>6</sup></td><td>2.7 × 10<sup>7</sup></td><td>1.6 × 10<sup>4</sup></td><td>1</td></tr><tr><td>LRa11</td><td>-18.2</td><td>50.9</td><td>34.3</td><td>4078</td><td>1.1 × 10<sup>6</sup></td><td>7.9 × 10<sup>6</sup></td><td>3.1 × 10<sup>7</sup></td><td>7.7 × 10<sup>3</sup></td><td>1</td></tr><tr><td>Lra12</td><td>1.0</td><td>52.0</td><td>36.9</td><td>142</td><td>7.1 × 10<sup>6</sup></td><td>1.8 × 10<sup>7</sup></td><td>5.3 × 10<sup>7</sup></td><td>3.7 × 10<sup>5</sup></td><td>3</td></tr><tr><td>Lra13</td><td>-3.5</td><td>38.4</td><td>24.2</td><td>99</td><td>2.3 × 10<sup>6</sup></td><td>4.6 × 10<sup>6</sup></td><td>1.7 × 10<sup>7</sup></td><td>1.7 × 10<sup>5</sup></td><td>3</td></tr><tr><td>Lra14</td><td>8.5</td><td>47.7</td><td>34.3</td><td>41</td><td>6.3 × 10<sup>6</sup></td><td>1.1 × 10<sup>7</sup></td><td>3.1 × 10<sup>7</sup></td><td>7.5 × 10<sup>5</sup></td><td>3</td></tr><tr><td>Lra15</td><td>6.8</td><td>49.4</td><td>35.4</td><td>57</td><td>6.8 × 10<sup>6</sup></td><td>1.3 × 10<sup>7</sup></td><td>3.8 × 10<sup>7</sup></td><td>6.6 × 10<sup>5</sup></td><td>3</td></tr><tr><td>HRa1</td><td>14.5</td><td>52.0</td><td>39.5</td><td>14</td><td>2.2 × 10<sup>7</sup></td><td>3.4 × 10<sup>7</sup></td><td>6.8 × 10<sup>7</sup></td><td>4.8 × 10<sup>6</sup></td><td>2</td></tr><tr><td>HRa2</td><td>-5.9</td><td>51.2</td><td>35.7</td><td>118</td><td>1.5 × 10<sup>7</sup></td><td>3.9 × 10<sup>7</sup></td><td>1.0 × 10<sup>8</sup></td><td>8.5 × 10<sup>5</sup></td><td>2</td></tr><tr><td>HRa3</td><td>22.4</td><td>53.1</td><td>42.2</td><td>7</td><td>2.5 × 10<sup>7</sup></td><td>3.3 × 10<sup>7</sup></td><td>5.9 × 10<sup>7</sup></td><td>7.9 × 10<sup>6</sup></td><td>2</td></tr><tr><td>HRa4</td><td>3.9</td><td>52.4</td><td>37.1</td><td>42</td><td>2.0 × 10<sup>7</sup></td><td>3.8 × 10<sup>7</sup></td><td>9.0 × 10<sup>7</sup></td><td>2.2 × 10<sup>6</sup></td><td>2</td></tr><tr><td>syr1m</td><td>-15.8</td><td>52.3</td><td>34.8</td><td>3106</td><td>2.1 × 10<sup>6</sup></td><td>1.1 × 10<sup>7</sup></td><td>4.4 × 10<sup>7</sup></td><td>1.4 × 10<sup>4</sup></td><td>1</td></tr><tr><td>syr2m</td><td>-15.6</td><td>52.3</td><td>36.0</td><td>3000</td><td>2.1 × 10<sup>6</sup></td><td>1.3 × 10<sup>7</sup></td><td>4.4 × 10<sup>7</sup></td><td>1.5 × 10<sup>4</sup></td><td>1</td></tr><tr><td>syr3m</td><td>-17.0</td><td>52.9</td><td>42.1</td><td>3984</td><td>2.0 × 10<sup>6</sup></td><td>2.2 × 10<sup>7</sup></td><td>4.7 × 10<sup>7</sup></td><td>1.2 × 10<sup>4</sup></td><td>1</td></tr><tr><td>j1m16</td><td>-1.7</td><td>47.9</td><td>34.3</td><td>236</td><td>1.7 × 10<sup>5</sup></td><td>5.3 × 10<sup>5</sup></td><td>1.6 × 10<sup>6</sup></td><td>6.8 × 10<sup>3</sup></td><td>1</td></tr><tr><td>j1m17</td><td>-1.8</td><td>47.9</td><td>33.8</td><td>241</td><td>1.7 × 10<sup>5</sup></td><td>5.1 × 10<sup>5</sup></td><td>1.6 × 10<sup>6</sup></td><td>6.6 × 10<sup>3</sup></td><td>1</td></tr><tr><td>j1m20</td><td>-15.8</td><td>47.7</td><td>32.5</td><td>2281</td><td>9.7 × 10<sup>4</sup></td><td>5.8 × 10<sup>5</sup></td><td>2.0 × 10<sup>6</sup></td><td>888</td><td>1</td></tr><tr><td>j1m22</td><td>7.6</td><td>48.3</td><td>35.3</td><td>64</td><td>2.4 × 10<sup>5</sup></td><td>4.8 × 10<sup>5</sup></td><td>1.3 × 10<sup>6</sup></td><td>2.1 × 10<sup>4</sup></td><td>1</td></tr><tr><td>j1m24</td><td>17.1</td><td>48.4</td><td>37.1</td><td>19</td><td>2.9 × 10<sup>5</sup></td><td>4.3 × 10<sup>5</sup></td><td>1.1 × 10<sup>6</sup></td><td>5.4 × 10<sup>4</sup></td><td>1</td></tr></tbody></table>
<figure><img src="image_63.png" /><figcaption>**Fig. 2.** Convective pattern as a function of the Rayleigh numbers at the top and bottom boundaries for 3D experiments in a cartesian box. Only experiments allowing to determine the convective pattern have been considered. Numbers in italic indicate the viscosity ratio. Regimes in the constant viscosity cases are indicated in the top and right bars. Open symbols designate experiments run with free surfaces (mostly numerical simulations) and solid symbols designate experiments run with solid boundaries. Black: our new data; light blue: Trompert and Hansen (1998); red: Tackley et al. (1994), Tackley (1996) and Ratcliff et al. (1997); green: Weinstein and Christensen (1991); dark blue: Guillou and Jaupart (1995); orange: Schaeffer and Manga (2001). The dashed grey border between the steady and time-dependent whole layer regimes is after White (1988). The thick dashed grey line shows the theoretical limit of Solomatov (1995) between the low-viscosity ratio regime and the transitional regime. The thick solid grey line between the whole layer and the stagnant lid regimes is after Solomatov (1995); the orange line is from the steady state calculations of Ogawa et al. (1991); the dashed thick purple line from Nataf and Richter (1982) and Richter et al. (1983); and the green thick line after Giannandrea and Christensen (1993). The limit between no convection and convection is drawn after Stengel et al. (1982). (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 2.2.4. Intermediate viscosity ratios and high Ra

In 3D, at high Ra but intermediate viscosity ratios (solid black disks on Fig. 2, typically for $\gamma > 50$), convective instabilities are observed on several time- and space-scales. Fig. 4a presents a typical vertical cross-section of the tank: a cluster of hot instabilities (plumes) are surrounded by cold downwellings. 3D scanning of the tank further shows that the latter constitute the sides of polygonal cells. Their associated temperature anomaly is imaged at least through half the tank depth, while the velocity field reveals a structure of sinking sheets all the way to the bottom of the tank (Fig. 4b). The time evolution (Fig. 5) of the TBLS reveals that the cold downwellings are quite stable through time. In the plume cluster, the plume closest to the downwelling never disappears completely, even though its position can meander through time. By comparison, the plumes inside the cluster are much more time-dependent. In this range of parameters, Fig. 2 shows that the cold downwellings (characterized by Ra_cold) are in the “cell” regime while the hot upwellings (characterized by Ra_hot) are in the “plume” regime. This is consistent with the co-existence of the two scales of convection just described. When both TBLS are in the “plumes” parameter range (black stars on Fig. 2), no cells are observed, but cold and hot plumes have different characteristic time-scales. This phenomenon was already documented by Schaeffer and Manga (2001) (orange stars on Fig. 2). In this regime, they showed that both time-scales were found in the hot TBL, and that cold downwellings could trigger hot upwellings nearby their impact in the hot TBL.

For $\gamma > 1000$, close up on the impact of the hot plumes under the cold TBL reveals a third scale of motions, much more localized

(Fig. 6): it appears that an impinging hot plume displaces the lower part of the cold TBL, generating a ring of cold material around the impact. The velocity field shows that the material in this cold ring is sinking but the recirculation stays very localized and does not reach the bottom of the tank.

## 2.3. Characteristics of multi-scale convection

Fig. 7 presents the vertical temperature structure recorded on the thermocouple probe as a function of time for two different viscosity ratios, and Fig. 8 the corresponding time-averaged thermal profiles. For low $\gamma$ (top Fig. 7), temperature fluctuations, with typical recurrence time of several minutes, are recorded in both TBLS. Direct visualization shows that each of them corresponds to the release of a thermal plume. Moreover the maximum amplitude of the associated temperature fluctuations are recorded at the top and bottom of the hot and cold TBLS, respectively. For intermediate $\gamma$, the temperature structure is asymmetric, with the top TBL thicker than the bottom one (Fig. 8; Richter et al., 1983; White, 1988; Manga et al., 2001). The time evolution shows two different characteristic time periods (bottom Fig. 7). Temperature fluctuations, similar to the ones just described are recorded on top of a longer time period modulation: the thermocouple probe records hot and cold plumes but also passes from preferentially hot areas to preferentially cold areas. These areas are correlated over the total height of the tank (Fig. 8 right), which means that the probe passes from one large-scale cell to the next. This meandering of the cells is responsible for the large amplitude of the temperature
<figure><img src="image_64.png" /><figcaption>**Fig. 3.** Experiment at low Rayleigh number and intermediate viscosity ratio (j1m20: γ = 888) (a) time evolution of the thermal structure in the cold TBL: the light intensity along the line of pixels indicated by the upper red line on (b) is followed through time; (b) snapshot at time t = 180 s of a vertical cross-section in the tank taken 10 cm from the front lateral wall: only one hot upwelling is observed within each cold cell, (c) light intensity along the line of pixels indicated by the lower red line on (b) is followed through time. The downwelling appears more stable than the upwellings, which meander or even disappear. The clear straight line labelled 'T' at position x = 115 mm is the thermocouple probe. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>**Fig. 4.** Experiment at high Rayleigh number and intermediate viscosity ratio (IRa12: Ra<sub>hot</sub> = 5.3 × 10<sup>7</sup>; γ = 142): Top: snapshot of a vertical cross-section taken 8 cm from the front lateral wall. The bright lines are isotherms, except the clear straight line labelled 'T' at position x = 110 mm which is the thermocouple probe. Bottom: Corresponding velocity field. Several hot plumes are seen while only one cold downwelling is observed. This situation resembles what is observed on Earth for example in the Pacific "box".</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>**Fig. 5.** (a) Time evolution of the thermal structure in the cold TBL of Fig. 4; (b) time evolution of the thermal structure in the hot TBL. The hot plume closest from the cold downwelling shows a strong continuous signature but its position fluctuates through time.</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>**Fig. 6.** Close up on the thermal (bright white lines) and velocity structures (coloured arrows) after a plume impact under the cold boundary (syr3m: γ = 3984; Ra<sub>hot</sub> = 4.7 × 10<sup>7</sup>). The two bottom isotherms of the cold TBL ("B" and "C") appear distorted and the hot plume puddle (delimited here by isotherm "A") is surrounded by a ring of colder material. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

anomalies (Fig. 8, right). But it is a chaotic process and we were not able to find any scaling for their duration (bottom Fig. 7).

On the other hand, both the temperature difference across the hot TBL and the amplitude of the hot temperature fluctuations depends on a temperature scale ΔT<sub>r</sub>, characteristic of strongly temperature-dependent viscosity fluids (Morris and Canright, 1984; Davaille and Jaupart, 1993; Solomatov, 1995; Manga et al., 2001):

$$ \Delta T_r = \frac{\eta}{\frac{d\eta}{dT}}(T). \qquad (2) $$

Fig. 9 shows that the temperature difference across the hot TBL scales as ΔT<sub>hot</sub> = 1.2(±0.2)ΔT<sub>r</sub>, as already reported by Solomatov
<figure><img src="image_68.png" /><figcaption>**Fig. 7.** Temperature as a function of dimensionless depth and time recorded on the vertical thermocouple probe. Time has been dimensionized by the convective time-scale $t_c = (H^2/\kappa)$. ($Ra_c/Ra_{hot}$)²¹/³. The temperature is colour-scaled in °C. Top: small viscosity ratio (LRa6: $Ra_{hot}$ = $1.8 \times 10^7$; $\gamma$ = 61). Bottom: intermediate viscosity ratio (LRa11: $Ra_{hot}$ = $3.1 \times 10^7$; $\gamma$ = 4078).</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>**Fig. 8.** Thermal structure recorded by the thermocouple probe as a function of the dimensionless depth for the experiment corresponding to the bottom Fig. 7 (LRa11: $Ra_{hot}$ = $3.1 \times 10^7$; $\gamma$ = 4078). Left: time-averaged temperature at each depth; right: maximum temperature difference recorded at each depth. The colours indicate the time range used to calculate the temperature structure: green for the whole duration, red for a hot area (yellow bands on previous figure), blue for a cold area. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

(1995) and Manga et al. (2001), while the temperature difference across the cold TBL increases with increasing viscosity ratio. The amplitudes of the temperature fluctuations follow the same trend. This is different from the regime of convection below a stagnant lid, where the temperature fluctuations amplitude recorded in the cold TBL scales also as $\Delta T_r$ (Fig. 9).

Meanwhile, the wavelength of plumes and cells depends on the local properties (and Ra) of each TBL. Fig. 10 shows that for $Ra > 10^5$, the minimum spacing $\lambda$ between two instabilities decreases with increasing Ra as $\lambda/H \sim Ra^{1/3}$. This means that the spacing is proportional to the TBL thickness (i.e. Howard, 1964). For $Ra < 10^5$, instabilities spacing is proportional to the thickness of the fluid layer, as already reported in the literature for rigid top and bottom boundaries.

### 3. At least three scales of convection in the Earth's mantle

To apply our experiments to the Earth, we therefore suppose that the parameters' ranges defining the different regimes seen with rigid boundaries are about the same with free boundaries. This assumption is seen to be true at least for low Rayleigh numbers

<figure><img src="image_70.png" /><figcaption>**Fig. 9.** Maximum temperature difference recorded in the cold downwellings (in blue) and in the hot upwellings (in red) as a function of the viscosity ratio. The diamonds are our new data, and the squares the data obtained by Davaille and Jaupart (1993) for quasi-steady convection under a stagnant lid. The green disks show the temperature differences measured across the hot TBL. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>**Fig. 10.** Minimum spacing $\lambda$ between rising hot currents (in red) and between cold currents (in blue) as a function of, respectively, the local Rayleigh numbers across the hot TBL $Ra_h = \frac{\alpha \rho g (T_{hot} - T_{min}) H^3}{K_y \nu_{cold}}$ and across the cold TBL $Ra_c = \frac{\alpha \rho g (T_{min} - T_{cold}) H^3}{K_y \nu_{cold}}$. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

and for the lid regime (Section 2.2). So, given the existence on Earth of a mobile lid (lithospheric plates) and of hot plumes creating hot-spots, a situation close to the laboratory experiment seen in Fig. 4, we can determine where the Earth's mantle stands in our regime diagram (Fig. 2):

#### 3.1. Subduction and large-scale circulation

The existence of subduction delimiting large cells (Fig. 1) points towards a cellular “cold” scale of convection; moreover, plate tectonics reconstructions shows that this cellular geometry evolves through geological times (e.g. Collins, 2003; Muller et al., 2008) and is therefore time-dependent. Hence the effective Rayleigh number $Ra_{cold}$ should range between $10^3$ and $10^5$ (Fig. 2). Since the temperature difference across the mantle top TBL (i.e. the lithosphere) is 1300 °C, taking $\alpha = 3 \times 10^{-4} \text{ K}^{-1}$, $H = 2900 \text{ km}$ and $\kappa = 10^{-6} \text{ m}^2 \text{ s}^{-1}$, will lead to a top lithosphere viscosity $\eta_{cold}$ ranging between $3 \times 10^{25}$ and $3 \times 10^{23}$ Pa s. This gives an estimate of the maximum effective viscosity of the lithosphere at the length-
and time-scales of mobile lid convection. Taking an average mantle viscosity of $10^{22}$ Pa s (e.g. Hager et al., 1985), we see that the lithosphere then acts as it was at most not more than ~30 to 3 000 times more viscous than the bulk mantle, or 300 to 30 000 times more viscous than an upper mantle of viscosity of $10^{21}$ Pa s.

These values are consistent with estimates from recent works on subduction using a wide range of different approaches. Matching both bathymetry and geoid in the vicinity of coherent subducting slabs and taking into account the asymmetrical nature of subduction, Billen et al. (2003) found a slab effective viscosity around $10^{23}$ Pa s. On the other hand, the comparison of regional models with subduction laboratory and numerical experiments reveal that the range of observed upper mantle slab shapes, the existence of both retreating slabs and advancing slabs, the trench curvature, and the back-arc deformations, require slabs to be only 150–500 times more viscous than the surrounding mantle (Schellart et al., 2007; Funiciello et al., 2008; Schellart, 2008; Ribe, 2010; Stegman et al., 2010). Moreover, since bending of slabs resists tectonic plate motions, the occurrence of the latter constrains the bending slab viscosity to be no more than 300 times more viscous than the upper mantle (Wu et al., 2008).

Our experiments are done with rigid top and bottom boundaries and the aspect ratio of the cold large-scale circulation ranges between 1 and 2 (Fig. 3 and 4). In the Earth's case of a stress-free upper boundary and mobile-lid regime, Weinstein and Christensen (1991) and Giannandrea and Christensen (1993) showed that the wavelength could increase up to 3.5 times the fluid layer depth. For a mantle thickness of 2900 km, it could therefore produce 10,000 km-long cells like the Pacific plate. But it would imply the presence of at least 6 cells to pave the surface of the Earth, which is different from the two mantle “boxes” (or degree 2) that we see nowadays. This discrepancy could be due to the fact that the laboratory experiments were not run in a spherical geometry. However, numerical studies run in the latter, although at lower Ra, show also the existence of more than two cells (Schubert et al., 2001; Ratcliff et al., 1997; Stemmer et al., 2006; Zhong et al., 2008). So the discrepancy may be more likely a consequence of the presence of continents on Earth.

## 3.2. Hot mantle plumes

One of the main results of our fluid mechanics study is that if the Rayleigh number associated with the lower (hot temperature) viscosity $Ra_{hot}$ is typically greater than $10^7$ (Fig. 2), several hot plumes can develop within each large-scale cold cell (Fig. 4), and the distance separating the hot plumes is proportional to the hot TBL thickness (Fig. 10). When we take the lower mantle parameter values ($\alpha = 2 \times 10^{-4}$ K$^{-1}$, $H = 2900$ km; $\kappa = 10^{-6}$ m$^2$ s$^{-1}$) and $\eta = 10^{22}$ Pa s), we obtain $Ra_{hot} \sim 4 \times 10^7$, which is indeed in the domain of multi-scale convection in Fig. 2. The situation seen on Earth with hot spots mainly confined in “boxes” delimited by subduction of cold plates is therefore well explained by thermal convection in a mantle with strongly temperature-dependent viscosity.

Fig. 10 shows that the hot plumes spacing would be ~0.2H or about 600 km. This could probably be increased up to 1200 km for a hot free slip boundary more relevant to the Earth case. But it would produce a much greater number of hot plumes than what is observed on Earth. Moreover, according to Eq. (2) and Fig. 9, the typical temperature difference across the hot TBL should be ~100 – 300 °C, which is much smaller than the temperature difference inferred across the “D” layer. Both observations point towards the presence of compositional heterogeneities at the bottom of the mantle. Denser material at the bottom of the mantle would increase the temperature difference across the TBL (Albers and Christensen, 1996; Farnetani, 1997; Davaille et al., 2002). Plumes could still be generated but also thermochemical domes and piles, with a larger spacing (Davaille et al., 2005).

According to Fig. 4 and Fig. 5, for this range of $Ra_{hot}$, the bottom zone influenced by subduction and predominantly horizontal flow is as large as the zone of predominantly upwelling flow. This situation had already been observed by Gonnermann et al. (2004) in experiments with a forced large-scale flow. We see here that it arises naturally in mobile-lid convection at high Ra. Moreover, we see also that the hot plumes generated at the border between the predominantly downwelling and upwelling regions are more long-lived than the plumes generated well within the hot area (Fig. 5). This explains well the observation whereby the long-lived hot spots (ex: Hawaii, Louisville, Easter, Tristan, Kerguelen) are usually on the edges of the two broad low-velocity zones on the CMB (e.g. Courtillot et al., 2003; Torsvik et al., 2006; Thorne et al., 2004). It is also to be noted that their long term existence does not imply their long term spatial fixity : in response to the “mantle wind” in the tank, both the plume origin in the hot TBL and its impact under the cold TBL can meander through time (Fig. 5). The typical meandering time- and length-scales are comparable to the corresponding scales in the hot TBL (Fig. 5). For the Earth, this could produce some 100–1000 km meandering on a 1–100 Ma time-scale. This could therefore explain kinks in hot spot tracks such as Hawaii (e.g. Tarduno et al., 2009).

Our experimental set up also allows us to follow the trajectories of particles (Fig. 11). We see that hot plumes bring to surface sampling only material issued from the hot TBL. Mantle material outside of the plume conduit are just recirculating within the mantle. So if hot plumes do not sample the bulk mantle, their presence nonetheless results in overall mantle stirring.

## 3.3. Heat and mass transfer in the mantle

Our results show that the large-scale pattern of convection in the Earth's mantle, with cold sheets of wide extent as preferred mode of downwellings and “clusters” of plumes as preferred mode of upwellings, results from thermal convection in a material with strongly temperature-dependent viscosity. This strengthens the case already made by previous works at lower Ra (e.g. Busse and Frick, 1985; White, 1988; Christensen and Harder, 1991) or with forced large-scale circulation (e.g. Jellinek et al., 2003; Gonnermann et al., 2004; Schubert et al., 2004).

Compared to the mantle, the laboratory experiments are done in much simpler conditions: cartesian geometry, no pressure- (i.e. depth-) dependence of material properties, no internal heating, no density heterogeneities of chemical origin. Numerical studies have shown that cold sheet and hot plume patterns are also favoured by spherical geometry (e.g. Bercovici et al., 1989; Tackley, 1996; Ratcliff et al., 1997; Schubert et al., 2001), internal heating (e.g. Sotin and Labrosse, 1999; Schubert et al., 2001; Labrosse, 2002; Moore, 2008) and the presence of a low-viscosity zone in the upper mantle (i.e. the pressure-dependence of viscosity; Cserepes and Christensen, 1990). There is also no secular heating or cooling in our tank (steady-sate conditions), and both the heat transfer and the mass transfer is symmetric between the top and the bottom of the tank. This is not the case in the mantle where the mass flux at subduction zones is much bigger than the mass flux due to hot plumes (e.g. Davies, 1988; Sleep, 1990). This asymmetry in mass fluxes is coherent with the secular cooling of the planet, whereby cooling by mobile-lid convection currently removes more efficiently the heat, than the latter is produced by internal heating (radioactive decay) and heating of the mantle by the core.

# 4. On the nature of the asthenosphere: the graveyard of hot plumes and cold small-scale instabilities

In terms of mantle convection, the lithosphere constitutes the upper cold thermal boundary layer (TBL) of the mantle, while the
<figure><img src="image_72.png" /><figcaption>**Fig. 11.** Particle trajectories obtained by taking the maximum light intensity recorded over 200 images (i.e. 2000 s) at each pixel (same experiment IRa12 as Fig. 4 and 5). We see that only the particles which have gone through the hot TBL will be sampled by the hot pancake under the surface (solid line). The other ones will recirculate (dashed line).</figcaption></figure>

asthenosphere is the region right underneath it. Studies of the mantle radial viscosity structure associate the asthenosphere with a zone with a viscosity (e.g. Hager et al., 1985; Forte and Mitrovica, 1996) lower than both the colder, more viscous lithosphere, and the mantle below 300 km. Seismology shows that the asthenosphere under the Pacific ocean also corresponds to a low seismic velocity zone (Ekstrom and Dziewonski, 1998; Ritsema and Allen, 2003; Ritzwoller et al., 2004; Maggi et al., 2006). One important question is whether the asthenosphere characteristics (low seismic velocity and low viscosity) could be the signatures of mantle flow or if they are due exclusively to the pressure-dependence of mantle mineral properties. Indeed, both a higher water content (Hirth, 2003), and the presence of partial melting (Kohlstedt et al., 1995) reduce seismic velocity and mantle viscosity. In terms of mantle convection, the asthenosphere is the region right underneath the cold TBL. So is there any fluid mechanics reason why this region should be different from the rest of the convective mantle?

As shown in our experiments, the region under the cold TBL is where all hot plumes finish their lives as well as generate small-scale cold localized instabilities (Figs. 12 and 13). This laboratory asthenosphere therefore presents an horizontally-averaged temperature hotter and an horizontally-averaged viscosity lower than both our experimental mantle at mid-depth and the cold TBL. This is also seen in numerical experiments (e.g. Jarvis and Peltier, 1982). Hence, this situation is reminiscent of the mantle asthenosphere, which would act in fact as the “graveyard” of all the hot plumes developing in the system. However, the experiments further show that such an asthenosphere should be patchy, since constituted of cold small-scale instabilities as well as hotter plumes pancakes. Recent tomographic models show that this is indeed the case, with most of the hot spots showing up as slow pancakes in the asthenosphere (Ritsema et al., 1999; Silveira and Stutzmann, 2002; Montelli et al., 2004; Debayle et al., 2005; Maggi et al., 2006).

The temperature anomalies associated with the cold small-scale instabilities should scale also as $\Delta T_r$ (Morris and Canright, 1984; Davaille and Jaupart, 1993; Solomatov and Moresi, 2000). This is well verified in our experiments (Androvandi, 2009). For mantle creep rheologies and a temperature difference of 1300 °C across the whole lithosphere, this would give small scale instabilities 50–200 °C colder than the ambient mantle, and issuing from the bottom 10–30 km of the lithosphere (Davaille and Jaupart, 1994; Singh et al., 2011). We propose that such cold rings could be responsible for the anomalously thin crust and/or depleted material found nearby hotspots tracks such as Kerguelen (Singh et al., 2011).

## 5. Conclusions

We have studied thermal convection at high Rayleigh number in a strongly temperature-dependent viscosity fluid in the

<figure><img src="image_73.png" /><figcaption>**Fig. 12.** Left: Hot instabilities in sugar syrup cooled from above and heated from below (syr1m: $R_{dhot} = 4.4 \times 10^7$; $\gamma = 3106$). Temperature field on a vertical plane section taken 8 cm from the front lateral wall. Beneath the cold TBL there are several hot pockets of fluid separated from one another by cold instabilities (seen in Fig. 13). These pockets come from older hot instabilities and are now disconnected from the hot TBL. We can also see the birth of new hot plumes. Right: Area with temperature above isotherm “A” as a function of depth. A hot temperature zone appears right below the cold TBL. Because of the temperature-dependence of viscosity, it also corresponds to a low viscosity zone.</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>**Fig. 13.** Oblique view of an horizontal plane just beneath the cold TBL of Fig. 12 (i.e. our laboratory “asthenosphere”). Plumes impacts appear as circular pancakes. They are separated from one another by cold downwelling rings.</figcaption></figure>

sluggish-lid regime (intermediate viscosity ratios), using new laboratory techniques which gave us access to both the temperature and velocity structure of the flow. We have shown that:

- in this parameter range, several scales of convection co-exist. In particular, cold downwellings delimiting cells can encase several hot upwellings, a situation reminiscent from what is observed on Earth where the Pacific subduction belt surrounds several deep seated hotspots and superswells.

- for this geometry to occur on Earth requires the lithosphere to be at most 30–3000 times more viscous than the bulk mantle, in agreement with estimates from recent works using a range of completely different approaches. It further means that there is a mechanism to weaken the cold top of the lithosphere in order to sufficiently reduce its effective viscosity.

- However, this analysis predicts typical cold cells extent around 10,000 km, which is compatible with the Pacific but would imply the presence of at least 6 cells to pave the surface of the Earth, which is different from the two boxes that we see nowadays. Convection in a strongly temperature-dependent viscosity fluid is therefore not sufficient to explain the degree 2 observed in mantle convection (Dziewonski et al., 2010). The latter is more likely a consequence of the presence of contaminants on Earth (Guillou and Jaupart, 1995; Li and Zhong, 2009). However, it would be interesting to see if the system could travel between several configurations (multiple cells, two boxes, one box).

- Hot spots plumes are stirring the bulk mantle and sampling the thermal boundary layer from which they originate. However the hot upwellings produced in those experiments present only the classical mushroom-shape morphology. Density heterogeneities of compositional origin seem required to explain superswells as well as hot spots.
- the asthenosphere is the graveyard of hot instabilities. Upon impact under the lithosphere, a hot plume will develop a hot puddle, but also a cold sheet of material around it.

The next step is now to run experiments with more complex rheology to achieve true plate tectonics and determine wether the simplified picture given in this paper still holds.

## Acknowledgements

The laboratory experiments were completed thanks to Christian Borget and Rafael Pidoux, who help to move and re-build them very quickly. This work has been supported by the DyETI program of CNRS-INSU and the French ANR projects BeGDY and PTECTO. We are also grateful to Nathanael Schaeffer and Michael Manga to have sent us their data points. The manuscript was improved, thanks to W.P. Schellart and an anonymous reviewer.

## Appendix A. Physical properties of the fluids

The physical properties of the fluid were measured in the laboratory, except for the specific heat which value is taken from the manufacturer (Table 3). The fluid rheology was studied between -20 and 65 °C using a rheometer RS600 from ThermoRheo (fluid 1 and 2) and a Physica MCR 501 from Anton Paar (fluid 3), using cylindrical geometries. The three fluids were found to be newtonian on the whole range of temperature and shear rates. Their viscosity can be modeled by:

$$ \eta = \exp(A_0 + A_1 \cdot T + A_2 \cdot T^2) \quad (A.1) $$

where the coefficients $A_0$, $A_1$, $A_2$ are given in Table 2 and T is the temperature in °C. The fluids densities were measured as a function of temperature using a Anton Paar DMA 5000. The thermal expansion coefficients and the densities at 20 °C are given in Table 3, as well as the thermal diffusivities.

<table><caption>**Table 2**
Coefficients of the viscosity law.</caption><thead><tr><th></th><th>A<sub>0</sub></th><th>A<sub>1</sub></th><th>A<sub>2</sub></th></tr></thead><tbody><tr><td>Sugar 1</td><td>4.7709</td><td>-0.1432</td><td>0.0007</td></tr><tr><td>Sugar 2</td><td>1.971</td><td>-0.1107</td><td>0.0006</td></tr><tr><td>Sugar 3</td><td>3.5263</td><td>-0.1342</td><td>0.0007</td></tr></tbody></table>

<table><caption>**Table 3**
Properties of the fluids at 20 °C.</caption><thead><tr><th></th><th>α (K<sup>-1</sup>)</th><th>ρ (kg m<sup>-3</sup>)</th><th>κ (m<sup>2</sup> s<sup>-1</sup>)</th><th>η (20:60 °C)(Pa s)</th></tr></thead><tbody><tr><td>Sugar 1</td><td>5.73 × 10<sup>-4</sup></td><td>1399</td><td>1.16 × 10<sup>-7</sup></td><td>8.91:0.27</td></tr><tr><td>Sugar 2</td><td>5.70 × 10<sup>-4</sup></td><td>1384</td><td>1.28 × 10<sup>-7</sup></td><td>1.0:0.08</td></tr><tr><td>Sugar 3</td><td>5.95 × 10<sup>-4</sup></td><td>1389</td><td>1.20 × 10<sup>-7</sup></td><td>3.07:0.13</td></tr></tbody></table>

## References

Albers, M., Christensen, U.R., 1996. The excess temperature of plumes rising from the core-mantle boundary. Geophys. Res. Lett. 23, 3567–3570.

Androvandi, S., 2009. Convection multi-chelles dans un Manteau dont la viscosité dépend de la température, Ph.D. Thesis, Institut de Physique du Globe de Paris, France.

Bercovici, D., 1996. Plate generation in a simple model of lithosphere-mantle flow with dynamic self-lubrication. Earth Planet. Sci. Lett. 144, 41–51.

Bercovici, D., 1998. Generation of plate tectonics from lithosphere mantle flow and void-volatile self-lubrication. Earth Planet. Sci. Lett. 154, 139–151.

Bercovici, D., 2003. The generation of plate tectonics from mantle convection. Earth Planet. Sci. Lett. (Frontiers) 205, 107–121.

Bercovici, D., Schubert, G., Glatzmaier, G.A., Zebib, A., 1989. Three-dimensional thermal convection in a spherical shell. J. Fluid Mech. 206, 75–104.

Billen, M.I., Gurnis, M., Simons, M., 2003. Multiscale dynamics of the Tonga Kermadec subduction zone. Geophys. J. Int. 153, 359–388.

Busse, F.H., Frick, H., 1985. Square-pattern convection in fluids with strongly temperature-dependent viscosity. J. Fluid Mech. 150, 451–465.

Christensen, U., 1984a. Heat transport by variable viscosity convection and implications for the Earth's thermal evolution. Phys. Earth Planet. Int. 35, 264–282.

Christensen, U.R., 1984b. Convection with pressure and temperature-dependent non-Newtonian rheology. Geophys. J. R. Astr. Soc. 77, 343–384.

Christensen, U.R., 1985. Heat transport by variable viscosity convection II: pressure influence, non-Newtonian rheology and decaying heat sources. Phys. Earth Planet. Int. 37, 183–205.

Christensen, U.R., Harder, H., 1991. 3-D Convection with variable viscosity. Geophys. J. Int. 104, 213–220.

Collins, W.J., 2003. Slab pull, mantle convection, and Pangaean assembly and dispersal. Earth Planet. Sci. Lett. 205, 225–237.

Courtillot, V., Davaille, A., Besse, J., Stock, J., 2003. Three distinct types of hotspots in the Earth's mantle. Earth Planet. Sci. Lett. 205, 295–308.

Cserepes, L., Christensen, U.R., 1990. Three-dimensional convection under drifting plates. Geophys. Res. Lett. 17, 1497–1500.

Davaille, A., Jaupart, C., 1993. Transient high-Rayleigh number thermal convection with large viscosity variations. J. Fluid Mech. 253, 141–166.

Davaille, A., Jaupart, C., 1994. Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle. J. Geophys. Res. 99, 19853–19866.

Davaille, A., Girard, F., Le Bars, M., 2002. How to anchor plumes in a convecting mantle? Earth Planet. Sci. Lett. 203, 621–634.

Davaille, A., Stutzmann, E., Silveira, G., Besse, J., Courtillot, V., 2005. Convective patterns under the Indo-Atlantic? box? Earth Planet. Sci. Lett. 239, 233–252.

Davaille, A., Limare, A., 2007. Laboratory Studies of Mantle Convection. In: Bercovici, D., Schubert, G. (Eds.), Treatise of Geophysics, vol. 7. Elsevier, pp. 89–165.

Davies, G.F., 1988. Ocean bathymetry and mantle convection, 1. Large-scale flow and hotspots. J. Geophys. Res. 93, 10467–10480.

Davies, J.H., 2005. Steady plumes produced by downwellings in Earth-like vigorous spherical whole mantle convection models. Geochem. Geophys. Geosyst. 6, Q12001. doi:10.1029/2005GC001042.

Davies, D.R., Davies, J.H., 2009. Thermally-driven mantle plumes reconcile multiple hot-spot observations. Earth Planet. Sci. Lett. 278, 50–54.

Debayle, E., Kennett, B., Priestley, K., 2005. Global azimuthal seismic anisotropy and the unique plate-motion deformation of Australia. Nature 433, 509–512.

Dziewonski, A.M., Lekic, V., Romanowicz, B., 2010. Mantle anchor structure: an argument for bottom up tectonics. Earth Planet. Sci. Lett. 299, 69–79.

Ekstrom, G., Dziewonski, A., 1998. The unique anisotropy of the Pacific upper mantle. Nature 394, 168–172.

Farnetani, C.G., 1997. Excess temperature of mantle plumes: the role of chemical stratification across D. Geophys. Res. Lett. 24, 1583–1586.

Forte, A.M., Mitrovica, J.X., 1996. New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post-glacial rebound data. Geophys. Res. Lett. 23, 1147–1150.

Funiciello, F., Faccenna, C., Heuret, A., Lallemand, S., Di Giuseppe, E., Becker, T.W., 2008. Trench migration, net rotation and slab mantle coupling. Earth Planet. Sci. Lett. 271, 233–240.

Giannandrea, E., Christensen, U.R., 1993. Variable viscosity convection experiments with a stress-free upper boundary and implications for the heat transport in the Earths mantle. Phys. Earth Planet. Sci. 78, 139–152.

Gonnermann, H.M., Jellinek, A.M., Richards, M.A., Manga, M., 2004. Modulation of mantle plumes and heat flow at the core mantle boundary by plate-scale flow: results from laboratory experiments. Earth Planet. Sci. Lett. 226, 53–67.

Grand, S.P., 2002. Mantle shear-wave tomography and the fate of subducted slabs. Phil. Trans. R. Soc. London Ser. A 360, 2475–2491.

Grand, S.P., van der Hilst, R.D., Widiyantoro, S., 1997. Global seismic tomography: A snapshot of convection in the earth. GSA Today 4, 1–7.

Guillou, L., Jaupart, C., 1995. On the effect of continents on mantle convection. J. Geophys. Res. 100, 24217–24238.

Gurnis, M., Zhong, S., Toth, J., 2000. On the competing roles of fault reactivation and brittle failure in generating plate tectonics from mantle convection. In: Richards, M.A., Gordon, R., Van der Hilst, R. (Eds.), The History and Dynamics of Global Plate Motions, AGU Geophysical Monograph 21, pp. 73–94.

Hager, B.H., Clayton, R.W., Richards, M.A., Comer, R.P., Dziewonski, A.M., 1985. Lower mantle heterogeneity, dynamic topography and the geoid. Nature 313, 541–545.

Hansen, U., Yuen, D.A., 1993. High Rayleigh number regime of temperature-dependent viscosity convection and the Earth's early thermal history. Geophys. Res. Lett. 20, 2191–2194.

Hirth, G., 2003. Laboratory constraints on the rheology of the upper mantle. In: Karato, S.-I., Wenk, H.-R. (Eds.), Plastic Deformation of Minerals and Rocks, Rev. Mineral. Geochem., vol. 51.

Houseman, G.A., 1990. The thermal structure of mantle plumes: axisymmetric or triple-junction? Geophys. J. Int. 102, 25–43.

Jarvis, G., Peltier, W.R., 1982. Mantle convection as a boundary layer phenomenon. Geophys. J. R. Astr. Soc. 68, 389–427.

Jaupart, C., Parsons, B., 1985. Convective instabilities in a variable viscosity fluid cooled from above. Phys. Earth Planet. Int. 39, 32–41.
Jellinek, A.M., Gonnermann, H.M., Richards, M.A., 2003. Plume capture by divergent plate motions: Implications for the distribution of hotspots, geochemistry of mid-ocean ridge basalts, and heat flux at the core-mantle boundary. Earth Planet. Sci. Lett. 205, 367–378.

Kameyama, M., Ogawa, M., 2000. Transitions in thermal convection with strongly temperature-dependent viscosity in a wide box. Earth Planet. Sci. Lett. 180, 355–367.

Kohlstedt, D.L., Evans, B., Mackwell, S.J., 1995. Strength of the lithosphere: constraints imposed by laboratory experiments. J. Geophys. Res. 100, 17,587–17,602.

Krishnamurti, R., 1970. On the transition to turbulent convection. Part 2. Transition to time-dependent flow. J. Fluid Mech. 42, 309–320.

Labrosse, L., 2002. Hotspots, mantle plumes and core heat loss. Earth Planet. Sci. Lett. 199, 147–156.

Li, Z.-X., Zhong, S., 2009. Supercontinent-superplume coupling, true polar wander and plume mobility: plates dominance in whole -mantle tectonics. Phys. Earth Planet. Int. 176, 143–156.

Maggi, A., Debayle, E., Priestley, K., Barruol, G., 2006. Multi-mode surface waveform tomography of the Pacific Ocean: a closer look at lithospheric cooling. Geophys. J. Int. 166, 1384–1397.

Manga, M., Weeraratne, D., 1999. Experimental study of non-Boussinesq Rayleigh-Bénard convection at high Rayleigh and Prandtl numbers. Phys. Fluids 10, 2969–2976.

Manga, M., Weeraratne, D., Morris, S.J.S., 2001. Boundary-layer thickness and instabilities in Bénard convection of a liquid with a temperature-dependent viscosity. Phys. Fluids 13, 802–805.

Masters, G., Jordan, T.H., Silver, P.G., Gilbert, F., 1982. Aspherical Earth structure from fundamental spheroidal-mode data. Nature 298, 609–613.

McNamara, A.K., Zhong, S., 2005. Thermochemical structures beneath Africa and the Pacific Ocean. Nature 437, 1136–1139.

Molnar, P., Stock, J., 1987. Relative motions of hotspots in the Pacific, Atlantic and Indian Oceans since late Cretaceous time. Nature 327, 587–591.

Montelli, R., Nolet, G., Dahlen, F.A., Masters, G., Engdahl, E.R., Hung, S.-H., 2004. Finite-frequency tomography reveals a variety of plumes in the mantle. Science 303, 338–343.

Moore, W.B., 2008. Heat transport in a convecting layer heated from within and below. J. Geophys. Res. 113, B11407. doi:10.1029/2006JB004778.

Moresi, L.N., Solomatov, V.S., 1995. Numerical investigation of 2D convection with extremely large viscosity variations. Phys. Fluids 7, 2154–2162.

Morris, S., Canright, D.R., 1984. A boundary-layer analysis of Benard convection with strongly temperature-dependent viscosity. Phys. Earth Planet. Inter. 36, 355–373.

Muller, R.D., Sdrolias, M., Gaina, C., Roest, W.R., 2008. Age, spreading rates, and spreading asymmetry of the worlds ocean crust. Geochem. Geophys. Geosyst. 9, Q04006. doi:10.1029/2007GC001743.

Nataf, H.C., Richter, F.M., 1982. Convection experiments in fluids with highly temperature-dependent viscosity and the thermal evolution of the planets. Phys. Earth Planet. Int. 29, 320–329.

Ogawa, M., Schubert, G., Zebib, A., 1991. Numerical simulations of three-dimensional thermal convection in a fluid with strongly temperature-dependent viscosity. J. Fluid Mech. 233, 299–328.

Olson, P., Schubert, G., Anderson, C., Goldman, P., 1988. Plume formation and lithosphere erosion: a comparison of laboratory and numerical experiments. J. Geophys. Res. 93, 15065–15084.

Parmentier, E.M., Sotin, C., 2000. Three-dimensional numerical experiments on thermal convection in a very viscous fluid: Implications for the dynamics of a thermal boundary layer at high Rayleigh number. Phys. Fluids 12, 609–617.

Ratcliff, J.T., Tackley, P.J., Schubert, G., Zebib, A., 1997. Transitions in thermal convection with strongly variable viscosity. Phys. Earth Planet. Inter. 102, 201–212.

Ribe, N.M., 2010. Bending mechanics and mode selection in free subduction: a thin? sheet analysis. Geophys. J. Int. 180, 559–576.

Ricard, Y., Bercovici, D., Schubert, G., 2001. A two-phase model for compaction and damage, 2. Applications to compaction, deformation and the role of interfacial surface tension. J. Geophys. Res. 106, 8907–8924.

Richter, F.M., Nataf, H.-C., Daly, S.F., 1983. Heat transfer and horizontally averaged temperature of convection with large viscosity variations. J. Fluid Mech. 129, 173–192.

Ritsema, J., van Heijst, H., Woodhouse, J., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland. Science 286, 1925–1928.

Ritsema, J., Allen, R.M., 2003. The elusive mantle plume. Earth Planet. Sci. Lett. 207, 1–12.

Ritzwoller, M., Shapiro, N., Zhong, S.-J., 2004. Cooling history of the Pacific lithosphere. Earth planet. Sci. Lett. 226, 69–84.

Schaeffer, N., Manga, M., 2001. Interaction of rising and sinking mantle plumes. Geophys. Res. Lett. 21, 765–768.

Schellart, W.P., Freeman, J., Stegman, D.R., Moresi, L., May, D., 2007. Evolution and diversity of subduction zones controlled by slab width. Nature 446, 308–311.

Schellart, W.P., 2008. Kinematics and flow patterns in deep mantle and upper mantle subduction models: influence of the mantle depth and slab to mantle viscosity ratio. Geochem. Geophys. Geosyst. 9, Q03014. doi:10.1029/2007GC001656.

Schmalzl, J., Breuer, M., Hansen, U., 2002. The influence of the Prandtl number on the style of vigorous thermal convection. Geophys. Astrophys. Fluid Dynamics 96 (5), 381–403.

Schubert, G., Turcotte, D.L., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press, Cambridge, England.

Schubert, G., Masters, G., Olson, P., Tackley, P., 2004. Superplumes or plume clusters? Phys. Earth Planet. Int. 146, 147–162.

Silveira, G., Stutzmann, E., 2002. Anisotropic tomography of the Atlantic Ocean. Phys. Earth Planet. Inter. 132, 237–248.

Singh, S., Carton, H., Chauhan, A.S., Androvandi, S., Davaille, A., Dyment, J., Cannat, M., Hananto, N.D., 2011. Extremely thin crust in the Indian Ocean possibly resulting from Plume-Ridge Interaction. Geophys. J. Int. 184, 29–42.

Sleep, N.H., 1990. Hotspots and mantle plumes: some phenomenology. J. Geophys. Res. 95, 6715–6736.

Solomatov, S., 1995. Scaling of temperature- and stress-dependent viscosity convection. Phys. Fluids 7, 266–274.

Solomatov, V.S., Moresi, L.N., 2000. Scaling of time-dependent stagnant lid convection: application to small-scale convection on Earth and other terrestrial planets. J. Geophys. Res. 105, 21785–21817.

Sotin, C., Labrosse, S., 1999. Three-dimensional Thermal convection of an isoviscous, infinite-Prandtl-number fluid heated from within and from below: Applications to heat transfer in planetary mantles. Phys. Earth Planet. Int. 112, 171–190.

Stein, C., Schmalzl, J., Hansen, U., 2004. The effect of rheological parameters on plate behaviour in a self-consistent model of mantle convection. Phys. Earth Planet. Int. 142, 225–255.

Stegman, D.R., Farrington, R., Capitanio, F.A., Schellart, W.P., 2010. A regime diagram for subduction styles from 3-D numerical models of free subduction. Tectonophysics 483, 29–45.

Stemmer, K., Harder, H., Hansen, U., 2006. A new method to simulate convection with strongly temperature-dependent and pressure-dependent viscosity in a spherical shell: applications to the Earth's mantle. Phys. Earth Planet Inter. 157, 223–249.

Stengel, K.C., Oliver, D.S., Booker, J.R., 1982. Onset of convection in a variable-viscosity fluid. J. Fluid Mech. 120, 411–431.

Tackley, P.J., Stevenson, D.J., Glatzmaier, G.A., Schubert, G., 1994. Effects of multiple phase transitions in a 3-dimensional spherical model of convection in Earth's mantle. J. Geophys. Res. 99, 15887–15901.

Tackley, P.J., 1996a. Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles. J. Geophys. Res. 101, 3311–3332.

Tackley, P.J., 1996b. On the ability of phase transitions and viscosity layering to induce long-wavelength heterogeneity in the mantle. Geophys. Res. Lett. 23, 1985–1988.

Tackley, P.J., 2000a. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 1. Pseudoplastic yielding. Geochem. Geophys. Geosyst. 1, 2000GC000043.

Tackley, P.J., 2000b. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 2. Strain weakening and asthenosphere. Geochem. Geophys. Geosyst. 1, 2000GC000044.

Tackley, P.J., 2000c. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations. Part 2. Strain weakening and asthenosphere. Geochem. Geophys. Geosyst. 1, 2000GC000043.

Tarduno, J., Bunge, H.P., Sleep, N., Hansen, U., 2009. The bent hawaiian-emperor hotspot track: inheriting the mantle wind. Science 324, 50–53.

Thorne, M.S., Garnero, E.J., Grand, S.P., 2004. Geographic correlation between hot spots and deep mantle lateral shear-wave velocity gradients. Phys. Earth Planet Int. 146, 47–63.

Torsvik, T.H., Smethurst, M.A., Burke, K., Steinberger, B., 2006. Large igneous provinces generated from the margins of the large low-velocity provinces in the deep mantle. Geophys. J. Int. 167, 1447–1460.

Travis, B., Olson, P., Schubert, G., 1990. The transition from two-dimensional to three-dimensional planforms in infinite-Prandtl-number thermal convection. J. Fluid Mech. 216, 71–91.

Trompert, R.A., Hansen, U., 1998. On the Rayleigh number dependence of convection with a strongly Temperature-dependent viscosity. Phys. Fluids 10, 351–360.

Weeraratne, D., Manga, M., 1998. Transitions in the style of mantle convection at high Rayleigh numbers. Earth Planet. Sci. Lett. 160, 563–568.

Weinstein, S.A., Olson, P.L., 1989. The proximity of hotspots to convergent and divergent plate boundaries. Geophys. Res. Lett. 16, 433–436.

Weinstein, S.A., Christensen, U.R., 1991. Convection planforms in a fluid with a temperature-dependent viscosity beneath a stress-free upper boundary. Geophys. Res. Lett. 18, 2035–2038.

White, D.B., 1988. The planforms and onset of convection with a temperature-dependent viscosity fluid. J. Fluid Mech. 191, 247–286.

Wu, B., Conrad, C.P., Heuret, A., Lithgow-Bertelloni, A., Lallemand, S., 2008. Reconciling strong slab pull and weak plate bending: the plate motion constraint on the strength of mantle slabs. Earth Planet. Sci. Lett. 272, 412–421.

Yoshida, M., Kageyama, A., 2006. Low-degree mantle convection with strongly temperature- and depth-dependent viscosity in a three-dimensional spherical shell. J. Geophys. Res. 111, B03412. doi:10.1029/2005JB003905.

Zhong, S., McNamara, A., Tan, E., Moresi, L., Gurnis, M., 2008. A benchmark study on mantle convection in a 3-D spherical shell unsing CitcomS. Geochem. Geophys. Geosyst. 9, 2008GC002048.
<figure><img src="image_75.png" /></figure>

# Optical method for measuring bed topography and flow depth in an experimental flume

A. Limare¹, M. Tal¹,*, M. D. Reitz², E. Lajeunesse¹, and F. Métivier¹

¹Equipe de Dynamique des Fluides Géologiques, Institut de Physique du Globe de Paris (IPGP, UMR 7154, CNRS, Univ. Paris Diderot), 1 rue Jussieu, 75238 Paris Cedex 05, France

²Department of Physics and Astronomy, University of Pennsylvania, 209 S. 33rd St., Philadelphia, PA 19104, USA

*now at: Univ. Aix-Marseille, CEREGE UMR 6635, Europôle de l'Arbois, BP 80, 13545 Aix-en-Provence cedex 04, France

Received: 1 March 2011 – Published in Solid Earth Discuss.: 15 March 2011

Revised: 17 June 2011 – Accepted: 7 July 2011 – Published: 28 July 2011

**Abstract.** We describe an optical method known as moiré for acquiring quasi-simultaneous measurements of bed topography and flow depth in laboratory experiments. The moiré method is based on projecting a fringe pattern (grating) on the bed and analyzing the deformation of the pattern caused by the topography with respect to a reference plane. The height of the object is encoded in the phase of the pattern and can be retrieved either through Fourier transform or phase shifting algorithms. The methodology enables image-based non-contact measurements over a continuous surface at very high spatial and temporal resolutions. We use a commercial software package of a moiré method called Light3D to map bed topography and flow depth in an experimental braided channel and demonstrate how the method can be used to characterize a full range of statistics not previously possible.

## 1 Introduction

Laboratory experiments are an important tool in the study of morphodynamics and have been used to investigate a wide range of landforms ranging from submarine channels (Métivier et al., 2005) to the response of drainage basins to tectonic uplift (Lague et al., 2003; Turowski et al., 2006). Analog models offer an opportunity to examine interactions underlying naturally complex systems in a well-controlled environment, and to observe how systems evolve over much longer timescales than possible in the field. In fluvial morphodynamics, experiments have been used successfully to

reproduce and visualize sediment-flow interactions and bed morphology (Lajeunesse et al., 2010a) and have helped improve our understanding of self-formed channels with mobile beds and stable banks (Armstrong, 2003; Malverti et al., 2008), braiding kinematics (Métivier and Meunier, 2003), bedforms (e.g., Coleman and Eling, 2000), seepage channels (Ni and Crave, 2006), alluvial fans (Reitz et al., 2010), knickpoints (Malverti et al., 2007), and alluvial bed morphodynamics (Coleman and Eling, 2000; Devauchelle et al., 2010a,b). (For a more complete review, see Lajeunesse et al., 2010b; Paola, 2010).

Several experimental techniques have been developed specifically for data acquisition in experiments with the objective of monitoring both bed topography and flow characteristics. Classical approaches to the measurement of bed topography in these types of experiments include moving point or line probes along the experimental surface such as single beam time-of-flight laser profiler mounted on a controlled carriage (Lague et al., 2003), or scanning a stripe laser over the experimental surface (Rice et al., 1988; Darboux and Huang, 2003; Huang et al., 2007a). These techniques often require interrupting the experiment (i.e., shutting off the flow) at regular intervals and moving the probe, often very slowly, over the temporarily frozen surface. A faster and continuous alternative is based on monitoring the experiment with a camera set to capture time-lapse images of the evolving surfaces. This method allows reconstruction of the topography using stereoscopy (Chandler et al., 2001).

To map flow depths, depth-from-luminosity measurements have been successfully acquired over both fixed and evolving topography. Aureli et al. (2008) mapped the evolving flow depth over a flat, transparent floor, illuminated from below, by adding blue dye to the water and monitoring from above

<figure><img src="image_76.png" /></figure>

Correspondence to: A. Limare (limare@ipgp.fr)
the distribution of light attenuation. Tal and Paola (2007) used rhodamine dye to estimate the distribution of water depth in braided and single-thread channels. Both methods involve using flow with a constant dye concentration and calibrating the depth dependence of the luminosity as recorded by camera images. These methods are closely related to the common use of fluorescent dye to monitor ground-water flow patterns (Weiler and Flühler, 2004; Huang et al., 2007b).

While independent methods for measuring bed topography and flow depth exist, simultaneous measurement of bed topography and flow depth remains a challenge. Huang et al. (2010) proposed a novel combination of methods to study the morphodynamics of sediment beds under the action of flowing water. In this method, a laser stripe is scanned back and forth across the experimental surface while a fluorescent dye is mixed with the flowing water, allowing simultaneous mapping of both bed topography and water depth at millimetric resolution and sub-millimetric accuracy (0.6 mm).

In this paper we describe a method for acquiring quasi-simultaneous measurements of bed topography and flow depth at very high resolution using an optical method known as moiré. This method requires the interrupting of the flow and the acquisition of two sets of topographies with and without flow. We begin by describing the different steps in applying this version of the moiré method to reconstruct the topography of a dry bed. We then discuss the case of an experimental alluvial plain with water flowing over the bed. Finally, we apply this method to a microscale braided channel experiment and discuss the quality of the results obtained.

The moiré method has been successfully applied in metrology studies (Chiang, 1983), industrial inspection (Sansoni et al., 2000), human body mapping for medical diagnosis (Halioua and Liu, 1989; Kozłowski and Serra, 1997) and art inspection (Brémand et al., 2007).

In geomorphology, simplified versions of the moiré method have been previously used in laboratory physical experiments to study granular flows (Pouliquen and Forterre, 2002), incision dynamics of subaqueous channels (Lancien et al., 2005), and micro-scale braided rivers (Meunier, 2004). In these studies the method was implemented by projecting a single grid onto the evolving surface and analyzing how it is deformed (captured by time-lapse images) using a Fourier transform. This approach has the advantage of being low cost (minimum materials required include only a digital camera and an overhead or slide projector) and is relatively simple to implement. However, it offers limited resolution and a limited field of view, and therefore limited object dimensions.

Here we present a more elaborate implementation of the moiré method to two fluvial geomorphology studies. This method was developed at the Laboratory of Solid State Mechanics of the University of Poitiers, France (Breque et al., 2004). This commercial software was already used for the experimental investigation of bedforms dynamics in a configuration where the sediment supply was limited (Dreano et al., 2010). It uses a phase shifting method for phase calculation,

robust phase unwrapping combined with gray coding and a calibration procedure that allows assessing the geometric parameters. The method proves relatively easy to use and to calibrate. The number of projected images can be varied so one can switch between high resolution and fast acquisition time according to the requirements of the experiments.

Braiding has always been easy to reproduce in the lab and many experiments have been developed since the fifties (Ashmore, 1982; Federici and Paola, 2003; Schumm, 1977; Schumm et al., 1987). Yet the existence of multiple highly dynamic threads poses challenging issues. The technique we present here enables us to address two key issues in braided river studies: topographic evolution of the bed and shear stress distribution.

The novelty here is that we are able to get the instantaneous depth everywhere on the braid plane. Thus, for the first time, we are able to compute PDFs of flow depth and probably shear stress distributions. To our knowledge there exist only two published PDFs for natural braided streams (Paola, 1996; Métivier et al., 2011) and one for experimental braided rivers (Tal and Paola, 2007). Both of these records are rough and limited, precisely because of the acquisition techniques. Yet these studies and both theoretical and numerical considerations (Paola, 1996, 2001) clearly point to the key importance of such measurements in order to understand the physics of sediment transport in a braided river.

# 2 Presentation of the moiré method

## 2.1 General considerations

The moiré phenomenon can be observed when superimposing two (or more) periodic or quasiperiodic structures; the so-called moiré method is not in fact one method, but rather a group of methods. Depending upon how or from what type of geometric patterns the moiré fringes are created, they assume different physical signification. The method of moiré contouring is a well known technique for 3-D shape measurement (Patorski, 1993): the moiré fringes are created by the superposition of two periodic alternating line-space patterns, i.e., linear gratings. One of these gratings is to be distorted by an object whose deformation or shape is represented by the moiré fringes. The main merits of the moiré topography are:

1. Measurements are whole-field, an advantage over point-by-point or line methods.

2. Data acquisition and processing is relatively fast so investigations can be done in quasi real-time and the process automated.

3. The resolution may be varied.
4. Differential as well as absolute measurements are possible so that the influence of system geometry and optical aberrations may be removed.

The output of a moiré fringe technique is a fringe pattern which must be analysed so that the results can be presented in a numerical or graphical form. In the process of fringe pattern analysis, the only measurable quantity is the intensity that has a cosinusoidal profile:

$$I(x, y) = a(x, y) + b(x, y)\cos[\phi(x, y)] \quad (1)$$

where $I(x, y)$ is the measured intensity, $\phi(x, y)$ is the phase of the signal, $a(x, y)$ is the background intensity and $b(x, y)$ is the local contrast; both $a(x, y)$ and $b(x, y)$ being proportional to the surface reflectivity. The phase $\phi(x, y)$ contains the desired out-of-plane information and should be extracted from Eq. (1). The main problems related to the phase evaluation are:

1. The phase is screened by two other functions: the background intensity $a(x, y)$ and the local contrast $b(x, y)$.

2. Due to the periodicity of cosinusoidal intensity profile, the phase is only determined modulo $2\pi$.

3. Due to the even character of the cosine function, the sign of $\phi(x, y)$ cannot be extracted from a single measurement, so that one cannot distinguish between depression and elevation without a priori knowledge.

If the moiré method usually supposes the superposition of two or more periodic structures, in this work we used a cosine fringe pattern projection. Simple Fourier analysis of the fringe pattern in the presence of the spatial carrier frequency, for example in the x direction, gives:

$$\phi'(x, y) = 2\pi f_0 x + \phi(x, y) \quad (2)$$

The phase to be measured, $\phi(x, y)$, modulates the carrier fringe pattern and this modulation is indicated by changes in the position and spatial frequency of the fringes.

Experimentally, the method consists of projecting a grating of parallel lines onto an object, in our case a stream table (Fig. 1). The topography of the object will deform the grating (Fig. 2b) and this deformation is calculated with respect to a reference plane which in our case is a flat bed (Fig. 2a). The height of the object is encoded in the phase of the pattern and can be retrieved either through the Fourier transform method or phase shifting algorithms.

## 2.2 Fourier-transform profilometry

If $a(x, y)$, $b(x, y)$ and $\phi(x, y)$ are slowly varying functions compared to the variation introduced by the carrier frequency $f_0$, simple Fourier spectrum analysis can be used to obtain

<figure><img src="image_77.png" /><figcaption>Fig. 1. Experimental setup.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Fig. 2. Fringe pattern on (a) flat surface, (b) relief.</figcaption></figure>

the phase $\phi(x, y)$ (Takeda et al., 1982). The cosine function in Eq. (1) can be rewritten as a sum of complex exponentials:

$$I(x, y) = a(x, y) + c(x, y) \exp (j2\pi f_0 x) \quad (3) \\ + c^*(x, y) \exp (-j2\pi f_0 x)$$

with

$$c(x, y) = \frac{1}{2} b(x, y) \exp [j\phi(x, y)] \quad (4)$$

After one-dimensional Fourier transform in the x direction, the following holds:

$$F(I) = A(f, y) + C(f - f_0, y) + C^*(f + f_0, y) \quad (5)$$

where $A$, $C$ and $C^*$ are complex Fourier amplitudes and $f$ is the spatial frequency in the x direction. The spectrum term $C(f - f_0, y)$ can be isolated using a filter centered at $f_0$ and the carrier frequency can be removed by shifting $C(f - f_0, y)$ on the frequency axis toward the origin by $f_0$. The unwanted background variation $a(x, y)$ has been filtered out at this stage. Applying the inverse Fourier transform to
the C(f, y), c(x, y) is obtained and the phase function may be calculated:

$$ \phi(x, y) = \arctan \frac{\mathrm{Re}[c(x, y)]}{\mathrm{Im}[c(x, y)]} \qquad (6) $$

The main advantage of the FFT method is that phase values are extracted only from a single image. The method is mainly suitable for real-time applications where speed is more important than accuracy. The various sources of error such as the influence of noise, fringe discontinuities, and choice of filtering window are discussed in (Takeda, 1990; Kozłowski and Serra, 1999) and a few tracks are given to minimize them.

## 2.3 Phase-shifting method

The principle of the phase-shifting method consists of generating a series of images where the projected fringe pattern or grating is shifted in its plane (Schmit and Creath, 1995). A minimum of three measurements is necessary to calculate the phase since there are three unknowns in the fringe pattern (1). When the phase shifts are chosen such that N measurements are equally spaced over one period: i.e., $\delta i = i2\pi/N$ where i = 1N, the phase is given by the equation:

$$ \phi(x, y) = \arctan \frac{\sum I_i(x, y) \sin \delta_i}{\sum I_i(x, y) \cos \delta_i} \qquad (7) $$

Because the phase is continuously distributed within its range of non-ambiguity (0–2π), the height resolution is limited only by the light level quantization and noise (an important source of noise can come from imperfect phase shifting if the grid pitch is not constant over the surface). On the other hand, due to the nature of arctan calculations, the equations for phase calculations (6) and (7) are determined only modulo 2π. This fact strongly reduces the height range, and an unwrapping process is necessary to reconstruct the surface map.

## 2.4 Phase unwrapping

Unwrapping uses an integration technique that sums up the phases to remove jumps between adjacent pixels larger than $\pi$. When the phase difference between adjacent pixels is greater than $\pi$, a multiple of 2π is added or subtracted to make the difference less than $\pi$. Unwrapping becomes more difficult when the absolute phase difference between adjacent pixels is greater than $\pi$. These discontinuities may be introduced by high-frequency, high-amplitude noise and discontinuous phase jumps. Several algorithms were proposed with different degrees of robustness for phase unwrapping as reviewed in the book of Ghiglia and Pritt (1998). Two types of approaches of 2-D phase unwrapping exists: one is based on path-following methods or local methods, while the other is based on minimum-norm methods or global methods.

For objects presenting marked discontinuities of the surface, a technique using gray-code structured light can be

employed to remove phase unambiguities at the expense of reduced resolution (Sansoni et al., 1997). It has been shown that this technique, also called temporal phase unwrapping reduces the computation time and improves the unwrapping reliability (Huntley and Saldner, 1997).

In the case of multiple-shot combination of gray-code and phase shifting, both extended measuring range and high resolution can be attained, as in the commercialized methods described by Sansoni et al. (2003) and Breque et al. (2004). These methods require the projection of 8 and respectively 16 computer-generated fringe patterns, which is time consuming. The most important drawback of the multiple-shot methods such as gray-code methods, phase shifting methods, and their combination, is their requirement for a series of images to be recorded at different moments of time. This feature ultimately imposes a limit on the speed of data capture and analysis.

For applications requiring real-time or dynamic measurements, a different projection scheme was proposed: unwrapping using local wave algorithm (Liebling et al., 2004) from a single image (Cochard and Ancey, 2008).

But phase unwrapping can be avoided, as shown by Kozłowski and Serra (1997). Their technique, called phase locked loop, uses two images with a projected grid shifted in phase by half of its period.

In general, the coding method is a trade-off between the resolution of the measurement, the characteristics of the target shape, the whole dimensions of the surface, and the time available for the projection.

## 2.5 Phase-to-height conversion

Once the phase map has been calculated, the conversion from phase to actual object height needs to be performed. In the case of parallel projected and observed beams, the relief Z is proportional to the phase difference between the reference grating and the grating deformed by the object (Brémand, 1994):

$$ Z = \frac{p}{2\pi} \frac{\phi(x, y)}{\tan \alpha} \qquad (8) $$

where p is the projected grating spacing. The simplified assumption of parallel beams is valid only if:

1. The projection and observation distances are the same ($h_p = h_o = h$);

2. The relief is small compared to the observation distance ($z \ll h$);

3. The size of object is small compared to the distance between the light source and the camera ($x \ll d$).

Moiré methods suppose that the projection grid on the reference plane has a constant frequency and that the light projected onto the object is parallel. These requirements hold true for small objects and for a light projector located far
from it. In practice, in order to obtain a larger field of view, diverging illumination needs to be used. Also, for large objects the optics of the camera could introduce significant distortion, but this can be corrected through image processing. The limiting requirements of constant frequency grid and parallel beams condition can be overcome by the use of a new method (Breque et al., 2004). The method is based on a pinhole model (Cardenas-Garcia et al., 1995), generally used in the principles of stereovision. Because of the diverging nature of the illumination, this situation requires different, more complicated analysis than Eq. (8). In this general case (Eq. 10), the relief is calculated with the phase $\phi(x, y)$, and the geometric parameters ($h_p, d, h_o, f_p$ and $P_r$); $h_p, d, h_o$ are defined in Fig. 1; $P_r/f_p$ is the projector magnification, $f_p$ and $P_r$ are the focal length of the projector and the grating pitch at the output lens of the projector respectively:

$$Z = \frac{h_p h_o [(2\pi h_p - d\phi P_r/f_p) i \gamma_{CCD} - P_r/f_p (d^2 + h_p^2)]}{h_o [(2\pi h_p - d\phi P_r/f_p)d + P_r/f_p (d^2 + h_p^2)\phi] - h_p (2\pi h_p - d\phi P_r/f_p) i \gamma_{CCD}} \quad (9)$$

The projected light being non-parallel, the coordinates X, Y of each point need to be corrected as a function of pixel index $(i, j)$, where $\gamma_{CCD}$ is the magnification of the CCD entrance lens in mm/pixel:

$$X = \frac{(Z + h_o)}{h_o} \gamma_{CCD} i \quad Y = \frac{(Z + h_o)}{h_o} \gamma_{CCD} j \quad (10)$$

The determination of these parameters is crucial for the accuracy of the final results.

## 2.6 Calibration

If the simplifying assumptions of parallel illumination and constant frequency apply, Eq. (8) holds true, and relief is proportional to the phase. The proportionality constant is calculated either by measuring the two unknown parameters (the angle between projection and observation axis ($\alpha$) and the grating pitch on the reference surface ($p$)) or by a simple calibration using an object of known height. As discussed in the previous paragraph, in the general case of diverging illumination the relief depends on five geometrical parameters and a more complex calibration needs to be implemented. Several calibration procedures have been proposed to measure relief by means of triangulation. Guangjun and Liqun (2000) proposed a combination of vertical and horizontal gratings projected onto a reference plane. San-soni et al. (2000) used a calibration based on an iteration algorithm and the measurement of a rectangular prism shape object of known height. This method is valid in the case of equal projection and observation height ($h_p = h_o$). For the general case (diverging illumination and $h_p \neq h_o$) the method proposed by Breque et al. (2004) uses a rotating reference plane. More details about this method are given in the experimental section (Sect. 3.1) as this is the method we implemented for our study.

## 2.7 The influence of water on the demodulation

Because of light refraction at the air/water interface, bed topography reconstructed by any method based on grid deformation will be distorted relative to the true topography when water is present above the surface. By measuring bed topography with and without the flow, water depth can thus be calculated. As for the setup proposed by Huang et al. (2010), we assumed that the flow is locally uniform (free surface parallel to the bed surface). This assumption was verified in similar experimental conditions in the work of Malverti et al. (2008): water velocities obtained experimentally compared rigorously with predicted values from laminar flow Navier-Stokes equations.

Light produces a bright spot when hitting the surface. From the point of view of the camera, the bright spot is observed in the direction corresponding to the angle $i_c$ with respect to the normal (Fig. 3):

$$\frac{\sin i_p}{\sin r_p} = \frac{\sin i_c}{\sin r_c} = n_{\text{water}} \quad (11)$$

where $n_{water}$ is the refractive index of water. The difference in elevation ($\Delta h$) between the apparent location of the bright spot ($h$) and the true depth ($H$) is obtained by trigonometry, then the true depth can be deduced:

$$H = \frac{1}{1 - \frac{\tan r_p + \tan r_c}{\tan i_p + \tan i_c}} \Delta h \quad (12)$$

The difference between topographies obtained with and without flow will give a distorted bathymetry ($\Delta h$) that is always shallower than the true bathymetry ($H$). Incidence angles can be calculated in any point by simple trigonometry: $i_p, i_c = f(x, y)$, and the correction factor in Eq. (12) can be applied to any point on the surface.

## 3 Experimental results

Figure 1 is a schematic of the experimental setup used to reproduce microscale braided channels. Details of experiments similar to the ones described here can be found in Métivier and Meunier (2003). The experiments were conducted at the IPGP experimental facility in a microscale flume made of plexiglass (1.5 m long and 0.75 m wide) with a fully adjustable bed (height and tilt). The initial slope of the bed was set to approximately 0.05. The sediment in these experiments was composed of glass beads with a D50 of 250 micron and a density of 2500 kg m⁻³. The initial condition for each experiment was a flat bed with a straight channel (1 cm deep x 2 cm wide) carved down the middle of the bed. The flow was laminar and a fully braided morphology spontaneously developed and reached steady state (sediment input equals sediment output) after several hours. Sediment and water were fed continuously at a constant rate at the upstream end of the flume. Sediment transport was measured at the outlet
<figure><img src="image_79.png" /><figcaption>Fig. 3. Light refraction into a water layer.</figcaption></figure>

using a constant head tank which rests on a balance. Flow discharge was typically $2.5 \times 10^{-5} \text{ m}^3 \text{ s}^{-1}$ and sediment discharge was typically around $8 \times 10^{-8} \text{ m}^3 \text{ s}^{-1}$. Pieces of plexiglass were placed at an angle along the sidewalls in order to deflect the flow and discourage it from sticking to the sidewalls. Glare (specular reflection) produced by the free surface of water was reduced by using two cross axis polarizers placed in front of the projector and the camera respectively.

## 3.1 Light3D description, calibration and verification experiments

The algorithm we used was developed at the Laboratory of Solid State Mechanics of the University of Poitiers, France and is available as a commercial, software package called Light3D. The software allows the user to choose the number of phase-shifted projected patterns (3, 8, 16 or 32) and phase unwrapping can be performed with or without 8 gray code images. The encoded patterns are sent from the computer to a Sanyo PLV-Z5 video-projector. A black and white, 8-bit, $1280 \times 1024$ pixels μeye Stemmer Imaging CCD camera was used for image acquisition. Since light beams are not parallel, the introduced phase shifting is not constant across the whole image. The algorithm implemented into the software recomputes the true phase shifts using Fast Fourier Transform (FFT). Therefore, the need to have a constant grid pitch over the entire image is also overcome. This requirement is essential when using the Fourier-transform profilometry as proposed by Takeda et al. (1982) (Sect. 2.2).

Inside of the rectangular image, the object can be of any shape. The object mask is determined as follows: if the intensity of the fringe pattern is the same, this indicates that it is an outside point, so the mask is set to zero and is set to 255 for an inside point if the intensities are different. The phase

is obtained modulo $2\pi$ for each internal point (i.e. when mask = 255). Then the phase needs to be unwrapped, but because of the noise or discontinuities, a second mask step is applied, according to a modified Bone algorithm (Brémand, 1994). This method adopts a “quality” measure to guide the unwrapping path by calculating the second order partial derivatives of the phase map and comparing them to a threshold value to detect the discontinuities. A point of discontinuity is then set to zero. When the mask procedure is finished, the phase unwrapping algorithm can be applied. Each point having a mask value of 255 should give only one fringe number. If this is not the case, a new threshold is applied until the discontinuity disappears. The image is scanned from one corner, the phase is unwrapped only if at least one of its 8 neighbours is treated. After the first pass, it is possible that several points could not be unwrapped. Then a second pass is performed from the opposite corner, and the cycle is repeated until all the points are analysed.

The phase-shifting technique of generating phase maps combined with a “quality” based unwrapping algorithm makes this method suitable for the measurement of smooth objects with a very high resolution of measurement. When both extended height measuring range and high resolution are required, the software can generate and use 8 gray code images, at the expense of longer acquisition times.

The software allows the unwrapping of absolute phase maps as well as the unwrapping of the difference between two phase maps. The first type of calculation will produce a relief expressed with respect to the reference plane and the second one the relief deformation with respect to an initial state.

For the results shown here the phase field was calculated with 8 phase-shifted images and robust unwrapping using a series of 8 gray code images. Prior to each experiment, a calibration procedure that uses a rotating plane at two positions separated by a given angle was performed to determine the geometric parameters. Two experimental measurements of phase fields were then carried out by means of phase shifting. The first one corresponds to the reference plane $\phi_{\text{ref}}$ (position 1 of the calibration board in Fig. 1) and the second one to the same reference plane rotated by an angle $\beta$ $\phi_{\text{cal}}$ (position 2 of the calibration board in Fig. 1). The two phase fields obtained for the two positions of the plane are interpolated by a least square fit considering some (40) different lines distributed into the image, according to the algorithm developed by Breque et al. (2000), and given that the relief of the rotating plane is known in every point. The rotation angle $\beta$ was measured by using a digital slope meter. To begin the unwrapping process, it is necessary to locate the initial fringe, which allows setting the origin O of the working reference. By projecting a single horizontal line onto the reference plane, the alignment of the projector/camera system is adjusted so that the projected line corresponds to the rotation axis of the reference plane during calibration. The zero fringe is therefore the same in the reference and calibration
<table><caption>**Table 1.** Symbols and notations used in this paper</caption><thead><tr><th>Symbol</th><th>Description</th><th>Dimension</th></tr></thead><tbody><tr><td>I(x, y)</td><td>intensity of fringe pattern</td><td>[counts]</td></tr><tr><td>a(x, y)</td><td>background intensity</td><td>[counts]</td></tr><tr><td>b(x, y)</td><td>local contrast</td><td>[counts]</td></tr><tr><td>c(x, y), c*(x, y)</td><td>complex contrast coefficients</td><td>[counts]</td></tr><tr><td>φ(x, y)</td><td>phase</td><td>[rad]</td></tr><tr><td>φ'(x, y)</td><td>phase with carrier pattern</td><td>[rad]</td></tr><tr><td>f<sub>0</sub></td><td>carrier frequency</td><td>[L]<sup>-1</sup></td></tr><tr><td>j</td><td>√-1</td><td>[0]</td></tr><tr><td>F(I)</td><td>fourier transform of intensity</td><td>[L]</td></tr><tr><td>A</td><td>fourier transform of background intensity</td><td>[L]</td></tr><tr><td>C, C*</td><td>fourier transform of contrast coefficients</td><td>[L]</td></tr><tr><td>N</td><td>number of phase shifts</td><td>[0]</td></tr><tr><td>δ<sub>i</sub> = i2π/N</td><td>phase shift increment</td><td>[rad]</td></tr><tr><td>I<sub>i</sub>(x, y)</td><td>intensity of fringe pattern for each phase shift</td><td>[counts]</td></tr><tr><td>Z</td><td>relief</td><td>[L]</td></tr><tr><td>p</td><td>projected grating spacing</td><td>[L]</td></tr><tr><td>h<sub>p</sub></td><td>projection distance</td><td>[L]</td></tr><tr><td>h<sub>o</sub></td><td>observation distance</td><td>[L]</td></tr><tr><td>d</td><td>distance between the light source and camera</td><td>[L]</td></tr><tr><td>f<sub>p</sub></td><td>projector focal length</td><td>[L]</td></tr><tr><td>P<sub>r</sub></td><td>grating pitch at the output lens</td><td>[L]</td></tr><tr><td>γ<sub>CCD</sub></td><td>magnification of CCD entrance length</td><td>[L] / pixel</td></tr><tr><td>i<sub>p</sub></td><td>projection incidence angle</td><td>[rad]</td></tr><tr><td>i<sub>c</sub></td><td>camera incidence angle</td><td>[rad]</td></tr><tr><td>r<sub>p</sub></td><td>projection refraction angle</td><td>[rad]</td></tr><tr><td>r<sub>c</sub></td><td>camera refraction angle</td><td>[rad]</td></tr><tr><td>Δh</td><td>difference in elevation between dry and wet relief</td><td>[L]</td></tr><tr><td>H</td><td>true water depth</td><td>[L]</td></tr><tr><td>n<sub>w</sub></td><td>water refraction index</td><td>[0]</td></tr></tbody></table>

planes. This simple calibration allows the determination of the geometrical parameters ($h_p$, $d$, $h_o$, and $P_r/f_p$). These geometrical parameters are difficult to measure manually and accurately. The precision of the calibration depends on the quality of the calculated phase fields: $\phi_{ref}$ and $\phi_{cal}$. A flat, light color, highly diffusive board was used in the calibration procedure. The results of the calibration were the following: $h_p = 2050$ mm, $d = 1857$ mm, $h_o = 1957$ mm, $P_r/f_p = 0.0051$, for a $\gamma_{CCD} = 0.78$ mm/pixel.

In order to evaluate the precision and resolution of the method, two preliminary experiments were performed. Firstly we analysed the distribution of residual noise of a perfectly flat surface and secondly we discussed the results obtained by demodulation on a object of a known height (a reference scale). In order to obtain a perfectly flat surface we used milk, which has the advantage of being not only flat but opaque and of light colour. Figure 4 shows the results obtained as a difference of 2 levels of milk. One can also see the bubbles that existed on the first surface (black marks) and on the second surface (white marks). Except for the bubbles zone, the surface was perfectly flat, with a “residual”

<figure><img src="image_80.png" /><figcaption>**Fig. 4.** Relative topography between two milk level surfaces</figcaption></figure>
<figure><img src="image_81.png" /><figcaption>Fig. 5. (a) Topography of a stairs object. (b) Cross-section through the object obtained at a length = 230 mm.</figcaption></figure>

modulation still visible. The milk depth level measured on the flat zone has a value of 2.42 mm with a standard deviation remarkably low (0.02 mm). This precision is comparable to the one reported by using a combination of stereo camera system with structured light illumination (ATOS system by GOM, France), used in the work of Turowski et al. (2006) (20 µm), and a real improvement with respect to the precision reported by Huang et al. (2010) (0.6 mm).

Figure 5a shows the relief obtained with our reference scale. One can notice the shadow zone behind the object, and the black zones on the top of it where a slopemeter is mounted on. The first stair is 10mm high and the 5 others are 2 mm high. A cross section passing at a length equal to 230 mm is shown in Fig. 5b. As one can notice, the stairs are well resolved. The uncertainty measured on each stair and on the flume board is about 0.1 mm and is mainly due to the non-uniformity of the flume board and of the white painting applied on the scale. The uncertainty on the water depth is the same as on the topography, provided that the water refraction is taken into consideration properly (Eq. 12). We verified this statement by introducing our scale into a shallow tank filled with water and by calculating the water depth from the difference between the topography obtained with and without water.

## 3.2 Experimental braided channel: topography and water depth

Bed topography was measured one minute after the flow was turned off (to ensure the bed was uniformly drained). The impact of stopping and starting the flow on bed topography was negligible. Water depth was calculated by differencing the topography obtained with the flow on immediately before it was shut off and with the flow off, taking into account the water refraction as described in Sect. 2.7. This is facilitated in our case by the fact that our experiments are conducted under laminar conditions and the slopes are small without steep variations, so that water free surface can be considered uniform and parallel to bed surface at the scale of moire fringe.

Bed topography and flow depth measured during an experiment are shown in Fig. 6. Dry bed topography is between

<figure><img src="image_82.png" /><figcaption>Fig. 6. (a) Bed topography obtained by switching of water flow and (b) water depth.</figcaption></figure>

−4 and 4 mm, while maximum water depth is approximately 3 mm (flow depth increases in the negative direction). Figure 6 shows bed topography and flow depth in cross-section at 530 mm downstream from the top of the study reach. Several occupied and abandoned channels are visible along this transect.

We excluded data from approximately 25 mm on either side of the flume where the roughness elements were located. A PDF of flow depth between −2.0 and 2.0 mm was calculated for the study area (red curve in Fig. 8). The average value of the water depth and the standard deviation were −0.24 mm and 0.43, respectively. If we consider positive flow depth noise and if we consider noise distribution as being symmetrical with respect to zero, we obtain the green curve in Fig. 8, which has a standard deviation of 0.28 mm. The difference between total depth distribution and noise results in the true depth distribution shown in blue in Fig. 8. The mean and standard deviation of this distribution are −0.49 and 0.41 mm, respectively. The area of the flume occupied by flow after correcting for noise is 49%.

In a separate set of microscale braided river experiments, the focus was on the dynamically evolving topography, and particularly on channel characteristics and dynamics. The study demonstrates how a moiré method can be used to characterize a full range of statistics beyond just topography and depth. Calibration proceeded as above, with slightly different camera distances resulting in 0.93 mm/pixel resolution. Data were again collected via a series of 8 phase-shifted images, and deconstructed with Light3D software. Topography data for the measurements shown here were collected after the system had evolved from a flat bed to the dynamic equilibrium condition of equal input/output sediment flux, after one hour of steady water and sediment input. The 60 cm by 60 cm focus region of measurements was located in the downstream half of the system, to avoid potential inlet effects. The bulk relief values of this area are characterized in Fig. 9a, generally within a range of −5 mm incision and +2 mm accretion. An initial scan with water running was followed by a scan
<figure><img src="image_83.png" /><figcaption>**Fig. 7.** Transects into wet, dry and water depth topographies at length value of 530 mm.</figcaption></figure>

<figure><img src="image_84.png" /><figcaption>**Fig. 8.** PDFs of water depth distribution (red curve) with the noise (green curve) removed (blue curve).</figcaption></figure>

taken two minutes after water flow had been stopped, and water depths were determined following the procedure outlined in Sect. 2.7 (Fig. 9b). Water depth values peak at a value of about 0.5 mm, and vary down to depths of 2–3 mm.

<figure><img src="image_85.png" /><figcaption>**Fig. 9.** Histograms of (a) relief values of developed braided river topography and (b) water depth values within channels. Measurements are with respect to an initial flat bed.</figcaption></figure>

In this study data were also collected on various aspects of the channels. In order to isolate the channelized areas in a systematic manner, small, noisy water depth values were eliminated by setting a threshold on the water depth data set to only include depths up to the third statistical quartile of depth values (see Fig. 10a for a color-depth image of channelized water depths as isolated by this threshold). Having used this method to isolate the channels, channel characteristics such as number, width, length, and spacing were measured. The within-channel slope profiles could also be tracked; see Fig. 9b for profiles of the five longest continuous threads of the study region, where it was possible to resolve sub-millimeter-scale features on the channel floors. Additionally, within-channel slopes at the banks were measured perpendicular to flow, as defined over a given distance inward of channel boundaries (Fig. 9c). Measurements such as these can improve our understanding of channel geometry organization, especially when combined with information about local Shields stress and bank migration rate, both of which can be extracted from a time-sequence of relief and water depth data.

Finally, by differencing successive topographic scans (Fig. 11), we can produce a series of maps of erosional and depositional features as they evolve through time. With the high temporal resolution feasible with this method of measurement, we can gain new insight into the dynamic evolution of braided rivers.
<figure><img src="image_86.png" /></figure>

**Fig. 10.** (a) Planview of the analyzed section of the braided river system. Colors correspond to water depth values. The direction of water flow is from the top to the bottom of the image; channelized threads are visible as vertical lines. (b) Within-channel elevation profiles of the five longest continuous channel threads. Profiles are plotted from an x-value corresponding to the distance down Sect. (a) at which the channel begins. (c) Histogram of cross-stream slope values within the channels, at the channel banks (outer 3 mm). Channel boundaries are defined using the wet fraction 3rd-quartile threshold discussed in the text.

<figure><img src="image_87.png" /></figure>

**Fig. 11.** (a) Topography scan of a developed braided river system, with flow direction from top to bottom of the image, and (b) scan repeated five minutes later. By subtracting the first scan from the second, we can produce (c) a map of the erosional and depositional features of the landscape.

## 4 Conclusions

The moiré method provides a way of characterizing evolving experimental systems at very high spatial and temporal resolutions. This detailed description is especially useful for systems as complex and ephemeral as braided rivers, for which patterns of channel organization and feature migration are sensitive to local topographic conditions that vary widely. With high resolution data, dynamic features can be characterized with respect to local constraints, and in this way we

can learn about processes that operate at small scales. With significant statistics, we can also identify bulk system patterns, and either track their evolution through time or test their sensitivity to changing boundary conditions. Then, with information about both scales of activity, we can hopefully connect the evolution driven by small-scale processes to the emergence of observed large-scale patterns. The ease and resolution of this method of data collection will enable us to describe the evolution and organization of braided river systems with more robust statistics than have yet been reported. With this description, we will hopefully be able to parse the relevant processes and timescales that drive braided river evolution.

A certain number of requirements have to be fulfilled in order to obtain the best raw images as input data regardless of which method is used for phase measurement. Moiré methods work well for light color, relatively uniform, diffusive substrate. As in any triangulation method, shadow zones may appear, where all the information is lost. In order to image those zones, the projection and observation angles should be changed. The method is relatively insensitive to stray light since the information is encoded into the phase. Background light should, if possible, be minimized, in order to take full advantage of the dynamic range of elementary detectors of the CCD camera.

Under optimised experimental conditions we report here uncertainties on both topography and bathymetry varying between 20 to 100 µm with sub-millimeter spatial resolution. Under dynamic conditions, due to the fact that this is a multiple-shot method, this uncertainty increases, since the different images are recorded at different moments in time. For time evolving systems, such as the examples shown here, the noise on the water depth distribution increases to about 300 µm, while the uncertainty on the topography stays unchanged, since it is measured without flow, in a frozen condition.

*Acknowledgements.* We thank Y. Gamblin and A. Vieira e Silva of the Institut de Physique du Globe de Paris for their technical assistance in designing and realizing the experimental set-up.

Edited by: D. R. Gröcke

cnrs
INSU
Institut national des sciences de l'Univers

The publication of this article is financed by CNRS-INSU.
## References

Armstrong, L.: Bank erosion and sediment transport in a microscale straight river, Ph.D. thesis, Université Paris7-Denis Diderot, France, 2003.

Ashmore, P.: Laboratory modelling of gravel bed braided stream morphology, Earth. Surf. Process. Landforms, 7, 201–225, 1982.

Aureli, F., Maranzoni, A., Mignosa, P., and C., Z.: Dam-break flows: acquisition of experimental data through an imaging technique and 2-D numerical modeling, J. Hydraul. Eng., 134, 10891101, doi:10.1061/(ASCE)0733-9429(2008)134:8(1089), 2008.

Brémand, F.: A phase unwrapping technique for object relief determination, Optics and Lasers in Engineering, 21, 49–60, doi:10.1016/0143-8166(94)90059-0, 1994.

Brémand, F., Doumalin, P., Dupré, J., Hesser, F., and Valle, V.: Optical Techniques for Relief Study of Mona Lisas Wooden Support, Chap. in: Experimental Analysis of Nano and Engineering Materials and Structures, Springer, Netherlands, 2007.

Breque, C., Brémand, F., and Dupré, J.: Characterization of biological materials by means of optical methods of measurement, in: Second International Conference on Experimental Mechanics, edited by Chau, F. and Quan, C., 463–468 pp., 2000.

Breque, C., Dupré, J., and Brémand, F.: Calibration of a system of projection moiré for relief measuring: biomechanical applications, Optics and Lasers in Engineering, 41, 241–260, doi:10.1016/S0143-8166(02)00198-7, 2004.

Cardenas-Garcia, J., Yao, H., and Zheng, S.: 3-D reconstruction of objects using stereo imaging, Opt. Lasers Eng., 22, 193–213, 1995.

Chandler, J. H., Shiono, K., Rameshwaren, P., and Lane, S. N.: Measuring flume surfaces for hydraulics research using a Kodak DCS460, Photogramm. Rec., 17, 39–61, 2001.

Chiang, F.: Moiré Methods of Strain Analysis, Chap. 6 in: Manual on Experimental Stress Analysis, edited by: Kobayashi, A. S., Soc. for Exp. Stress Anal., Brookfield Center, CT, 51–69 pp., 1983.

Cochard, S. and Ancey, C.: Tracking the free surface of time-dependent flows: image processing for the dam-break problem, Exp. Fluids, 44, 59–71, 2008.

Coleman, S. and Eling, B.: Sand wavelets in laminar open-channel flows, J. Hydraul. Res., 38, 331–338, 2000.

Darboux, F. and Huang, C.: An instantaneous-profile laser scanner to measure soil surface microtopography, Soil Sci. Soc. Am. J., 67, 92–99, 2003.

Devauchelle, O., Malverti, L., Lajeunesse, E., Josserand, C., Lagree, P. Y., and Métivier, F.: Rhomboid beach pattern: A laboratory investigation, J. Geophys. Res., 115, F02017, doi:10.1029/2009JF001471, 2010a.

Devauchelle, O., Malverti, L., Lajeunesse, E., Josserand, C., Lagree, P. Y., and Thu-Lam, K. N.: Stability of bedforms in laminar flows with free-surface: from bars to ripples, J. Fluid Mech., 642, 329–348, 2010b.

Dreano, J., Valance, A., Lague, D., and Cassar, C.: Experimental study of sediment flux controls on transient and steady-state dynamics of bedforms in supply-limited conditions, Earth Surf. Proc. Land., 35, 1730–1743, 2010.

Federici, B. and Paola, C.: Dynamics of channel bifurcations in noncohesive sediments, Water Resour. Res., 39, 1162, doi:10.1029/2005JF000306, 2003.

Ghiglia, D. C. and Pritt, M. D.: Two-dimensional Phase Unwrapping: theory, algorithms, and software, Wiley Science, 1998.

Guangjun, Z. and Liqun, M.: Modeling and calibration of grid structured light based 3-D vision inspection, J. Manuf. Sci. Eng., 122, 734–738, 2000.

Halioua, M. and Liu, H.-C.: Optical three-dimensional sensing by phase measuring profilometry, Optics and Lasers in Engineering, 11, 185–215, doi:10.1016/0143-8166(89)90031-6, 1989.

Huang, A., Capart, H., Chen, R., and Huang, M.: Laser scanning technique for the acquisition of digital elevation models of debris flow material, in: Fourth Int Conf on Debris-Flow Hazard Mitigation, Chengdu, China, September 2007, 539–546 pp., 2007a.

Huang, A., Capart, H., Chen, R., and Huang, M.: Experimental analysis of the seepage failure of a sand slope, in: Fourth Int Conf on Debris-Flow Hazard Mitigation, Chengdu, China, September 2007, 259267, 2007b.

Huang, M., Huang, A., and Capart, H.: Joint mapping of bed elevation and flow depth in microscale morphodynamics experiments, Experiments in Fluids, 1–14 pp., doi:10.1007/s00348-010-0858-4, 2010.

Huntley, J. and Saldner, H.: Shape measurement by temporal phase unwrapping: comparison of unwrapping algorithms, Meas. Sci. Technol., 8, 986–992, 1997.

Kozłowski, J. and Serra, G.: New modified phase locked loop method for fringe pattern demodulation, Opt. Eng., 36, 2025–2030, 1997.

Kozłowski, J. and Serra, G.: Analysis of the complex phase error introduced by the application of the Fourier transform method, J. of Modern Opt., 46, 957–971, 1999.

Lague, D., Crave, A., and Davy, P.: Laboratory experiments simulating the geomorphic response to tectonic uplift, J. Geophys. Res., 108, 2008, doi:10.1029/2002JB001785, 2003.

Lajeunesse, E., Malverti, L., and Charru, F.: Bed load transport in turbulent flow at the grain scale: Experiments and modeling, J. Geophys. Res., 115, F04001, doi:10.1029/2009JF001628, 2010a.

Lajeunesse, E., Malverti, L., Lancien, P., Armstrong, L., Métivier, F., Coleman, S., Smith, C., Davies, T., Cantelli, A., and Parker, G.: Fluvial and Submarine Morphodynamics of Laminar and Near-Laminar Flows: a synthesis, Sedimentology, 57, 1–26, 2010b.

Lancien, P., Métivier, F., Lajeunesse, E., and Cacas, M.: Incision dynamics and shear stress measurements in submarine channels experiments, in: River, Coastal and Estuarine Morphodynamics 2005, edited by: Parker, G. and Garcia, M. H., Taylor and Francis Group, London, 527–533 pp., 2005.

Liebling, M., Blu, T., and Unser, M.: Complex-wave retrival from a single off-axis hologram, J. Opt. Soc. Am. A, 21, 367–377, 2004.

Malverti, L., Lajeunesse, E., and Métivier, F.: Experimental investigation of the response of an alluvial river to a vertical offset of its bed, in: 5th IAHR Symposium on River, Coastal and Estuarine Morphodynamics, edited by: Dohmen-Janssen, C. and Husche, S., 179–184 pp., 2007.

Malverti, L., Lajeunesse, E., and Métivier, F.: Small is beautifull: upscaling microscale experimental results to the size of natural rivers, J. Geophys. Res., 113, F04004, doi:10.1029/2007JF000974, 2008.

Métivier, F. and Meunier, P.: Input and Output mass flux correlations in an experimental braided stream. Implications on the
dynamics of bed load transport, J. Hydrol., 271, 22–38, 2003.

Métivier, F., Lajeunesse, E., and Cacas, M.: Submarine Canyons in the Bathtub, Journal of Sed. Res., 75, 6–11, 2005.

Métivier, F., Narteau, C., Lajeunesse, E., Devauchelle, O., Liu, Y., and Ye, B.: A new simple integral technique to analyze bedload transport data, in: 7th IAHR Symposium on River, Coastal and Estuarine Morphodynamics, 179–184 pp., 2011.

Meunier, P.: Braided rivers dynamics, Ph.D. thesis, Institut de Physique du Globe de Paris, 2004.

Ni, W.-J. and Crave, H.: Groundwater drainage and recharge by networks of irregular channels, J. Geophys. Res., 111, F02014, doi:10.1029/2005JF000410, 2006.

Paola, C.: Incoherent Structure: Turbulence as a Metaphor for Stream Braiding, in: Coherent Flow Structures in Open Channels, edited by: Ashworth, P., Bennett, S., Best, J., and McLeland, S., 705–723 pp., John Wiley & Sons Ltd, 1996.

Paola, C.: Modelling stream braiding over a range of scales, in: Gravel Bed Rivers, edited by: Mosley, M., 11–46 pp., 2001.

Paola, C.: Quantitative models of sedimentary basin filling, Sedimentology, 47, 121–178, 2010.

Patorski, K.: Handbook of the moiré fringe technique, Chap. 10 and 13, Elsevier Science Publishers, 1993.

Pouliquen, O. and Forterre, Y.: Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane, J. Fluid Mech., 453, 131–151, 2002.

Reitz, M. D., Jerolmack, D. J., and Swenson, J. B.: Flooding and flow path selection on alluvial fans and deltas, Geophys. Res. Lett., 37, L06401, doi:10.1029/2009GL041985, 2010.

Rice, C., Wilson, B. N., and Appleman, M.: An instantaneous-profile laser scanner to measure soil surface microtopography, Comput. Electron. Agric., 3, 97–107, 1988.

Sansoni, G., Corini, S., Lazzari, S., Rodella, R., and Docchio, F.: Three-dimensional imaging based on Gray-code light projection: characterization of the measuring algorithm and development of a measuring system for industrial applications, Appl. Opt., 36, 4463–4472, 1997.

Sansoni, G., Carocci, M., and Rodella, R.: Calibration and performance evaluation of a 3-D imaging sensor based on the projection of structured light, IEEE Trans. Instrum. Meas., 49, 628–635, 2000.

Sansoni, G., Patrioli, A., and Docchio, F.: OPL-3D: A novel, portable optical digitizer for fast acquisition of free-form surfaces, Review of Scientific Instruments, 74, 2593–2603, doi:10.1063/1.1561602, 2003.

Schmit, J. and Creath, K.: Extended averaging technique for derivation of error compensating algorithms in phase shifting interferometry, Appl. Opt., 34, 3610–3619, 1995.

Schumm, S. A.: The fluvial system, John Wiley & Sons Inc., 1977.

Schumm, S. A., Morley, M. P., and Weaver, W. E.: Experimental fluvial geomorphology, Wiley-Interscience, New York, 1987.

Takeda, M.: Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: An overview, Industrial Metrology, 1, 79–99, 1990.

Takeda, M., Ina, H., and Kobayashi, S.: Fourier-transform method of fringe-pattern analysis for computer based topography and interferometry, J. Opt. Soc. Am., 72, 156–160, 1982.

Tal, M. and Paola, C.: Dynamic single-thread channels maintained by the interaction of flow and vegetation, Geology, 35, 347–350, 2007.

Turowski, J. M., Lague, D., Crave, A., and Hovius, N.: Experimental channel response to tectonic uplift, J. Geophys. Res., 111, F03008, doi:10.1029/2005JF000306, 2006.

Weiler, M. and Flühler, H.: Inferring flow types from dye patterns in macroporous soils, Geoderma, 120, 137–153, 2004.
# Microwave-heating laboratory experiments for planetary mantle convection

A. Limare<fnref n="1" />,†, K. Vilella<fnref n="1" />, E. Di Giuseppe<fnref n="1" />,<fnref n="3" />, C. G. Farnetani<fnref n="1" />, E. Kaminski<fnref n="1" />, E. Surducan<fnref n="2" />, V. Surducan<fnref n="2" />, C. Neamtu<fnref n="2" />, L. Fourel<fnref n="1" /> and C. Jaupart<fnref n="1" />

(Received 19 December 2014; revised 20 April 2015; accepted 18 June 2015)

Thermal evolution of telluric planets is mainly controlled by secular cooling and internal heating due to the decay of radioactive isotopes, two processes that are equivalent from the standpoint of convection dynamics. In a fluid cooled from above and volumetrically heated, convection is dominated by instabilities of the top boundary layer and the interior thermal structure is non-isentropic. Here we present innovative laboratory experiments where microwave radiation is used to generate uniform internal heat in fluids at high Prandtl number (>300) and high Rayleigh–Roberts number (ranging from $10^4$ to $10^7$), appropriate for planetary mantle convection. Non-invasive techniques are employed to determine both temperature and velocity fields. We successfully validate the experimental results by conducting numerical simulations in three-dimensional Cartesian geometry that reproduce the experimental conditions. Scaling laws relating key characteristics of the thermal boundary layer, namely its thickness and temperature drop, to the Rayleigh–Roberts number have been established for both rigid and free-slip boundary conditions. A robust conclusion is that for rigid boundary conditions the internal temperature is significantly higher than for free-slip boundary conditions. Our scaling laws, coupled with plausible physical parameters entering the Rayleigh–Roberts number, enable us to calculate the mantle potential temperature for the Earth and Venus, two telluric planets with different mechanical boundary conditions at their surface.

**Key words:** boundary layers, geophysical and geological flows, mantle convection

## 1. Introduction

The secular evolution of the Earth through geological time is best tackled from a thermal perspective because geological events are driven by internal heat dissipation. The present-day energy budget is such that the rate of heat loss is approximately twice as large as the amount of heat released by the radioactive decay of long-lived radioactive isotopes ($^{238}$U, $^{235}$U, $^{232}$Th and $^{40}$K), implying that the Earth is currently

<footnote n="1">¹Institut de Physique du Globe de Paris, Université Paris Diderot, CNRS, 1 rue Jussieu,
75238 Paris, France</footnote>

<footnote n="2">²National Institute for Research and Development of Isotopic and Molecular Technologies,
67-103 Donath St., 400293 Cluj-Napoca, Romania</footnote>

<footnote n="3">³CEMEF, MINES ParisTech, CNRS, CS 10207, 06904 Sophia Antipolis, France</footnote>

<footnote>†Email address for correspondence: limare@ipgp.fr</footnote>