Convective transport in crustal magma bodies

Journal of Volcanology and Geothermal Research, 19 ( 1 9 8 3 ) 4 5 - - 7 2

45

Elsevier Science P u b l i s h e r s B.V., A m s t e r d a m - - P r i n t e d in T h e N e t h e r l a n d s

CONVECTIVE T R A N S P O R T IN C P U S T A L MAGMA BODIES

H.C. H A R D E E

Sandia National Laboratories, Geophysics Division 1541, Albuquerque, NM 87185 (U.S.A. ) (Received April 5, 1 9 8 2 ; reversed and a c c e p t e d N o v e m b e r 21, 1 9 8 2 )

ABSTRACT H a r d e e , H.C., 1 9 8 3 . C o n v e c t i v e t r a n s p o r t in crustal m a g m a bodies. J. Volcanol. G e o t h e r m . Res., 19: 4 5 - - 7 2 . C o n v e c t i v e t r a n s p o r t in enclosures is reviewed as it relates to m a g m a bodies. Convective effects are discussed t h a t are p e c u l i a r t o t h e rheological c h a r a c t e r of m a g m a . T h e large P r a n d t l n u m b e r o f m a g m a delays o n s e t of t u r b u l e n c e and causes velocity disturb a n c e s t o e x t e n d far b e y o n d t h e t h e r m a l b o u n d a r y layer and well i n t o t h e fluid core. A b o u n d a r y layer analysis is used to s h o w t h a t c o n v e c t i v e t r a n s p o r t in e n c l o s e d m a g m a b o d i e s p r o d u c e s an u p w a r d h e a t flux w h i c h can lead to m i g r a t i o n of m a g m a by u p m e l t ing. Multicellular c o n v e c t i o n is discussed and s h o w n to p r o d u c e u p f l o w at the walls in cases w h e r e extensive cooling occurs at t h e t o p of t h e m a g m a c h a m b e r . Multicellular conv e c t i o n can lead t o b i m o d a l m a g m a s y s t e m s and can also p r o d u c e vertical z o n a t i o n and layering. E x a m p l e s of oscillatory c o n v e c t i v e c i r c u l a t i o n are given and this p h e n o m e n a is discussed as a possible m e c h a n i s m for e r u p t i v e cycles and for vertical z o n a t i o n in m a g m a c h a m b e r s . N o n - N e w t o n i a n r h e o l o g y is also e x a m i n e d as it applies to c o n v e c t i o n in m a g m a .

INTRODUCTION

Models of convective transport for magma in the lithosphere or crust usually deal with ascent of magma through the lithosphere, convection within the magma body, or interaction of the magma with a hydrothermal zone. Although the emphasis of this paper is on convective transport in magma bodies in the lithosphere, models for the formation and ascent of magma in the asthenosphere are important in understanding the character and configuration of magma as it enters the lithosphere. A number of models have been proposed for the formation of magma bodies in the asthenosphere and for the ascent of these bodies through the asthenosphere and lithosphere (Yoder, 1976; Oxburgh, 1980; Shaw, 1980; Spera, 1980a). Most models of magmatic ascent are based on the assumption of a deep heat source which produces rising diapirs or columns of molten magma in the asthenosphere (Fedotov, 1975, 1976, 1979; Marsh, 1978). In some models the magma continues to rise as a plume or column (Fedotov, 1975) while in other cases Sandia Laboratories Report SAND 82-0907J.

0377-0273/83]$03.00

© 1983 Elsevier Science Publishers B.V.

46 the column is assumed to pinch off into an isolated spheroidal body ~t rising magma (Whitehead and Luther, 1975; Fedotov, 1976; Marsh, 1979: Williams and McBirney, 1979). Other variations on this general model exist. For instance, Turcotte and Ahem (1978) proposed a porous flow model in which magma ascends in the asthenosphere through a distributed collection of conduits. Once a diapiric magma body reaches the lithosphere, the magma can continue to ascend in a similar diapiric shape although the mechanism of propagation may be altered. In the asthenosphere, buoyant viscous fluid motion may be the principal mechanism of convective diapir ascent whereas in the lithosphere, other mechanisms may be dominant such as melting, fracturing, or stoping (Yoder, 1976; Shaw, 1980). Magma may also be transported through the lithosphere by a process of forced convection (intrusion) through small conduits. Again, fracturing, possibly coupled with some melting, may be a dominant mechanism in this forced convective transport process. Even if fracturing is a dominant factor in magma ascent in the lithosphere, convective transport is still a subject of major importance in understanding magma bodies in the crust. Natural, buoyancy
47 on the composition of silicate magmatic liquids during ascent of the magma. Shaw (1974) analyzed heat and mass transfer in magma chambers using conventional Newtonian-liquid boundary layer models. Shaw also examined the effect of chemical exchange with convective boundary layers. Spera et al. (1981, 1982) have extended boundary layer analysis to include the case in which the viscosity varies by several orders of magnitude across the boundary layer. Bartlett (1969) has examined the combined effects of convection and crystal settling in solidifying magma bodies. Bhattacharji (1974, 1979) has investigated the use of scale models for natural convection simulation of magmatic systems. Many models of magmatic ascent assume a stoping or fracturing process, frequently in combination with convection and melting (Weertman, 1971; Anderson and Grew, 1977; Shaw, 1980; Spera, 1980a). Convective circulation in the diapir is still important in these models since the internal convective transport process may establish thermal gradients ahead of the ascending melt which affect and in some cases control the fracturing process. Also, the convective circulation may control the rate of removal of stoped or fractured material at the top of an ascending magma b o d y (McBirney, 1959). The analysis of convective transport in enclosed fluid bodies usually involves numerical methods such as finite difference or finite element techniques or analytical methods, such as perturbation or boundary layer techniques. Each technique has its own advantages and disadvantages. Finitedifference numerical techniques work well for transient convection problems b u t have difficulty with irregular or varying boundaries (Brown, 1967: Larson et al., 1977). Finite element techniques handle irregular boundaries very well b u t are best suited for steady state convective transport problems (Gartling, 1976; Larson et al., 1977). Numerical techniques have considerable difficulty in handling convective transport problems at large Rayleigh number for fluids with large Prandtl number and unfortunately this includes most magmas. This difficulty does not seem to be caused by any numerical or physical instability associated with large Prandtl number. Rather it appears to be due to the vastly different grid size required for the temperature field and velocity field. Analytical perturbation techniques can be used on enclosed convection problems when the geometry is relatively simple and the Rayleigh number is small (Shaposhnikov, 1952; Batchelor, 1954; Mack and Hardee, 1968; Mack and Bishop, 1968; Hardee, 1974). At high values of Rayleigh number, boundary layer techniques give good results and can be used with a variety of smoothly varying boundary shapes (Hardee, 1974). By combining boundary layer techniques and perturbation techniques, virtually the entire Rayleigh n u m b e r range of possible interest can be covered (Brown, 1967~ Hardee, 1974). Another type of analytical technique known as the "Parameterized Approach" is very useful for studying the cooling times of magma bodies (Spera, 1980b). The "Parameterized Approach" assumes that convection in the magma chamber is sufficient to keep the entire chamber at

