Earth and Planetary Science Letters, 88 (1988) 3746 Elsevier Science Publishers B.V., Amsterdam  Printed in The Netherlands
37
[6]
Dynamical geochemistry of the Hawaiian plume Neil M. Ribe Department of Geology and Geophysics, Yale University, Box 6666, New Haven, CT 06511 (U.S.A.)
Received June 15, 1987; revised version received January 18, 1988 This paper presents a simple dynamical model for melting and trace element distribution in the Hawaiian mantle plume. I model the plume as a partially molten stagnation point flow against the oceanic lithosphere, and obtain solutions for the temperature, melt migration rate, and trace element concentration within it. Trace dement concentrations in the melt exceed simple batch melting predictions by up to 70%. The magnitude of this effect depends strongly on the solidmelt partition coefficient K. Trace elements with different K therefore experience a "dynamical fractionation" within the plume, and incompatible trace element ratios such as La/Ce always exceed the batch melting predictions. I suggest a simple model for plumelithosphere interaction in which melts from these two sources mix in proportions determined by thermodynamic constraints. The model can explain the composition of basalts from Haleakala if the degree of melting of the lithosphere F l decreases with time from roughly 10% for tholeiites to 2% for alkalic basalts. These values are considerably higher than previous estimates of F l < 1%, and imply correspondingly smaller and more realistic values (  10 kin) for the thickness of the melted part of the lithosphere. Partial melting of additional depleted sources such as the asthenosphere is therefore not required by the Haleakala data. Estimates of F~ are highly sensitive to the values chosen for the partition coefficients, however, and should therefore be interpreted with caution.
1. Introduction T h e H a w a i i a n I s l a n d s are b y far the m o s t s t u d i e d v o l c a n o e s in the world, a n d are the m o s t i m p o r t a n t e x a m p l e of oceanic i n t r a p l a t e volcanism. T h e H a w a i i a n I s l a n d s are g e n e r a l l y t h o u g h t to have b e e n p r o d u c e d b y m o t i o n o f the Pacific p l a t e over a n u p w e l l i n g o r " p l u m e " w h i c h rises f r o m a fixed source (a " h o t s p o t " ) d e e p in the m a n t l e [1,2]. T h e u p w e l l i n g m a t e r i a l u n d e r g o e s pressurerelease m e l t i n g d u r i n g its ascent, a n d the m e l t thus p r o d u c e d migrates u p w a r d to s u p p l y v o l c a n i s m at the surface. Because melt p r o d u c t i o n a n d m i g r a t i o n within the H a w a i i a n p l u m e c a n n o t b e o b s e r v e d directly, they are p o o r l y u n d e r s t o o d . R e c e n t geochemical d a t a f r o m H a w a i i a n volcanoes, however, p r o v i d e i m p o r t a n t new c o n s t r a i n t s on these processes. T h e a i m of this p a p e r is to use these d a t a to e v a l u a t e a simple p h y s i c a l m o d e l for m e l t i n g within the H a w a i i a n plume. N o m a t t e r where v o l c a n i s m occurs, a n a c c e p t a b l e melting m o d e l m u s t b e a b l e to a c c o u n t for 0012821X/88/$03.50
© 1988 Elsevier Science Publishers B.V.
the c h e m i s t r y of the lavas o b s e r v e d at the surface. I n the case of H a w a i i , however, this r e q u i r e m e n t is especially difficult to satisfy, b e c a u s e the c o m p o s i tion of the lavas e r u p t e d f r o m a n y given v o l c a n o changes w i t h time. M a c d o n a l d a n d K a t s u r a [3] d e s c r i b e d this t r e n d as a sequence of distinct stages, b e g i n n i n g w i t h a tholeiitic s h i e l d  b u i l d i n g stage which f o r m s a b o u t 99% of the volcano, followed b y p o s t  s h i e l d a n d p o s t  e r o s i o n a l stages d u r i n g w h i c h smaller a m o u n t s of alkali b a s a l t s are e r u p t e d . M o r e r e c e n t w o r k has s h o w n t h a t the trace e l e m e n t a n d i s o t o p i c c o m p o s i t i o n s of H a w a i i a n lavas also v a r y s y s t e m a t i c a l l y t h r o u g h time, T h e s e v a r i a t i o n s h a v e b e e n m o s t t h o r o u g h l y docum e n t e d b y C h e n a n d F r e y [4,5], w h o m e a s u r e d the trace e l e m e n t a n d i s o t o p i c c o m p o s i t i o n s of a large number of stratigraphically controlled samples f r o m H a l e a k a l a v o l c a n o o n E a s t Maui. Similar studies h a v e also b e e n c a r r i e d o u t on M a n n a Loa, M a n n a Kea, K o h a l a a n d East M o l o k a i v o l c a n o e s (see [5] for references). T h e d a t a of C h e n a n d F r e y [4,5] exhibit three
38
essential features: (a) incompatible trace element ratios (e.g., L a / C e ) increase with time; (b) 87Sr/86Sr decreases with time while 143Nd/]44Nd increases; and (c) the Sr and Nd isotopic ratios are inversely correlated with their respective p a r e n t / daughter ratios. Chen and Frey [4,5] interpreted their data as requiring mixing of melts derived from an enriched source with melts derived from a MORBlike source. They identified the enriched source as the mantle plume itself, and the MORB source as the lithosphere a n d / o r the asthenosphere. Chen and Frey calculated the trace element and isotopic composition of melts from these sources using simple batch melting relationships. They then determined the degree of melting of each source and the mixing proportions which provided the best fit to the data. The data are consistent with a variety of choices for these parameters, and the model is therefore not unique. All the successful models discussed by Chen and Frey [4,5], however, require that the degree of melting of the MORB source be low (0.21.0%) in order to explain the high incompatible element ratios of the alkali basalts. The mixing model of Chen and Frey [4,5] provides valuable insight into the processes responsible for Hawaiian volcanism. Like all simple mixing models, however, it ignores the mechanism by which mixing takes place. A real mixing process must not only satisfy a mass balance, but must conserve momentum and energy as well. It is not obvious that a simple mixing model can meet these requirements. A number of questions come to mind in this regard. Are batch melting relationships appropriate for estimating the trace element composition of melts produced by dynamic (nonstatic) melting? What determines the rate at which a mantle plume can supply melt to the base of the lithosphere? How much melt production occurs within the lithosphere itself in response to heating by the plume? In the following sections I investigate a simple physical model in an attempt to answer these questions. 2. The model
The model to be investigated is shown in Fig. 1. Unmelted mantle material ascends at velocity W0 and begins melting at a depth z =  d , where the
Tz =  Nb(TT I ) / d W= Uz = 0
z=0
r,u~U
W=Wo, U = Wo r / 2 d r=To, ~=0
z=d
Fig. 1. M o d e l for m e l t i n g in the H a w a i i a n plume. A s c e n d i n g m a n t l e m a t e r i a l begins m e l t i n g at a d e p t h z =  d where it crosses its solidus. D e p t h z = 0 is the base of the oceanic lithosphere. P a r t i a l l y m o l t e n region b e t w e e n z =  d a n d z = 0 is a d e f o r m a b l e p o r o u s m e d i u m or " m a t r i x " s a t u r a t e d w i t h less dense melt. Streamlines for the m a t r i x (solid lines) a n d the m e l t (dashed lines) are shown. T = temperature, q, = porosity, (U, W ) = (radial, vertical) m a t r i x velocity, (u, w ) = melt velocity, N B = Biot n u m b e r (see text).
solidus temperature is TO. The material between z =  d and the base of the lithosphere z = 0 is partially molten, and melting is occurring everywhere within this region. Following Sleep [6] and McKenzie [7], I treat the partially molten region as a deformable porous medium, consisting of an unmelted viscous "matrix" saturated with a less dense melt. The flow of the matrix is shown by the solid lines (streamlines) in Fig. 1. In calculating the matrix flow, I assume that the tangential stress at the base of the lithosphere is zero; similar results are obtained if this surface is rigid. The dashed lines are typical streamlines for the melt, which migrates through the matrix because of its buoyancy. Note that the melt streamlines intersect the base of the lithosphere, indicating that melt is extracted from the plume. The thermal boundary condition imposed at the base of the lithosphere contains a dimensionless parameter Nr~, the Biot number, which determines how much the lithosphere can be heated by the plume. If N a = 0, the lithosphere is a perfect insulator, and the heat flux into it is zero. If N B ~ ~ , however, the lithosphere behaves as a perfectly conducting thermal reservoir at constant temperature Tr The real situation, of course, is between these two extremes. The governing equations for the model are those of McKenzie [7], with a few significant differences
39 to be explained below. For steady flow, these equations are: v.
[(1  , ) u ]
=  r/p
(1)
= r/p
(2) (3)
u = u + kgAo 2
~f*
k = Ca 2,3
(4)
17s(V2 71)(UzWr)=g(Ap*+par)r
(5)
[(1*)U
+*u]'vT=
L r x v 2 r  p~p
(6)
F = p[(1  * ) U + ,u]" v F
(7)
F=A(T
(8)
To) + B ( z + d )
[K(IO)U+*u]Vc
=
1~) ( K  rc
(9)
P
The dependent variables in (1) to (9) are the matrix velocity U, the melt velocity u, the porosity (volume fraction of melt present) , , the temperature T, the rate of mass transfer from the matrix to the melt F, the fraction F of the rpatrix which is melted, and the trace element concentration c. Subscripts z and r denote partial differentiation. Equation (1) is a statement of conservation of mass for the matrix, and requires the net flux of matrix from a given volume (left side) to be equal to the rate at which matrix is destroyed within the volume by melting (right side). Equation (2) is a similar statement of conservation of mass for the melt. Equations (1) and (2) are obtained from equations (A4) and (A3) of McKenzie [7] by applying the Boussinesq approximation, in which the matrix density Ps and the melt density Pr are set equal to the same value p everywhere except in the buoyancy term g(Ps  Pe) = gAp [8]. Equation (3) is Darcy's Law, and states that melt migrates through the matrix at a rate proportional to its bouyancy. The quantity ~/f is the viscosity of the melt, k is the permeability of the matrix, and 2 is a vertical unit vector. Melt migration can also be driven by the nonhydrostatic pressure gradient associated with the flow of the matrix [6,7]. In a mantle plume, however, this pressure gradient is much smaller than the buoyancy gAp [8], and has therefore been omitted in (3). Equation (4) is a relation between the permeability and the poros
ity; C is a dimensionless constant and a is the average grain size of the matrix. Equation (5) is a vorticity equation which states that flow of the matrix is driven by buoyancy associated with radial gradients in porosity and temperature. The quantity 7/S is the viscosity of the matrix. Equation (6) expresses conservation of energy, and requires that advection of temperature gradients by the flow (left side) be balanced by thermal diffusion (first term, right side) and latent heat absorption due to melting (second term). Equation (6) is obtained from equation (A39) of McKenzie [7] by ignoring viscous dissipation and adiabatic decompression, and assuming that the melt and the matrix have the same thermal diffusivity K and heat capacity Cp. The quantity L is equal to TAS, where AS is the entropy change due to melting [7], and is assumed constant. Equation (7) states that the rate of melting is proportional to the rate at which material is being advected to a region where the intrinsic melt fraction F is greater. A derivation of the onedimensional analog of (7) is given by Ribe [9]. Equation (8) is a simplified equilibrium melting relation, which states that the melt fraction F is a linear function of temperature and pressure (depth). Finally, equation (9) expresses conservation of mass for a trace element in the melt, which requires that advection of trace element concentration gradients by the flow (left side) be balanced by transfer of trace element from the matrix to the melt (right side). Equation (9) is obtained from equation (A42) of McKenzie [7] by ignoring trace element diffusion and assuming that the matrixmelt partition coefficient K is a constant. The solution of equations (1) to (9) for the model shown in Fig. 1 is discussed in Appendix 1. This solution has the form of a "stagnation point flow" in which the radial velocities U and u are proportional to r, while all other dependent variables are functions of the depth z alone. I choose this form for the solution because of its simplicity, and because the streamlines match those at the top of an upwelling hot plume in buoyancydriven convection [10]. The solution of equations (1) to (9) proceeds sequentially: equations (1) to (8) are first solved simultaneously, and the results are then used to determine the trace element concentration c from (9). The solutions obtained in this way are discussed in the next two sections.
40
3. Trace element concentration in the plume
F
Equations (1) to (9) show that the concentration of a trace element in the plume depends in a complicated way on the flow and thermal structure. The best way to understand such a complex situation is to begin with a simple case. In this section, I assume that thermal interaction between the plume and the lithosphere can be neglected. A typical plume flow is thermally "fast", i.e., its characteristic Peclet number Pe = Wod/K is much greater than unity. This means that the lithosphere can only influence the thermal structure of the plume in a thin boundary layer at the plumelithosphere interface. The solution for this case is obtained from (A11) by ignoring the last term, which represents the thermal boundary layer at the base of the lithosphere. We then find that the temperature gradient within the plume is Constant, viz.: dT dz
LB 1 + Cp
(10)
This negative temperature gradient is due solely to absorption of latent heat during pressurerelease melting (recall that adiabatic decompression has been neglected). The melt fraction F varies linearly from zero at the bottom of the melting zone to a value Fp at the top, where Fp= Bd 1 + ~
(11)
The solution for the trace element concentration ¢(z) as a function of depth is obtained in Appendix 1. Fig. 2 shows the concentration Cp at the top of the plume z = 0 , as a function of the melt fraction Fp. The concentration is normalized by the concentration c b predicted for batch melting, viz.: cb =
CO K + Fp(1  K )
(12)
where c o is the concentration in the unmelted source. Fig. 2 thus shows how much the actual trace element concentration differs from the batch melting prediction. Results are shown for three values of the partition coefficient K, and for two values of a dimensionless parameter R defined by: R
Ca2gAp nfW0
(13)
1.6
t
,I /
///
/ I " 0.01
" t/ tall I / t~
/
l/
o l
." 0.05
........
0.10
0.15
.......
OI 0.20
0.25
F0 Fig. 2. Trace element concentration Cp at the top of the plume z = 0 as a function of the degree of partial melting Fp. Concentrations are normalized by the value c b predicted by the batch melting relation (12). Results are shown for three values of the partition coefficient K, and two values of the dimensionless buoyancy R = 1 0 3 (solid lines) and 10 5 (dashed lines). Quantity R, defined by (13), measures the importance of melt migration.
The parameter R is proportional to the buoyancy gAp of the melt, and measures the intensity of melt migration within the plume. A value of R = 0 corresponds to no melt migration, i.e., pure batch melting. In a mantle plume, however, buoyancydriven melt migration is important, and R 103105. The solid lines in Fig. 2 are for R = 103, and the dashed lines are for R = 105. The actual trace element concentration exceeds the batch melting prediction by up to 70%. The discrepancy is greater for more incompatible elements (smaller values of K ) and for faster melt migration (greater values of R). These results show that the batch melting relationship (12) does not provide good estimates of trace element concentrations in mantle plume melts. Fig. 2 refers only to absolute trace element concentrations. Melt source characteristics, however, are generally inferred from trace element concentration ratios. This procedure is based on the assumption that concentration ratios of highly incompatible elements can only be changed by very small degrees of partial melting. For the relatively high degrees of partial melting ( 0.050.3) that produce most common volcanic rocks, therefore, the incompatible element ratio in the rock should be similar to that in the source.
41
1.6
' / /' ~ '
'
'
'
1.4 c.! 5
c2
1.0 0.0
~
0.05
'
'
,
0.10
0.15
0.20
0.25
Fp
Fig. 3. Concentration ratio of two hypothetical trace elements at the top of the plume z = 0 as a function of the degree of melting Fp, for three values of the dimensionless buoyancy R. Value R = 0 corresponds to pure batch melting (no melt migration). Trace elements have partition coefficients K a = 0.001 and K 2 = 0.01.
This assumption, however, is based on very simple models of the melting process (batch or perfect fractional melting), and may not be valid in a more realistic dynamical melting situation. Fig. 3 shows the concentration ratio c l / c 2 of two hypothetical incompatible elements in the'melt exiting the plume as a function of the melt fraction Fp, obtained from the curves shown in Fig. 2. These elements have partition coefficients K 1 = 0.001 and K 2 = 0.01, and the ratio of their concentrations in the source is unity. The curve labeled R = 0 shows the concentration ratio predicted by the batch melting relation (12). The curves for R = 103 and R = 105 show the corresponding resuits when melt migration is allowed. These curves differ substantially from the batch melting prediction, indicating that dynamical melting fractionates incompatible trace elements more than batch melting does. This "dynamical fractionation" effect is a consequence of the threedimensionality of the flow, as discussed below in section 5.
4. Plumelithosphere interaction The results of the previous section were obtained by temporarily ignoring the fact that the hot plume can melt the oceanic lithosphere. In this section, I investigate how such melting affects the trace element and isotopic composition of plume
generated melts. The derivation comprises three steps: (1) determine the melt production rates qp and ql in the plume and the lithosphere from thermodynamic constraints; (2) determine the trace element concentrations Cp and e 1 in these melts; (3) estimate the isotopic and trace element composition of the erupted melt by mixing the melts from the plume and the lithosphere in proportion to their production rates. One source of energy for melting the lithosphere is the conductive heat flux from the plume, which can be calculated directly from the full solution (A11) for the temperature. The other important energy source is the excess heat content of the melt extracted from the plume, which is proportional to the difference between the temperature T(0) at the base of the lithosphere and the surface eruption temperature Te. The total energy flux (J m 2 s 1) available to melt the lithosphere is therefore: O=oCp(xTz(O)+qp[T(O)
 Tel }
(14)
where the quantity qp = RW0dp3(0) is the volume flux (m s 1) of melt from the plume. The heat flux (14) must supply both the energy required to raise cold lithospheric material to its solidus and the subsequent latent heat of melting. The rate at which lithospheric material melts is therefore: = ql
Q p [ C p ( T ~ o l  Tuth) + L]
(15)
where T~oI is the solidus temperature and THth is the temperature of the lithosphere before it is heated. Both of these temperatures depend on the depth at which the melting occurs. For simplicity, I use the values appropriate at the base of the lithosphere, i.e., TsoI = T o  B d / A and Tuth = T 1. Now that we know the melt fluxes qp and q~ from the plume and the lithosphere, the next step is to determine the associated trace element concentrations. The trace element concentration at the top of the plume, Cp, is determined from equations (A11) to (A14). The value of Cp obtained in this way depends primarily on the degree of melting F ( 0 )  Fp in the plume. Reasonable values for Fp are 0.050.1 for alkali basalts to 0.20.3 for tholeiitic basalts [11]. The trace element concentration c 1 of melts from the lithosphere is determined from the batch melting equa
42
tion (12) as a function of the average degree of lithospheric melting F 1. Note that F 1 is a free parameter which is independent of the thermodynamically determined melt production rate q~. This is not a very satisfactory procedure, but must suffice until more is known about the dynamics of melting and melt transport in the lithosphere. The final step is to determine the isotopic and trace element composition of the mixture (plume melt plus lithospheric melt). The isotopic ratio is given by: (16)
qpCprp + qlclrl qpCp + qlC l
,
,
I
,
,
i
008
i
i
I
• aTSr/aSS
r=
.7
0.04
O.OZ
"
0.7033 DL T 0 3 1
0.0 005
,
•
•
•
•
, 010
L
,
0.15
t 0.20
0.25
0.:50
%
where rp and r 1 are the isotopic ratios in the plume and the lithosphere, respectively, for the element in question. The concentration ratio of two trace elements 1 and 2 in the mixture is:
( qpCp +
O,lO
(17)
qlCl)l
Fig. 4. 87Sr/86Sr of the mixture (plume melt plus lithospheric melt) as a function of the degree of partial melting in the plume (Fp) and the lithosphere (F1). Solid symbols indicate typical values of 87Sr/86Sr for Honomanu tholeiite (squares), Kula alkalic basalt (triangles), and Hana alkalic basalt (circles) from Haleakala volcano (Table 2). Source compositions and partition coefficients assumed are given in Table 1.
(qpCp + qlCl)2
To obtain the solutions presented below, I have assumed the following parameter values: g = 9.8 m s 2, a = 103 m, C = 1 0 3, A p = 5 0 0 kg m 3, ~Tf= 1.0 Pa s, p = 3300 kg m 3, i¢ = 1.0 × 10 6 m 2 s l, Cp = 1000 J kg 1 ° C  1 , L = 3.3 × 10 5 J kg 1, A = 2 . 5 x 1 0 3 ° c  t , B = 9 . 7 × 1 0 6 m 1, T~= 1200 o C, Te = 1250 o C, and N B = 1.0. The depth d of the melt zone is a free parameter which is varied to obtain a range of values of Fp; smaller values of d yield smaller values of F_. For a given value of d, W0 = 3 x 1014 d m s  ~ and TO= T~ + B ( d + h)/A, where T~ = 1150°C is the solidus temperature at zero pressure and h = 6 x 10 4 m is thickness of the lithosphere. I show model results for three ratios: 8VSr/86Sr, 143Nd/144Nd, and L a / C e . The source concentra
tions and partition coefficients for Sr, Nd, La and Ce and the source values of the isotopic ratios are given in Table 1. The partition coefficients are weighted averages of the individual mineralmelt partition coefficients, as discussed by Chen and Frey [4,5]. Fig. 4 shows the 87Sr/86Sr ratio of the mixture as a function of the degrees of melting Fp and F 1 of the plume and the lithosphere, respectively. Figs. 5 and 6 are similar plots for a43Nd/144Nd and L a / C e . The solid symbols show typical values of these ratios for H o n o m a n u tholeiite (squares), Kula alkalic basalt (triangles), and H a n a alkalic basalt (circles) from Haleakala volcano [4]. These values are given in Table 2 for reference. Figs. 4  6 show that the composition of the mixture is primarily controlled by the degree of melting F 1 of the lithosphere. This is because F1
TABLE 1 Source compositions and partition coefficients [4,5]
TABLE 2
Element or ratio
Value in plume source
Value in lithosphere
Partition coefficient
Isotopic and trace element ratios of three typical Haleakala basalts [4]
Sr Nd La Ce
23.7 ppm 1.2 ppm 0.71 ppm 1.9 ppm
13.2 0.86 0.31 0.95
0.008 0.03 0.0025 0.005
Ratio
Honomanu tholeiite C122
Kula alkalic basalt H654
Hana alkalic basalt H6511
87Sr/SrSr 143N d / l a n Nd La/Ce
0.70387 0.51291 0.378
0.70334 0.51302 0.446
0.70314 0.51310 0.468
875r/86Sr 143Nd/a~Nd
0.7047 0.51265
ppm ppm ppm ppm
0.7023 0.51330
43 0.I0
.
.
.
.
.
,
,
0.08
0.08
0.06
F 0.040"06 ~
F i
0.04
~
~La
0 02
k
~
J
ICe=0 415
0021 o4°:
0.°,302
0475
0.0 0.05
0.10
0.15
0.20
0.25
0.30
Fp Fig. 5. Same as Fig. 4, but for ]43Nd/144Nd.
O.lO
0.08 / 006
F, 0.04
0 15
Fp
0.20
0.25
0.50
however, suggests a plausible explanation for the discrepancy: the large uncertainty in the values of the partition coefficients. The partition coefficients for La and Ce are known to be small (<< 1), but their magnitudes are poorly known (F.A. Frey, personal communication). The weighted averages for KLa and Kce given in Table 1 are therefore highly uncertain, and it is advisable to test the sensitivity of the model to changes in these values. Fig. 7 is identical to Fig. 6, except that the partition coefficients for La and Ce are assumed to be KLa = 0.001 and Kce = 0.01. These values increase the incompatibility of La relative to Ce by a factor of 5 over the values given in Table 1. The model results are now quite different: the values of F 1 predicted by Fig. 7 are much higher than those predicted by Fig. 6, and are similar to those required by the Sr and N d isotopic data (Figs. 4 and 5). A further prediction of Fig. 7 is that the degree of melting in the plume Fp must be about 2530% for the H o n o m a n u tholeiite, which is quite reasonable on petrological grounds [11]. 5. Discussion
Lo/Ce=0575
0.02
0.s95
0.41s~ • 0455_.I _ • _
0.05
0.10
Fig. 7. Same as Fig. 6, but for partition coefficients KLa = 0.001 and Kce = 0.01.
is generally smaller than Fp, and trace element concentration varies most rapidly as a function of F when F is small. Moreover, in each of Figs. 4  6 the required value of F 1 decreases from the Honomanu tholeiite to the Kula alkalic basalt to the Hana alkalic basalt. Because the alkalic basalts erupted later than the tholeiites, this result suggests that F 1 decreases with time. The L a / C e data, however, suggests much smaller values of F 1 than the Sr or Nd isotopic data. For example, a value of F 1 1  2 % is required to explain the L a / C e data for the Honomanu tholeiite, but the Sr and Nd isotopic data require F j  510%. A similar discrepancy is apparent for the Kula and Hana alkalic basalts. The model thus appears to be fundamentally in error. Further consideration,
o.o
0.0 0.05
,

,
010
,
,
0.15
,
,
0.20
,
",',
0.25
0.30
F0 Fig. 6. Same as Fig. 4, but for La/Ce. The two lines labelled 0.375 are part of the same contour.
The principal aim of this paper has been to combine geochemical and dynamical modeling to better understand the melting processes in the Hawaiian plume. The most robust result of this work is that trace element concentrations produced by realistic dynamical melting models differ substantially from simple batch melting predic
44 tions. This result may at first sight seem surprising, especially when one reflects that the model assumes perfect trace element equilibrium between the melt and the matrix. This behavior, however, may be easily understood by referring to the conservation law (9). The lefthand side of (9), which represents advection of trace element by the flow, is the product of an effective flow velocity and the concentration gradient. Because melt migration increases the effective vertical velocity, the vertical concentration gradient must be smaller (i.e., less negative) to maintain overall conservation of the trace element. Melt migration therefore causes the trace element concentration to be greater than it would be if no melt migration were occurring (batch melting). Batch melting relationships for trace element concentration are only valid if there is no melt migration (R = 0) or the flow is steady and unidirectional [9]. Richter [12] has obtained solutions for trace element concentration in a simple unsteady unidirectional flow, and has shown that batch melting predictions are not valid in this case. The results for the plumelithosphere interaction model are less satisfactory than those for the plume model alone, primarily because no dynamical model for melting in the lithosphere is available. The interaction model does lead to several important conclusions, however. The first is that the model predictions are highly sensitive to the values chosen for the partition coefficients. Estimates of Fp and F 1 obtained from geochemical modeling should therefore be regarded with skepticism. Nevertheless, Fig. 7 shows that reasonable values of the partition coefficients can be found which give a good match to both isotopic and trace element data. It thus appears that the model can account for the observed compositions of Hawaiian basalts if the degree of melting of the lithosphere decreases with time, from roughly 10% for tholeiites to about 2% for alkalic basalts. Partial melting of other depleted sources such as the asthenosphere is not required to explain the data. Moreover, the fact that the model embodies thermodynamic constraints shows that the posited mixing process is thermodynamically viable. Chen and Frey [4,5] have suggested that the compositions of Hawaiian basalts require much smaller degrees of lithospheric melting ( F 1 < 1%) than those predicted by this model. The difference
is due primarily to the fact that they used partition coefficients for La and Ce similar to those in Table 1. Values of F l < 1% are unlikely on dynamical grounds, because very small amounts of melt are difficult to extract from the source rock. Richter and McKenzie [13], for example, found that melt extraction from a source containing < 1% melt occurs far too slowly to account for the observed temporal isotopic evolution of Hawaiian lavas. More importantly, small values of F~ require melting of unrealistically large volumes of oceanic lithosphere. An estimate of the thickness h.1 of the melted lithosphere can be obtained from the formula h m = qlro/UoFl, where ql is the rate of melt production in the lithosphere, r0 is a characteristic horizontal scale for the plume, and U0 is the horizontal velocity of the lithosphere relative to the plume. The model discussed in this paper gives a typical value for ql of 3 × 10 11 m s ~. With r0 = 100 km, U0 = 3 × 1 0  9 m s 1 (10 cm yr 1) and F 1 0.01, we find h m  100 km, an impossibly large value. For F l = 0.1, however, h m  10 km, which is much more reasonable. The main limitation of the model presented here is that the temperature within the plume varies only with depth, so that the degree of melting Fp is laterally uniform. In reality, however, Fp (and hence the melting depth d) should decrease from a m a x i m u m value at the center of the plume to smaller values at the margins. In this paper, I have modelled the lateral variation of Fp and d by treating them as variable parameters of a laterally uniform model. This variation in Fp is then related to the temporal evolution of Hawaiian lavas by assuming that alkalic basalts correspond to smaller values of Fp than tholeiitic basalts. To properly understand the role of lateral temperature variations, however, a fully threedimensional model is required. Work on such a model is currently underway (Ribe and Daly, in preparation).
Acknowledgements I am grateful to F.A. Frey, J. Longhi and an anonymous referee for helpful comments. This research has been supported by National Science Foundation grant E A R 8618236, and by an award from the Atlantic Richfield Junior Faculty Program in Science and Engineering. Acknowledge
45 ment is also made to the Donors of The Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research.
subject to the conditions T (  d ) = TO and Tz(0)  NB[T(0)  T l l / d to yield:
SB T=T o
A(l+s)(z+d)
+ (1 + N B ) S B d  NB(1 + S ) A ( T o  T1) I ( z ) (1 + S ) A [ 1 + N.I(0)]
Appendix 1
A stagnation point flow solution of equations (1) to (9) is obtained by assuming that u = U = WorH(z)/2d and that all other dependent variables are functions of z alone. Equations (1) to (9) can then be written: Wo F [(1  q))W]z + d(1  i f ) H = 
(Al)
(CoW + RWoeO3)z + ~q~H = F__
(A2)
Hzzz = 0
(A3)
(A11) where:
I(z) =
f~/aexpl(1 S)Pe +
(A121 P
L
( W + RWo~3)Tz = rT~z  )~pF
(A4)
r= o(W + RWo~3)(AT~+ B) ([(h + K(1  cp)] W + RWo~ 3 } cz
(A5) (KDr p
c
and Pe = Wod/K. An equation for the porosity is obtained by using (A7) and (A8) to eliminate H and W from (A2), yielding: ~ =
a(Ar +B) G + R~2(3  4q))
(A13)
Finally, the solution to (A6) is:
c0 {  ( l  K )
c=~exp
(A6) The quantity R is defined by (13), and a subscripted z denotes differentiation. Equation (A3) can be integrated subject to the boundary conditions Hz(0 ) = H (  d)  1 = 0 to yield:
H(z)=I+D
~1
(A7)
×fz
W = Wo[G(z)  R e 3]
(A8)
where: 2
G(z) = ~D 
(1  D ) z 2
D z3 3 d3
(A9)
U p o n using (A8) to rewrite (A4), we obtain:
rTzz=WoG[(l+ S)Tz + ~B]
(Ae0)
where S = LA/Cp. Equation (A10) can be solved
}
(A14) The unknown constant D which appears in (A9) is determined by requiring W(0)0, which by using (A8) can be written: R e ( 0 ) = ]D
where D is an undetermined constant. We now add (Al) and (A2), and integrate the result subjeet to the conditions e p (  d ) = W (  d ) Wo = 0 to obtain:
G(AT~ + B)
_ d G + (1K)O~7(Rq,3_ G) dz
(A15)
Condition (A15) can be satisfied using a simple iterative procedure, as follows: (1) guess D; (2) determine T(z) from (A11); (3) determine q)(z) by solving (A13), using a highorder numerical scheme such as the RungeKutta algorithm; (4) repeat this procedure until (A15) is satisfied, using Newton's method to assure convergence. Solutions obtained in this way converge to within a fraction of a percent in 510 iterations. References
1 J.T. Wilson,A possible origin of the Hawaiian islands, Can. J. Phys. 41, 863870, 1963. 2 W.J. Morgan, Plate motions and deep mantle convection, Geol. Soc. Am., Mem. 132, 722, 1972.
46 3 G.A. Macdonald and T. Katsura, Chemical composition of Hawaiian lavas, J. Petrol. 5, 82133, 1964. 4 C.Y. Chen and F.A. Frey, Origin of Hawaiian tholeiite and alkalic basalt, Nature 302, 785789, 1983. 5 C.Y. Chen and F.A. Frey, Trace element and isotopic geochemistry of lavas from Haleakala volcano, East Maui, Hawaii: implications for the origin of Hawaiian basalts, J. Geophys. Res. 90, 87438768, 1985. 6 N.H. Sleep, Segregation of magmas from a mostly crystalline mush, Geol. Soc. Am. Bull. 85, 12251232, 1974. 7 D.P. McKenzie, The generation and compaction of partially molten rock, J. Petrol. 25, 713765, 1984. 8 N.M. Ribe and M.D. Smooke, A stagnation point flow model for melt extraction from a mantle plume, J. Geophys. Res. 92, 64376443, 1987.
9 N.M. Ribe, The generation and composition of partial melts in the earth's mantle, Earth Planet. Sci. Lett. 73, 361376, 1985. 10 G.O. Roberts, Fast viscous Benard convection, Geophys. Astrophys. Fluid Dyn. 12, 235272, 1979. 11 B.O. Mysen and I. Kushiro, Compositional variations of coexisting phases with degree of melting of peridotite in the upper mantle, Am. Mineral. 62, 843856, 1977. 12 F.M. Richter, Simple models for trace element fractionation during melt segregation, Earth Planet. Sci. Lett. 77, 333344, 1986. 13 F.M. Richter and D.P. McKenzie, Dynamical models for melt segregation from a deformable matrix, J. Geol. 93, 729740, 1984.