a uniform but slowly decreasing temperature as the magma cools by solid. ification at the boundary and conduction of heat into the surrounding country rock. This assumption is similar to the classical heat transfer assumption of "negligible internal resistance" or "Newtonian cooling" (Schneider~ 1957). Other techniques for the study of convective transport in magma include laboratory simulation experiments and field measurements in lava lakes or above suspected shallow crustal magma bodies. Laboratory simulation experiments usually involve the use of waxes (Hardee and Larson, 1977) or other liquids IHardee, 1974; Whitehead and Luther, 1975; McBirney, 1980; Turner and Gustafson, 1981) to simulate a model magma system. It is very difficult in such laboratory experiments to simulate or model both the Prandtl number and Rayleigh number of magma properly, yet dimensional analysis (Chapman, 1967) demonstrates that these two parameters must be properly modeled in any natural convective transport simulation. In addition to the simple simulation experiments using waxes or other fluids, there are precise laboratory techniques, such as laser holographic interferometry {Schimmel, 1976~ Larson et al., 1977}, which can be used to visualize and analyze isotherms and thereby determine convective heat transfer coefficients for enclosed convecting fluids in the laboratory. This technique also suffers the same difficulties as numerical methods in being unable to handle large Prandtl number fluids. When large Prandtl number fluids are used with laser holographic interferometry, the isotherm fringes are compressed to such a degree that a determination can not be made of the temperature gradient and heat transfer rate at the wall. This problem with large Prandtl number can be overcome to some degree by using a laser anemometer instrument (Durst et al., 1976) to measure velocities in the fluid instead of temperature. This approach does not lead directly to a determination of the heat transfer coefficient; however, the experimental measurement of the velocity field can be used to verify heat transfer numerical codes which in turn can be used to determine heat transfer rates at large Prandtl number. The effects of large Prandtl number can also be treated by running laboratory convection experiments in the desired fluid such as molten lava. Laboratory convection experiments have been run in furnaces at atmospheric pressure using 10- to 20-liter samples of degassed molten basalt (Hardee, 1981; Hardee and Dunn, 1981). Recently, enclosed convection tests were run in a high pressure furnace using a l~-liter sample of water-saturated molten basalt at a pressure of 2.0 kbar and temperature of 1300°C which corresponds to typical conditions in a magma body at 7 km depth in the crust (Dunn et al., 1981: Dunn, 1982). Other experimental techniques such as field experiments in lava flows and in ponded molten lava (Shaw et al., 1968; Pinkerton and Sparks, 1978; Peck, 1978~ Hardee, 1979, 1980) or field measurements in fossilized or exhumed magma bodies (McBirney and Noyes, 1979; Norton and Taylor, 1979), can sometimes be used in studies of convective transport.

49 CONVECTIVE CIRCULATION AND THE EFFECT OF PRANDTL NUMBER A spherical body of magma is a simple geometry for the study of enclosed convective magma transport. Since the spherical body is cooled at its contact with the surrounding medium, magma at the boundary becomes denser and descends while magma at the center rises due to continuity of flow. At low values of Rayleigh number, the entire spherical magma body will be slowly circulating in a uniform torodial pattern (Pustovoit, 1958) as indicated in the typical case shown in Fig. 1. At larger Rayleigh mtmbers,

Fig. 1. Typical torodial circulation patterns or steamlines in a cooling spherical magma body at low Rayleigh number, R a = 104. thermal boundary layer flow will develop and thermal gradients and major velocity gradients will be concentrated near the wall. If magma were a typical fluid with a Prandtl number near unity, then at large values of Rayleigh number all the thermal and velocity changes would be concentrated in a thin boundary layer region near the wall as in Fig. 2a. Liquidus magmas and subliquidus magmas are large Prandtl number fluids and for practical purposes can be treated as if they had an infinite Prandtl number (Hardee, 1981). For typical liquidus and subliquidus magmas at large Rayleigh numbers (106 and greater), the thermal boundary layer will occupy a thin region near the wall and the corresponding velocity boundary layer will go from zero velocity at the wall to peak velocity at a position corresponding to the edge of the thermal boundary layer (Kuiken, 1968). The complete velocity boundary layer or region of appreciable velocity change, as shown in Fig. 2b,c, will extend throughout the enclosed region because of the large Prandtl number (Gebhart, 1961; Hardee, 1974). The assumption of an infinite Prandtl number fluid actually simplifies convective transport calculations, since certain cross-convection terms drop out of the Navier Stokes equations of fluid motion or the corresponding reduced form of the bound-

THERMAL BOUNDARY LAYER THICKNESS

i

j

, I'

DISTRIBUTION

ORDINARY FLUID Pr = 1

MAGMA Pr = 100 MODERATE Ra

MAGMA Pr 100 HIGH Ra =

Fig, 2. Typical velocity and t e m p e r a t u r e profiles for a c o n v e c t i n g fluid in a spherical enclosure at m o d e r a t e to large Rayleigh n u m b e r s , a. O r d i n a r y fluid w i t h P r = 1. b. Magma with P r ~ 100 at m o d e r a t e R a . c. Magma w i t h P r ~ 100 at large R a .

ary layer equations (Hardee, 1966, 1974; Lin and Chao, 1974). The assumption of infinite Prandtl number is valid when the Prandtl number exceeds the value 100 (Grober et al., 1961), and this is always the case for crustal

magmas (Hardee, 1981). TURBULENT CONVECTION If a magmatic body is large, the Rayleigh number, [pg(3(Tf--Tw)X ~/pa] being a function of length cubed, will also be large and the question of whether the convective circulation is laminar or turbulent becomes very important. For ordinary fluids, natural convective circulation will become turbulent when the Rayleigh number exceeds 109 (Rohsenow and Choi, 1961; Elder, 1965b). In enclosed bodies of fluid the onset of turbulence also appears to depend somewhat on the aspect ratio of the enclosure (Elder, 1965b). A large value of Prandtt number however tends to raise the Rayleigh number turbulent transition value. Hieber and Gebhart (1971) studied the problem of turbulent transition in vertical natural convection boundary layers at large Prandtl number. They found two distinct modes of instability. One of these modes arises from the coupling effect of the disturbance velocity and disturbance temperature via b u o y a n c y , and this mode is dominant and responsible for turbulent transition at large Prandtl number (Pr > 100). For a vertical isothermal surface, Hieber and Gebhart gave an approximate criteria for the onset of turbulence at large Prandtl number which in terms of a critical Rayleigh number becomes: Racritical = 7 Pr= × 104

(for Pr > 100)

(1)

51 Hieber and Gebhart noted that the coupled instability mode of eq. 1 occurred at smaller values of Rayleigh number than the uncoupled instability mode when Pr > 100. They further noted that at very large Prandtl number (Pr >~ 100) the boundary layer convective flow will always be turbulent because turbulent transition will occur before the full laminar boundary layer structure can develop. In other words, turbulent transition will occur in the thermal boundary layer region near the wall in Fig. 2c before the velocity profile outside this region can fully develop. There have been a number of theoretical studies of the effect of Prandtl number on turbulent transition (Elder, 1965b; Hieber and Gebhart, 1971: Gershuni et al., 1974; Golitsyn, 1979). These studies all show the trend for the critical Rayleigh number to increase with increasing Prandtl number. The variation with Prandtl number ranges from cases where the critical Rayleigh number is predicted to be proportional to the square root of Prandtl number to cases where the critical Rayleigh number is predicted to be proportional to the fourth power of Prandtl number (Elder, 1965b). The theoretical study by Hieber and Gebhart (1971} leading to eq. 1 appears to be the best current analysis. Experimental supporting data for any of these theoretical studies is lacking, however. Very little experimental data exists for transition to turbulent natural convection at large Prandtl numbers. Experiments by Elder (1965b) only considered fluids with Prandti number up to 7. The critical Rayleigh number for a basalt with a Prandtl number of 7.9 × 103 (Hardee, 1981) using eq. 1 becomes Racritica! = 4.5 × 10:'. This value is much greater than the usual value of Racritical = 109 for ordinary fluids with Prandtl number near unity. The Rayleigh number however, being a function of length cubed, tends to be large for typical magma bodies and can easily exceed the critical value in eq. 1. Turbulent natural convection in magma bodies is therefore a likely possibility. Although the critical Rayleigh number value for the onset of turbulence for both ordinary and large Prandtl number fluids was developed for vertical surfaces, the results do not differ much for the case o f the spherical magma body. For instance, Kreith (1970, 1973) notes that the transition Rayleigh number for the onset of turbulence for ordinary fluids in spherical cavities is also 109 just as it was for a vertical wall. The only difference is that the Rayleigh number is based on the diameter of the spherical cavity rather than the height of the vertical wall. For the case of a large Prandtl number fluid, such as magma, in a spherical cavity, eq. I should apply provided the Rayleigh number is based on the sphere diameter. Once turbulence occurs, the convective heat transfer to the wall increases because of the additional convective transport of the turbulence and the resulting steeper temperature gradient in the boundary layer region at the wall. In turbulent convective transport the usual Prandtl number is replaced by the turbulent Prandtl number P r t (Schlichting, 1960). The turbulent Prandtl number is the ratio of the eddy viscosity to the eddy diffusivity

52

and since these parameters are usually of the same order in turbulent convection (Schlichting, 1960; Kestin and Richardson, 1963), the turbulent Prandtl n u m b e r is close to unity regardless of the molecular values of viscosity and diffusivity. The assumption that Prt = 1 has been verified by experimental measurements (Ludwieg, 1956). The net result is that after fully developed turbulence occurs, the thermal and velocity boundary layers will be of equal thickness, as in Fig. 2a, even though prior to the onset of turbulence the velocity boundary layer in laminar convection was much thicker than the thermal boundary layer, as in Fig. 2b. CONVECTIVE UPMELTING

In a body with a simple form, such as the spherical one considered earlier, magma descends along the walls as it cools and warmer replacement magma rises at the center. The hottest magma is in contact with the wall at the top of the sphere while the coolest contact occurs at the bottom. The heat transfer to the wall varies around the sphere from top to bottom, being greater at the top of the sphere and less at the bottom. In some cases the heat transfer rate at the top of the sphere can be several times the heat transfer rate at the b o t t o m of the sphere. The variation in convective heat transfer rate around the cavity can be easily analyzed for the simple case of a cooling, spherical, convecting, magma body. If sufficient convection exists to produce boundary layer flow (Ra > 104) in the cavity, then the integral boundary layer equations can be used. The three-dimensional axisymmetric integral boundary layer equation for m o m e n t u m (Schlichting, 1960; Chapman, 1967} for a spherical cavity is: d --

~ 1 dr ~ ~ [du~ ~[ u2dy + - - ~[ u2dy = g~ sin O ; ~ d y - - v ~ - ~ y ] o

dx 0

r dx 0

b

(2)

and the corresponding integral boundary layer equation for energy is (Chapman, 1967): - d- f dx o

uC~dy+ -1- -dr r dx

/ o

uCdy = --a (d~_y)o

(3,

where x = RO is the distance along the surface in the direction of flow measured from the stagnation point and r = R sin 0 is the distance from a surface element to the axis of symmetry. The temperature (¢ = T -- Tf) and the velocity functions in the boundary layer are assumed to be of the form (Eckert and Drake, 1959):

53

This satisfies the b o u n d a r y conditions that the temperature is ¢w and the velocity is zero at the wall (y = 0) and the temperature is zero (4 = 0 or T = Tf) and the velocity is m a x i m u m (u = u~) at the edge of the thermal boundary layer (y = 5 ). For a large Prandtl n u m b e r fluid such as magma, the velocity reaches its m a x i m u m value at the edge of the thermal boundary layer (Kuiken, 1968) as noted earlier in Fig. 2. Substituting eqs. 4 and 5 into eqs. 2 and 3 and noting that the Prandtl number (Pr -- v/a) is large for magma leads to the result: d --

d~

cos 0 (u~5)

+ - -

sin 0

5 ul~

=

2

g~w

R

a

5 sin 0

v

u~

(6)

and d cos 0 dO (u15) + u15 sin 0

aR 15 -

(7)

6

Eliminating ul from eqs. 6 and 7 then gives: d6 --

d0

2 + -

3

cos 0 5

sin 0

R4 = 30

(8)

Ra 53 sin 0

where the Rayleigh n u m b e r is:

Ra -

g~gJw R3

(9)

p~

Solving eq. 8 for the thermal boundary layer thickness 5 then gives: If0° sinS/30d0] 1j4 5 _ (120)1z, Ra_V" R

(10)

sin 2j30

The local convective heat flux at the wall is:

q= hepw = - - K ( d_~y~) w

_ 2K~ w 5

(11)

The local convective heat transfer coefficient, h, can then be found by combining eqs. 10 and 11 to give:

hR K

sin 2'3 0

- 0.604 Ra I/4 0

If

for large Pr

(12)

d0i 1,4 sinS/30 J

0

This result is similar to the result found by Lin and Chao (1974) for the convective heat transfer coefficient on the outer surface of a sphere with natural convection in a large Prandtl n u m b e r fluid.

5.1

The average value of the convective heat transfer coefficient in the cavity is found by averaging eq. 12 over the spherical region or" hR_

4 /~

0.604 Ra-1

2zR 2(sin2a0lsin0 dO

J

K

o 4rrR ~

s i n s 30 dO

1:4

After integrating, this becomes:

BR K

-

0.604

R a -114

[{2/3) x/~ F(4/3) ] 34

?(11/6ij

or hR/K

=

0.595

(13)

Ra

This result compares well with data from Kreith (1970, 1973) for natural convection inside spheres. The local Nusselt number, defined as the ratio of the local convective heat transfer to the total convective heat transfer, is found by dividing eq. 12 by eq. 13 or: Nu =

sin: 30

1.015

(141

1

Figure 3, based on eq. 14, shows the variation of the local convective heat transfer around the wall of the spherical cavity. The average value of this Nusselt n u m b e r over the whole cavity is unity and is denoted by the dashed line in Fig. 3. The local convective heat transfer rate, denoted by the Z I

I

I

I

I

;

I 60

I 90

J 120

t 150

UJ

no

O

2.0

O.O 0

3

(TOP OF C A V I T Y )

180

( B O T T O M OF C A V I T Y )

POLAR ANGLE (degrees)

Fig. 3. Local Nusselt n u m b e r or d i m e n s i o n l e s s heat flux as a f u n c t i o n o f polar angle for b o u n d a r y layer f l o w inside a spherical cavity.

55 solid line in Fig. 3, is greater than the average value in the upper hemisphere (0 < 90 °) and less than the average value in the lower hemisphere (0 > 90°). In fact, the local convective heat transfer rate in the upper hemisphere is several times the rate near the b o t t o m of the spherical cavity (0 ~ 180°). This Nusselt number is inversely proportional to the thermal boundary layer thickness, 5, determined earlier in eq. 10. The b o t t o m of the sphere becomes a region of low heat transfer, partly because of the tendency of thermal boundary layers to thicken with length of flow down the wall ( R o h s e n o w and Choi, 1961), b u t largely because of the boundary layer flow at the b o t t o m of the sphere converges into a region of decreasing surface area causing the thermal boundary layer to thicken at a very rapid rate. The unusually thick boundary layer near the b o t t o m of the sphere results in very low heat transfer rates in this region as noted in Fig. 3. For these reasons, situations can occur where a spherical magma b o d y is melting wall rock at the upward boundary while magma is solidifying at the lower boundary (Ahem et al., 1981). Such p h e n o m e n a have been observed in simulated laboratory experiments (Hardee, 1974). A magma b o d y in this case could ascend a short distance through the crust b y this process of preferential melting and solidification. Unless the spherical magma b o d y has appreciable superheat, the upward movement of the b o d y by upmelting is limited to a distance of about one magma chamber diameter (Hardee, 1974). Calculations for cylindrical or diapiric magma bodies have also shown an upmelting tendency (Hardee and Larson, 1977). If a diapiric magma b o d y extended downward to a deep heat source or if heat is replenished by a new influx of fresh magma, then upmelting could be a significant factor in magma ascent. Such a magmatic b o d y would grow radially as well as in the upward direction by melting (Hardee and Larson, 1977). The shape of an upmelting magma b o d y varies considerably depending upon the source of heat. A non-replenished magma b o d y will tend toward an oblate spheroidal shape (Hardee, 1974). A replenished magma b o d y will grow in the shape of a vertical cylinder with vertical growth being five to ten times the rate of radial growth (Hardee and Larson, 1977). Magma chambers associated with caldera systems tend to be cylindrical in shape and those associated with small volume systems tend to be more conical (Spera and Crisp, 1981). VERTICAL ZONATION AND LAYERING Vertical zonation or layering in magma chambers is indirectly observed in singular or successive eruptions of compositionally zoned volcanic rocks and as fossilized evidence in solidified magma chambers (McBirney and Noyes, 1979; Irvine, 1979; McBirney, 1980). A prominent explanation for this observed zonation and layering is based on the convective process of double
56 explain observations of zonation and layering. Zonation in magma chambers is important, because it provides information for constructing models eft magma chamber formation and evolution (Spera and Crisp, 1981) as well as constructing models for understanding eruptive cycles. Double diffusive convection can result in a double boundary layer at the wall with one boundary layer flowing upward and one downward. This in turn can lead to stratification and layering in the upper part of a convecting enclosure. Double diffusive convection has been studied in several classical applications: thermohaline circulation in the ocean, sea water intrusion in lakes, solidification in metal castings, surface catalyzed chemical reactions, solution mining of salt cavities for crude-oil storage, melting of ice in sea water and melting of nuclear reactor containment vessels (Jensen, 1971: Turner, 1974; Nilson and Baer, 1981). Shaw (1974) analyzed the case of absorption of water from walls into convecting boundary layers of magma chambers. He concluded that chemical diffusion could produce a counterb u o y a n t boundary layer which would flow in a direction opposite to the thermal boundary layer and this could in turn cause the upper part of a solidifying magma chamber to become chemically zoned. More recently, McBirney (1980) and Turner (1980) have suggested that crystallization at the boundary of a convecting magma chamber can result in double-diffusive convection which will also lead to zonation in the upper part of a magma chamber. Nilson and Baer (1981) concluded that c o u n t e r b u o y a n t boundarylayer flow for large Prandtl number fluids like magma could occur in special cases. In particular, they showed that when the flow is dominated by the interior thermally driven boundary layer, a backflow region adjacent to the wall could take place at certain low values of relative buoyancy. Their analysis is somewhat clouded by their conclusion that the classical selfsimilar theory they used is not applicable in situations where strongly reversed counterflows are expected to occur. There are a number of other convective mechanisms, besides doublediffusive reversed-flow, that can explain observed layering and zonation in magma chambers. These mechanisms include reversed-flow in bicellular or multicellular circulation caused by extensive cooling at the top of a magma b o d y , periodic oscillatory convective circulation leading to reversed flow. and non-Newtonian b o u n d a r y layer effects. REVERSED FLOW AND MUTICELL CIRCULATION If the t o p of a spherical magma chamber is cooled extensively by contact with a two-phase geothermal zone or superconvecting zone (Dunn and Hardee, 1981), the negative b u o y a n c y due to cooling and crystallization will cause the fluid in the top of the chamber to descend vertically with replacement flow rising upward along the walls of the magma chamber. This reversed circulation can dominate the usual circulation caused by gradual cooling of the fluid due to heat loss at the side wall of the chamber.

57 Steady bicellular or multicellular flow, frequently observed in enclosed convection problems {Elder, 1965a; Bishop et al., 1966; Mack and Bishop, 1968; Powe et al., 1969, 1971; Szekely and Todd, 1971; Gartling, 1982), can also produce reversed flow at the wall. For instance, Bishop et al. (1966) observed that the flow at the cold outer wall of a spherical cavity was upward in the upper portion of the cavity. A perturbation analysis (Hardee, 1966, 1974; Mack and Hardee, 1968) can be used to illustrate a case of reversed convective flow. An idealized spherical magma chamber is assumed in which magma in the top of the chamber is cooled by contact with a near-surface convecting geothermal zone. Magma in the lower portion of the chamber is assumed to be heated by the injection of fresh magma. In the pseudo-steady case of a slowly convecting fluid, the cooling and heating conditions can be represented as volumetric internal heat generation (Pustovoit, 1958). The governing NavierStokes fluid-momentum equation and energy equation can be simplified using a technique described by Batchelor (1954) and others (Shaposhnikov, 1952; Pustovoit, 1958; Hardee, 1966, 1974; Mack and Bishop, 1968) in which the temperature and stream function are expanded in a Rayleigh number power series:

T = T o + R a TI+Ra2Tz + . . .

(15)

= Ra ~ I + Ra2~2 + • • •

(16)

Batchelor (1954) discussed the basis for this m e t h o d and noted that the m a x i m u m useful value for the Rayleigh number in this approach is around Ra = 1000. Assuming a spherical geometry with azimuthal symmetry, the dimensionless m o m e n t u m and energy equations for the zeroth-order ternperature and first-order stream function terms become (Hardee, 1974): -

-

OR 2

+---+ sin0 R OR R 2sin0 00

D4~1 = --

(

OTo

To--cos0 =0

~To )

R s i n : 0 ~OR + cos 0 sin 0 00 -

(17)

(18)

where the operator D 4 is: D4~ = D 2 (D2¢)

(19)

where D2~=

(02 1 O2 cot0 ~-~2 ~ R 2 ~02 -- - R 2

~) -~- ¢

(20)

The internal heat generation function, --cos 0, in eq. 17 is selected as a simple function which corresponds to cooling in the upper hemisphere and heating in the lower hemisphere of the spherical magma chamber. The chosen function, --cos 0, concentrates cooling and heating near the vertical

5~

centerline as it might be if cooling was due to c o n t a c t with a h y d r o t h e r m a l system at the top of a spherical magma cham ber and heating was due to injection o f magma at the b o t t o m of the spherical chamber. Many ot her variations o f cooling and heating can be generated by using related heat generation functions. For instance, the function --R cos 0 could be used t o p r o d u ce a horizontally stratified but vertically varying cooling and heating condition. Solution o f eq. 17 with the appropriate dimensionless b o u n d a r y condition, T i = 0 at R = 1, gives: 1 To = -~- (R E - - R ) c o s 0

(21)

This is substituted into eq. 18 to give: 1

D4@, -

R: sin: 0 cos 0

4

(22)

Solution o f this equation with t he appropriate dimensionless b o u n d a r y conditions, ~i = 3 @ i / O R = 0 at R = 1 and $i = 0 at R = 0, gives: , -

1 1152

(--2R 6 + 3R s ~ - R 3) sin 2 0 cos 0

(23)

The stream function from eq. 10 is then: Ra

-

1152

[Ra (R -- 1): (2R + 1)] sin: 0 cos 0

(24)

An examination o f eq. 24 easily reveals t hat the term in the brackets is always positive for all real values o f R, and t he cos 0 term is positive in the u pp er hemisphere (0 < ~/2) and negative in the lower hemisphere (0 > ~/2). The stream f u n ct i on ~ in eq. 24 is therefore negative in the u p p e r hemisphere and positive in the lower hemisphere indicating a bicellular circulation pattern with opposite r ot a t i on in the u p p e r and lower hemisphere as shown in Fig. 4. The direction o f r ot a t i on or circulation is made more obvious by determining the tangential velocity c o m p o n e n t v (Mack and Hardee, 1968): u-

(25)

R sin 0 aR F r o m eq. 24 : ~

0R

Ra -

384

[R3(R

--

1) (4R -- 1)] sin=0 c o s 0

(26)

and substituting this into eq. 25 gives: Ra v = -

384

[R(R

--

1) (4R -- 1)] sin 0 cos 0

(27)

59

2.5

J

O

1

Fig. 4. Bicellular convection pattern in a cooling spherical magma body showing reversed upflow at the wall in the upper hemisphere. The quantity in the brackets in eq. 27 is negative near the wall or boundary, where R is less than, but close to, unity. The function sin 0 cos 0 is positive in the upper hemisphere and negative in the lower hemisphere. The tangential velocity is therefore negative in the upper hemisphere or opposite to the normal clockwise direction of the angle 6 in the spherical coordinate system. This means t h a t the tangential velocity near the wall is counterclockwise or upward in the upper hemisphere. A similar argument shows that the tangential velocity near the wall is downward in the lower hemisphere. The circulation pattern in the upper hemisphere in this example rotates in a direction opposite to the usual expected sense. Examples of multicell circulation are observed in laboratory experiments and in numerical convection calculations. No precise theory is available to predict the onset of multicell circulation although the onset is generally found to be a function of Rayleigh number (Elder, 1965a). Multicell con-

60 vection is more likely to occur at large Rayleigh numbers. The larger the Rayleigh number, the more pronounced the multicell circulation becomes. Elder (1965a) explains that this dependence on Rayleigh n u m b e r occurs because the perturbation amplitudes, which lead to multicell circulation, increase with increasing Rayleigh number. The onset of multicell circulation is also a function of geometry, boundary and initial conditions and fluid properties (Liu et al., 1961; Elder, 1965a; Mack and Bishop, 1968). Powe et al. (1971) have speculated that certain multicell flow patterns are related to a Benard type instability which occurs when a cold boundary overlies a heated fluid region (Benard, 1900). This speculation is based on the observation that the Rayleigh number at the onset of multicell circulation is the same order as the Rayleigh number for the Benard instability in some c o m m o n enclosure geometries. The bicellular circulation pattern shown in Figure 4 developed when the heating in the lower hemisphere and cooling in the upper hemisphere were of equal magnitude. Use of a different temperature b o u n d a r y condition (Ti = 0 at R = 1) in eq. 21 would not change the convection result in eqs. 22--27. This bicellular circulation pattern results f~om the cooling/heating, cos 0 term in eq. 17 which corresponds to cooling in the upper part of the magma chamber by circulating ground water and heating in the lower part of the magma chamber by resupplied magma. The cos 0 term in eq. 17 leads to the presence of the first or R'- term in eq. 21 which in turn leads to the result on the right side of eq. 22 for the stream function. The temperature b o u n d a r y condition (Ti -- 0 at R = 1) used in eq. 21 only affects the second or R term in eq. 21 and this does n o t have an impact on the stream function eq. 22. Physically, this means that it is the cooling/heating ratio between the upper and lower hemispheres that produces the bicellular circulation pattern. If extensive cooling occurs at the top of the magma chamber and there is no resupplied heat at the b o t t o m of the chamber, then a single torodial circulation cell will develop similar to Figure 1. The direction of circulation will be reversed f~om the usual sense however and flow will be upward near the wall and downward along the vertical centerline of the chamber. This example can be analyzed similar to the previous example by changing the cos 0 term in eq. 17 to (1 + cos 0 ). This new term corresponds to extensive cooling by groundwater at the top of the chamber: 1+cos0 =2

at0 =0

and no cooling or heating at the b o t t o m of the chamber: 1 +cos0 =0

at0 =~

Repeating the solution procedure using this cooling term leads to the following result for tangential velocity: v =Ra

(3R ~ + 3 R - 1 ) + -

384

(4R - - 1 ) cos 0

( R - - 1 ) sin0

(28)

61 A study of this equation shows that the tangential velocity is upward near the wall. For instance, near the wall R < I b u t R ~ 1 and the tangential velocity is approximately: v ~

--Ra

[0.0119 + 0.00781 cos 0 ] (1 - - R ) sin 0

(29)

The tangential velocity is negative or upward even though the cos 0 changes sign between the upper and lower hemisphere. It can also be seen that the velocity is larger in the upper hemisphere near the region of intense cooling because the two terms in the bracket in eq. 29 are additative in the upper hemisphere where cos 0 is positive. Field examples of bicellular circulation and reversed circulation may have occurred at such places as the Skaergaard intrusion where the compositionally independent b o d y of the Upper Border Group might have developed because cooling was enhanced in the fractured basalt of the roof while cooling was restricted in the low permeability basement rocks (A.R. Mc Birney, pers. commun., 1982). The bicellular circulation pattern may also be related to bimodal eruptive patterns. OSCILLATING CONVECTIVE CIRCULATION PATTERNS

Another explanation for zonation or layering in magma bodies is reversed convective flow caused by transient or oscillatory convective circulation. Certain oscillatory circulation patterns can temporarily cause reversed-flow or up-flow near the top of spherical magma chambers. Cyclic circulation patterns in spherical enclosures were first observed experimentally by Bishop et al. (1966). Bishop et al. (1966) observed an unstable circulation pattern in a concentric spherical annulus enclosure. They referred to this unstable circulation pattern as "falling vortices flow". This flow pattern consisted of counter-rotating circulation cells which began forming in the upper part of the enclosure and spread downward, as illustrated in Fig. 5,

I

Fig. 5. Sequence of convective circulation in "falling vortices flow".

62

until the entire circulation pattern in the annulus had been disrupted. The original simple circulation cell then reformed and the process repeated in a periodic fashion. This instability was observed to begin at a gap Rayleigh number of R a * = 3600 where this Rayleigh number is based on a length equal to the width of the annulus. Liu et al. (1961) observed a similar p h e n o m e n o m for cylindrical annuli when the gap Rayleigh number reached 1600. This type of instability is probably related to the classical Benard cellular convection patterns which occur between horizontal flat plates when the gap Rayleigh number reaches 1700 (Benard, 1900; Rayleigh, 1916). It is also possible that unstable circulation patterns of this type could occur in a simple spherical cavity. Gartling (1982) has observed oscillatory reversed circulation in his finiteelement numerical solutions for natural convective circulation in cylindrical enclosures. Finite-difference numerical solutions by Larson (1982) for enclosed convecting magma bodies also showed an oscillatory circulation pattern which could lead to reversed flow. As an example, consider a cooling spherical magma b o d y that is convecting and losing heat to the wall. The hottest part of the spherical magma body is initially at the center of the sphere since the magma nearer the walls is cooled by contact with the walls. The hot magma at the center of the sphere is the most b u o y a n t magma, and it will rise to the top of the sphere and spread and circulate down the walls as indicated in the circulation pattern shown earlier in Fig. 1. After a while however, the magma in the non-circulating nodal region in Fig. 1 will in turn become the hottest magma in the sphere simply because it has n o t circulated and cooled like the other magma. Once the magma at the nodal region becomes the hottest and most b u o y a n t magma in the sphere, it will begin to rise and circulate around a newly created nodal region, which is shifted toward the wall as shown in Fig. 6. A similar process involving the movement of the maximum temperature node has been obNODAL POINT

a

b

c

Fig. 6. Transient convective circulation p a t t e r n or s t r e a m l i n e s e q u e n c e in a cooling spherical m a g m a b o d y .

63

served by Gartling (1982) for natural convection in cylindrical enclosures. This process will continue until the nodal region is so close to the wall that it cools rapidly by conductive heat transfer to the wall. The convection circulation pattern will then momentarily break up and reform into the original pattern of Fig. 6a. The process will then repeat at a rate governed by the circulation rate of magma around the sphere. When the circulation pattern changes from Fig. 6c to 6d, the circulation cell is momentarily disrupted and backflow at the wall might occur. This could result in material of a particular composition near the nodal point in Fig. 6c being swept up into the upper portion of the cavity. It is possible that these cyclic circulation patterns could result in some observable natural phenomena such as volcanic eruption cycles or changes in lava compositions during a volcanic eruption. Reversed flow accompanying the different circulation cycles could bring up magma of different composition from lower in the magma chamber to the top of the chamber where it could be erupted. This cyclic process could also release gases and lead to pressure oscillations which could result in cyclic eruption episodes. PROPERTY VARIATIONS

In analytical studies of natural convection, fluid properties are usually considered to be constant except for density, and density variations with temperature are considered only to the extent that such variations produce the b u o y a n c y term in the fluid m o m e n t u m equations. The density variation in the b u o y a n c y term is specified with an equation of state which is usually a Taylor series expansion in temperature with all terms neglected above the first order (Bird et al., 1966): ~P

P =Pr + ~'-TT r

( T - - Tr) + . . .

(30)

where Pr and Tr are the reference density and reference temperature. Frequently, this result is further simplified by replacing the density derivative with the coefficient of volumetric expansion, ~; ....

p

(31) P

to give the buoyancy term in the fluid momentum equations (Bird et al., 1966): P - - P r = p r f 3 g ( T - - Tr)

(32)

This approach works well with analytical solutions of simple convection problems when the reference temperature in eq. 32 is obvious. When the reference temperature is n o t obvious and when large temperature and density variations occur, analytical solutions based on the simplification

H4

of eq. 32 can lead to large errors [Gebhart, 19611. Sometimes these errors can be reduced by evaluating the fluid properties at an appropriate, weighted, reference temperature (Sparrow and Gregg, 1958; Gebhart, 1961). In other cases it may be necessary to resort to a numerical convection solution in which property variations, particularly density, are determined from an equation of state such as eq. 30 (Gartling, 1982). Property variations are extreme in the melt range. Density variations become large as crystals solidify or melt. If the crystals are small and entrained in the flow, the buoyancy term m the m o m e n t u m equations can be based on the density of the liquid/solid mixture. The solidification process also causes a large variation in viscosity with temperature. The specific heat changes significantly if this property is considered to account for the latent heat of solidification in the melt range (Carslaw and Jaeger, 1959). When lumped together in Lhe Rayleigh number, all these variations tend to partially compensate. Although dimensional analysis indicates that natural convection heat transfer is a function o f two dimensionless variables such as the Rayleigh number and the Prandtl number, Grober et al. (1969) noted that when the inertia terms are small the convection heat transfer, or Nusselt number, is a function of only the Rayleigh number:

[gp "'flCp.ST]

Nu = lira]

=f L

~K-

(33)

Here, p, f~, Cp and ~ are possibly _functions of temperature and ~ is considered to include the effect of entrained crystals. In the melt range, Cp and 13 become large in the numerator of the Rayleigh number and/x becomes large in the denominator. This tends to partially compensate the overall variation of Rayleigh number with temperature. In cases where solidification can not be treated as a macroscopic effect with entrained crystals, then the crystals and fluid must be treated as separate entities and a numerical solution approach is required. NON-NEWTONIAN EFFECTS

A Newtonian fluid is one which exhibits a direct proportionality between shear stress and shear rate in laminar flow (Skelland, 1967). Magma at a temperature well above the liquidus tends to behave like a Newtonian fluid (Hardee, 1981). In magmas at or below their liquidus, however, the presence of crystals in the liquid causes the magma to exhibit non-Newtonian behavior (Bird et al., 1966; Spera, 1980). In this case the magma may exhibit a finite shear stress at zero shear rate and the variation o f shear stress with shear rate is generally non-linear. As a first approximation, liquidus and subliquidus magmas are frequently considered to be Bingham fluids (Shaw et al., 1968). A Bingham fluid has a finite yield stress at zero-shear rate but otherwise the variation of shear stress with shear rate is linear. Once the initial yield stress is exceeded, a Bingham fluid behaves much like a New-

65 tonian fluid. Although the Bingham model is a reasonable assumption since m a n y subliquidus magmas do exhibit a finite initial yield stress (Murase and McBirney, 1973), this initial yield stress is n o t constant but varies with time (Murase and McBirney, 1973; McBirney and Noyes, 1979) and sheax rate history (Shaw et al., 1968; Hardee and Dunn, 1981). The initial yield stress of subliquidus magmas generally increases with time for a magma initially at rest (McBirney and Noyes, 1979; Hardee and Dunn, 1981} probably due to alignment and link-up of crystals in the melt. In convective transport analysis of a Bingham fluid, the initial convective start-up will be delayed until the thermal boundary layer thickness has grown by conduction to the point where buoyancy can overcome the initial yield stress (Hardee, 1981). Once convective circulation has started in an enclosed Bingham fluid, there will be two positions where the shear stress will be small or zero and will therefore be less than the Bingham yield stress. These positions are at the edge of the thermal boundary layer and in the core near the center of the enclosure where, as can be seen in Fig. 2b, c, the velocity gradient d r / d r goes to zero. Since the shear stress at these two positions is less than the Bingham yield stress, the fluid in these unsheared regions would move as a "solid plug" (Skelland, 1967). Data and analysis is lacking for natural convective circulation in Bingham fluids (Skelland, 1967), however, some speculations can be made concerning the expected flow. The solid plug, particularly the one at the edge of the thermal boundary layer, would certainly distort the circulation patterns t h a t would normally be expected. The solid plug at the edge of the thermal boundary layer would prevent fluid in the thermal boundary layer from circulating into the core after this boundary layer fluid had made a circuit past the wall of the enclosure. This could lead to a mechanism for stratification or layering and, if so, the layers would probably have a thickness on the order of the thermal boundary layer thickness. In a Bingham fluid the thermal boundary layer thickness at the initiation of flow is (Hardee, 1981): 5 = Oo/pgfJ

2

(34)

Substituting properties for a typical basaltic magma body (Hardee, 1981): p = 2700 kg/m 3 = 10 × 10 -s 1/°C T f - - T w = 200°C and the value of the Bingham yield stress given by Shaw et al. (1968} for molten basalt: o0 = 70 N/m 2 into eq. 34 gives: 5 = 26 cm

This layer thickness is in the range of field observations of layering (Irvine, 1979). This simple calculation merely suggests that Bingham fluid properties could lead to a mechanism for layering or stratification. A more detailed analysis of the cooling of a convecting Bingham fluid would be required to substantiate this layering mechanism. In order to avoid some of the analytical peculiarities of Bingham fluids, it is sometimes more convenient to represent magma as a power law fluid and use appropriate power law convection correlations (Yuen and Schubert, 1976; Hardee and Dunn, 1981). Subliquidus magmas can be fitted to power law fluid relations as well as to a Bingham relation (Hardee and Dunn, 1981). Power law fluids occur in many chemical engineering problems, and many convective transport studies have been made for these fluids (Metzner, 1956: Acrivos, 1960: Tien, 1967: Skelland, 1957) and it is possible to take advantage of these studies and apply them to the problem of convective transport in magma. SUMMARY

Convective transport in enclosures is a very useful subject for understanding the nature and evolution of magma bodies in the crust. Many analytical, numerical and experimental techniques have been used to study convective transport in magma. Magma is a peculiar fluid because of its large Prandtl number and its frequently non-Newtonian rheological character. The large Prandtl number of magma leads to a number of interesting effects. In boundary layer convective transport, the large Prandtl number causes the velocity boundary layer to be much thicker than the thermal boundary layer. In this situation, velocity variations may be spread throughout the magma chamber while thermal changes are restricted to regions near the wall. The Prandtl number has a minor effect on increasing heat transfer rates, and large values of Prandtl number tend to delay significantly the onset of turbulence and turbulent natural convection. For most convective transport calculations and heat transfer calculations, the Prandtl number of typical magmas can be assumed to be infinite and this simplifies the calculations. Convective transport in an enclosed magma chamber produces heat fluxes at the top of the chamber that are much greater than those at the b o t t o m of the chamber. This process can result in convective upmelting in which the top of the magma chamber melts as the b o t t o m solidifies. If the chamber is isolated from additional heat supply, the convective upmelting process will be restricted to magma migration distances on the order of a chamber diameter. If fracturing or stoping occurs at the top o f the chamber or if additional heat is supplied to the chamber, then greater upward magma migration distances are possible. Oscillatory convective circulation patterns sometimes occur in enclosed convecting fluid bodies and these oscillations are frequently periodic in

67

nature. Periodic convective oscillations could be related to observed phenomena such as eruptive cycles or changes in lava composition. Observations of vertical zonation and layering in magma chambers is sometimes attributed to double-diffusive boundary layer convection. Other simple convective phenomena can also explain observations of zonation and layering. Multicellular convective circulation may lead to reversed flow or upflow at the walls of an enclosed convecting b o d y . Upflow at the wails could bring wall material into the upper portion of a magma chamber where the material would then solidify into layers. Oscillatory convective circulation could also produce the same effect. The non-Newtonian properties of Bingham fluids can also lead to a possible mechanism for layering. NOMENCLATURE specific heat at c o n s t a n t pressure spherical o p e r a t o r function where D4¢ = D : ( D : o )

Cp

D4

and D2¢ =

[~

+ 2

f(x) g h K Nu Pr Pr t q r R Ra Ra* R a critical Racr, R a ~r

T Tf Tr Tw To T, T2 U U X

X y c~

r( 5

)

1 ~2 R 2 ~0 ~

cot0 R2

~1 V ~0

f u n c t i o n of x acceleration o f gravity convective heat transfer coefficient average convective heat transfer coefficient thermal c o n d u c t i v i t y Nusselt n u m b e r or dimensionless heat flux Prandtl n u m b e r t u r b u l e n t Prandtl n u m b e r convective heat flux distance f r o m surface element to the axis of s y m m e t r y radius in spherical c o o r d i n a t e system Rayleigh n u m b e r based on radius Rayleigh n u m b e r based on gap thickness Rayleigh n u m b e r at onset of turbulence Rayleigh n u m b e r at onset of convective instability temperature fluid or magma t e m p e r a t u r e reference t e m p e r a t u r e wall or b o u n d a r y t e m p e r a t u r e zero-order t e m p e r a t u r e p e r t u r b a t i o n term first-order t e m p e r a t u r e p e r t u r b a t i o n term second-order t e m p e r a t u r e p e r t u r b a t i o n t e r m velocity in direction o f b o u n d a r y layer flow velocity at edge of thermal b o u n d a r y layer flow tangential velocity c o m p o n e n t in 0 direction distance along surface in direction of flow characteristic length distance normal to surface or b o u n d a r y thermal diffusivity, K / p Cp coefficient of v o l u m e t r i c expansion gamma function thermal b o u n d a r y layer thickness

6~

~J

Pr O 0

~w

polar angle in spherical coordinate system viscosity kinematic viscosity, ~ ]p density reference density Bingham yield stress temperature difference, T - - Tf temperature difference, T w -- Tf stream function first-order stream function perturbation term second-order stream function perturbation term

REFERENCES Aerivos, A., 1960. A theoretical analysis of laminar natural convection heat transfer to non-Newtonian fluids. A.I.Ch.E.J., 6: 584--590. Ahem, J.L., Turcotte, D.L. and Oxburgh, E.R., 1981. On the upward migration of an intrusion. J. Geol., 89(4): 421--432. Anderson, O.L. and Grew, P.C., 1977. Stress corrosion theory of crack propagation with applications to geophysics. Rev. Geophys. Space Phys., 15 : 77--104. Batchelor, G.K., 1954. Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures. Q. Appl. Math., 12 : 209--233. Bartlett, R.W., 1969. Magma convection, temperature distribution, and differentiation. Am. J. Sci., 267: 1 0 6 7 - 1 0 8 2 . Benard, H., 1900. Tourbillons cellulaires dans une nappe liquide. Rev. Gen. Sci., 11: 1261--1271 and 1309--1328. Bhattacharji, S., 1974. Scale models for natural thermal convection in magmatic systems. Indian J. Earth Sci., 1(1): 73--83. Bhattacharji, S., 1979. Convective plumes, magma chamber development, penetrative magmatism and interplate rifting. Proc. Hawaii Symposium on Intraplate Volcanism and Submarine Volcanism, Hilo, Hawaii, July 16--22, p. 35. Bird, R.B. Stewart, W.E. and Lightfoot, E.N., 1966. Transport Phenomena. Wiley, New York, N.Y., pp. 10--14 and 299. Bishop, E.H., Mack, L.R. and Scanlan, J.A., 1966. Heat transfer by natural convection between concentric spheres. Int. J. Heat Mass Transfer, 9: 649--662. Brown, R.J., 1967. Natural Convection Heat Transfer Between Concentric Spheres. Ph.D. Dissertation, University of Texas, Austin, Texas. Carmichael, I.S.E., Nicholls, J., Spera, F.J., Wood, B.J. and Nelson, S.A., 1977. Hightemperature properties of silicate liquids: applications to the equilibration and ascent of basic magma. Philos. Trans. R. Soc. London, Ser. A, 286: 373--431. Carslaw, H.S. and Jaeger, J.C., 1959. Conduction of Heat in Solids, Oxford University Press, London, 289 pp. Chapman, A.J., 1967. Heat Transfer. The MacMillan Company, New York, N.Y., pp. 231--235,251--253 and 310--320. Chen, C.F. and Turner, J.S., 1980. Crystallization in a double-diffusive system. J. Geophys. Res., 85 : 2573--2593. Cheng, P., 1978. Heat transfer in geothermal systems. Advances in Heat Transfer, Vol. 14, Academic Press, New York, N.Y., pp. 1--100. Dunn, J.C., 1982, unpublished data. Also in Sandia/BES Geoscience Report, Second Quarter 1982, Sandia National Laboratories, Albuquerque, N.M. Dunn, J.C. and Hardee, H.C., 1981. Superconvecting geothermal zones. J. Volcanot. Geotherm. Res., 11: 189--201. Dunn, J.C., Carrigan, C.R. and Wemple, R.P., 1981. The influence of dissolved H20 on convective heat transfer in a magma. Trans. Am. Geophys. Union, 62: 1055.

69 Durst, F., Melling, A. and Whitelaw, J.H., 1976. Principles and Practice of Laser-Doppler Anemometry, Academic Press, London, 405 pp. Eckert, E.R.G. and Drake, R.M. Jr., 1959. Heat and Mass Transfer. McGraw-Hill, New York, N.Y., pp. 132 and 312--315. Elder, J.W., 1965a. Laminar free convection in a vertical slot. J. Fluid Mech., 23: 77--98. Elder, J.W., 1965b. Turbulent free convection in a vertical slot. J. Fluid Mech., 23: 99--111. Fedotov, S.A., 1975. Mechanisms of magma ascent and deep feeding channels of island arc volcanoes. Bull. Volcanol., 39(2): 1--14. Fedotov, S.A., 1976. On basic magmas ascent in the earth's crust and the mechanism of fissure basalt eruption. Proc. USSR Acad. Sci., 10: 5--23. Fedotov, S.A., 1979. On viscous heating of lava flows, diameter of asthenosphere magma columns, and velocities of magma ascent and magma differentiation beneath island arc volcanoes. Volcanol. Seismol., Acad. Sci. of USSR, 1: 1--15. Gartling, D.K., 1976. Finite element analysis of problems in convective heat transfer. Sandia Laboratories Report SAND75-0571. Gartling, D.K., 1982. A finite element analysis of volumetrically heated fluids in an axisymmetric enclosure. In: R.H. Gallagher, D. Norrie, J.T. Oden, O.C. Zienkiewicz (Editors), Finite Elements in Fluids. Wiley, New York, N.Y., pp. 253--267. Gebhart, B., 1961. Heat Transfer. McGraw-Hill, New York, N.Y., pp. 256--265. Gershuni, G.Z., Zhukhovitsky, E.M. and Yakimov, A.A., 1974. On stability of plane parallel convective motion due to internal heat sources. Int. J. Heat Mass Transfer, 17 : 717--726. Golitsyn, G.S., 1979. Simple theoretical and experimental study of convection with some geophysical applications and analogies. J. Fluid Mech., 95: 567--608. Grober, H., Erk, S. and Grigull, U., 1961. Fundamentals of Heat Transfer. McGraw-Hill, New York, N.Y., pp. 183 and 304. Hardee, H.C., 1966. Natural Convection Between Concentric Spheres at Low Rayleigh Numbers. Ph.D. Dissertation, University of Texas, Austin, Texas. Hardee, H.C., 1974. Natural convection in a spherical cavity with uniform internal heat generation. Sandia Laboratories Report SAND74-0089. Hardee, H.C., 1979. Heat transfer measurements in the 1977 Kilauea Lava Flow, Hawaii. J. Geophys. Res., 84(B13): 7485--7493. Hardee, H.C., 1980. Solidification in Kilauea Iki lava lake. J. Volcanol. Geotherm. Res., 7: 211--223. Hardee, H.C., 1981. Convective heat extraction from molten magma. J. Volcanol. Geotherm. Res., 10: 175--193. Hardee, H.C. and Dunn, J.C., 1981. Convective heat transfer in magmas near the liquidus. J. Volcanol. Geotherm. Res., 10: 195--207. Hardee, H.C. and Larson, D.W., 1977. The extraction of heat from magmas based on heat transfer mechanisms. J. Volcanol. Geotherm. Res., 2 : 1 1 3 - - 1 4 4 . Hieber, C.A. and Gebhart, B., 1971. Stability of vertical natural convection boundary layers: expansions at large Prandtl number. J. Fluid Mech., 49: 577--591. Irvine, T.N., 1979. Rocks whose composition is determined by crystal accumulation and sorting. In: H.S. Yoder Jr. (Editor), The Evolution of the Igneous Rocks. Princeton University Press, Princeton, N.J., pp. 245--299. Jensen, F.W., 1971. Total solution mechanism. Soc. Min. Eng. Trans. AIME, 250: 298-303. Kestin, J. and Richardson, P.D., 1963. Heat transfer across turbulent incompressible boundary layers. Int. J. Heat Mass Transfer, 6: 147--189. Kreith, F., 1970. Thermal design of high altitude balloons and instrument packages. J. Heat Transfer, 92 : 307--332. Kreith, F., 1973. Principles of Heat Transfer. Intex Educational Publishers, New York, N.Y., pp. 402--403.

70 Kuiken, H.K., 1968. An asymptotic solution for large Prandtl number free convection J. Eng. Math., 2(4): 355--371. Larson, D.W., 1982. Unpublished calculations. Sandia National Laboratories. Larson, D.W., Gartling, D.K. and Schimmel, W.P. Jr., 1977. Computational and experi~ mental methods for enclosed natural convection. Sandia Laboratories Report SAND770645. Lin, F.N. and Chao, B.T., 1974. Laminar Free Convection Over Two-Dimensional and Axisymmetric Bodies of Arbitrary Contour, ASME Transactions. J. Heat Transfer, 96: 435--442. Liu, C.Y., Mueller, W.K. and Landis, F., 1961. Natural convection heat transfer in long horizontal cylindrical annuli. International Developments in Heat Transfer, Part V, American Society of Mechanical Engineers, New York, N.Y., pp. 976--984. Lord Rayleigh (J.W. Strutt), 1916. On convection currents in a horizontal layer of fluid, when the higher temperature is on the underside. Philos. Mag., 32: 529--546. Ludwieg, H.,1956. Bestimmung des Verhaltnisses der Austauchkoeffizienten fiir Warme and Impuls bei turbulenten Grenzschichten. Z. Flugwiss., 4: 73--81. Mack, L.R. and Bishop, E.H., 1968. Natural convection between horizontal concentric cylinders for low Rayleigh numbers. Q. J. Mech. Appl. Math., 21(2): 223--241. Mack, L.R. and Hardee, H.C., 1968. Natural convection between concentric spheres at low Rayleigh numbers. Int. J. Heat Mass Transfer, 11: 387--396. Marsh, B.D., 1978. On the Cooling of Ascending Andesitic Magma, Philos. Trans. R. Soc. London, Ser. A, 288: 611--625. Marsh, B.D., 1979. Island-arc volcanism. Am. Sci., 67 : 161--172. Marsh, B.D. and Kantha, L.H., 1978. On the heat and mass transfer from an ascending magma. Earth Planet. Sci. Lett., 39 : 435--443. McBirney, A.R., 1959. Factors governing emplacement of volcanic necks. Am. J. Sci., 257: 431--448. MeBirney, A.R., 1980. Mixing and unmixing of magmas. J. Volcanol. Geotherm. Res., 7: 357--371. McBirney, A.R. and Noyes, R.M., 1979. Crystallization and layering of the Skaergaard intrusion. J. Petrol., 20(3): 487--554. Metzner, A.B., 1956. Non-Newtonian technology: fluid mechanics, mixing, and heat transfer. In: T.B. Drew and J.W. Hoopes Jr. (Editors), Advances in Chemical Engineering, Vol. 1. Academic Press, New York, N.Y., pp. 77--153. Murase, T. and McBirney, A.R., 1973. Properties of some common igneous rock and their melts at high temperatures. Geol. Soc. Am. Bull., 84: 3563--3592. Nilson, R.H. and Baer, M.R., 1981. Double-diffusive counterbuoyant boundary layer in laminar natural convection. Sandia National Laboratories Report, SAND80-2397. Norton, D. and Taylor, H.P., 1979. Quantitative simulation of the hydrothermal systems of crystallizing magmas on the basis of transport theory and oxygen isotope data: an analysis of the Skaergaard intrusion. J. Petrol., 20: 421--486. Oxburgh, E.R., 1980. Heat flow and magma genesis. In: R.B. Hargraves (Editor), Physics of Magmatic Processes. Princeton University Press, Princeton, N.J., pp. 161--199. Peck, D.L., 1978. Cooling and Vesiculation of Atae Lava Lake, Hawaii, U.S. Geol. Surv., Prof. Pap. 935-B. Pinkerton, H. and Sparks, R.S.J., 1978. Field measurements on the rheology of lava. Nature, 276: 383--385. Powe, R.E., Carley, C.T. and Bishop, E.H., 1969. Free convection flow patterns in cylindrical annuli. J. Heat Transfer, ASME, Ser. C, 91(3): 310--314. Powe, R.E., Carley, C.T. and Carruth, S.L., 1971. A numerical solution for natural convection in cylindrical annuli. J. Heat Transfer, 93: 210--220. Pustovoit, S.P., 1958. Transient thermal convection in a spherical cavity, J. Appl. Math. Mech., 22: 800--806. Rohsenow, W.M. and Choi, H., 1961. Heat, Mass and Momentum Transfer. Prentice-Hall, Englewood Cliffs, N.J., pp. 160--163,205.

71 Schimmel, W.P., 1976. An optical measurements laboratory for determining heat transfer coefficients. Sandia National Laboratories Report SAND76-0162, Albuquerque, N.M. Schliechting, H., 1960. Boundary Layer Theory. McGraw-Hill, New York, N.Y., pp. 256 and 492--499. Schneider, P.J., ]957. Conduction Heat Transfer. Addison-Wesley, Reading, Mass., pp. 230--233. Shaposhnikov, I., 1952. K, Theorii Slaboi Konvekstii (contribution to the theory of weak convection). Zh. Tekh. Fiz., 22: 826--828. Shaw, H.R., 1974. Diffusion of H20 in granitic liquids: Part I. Experimental data, Part II. Mass transfer in magma chambers, Geochemical Transport and Kinetics. Carnegie Inst, Washington Publ. 634, Washington, D.C., pp. 139--170. Shaw, H.R., 1980. The fracture mechanisms of magma transport from the mantle to the surface. In: R.B. Hargraves (Editor), Physics of Magmatic Processes. Princeton University Press, Princeton, N.J., pp. 201--264. Shaw, H.R., Wright, T.L., Peck, D.L. and Okamura, R., 1968. The viscosity of basaltic magma, an analysis of field measurements in Makaopuhi lava lake, Hawaii. Am. J. Sci., 266: 225--263. Skelland, A.H.P., 1967. Non-Newtonian Flow and Heat Transfer, Wiley, New York, N.Y., pp. 2--5, 73--102 and 413--421. Sparrow, E,M. and Gregg, J.L., 1958. The variable fluid-property problem in free convection. Trans. ASME, 80: 879--886. Spera, F.J., 1980a. Aspects of magma transport. In: R.B. Hargraves (Editor), Physics of Magmatic Processes. Princeton University Press, Princeton, N.J., pp. 265--323. Spera, F.J., 1980b. Thermal evolution of plutons: a parameterized approach. Science, 207: 299--301. Spera, F.J. and Crisp, J.A., 1981. Eruption volume periodicity, and caldera area: Relationships and inferences on development of compositional zonation in silicic magma chambers. J. Volcanol. Geotherm. Res., 11: 169--187. Spera, F.J., Yuen, D. and Kirschvink, S.J., 1981. A boundary layer approach for strongly dependent temperature viscosity. Trans. Am. Geophys. Union, 62: 1055. Spera, F.J., Yuen, D. and Kirschvink, S.J., 1982. Thermal boundary layer convection in silicic magma chambers: effect of temperature-dependent rheology and implications for thermogravitational chemical fractionation. J. Geophys. Res., 87(B10): 8755--8767. Szekely, J. and Todd, M.R., 1971. Natural convection in a rectangular cavity:Transient behavior and two-phase systems in laminar flow. Int. J. Heat Mass Transfer, 14: 4 6 7 482. Tien, C., 1967. Laminar natural convection heat transfer from vertical plate to power-law fluid. Appl. Sci. Res., 17 : 233--248. Turcotte, D.L. and Ahern, J.L., 1978. A porous flow model for magma migration in the asthenosphere. J. Geophys. Res., 83 (B2): 767--772. Turner, J.S., 1974. Double diffusive phenomena. Annu. Rev. Fluid Mech., 6: 37--56. Turner, J.S., ]980. A fluid dynamical model of differentiation and layering in magma chambers. Nature, 285: 213--215. Turner, J.S. and Gustafson, L.B., 198]. Fluid motions and compositional gradients produced by crystallization or melting at vertical boundaries. J. Volcanol. Geotherm. Res., 11: 93--125. Weertman, J., 1971. Theory of water-filled crevices in glaciers applied to vertical magma transport beneath ocean ridges. J. Geophys. Res., 76: 1171--1183. White, D.E. and Guffanti, M., 1979. Geothermal systems and their energy resources. Rev. Geophys. Space Phys., 17 : 887--902. Whitehead, J.A. Jr., and Luther, D.S., 1975. Dynamics of laboratory diapir and plume models. J. Geophys. Res., 80: 705--718. Williams, H. and McBirney, A.R., 1979. Volcanology. Freeman, Cooper and Company, San Francisco, Calif., p. 46.

72 Yuen, D.A. and Schubert, G., 1976. Mantle plumes: A boundary layer approach fo,Newtonian and non-Newtonian ternperature-dependent rheologies, d. Geophys. Res.. 81: 2499--2510. Yoder, H . S , 1976. Generation of Basaltic Magma, National Academy of Sciences. Washington, D.C,, pp, 1 --11.