Melt extraction from partially molten regions beneath mid-ocean ridges

Earth and Planetary Science Letters, 103 (1991) 69-78 Elsevier Science Publishers B.V., A m s t e r d a m

69

[vgl

Melt extraction from partially molten regions beneath mid-ocean ridges Li Yinting, R.O. M e i s s n e r , Fr. Theilen, X u e En Institut Jar Geophysik, Universit&'t Kiel, 2300 Kiel, FRG * Received August 22, 1990; revision accepted December 7, 1990

ABSTRACT A complete set of two-dimensional equations for conservation of mass, m o m e n t u m and energy of a two-phase system of melt in a deformable matrix is used to derive analytically a model for the vertical distribution of the melt fraction. Also, a set of asymptotic relations for parameters in the velocity field of the matrix and the melt in the neighbourhood of the symmetry axis of sea-floor spreading is obtained. These relations have been used to derive solutions appropriate to the melt extraction from the partially molten regions beneath mid-ocean ridges and the conditions for their existence. The results show the convergence of the melt streamlines towards the spreading center. It is caused by: (1) a horizontal gradient of the piezometric pressure associated with the deformations of the matrix, and (2) an upward increase of the porosity due to melting. For accepted values of the relevant parameters, a focusing effect of melt towards spreading centers exists, thereby providing an explanation for the mid-ocean " p a r a d o x " of a narrow region of ridge axis volcanism overlying a broad melt production zone.

1. Introduction Since McKenzie [1] and Scott and Stevenson [2] treated the partially molten region as a deformable porous medium, consisting of an unmelted viscous matrix saturated with a less dense melt, investigations concerning melt segregation and extraction have gained increasing attention. But in most of the publications about this topic the actual process of melting is not included in the models. Without an incorporation of melting into transport processes the equation of energy is not coupled with the equations for conservation of mass and momentum, and so is not included in the governing equations. The principal purpose of the present paper is to study the effect of melting on melt segregation. In considering the porous flow process for the melt and the viscous deformation process for the matrix we have to note that the melting of the parent rock may occur at the same time. In some cases which are of interest for geoscientists the

* Contribution no. 414 from Institut fi~r Geophysik, Kiel. 0012-821X/91/$03.50

© 1991 - Elsevier Science Publishers B.V.

rate of melting could be the rate-limiting factor in all processes of melt migration. As pointed out by Scott and Stevenson [3], this is true for the partially molten regions beneath spreading centers. In order to discuss the melt extraction beneath mid-ocean ridges, Spiegelman and McKenzie [4] have presented a simple corner flow model in which the porosity was assumed to be constant and melting was ignored. It has been concluded that matrix corner flow at ridges causes melt to be drawn to the ridge axis enabling the extraction of small melt fractions from a wide melting zone while showing a narrow zone of volcanism at the surface. F r o m their model, however, to pull melt out of a wide area requires an asthenospheric viscosity greater than 1 0 22 poise which is two orders of magnitude higher than presently thought. Due to this result it has been thought that the mid-ocean " p a r a d o x " of a narrow region of ridge axis volcanism overlying a broad melt production zone has not yet been explained by the corner flow model [5]. To explain this " p a r a d o x " , we have to solve the complete set of equations governing both melting and transport processes in the partially molten regions.

70

LI YINTING ET AL.

At first, a set of steady two-dimensional equations which include the conservation of energy is given in section 2. After reviewing previous results of two-dimensional upwelling flow and using some basic features of them, two assumptions are introduced, (1) that the temperature in the neighbourhood of the symmetry axis depends only on depth and (2) that the Peclet number is much greater than unity. As a result the vertical distribution of the melt fraction is found from the equation of energy in section 3. In section 4 an asymptotic analysis with second-order approximation is used to obtain a set of asymptotic relations among the parameters of the velocity field of the matrix and the melt. In section 5 we suggest a concept of "appropriate solutions" which could satisfy all constraints for testing models of melt migration beneath mid-ocean ridges. Here also the conditions for focusing melt towards the ridge axis are discussed. The physical aspect of results obtained in the present paper, such as the mechanism of melt extraction, the basis of the convergence of melt streamlines, the effects of the deformation of the matrix and the limitation of models with constant porosity suggested before are discussed in section 6.

where G is the density of the matrix. (c) Conservation of m o m e n t u m for the melt in X and Z directions:

2. Governing equations

0 --

Consider a partially molten layer of thickness L beneath the mid-ocean ridge, where the melt extraction processes can be modeled by a steady two-dimensional two-phase flow with melting, shown in Fig. 1. In this model, unmelted mantle material ascends at constant velocity Woo in the vicinity of the axis of symmetry until it reaches a depth Z = 0 where melting begins. In the melt zone between Z = 0 and Z = L the melt occupies a local volume fraction (porosity) ~ and moves with lateral velocity u and vertical velocity w relative to a fixed reference frame (X, Z}. The matrix flows with velocity (U, W), and its vertical velocity W is zero at Z = L. Then the governing equations can be written as follows: (a) Conservation of mass for the melt:

O ((Pu) + O (¢Pw) = F/Of

0X

z

Woo Fig. 1. Schematic diagram of a steady 2-D two-phase flow and

conceptual paths of melt (broken lines) and matrix (solid lines).

(b) Conservation of mass for the matrix: 3 3 X [ ( 1 - ~P)U] + -0-2 [(1 - (I)) W] =

-F/p s (2-2)

7Jf(I)2 (U -- U ) -- • 3___PP

K¢

O=-dgpfg-

(2-3)

3X

3P K¢ ( w - W)-c~3---Z

TJf (I)2

(2-4)

where ~f is the viscosity of melt, K¢ is the permeability, P is the pressure and g is the gravitational acceleration. (d) Conservation of m o m e n t u m for the matrix in X and Z directions:

3P T~f(I:)2 (U -- U ) - (1 - (I)) OX K¢

+n

+ 5-S] + (

= 0

+ - - o xoz (24)

0PZ TJf(I~2 ( W - - W ) - (1 - O ) O - (1 - ¢ ) p , g + ---~¢

(2-1)

where F is the mass transfer rate per unit volume from the matrix to melt and Pf is the density of the melt.

O2W + n -ff = 0

021411 + oz2 j +

(~+ 3}t 3 X O Z + ~ - ~ ] 1 (2-6)

MELT EXTRACTION FROM PARTIALLY MOLTEN REGIONS

71

where ~ is shear viscosity of the matrix and ~ is bulk viscosity of the matrix. In eqns. (2-3)-(2-6) the inertia terms are ignored because the Reynolds number is small for the slow viscous flow of both phases [3]. (e) Conservation of energy for the two-phase system: ~T

C p [ 10s(1 - ¢~) W + Of*w] OZ

+ cp], o~(1 - *

= K[ 02T

)u +

OT o f * . l . ~X

O2T

\-g2-5 + -~-Z ) - L Pr

(2-7)

[0s(1 - , ) w +

where T is temperature, Lp is the latent heat of melting. It is assumed that the melt and the matrix have the same specific heat at constant pressure Cp and thermal conductivity K. It is also assumed that viscous dissipation and adiabatic temperature gradient can be neglected [61. The permeability of the matrix K , can be expressed as a function of *: K. = D~"

(2-8)

where D = rZ/b, and ra is the grain size, b and n are dimensionless parameters. The value of n is between 2 and 3. According to Ribe [6] F can be F = [&(1 + [&(1

- * ) U + Of*u]

OF OF

- * ) W + Or*w] ~

by Li et al. [7] and Loper and Stacey [8] is that the lateral scale of vertical motion is narrower than that of the thermal anomaly. Their results show that in the vicinity of the symmetry axis the lateral distribution of temperature remains almost constant while the upward velocity has decreased by 10%. This feature leads to the assumption that the temperature is only dependent on Z in the vicinity of the symmetry axis while the upward velocity is also dependent on X. Due to T = T ( Z ) in the vicinity of X = 0, equation (2-7) becomes: dT

k d2T

or*w] ~-~ = o~ -d~

where k = K / ( C p ~ ) s ) is the thermal diffusivity. Another feature of the uprising plume is that its characteristic Peclet number Pe = LWoo/k is much greater than unity [9]. When Pe >> 1 the first term on the right-hand side of eqn. (3-1) can be neglected and then eqn. (3-1) becomes: [0s(1 - ¢ ) W +

dT

Or*w] d-Z

Lp + ~ppV = 0

where F is the mass fraction of matrix which has melted. To close the governing equations, a relationship between F and the other parameters is required and it will be given for a special case in section 3. The boundary conditions are discussed in section 5.

3. Vertical distribution of the porosity along the symmetry axis of an upwelling flow To provide some conditions for simplifying the complex set of equations given in section 2, we begin by reviewing the previous results about the thermal and mechanical structures of upwelling flows beneath mid-ocean ridges. A specific feature of uprising plumes presented

(3-2)

An experimental relation between the melt fraction F and temperature and pressure (depth) has been suggested by McKenzie and Bickle [9]. For small melt fractions, for example less than 30%, it can be approximately written as a linear relation:

F = A ( T - To) + BZ (2-9)

Lp

E r (3-1)

(3-3)

where A and B are constants and To = T(0) is the temperature at X = 0 and Z = 0. It is obvious that at the depth of the beginning of melting, Z = 0, the melt fraction F is zero. Substituting eqn. (3-3) into eqn. (2-9), we have dF r = [o~(1 - * ) w + oe,~w] dZ

= [Os(1-*)W+pf*wl(A dT +

(3-4)

Using eqn. (3-4) and eqn. (3-2), we have

[o.(1 - , ) w +

o,*w] - ~ + ~

=0 and then dT

Lp f Cp~

dT + B \ az

A-d~ + B

72

LI YINTING ETA[,.

and

dT dZ

It is obvious that D 3 is the melt fraction at 2~ = 1 ( Z = L). Equations (3-4) and (3-8) become

LpB Cp + LpA

(3-5)

F = J [p~(1 - * ) W + pfOw]

(3-10)

and

or

LpB Cp + LpA Z

T = TO

(3-6)

Substituting eqn. (3-6) into eqn. (3-3) we obtain

F = JZ

(3-7)

where the constant J B(1 + ALp/Cp) ~ is the slope of the vertical distribution of the melt fraction F and can be interpreted as representing the strength of melting. This quasi-onedimensional distribution of F could be used to derive an approximate vertical distribution of (I) in the vicinity of X = 0, after using the following relation between F and (I), given by Ribe [10], =

~' + ( 0 / ~ o ) ~ = F

(3-8)

where ~ff = ~lfWoo/(AogD). To derive this relation the Boussinesq approximation is applied and small porosity, for example < 20%, is assumed. In addition, it is required that melting occurs over vertical length scales exceeding about 100 m, a condition which is certainly satisfied beneath mid-ocean ridges. Let 2~ = Z / L and D 3 = JL, then

+ ((I)/(P0)" = D32 Figure 2 shows the vertical distributions of q) for different n, Woo and D3, when A 0 = O s - Of = 0.5 g / c m 3, ~ f = 1 0 poise and D = 1 0 5 c m 2 are used [6]. Values of D 3 are chosen to be 0.21 and 0.14, which express an average and a lower melt fraction at Z = L. An important fact is the upward increase of (P for every curve in Fig. 2 due to the upward increase of F.

4. Approximate solution in the vicinity of plume axis After (I)(Z) and T(Z) have section 3, eqns. (2-1)-(2-6) do rest of the complete set given on eqn. (3-10), eqns. (2-1) and

been determined in not couple with the in section 2. Based (2-2) become

(-11 where E = Pr/Ps, and

(1 - ~') ~-% ~ i + -~2 [(1 - ~ , ) w ] = - J [(1 - (I))W + e(I)w I

F = D32

(3-9)

(4-2)

Substituting eqns. (2-3) and (2-8) into eqn. (2-5), we have T~f (~91

[ 02U

oZu ]

0 5

%

D3=0"21~

3% /

t

( a2u + ax, a2w a---) = 0

+(f+~)lOX2

and substituting eqns. (2-4) and (2-8) into eqn. (2-6) yields

.// / ~.,'--%'~'L--------"~

- (1 -

////1O/o

//

'~o°~ ~"--- ~

(4-3)

1.

=2

e)[zxt,

g +

og, (r-

'qf(i)1-"(w- w) + n/5-2

{ 02W

+9

~ 2

0'5 10 Fig. 2. Vertical distribution of the porosity • for different n, Woo and D3. D3 is a dimensionless upward gradient of the melt fraction F.

[

-f- ~ - 1 -

n ~[

02U

~-]~--~---~"3t-

02W)

+ 5-U

O2W] aZ 2 ] =0

/

(4-4)

where p g a ( T - T~) is the buoyancy due to thermal density differences and p is the density of

MELTEXTRACTIONFROMPARTIALLYMOLTENREGIONS whole rock. T~o= T(X--* oo, Z ) is the ambient temperature at the same depth. In the present section a method of asymptotic analysis is used to obtain an approximate solution of eqns. (4-1)-(4-4) in the vicinity of the plume axis. This method has been succesfully used to calculate the structure of an upwelling flow in upper mantle [7] and in lower mantle [8], and has also been used to present a stagnation point flow model for melt extraction from a cylindrical mantle plume [11]. Beneath the ridge axis melting occurs within a broad zone that can extend up to 50-100 km to each side of the ridge axis [5]. In comparison with such a large lateral characteristic length a range of a few kilometers from the plume axis could be considered as the vicinity of the axis, in which we can write the Maclaurin expansions of W( X , Z ) , U( X , Z ) , w( X , Z ) and u( X , Z ) about X at X = 0, as follows:

q'( x'z):'Vo(Z) + ( -[-

.

.

o,I,

)oX +

vicinity of X = 0, then the following relations can be obtained: ~ d l + d ~ ( ~ b ° ) = 7J (1 - ~ ) a 0 + Jt~b o

(4-9)

da0 dO (1 - ~ ) c , + (1 - q~) --d-z - a ° d Z =-d(1-~)a

Cl-dl=2D

rh

(4-10)

o-J~qbb o

r/ d2Cl ]

* " - l [ ( ~ + fls) da2

(4-11)

d2a2 1 d3cl OZ--~ = ~ d Z ~3

(4-12)

a o - b0 nf

k

Ap dc I

1 ( O2xI*I g 2

47/t d2a0 \

aX2 )o

(4-13)

.

where ' t ' ( X , Z ) represents any value of W, w, U and u, and the subscript 0 represents the parameters at X = 0 . According to the character of sea-floor spreading, W and w must be even functions of X; U and u are odd functions of X. Thus, in the vicinity of X = 0, the approximations of W, w, U and u to the second-order are as follows:

W = a o + a2 X2

(4-5)

w = bo +

(4-6)

b2X 2

U = ClX

(4-7)

U = dlX

(4-8)

where 1(32W] a0=W0(Z)

73

a2 = ~ t 3X2 ]o

bo = wo( Z )

l(~2W I

From eqn. (4-12), we have

1 dc, + • ( Z -

a2- 2 dZ

(4-14)

where X is an unknown constant. The boundary condition a 2 = 0, d c l / d Z = 0 at Z = L has been used to obtain eqn. (4-14). The asymptotic relations eqns. (4-9)-(4-11) and (4-13)-(4-14) will be called a set of asymptotic equations which govern the velocity fields of the matrix and the melt in the vicinity of a symmetry axis.

5. Solutions of asymptotic equations and conditions for focussing melt towards mid-ocean ridge axes The best way to study the solution of eqns. (4-9)-(4-11) and (4-13)-(4-14) is to introduce the following dimensionless variables:

Z 2 =T'

Substituting eqns. (4-5)-(4-8) into eqns. (4-1)-(44) and noticing that the sum of the coefficients of X °, X and X 2 must disappear because eqns. (41) (4-4) must be satisfied at any point within the

L)

ao HCl H L dcl Y~ -- I4~oo ' Y2 =W~oo' Y-~- Woo d Z '

H 2a 2

Y4-

H 2L

b0

H

WoO , Y ° = ~ o o X, Y s - Woo , Y6=-~oodl

(5-1)

74

LI Y[NTING ET AI...

where H is a characteristic length scale for the horizontal change of the velocity field. After substituting eqn. (5-1) into eqns. (4-9)-(4-11) and (4-13)-(4-14), the following dimensionless groups of parameters are obtained:

Apg L2 D 1 - rlsWo~' D 2 -

1)( 1 - * ) ( 1

+ D22)

(5-3) From eqn. (4-11), we obtain ')(37Y'3 + ,sD s r ; o )

(5-4)

Substituting eqns. (5-1)-(5-4) into eqns. (4-9) and (4-10), a set of simultaneous first-order ordinary differential equations with three unknown functions Y1, Y2 and ~ can be given as 5 ' = A,Y~ - DsY 2 + A 3 Y 3 + A4 Y ° + A s

(5-5)

Y2' = Y3

(5-6)

Y3'= B1Y, + B, Y3 + & Y 4 ° + B5

(5-7)

where Y,' = d Y,/d Z, i = 1, 2, 3, and 1-*

D3 - D3e l - *

A3 = ~eDBD4D5 1 - *

(5-8) (5-9)

,n

1)

(5-10)

D2Z)

(5-11)

A4 = 2eD3D4Ds-f g ¢

( 2-

As

+

=

-eD,D3D4*"(1

o3 1o,]

eD3 * 2 1-*

n *' D3 2 * +~-

B4=-1-{(Z-1)(n~

-~'-

- E-dp

(5-13)

D3 - D 3 ~ 1 _ ~ ) (5-14)

,' +n~-(1-*)(1

+ D2Z') - * ' ( 1 + D2Z )

+o2(1 - *) - o,(1 - ,)(1

+ 022)]

3

(5-15) where * ' = d * / d Z'. Boundary conditions are: Y ~ = I , Y2=0 at 2 = 0 YI=0, Y3=0at 2=1

- ~D4Ds*~"-I~Y3 - 2 D 4 0 s * ~ " - 1 ~ Y ~ ' ( 2 - 1)

A1-

o

(5-2)

where "/1 = d T / d Z = - L p B ( C p + L p A ) 1 and Y2 is the geothermal gradient; ~,~ = ~"= ~/, which is a reasonable assumption as pointed out by Scott and Stevenson [3]. The dimensionless group of parameters D 1 in eqn. (5-2) characterizes the ratio of the first term and the last one in the brackets of the right-hand side of eqn. (4-13) and D 1 >> 1 as pointed out by Ribe and Smooke [11] and Morgan [5]. Accordingly, the last term can be neglected, and after using eqns. (5-1)-(5-2), eqn. (4-13) then becomes:

Y6= Y 2 - D 4 * ' "

[

1 A1 + - ~ - - - D 3

' D3 = J L ,

D~ Ds- L / ) 4 - r/fL 2' H

Y5 = Y] + DID4 *(n

14D4Ds,.

(5-12) B3-

P0l(Vl -- V2)

Ap

3 B1 -

(5-16)

The last condition Y3 = 0 at 2~ = 1 is derived from the assumption that the shear stress at Z = L is zero. For the three equations (5-5)-(5-7) there are four boundary conditions at two end-points and the surplus condition which is used simultaneously to determine the unknown constant yO in eqn. (5-5). Equations (5-5) (5-7) with boundary conditions (5-16) can be easily transformed into a set of algebraic simultaneous equations by a method of backward differences at m equidistant nodes. It was found that when m > 30, the solutions of this set of algebraic equations for different values of m remain almost the same with a relative error less than 10 3. Obviously, the solutions of the dimensionless equations (5-5) (5-7) depend on the dimensionless groups Di, i = 1 , . . . , 5 , and the dimensionless parameters n and e. Besides, the features of the solutions should be checked by geophysical observations. For the problem discussed in the present paper, the foXlowing conditions should be satisfied: (1) the streamlines of the matrix must be diverging to move the spreading oceanic plate away; (2) at the symmetry axis, the upward velocity of the matrix

75

MELT EXTRACTION FROM PARTIALLY MOLTEN REGIONS

should monotonously decrease from Woo at 2~ = 0 to zero at Z = 1 as Z increases from 0 to 1; (3) the horizontal gradient of the horizontal velocity of the matrix increases upwards; and (4) the maxim u m upward velocity of the matrix for every depth is at the symmetry axis and so a 2 < O. In addition, la21 should decrease upward so that the horizontal profiles of upward velocity of the matrix could become even smoother when approaching the boundary Z = 1. To discuss the extraction of melt from partially molten regions beneath mid-ocean ridges and to find the conditions for the convergence of the streamlines of the melt towards the symmetry axis another two constraints should be included in the present work. They are: (1) the Z-component of the velocity of the melt at X = 0 should be upward, i.e. w > 0 and its magnitude should be larger than that of the matrix at the same depth due to its smaller density and greater buoyancy; and (2) the X-component of the melt velocity in the adjacent regions of the axis should be directed towards the axis, i.e. d 1 < 0. Solutions of eqns. (5-5)-(5-7) with eqn. (5-16) by which the above-mentioned six conditions have been satisfied will be named the "solutions appropriate to the melt extraction beneath midocean ridges" in the present work, or, for convenience's sake, the "appropriate solutions". Obviously, the appropriate solutions can be obtained only for a set of special values of dimensionless groups D,, i = 1 . . . . . 5, and parameter n.

Y1

Y2

n:3 L = 6× 106cm

0.L

J :3.5 × l(J 8/cm "qs: 1- lO21poise 0.3

0.2

0.1

o

0,5

1.0

Fig. 4. Dimensionless horizontal gradient of X-component of the matrix velocity, as a function of Z for the same parameters as in Fig. 3.

Figures 3 - 7 show the vertical distributions of the dimensionless upward velocity of the matrix, Y1 = ao/Woo, and of the melt, 115= Do/Woo, the horizontal gradient of the X-component of the matrix velocity, 112 = Hcl/Woo, and of the melt velocity, Y6=Hdl/Woo, and the dimensionless parameter II4 = H2a2/Woo, respectively. It is very easy to see that the solutions given in Figs. 3 7 satisfy the above-mentioned conditions. They are appropriate solutions. The next step is to find the critical conditions for the existence of the appropriate solutions. A general discussion is difficult in mathematics. We have found some of them by numerical variations of parameters within limits comparable to natural

0 -

-

-1.0

-2.0 /

05

\

n:3

L = 6 × 106cm J - 3.5 × 113~/cm "qs: 1~ 1021poise

i

0.5

n=3

\

-3.0

\

/

L = 6 ×106c'm J =3.5 × 10-8/cm "qs= lx10

21

poise

-/,.0

2 1.0

Fig. 3. Dimensionless upward velocity of the matrix at X = 0, as a function of Z, the dimensionless distance from the depth of the beginning of melting, for n = 3, L = 40 kin, J = 3.5 × 10-8/cm, "% =1021 poise.

0.5 1.o Fig. 5. Dimensionless parameter relating to horizontal changes of upward velocity of the matrix, Y4 = H2a2/Woo, as a function of Z for n = 3 , L = 4 0 km, J = 3 . 5 × 1 0 - 8 / c m , ~/s=1021 poise.

76

LI YINTING ET AL.

Y5 I-

n=3 L=6xlO6cm J =3.5xlff8/cm

8.0

3.0 ~

~

6.0

~ L =4xlObcm

~qs= 1" 1022POise ~ls =6x10 21 TIs =3x10 21 2.0

3

0.5

x101,

1.0

Fig. 6. D i m e n s i o n l e s s u p w a r d velocity of the melt at X = 0, as a function of 2~ with different ~s for n = 3, J = 3.5 × 1 0 - 8 / c m , L = 40 km.

conditions. We do not believe that these conditions are complete but we think they are suitable for providing the conditions for focussing the melt towards the mid-ocean ridge. Before discussing these conditions, values or limits of the parameters presented in eqns. (5-5)-

Y6

-10

-10

2

-103 0

0.5

__ 2 1.0

Fig. 7. D i m e n s i o n l e s s h o r i z o n t a l g r a d i e n t of X - c o m p o n e n t of the melt velocity, as a function of 2~ with different ~l~ for the same p a r a m e t e r s as in Fig. 6.

J :3.5~ 108/cm

1:10°

1:102,

3:10 ,

lx 1022poise

Fig. 8. Critical values of n for the existence of " a p p r o p r i a t e solutions", as a function of ~t~ for L = 40 k m and 60 km, J = 3.5 × 10 ~ / c m . F o r the values of n a n d rl~ located in the h a t c h e d area no " a p p r o p r i a t e solutions'" can be found.

(5-15) should be selected. All of them, both dimensional and dimensionless, were divided into two groups. The first includes parameters which could be exactly estimated, such as Or = 2.8 g / c m 3, p~ = 3.3 g / c m 3, g = 103 c m / s 2, ~f = 10 poise, D = 10 5 c m 2, k = 1 0 2cm2/s, a=3×10 5 / K , Y2= 2 x 1 0 4 K / c m , C p = 0 . 2 4 c a l / g K, L p = 8 0 c a l / g , and H = 106 cm, which does not affect the results. The second group includes the rest of the parameters which could not be estimated exactly, but the limits of their changes are known, such as n between 2 and 3 for a partially molten m e d i u m of low porosity [3]; J is not obtained exactly due to the uncertainty of A and B, but an upper- and a lower-limit of J can be assumed as 3.5 x 10 S / c m and 2.0 X 10 8 / c m in order to make the m a x i m u m value of the melt fraction at the upper b o u n d a r y Z = L between 8% and 21% [5], ~s is 102°-1022 poise [11] and L is between 4 × 106 and 6 × 106 cm [11], and about Woo, the uprising velocity at X = 0, Z = 0, its m a x i m u m could be up to 3 x 10 6 c m / s [7]. Figures 8 - 1 0 show the conditions for the existence of the appropriate solutions. Two curves in Fig. 8 show critical values of n as a function of ~7~ for L = 40 km and L = 60 km when J = 3.5 x 1 0 - S / c m . W h e n n < nor, there are appropriate solutions. Parameter points { n, ~s } located in two regions with oblique lines in Fig. 8 correspond to n o n - a p p r o p r i a t e solutions. For example, if ~/~ > 1 × 10 20 poise for 2 < n _< 3 the appropriate solutions could be obtained, but if 7/s = 3 × 1019 poise,

MELT EXTRACTION FROM PARTIALLY MOLTEN REGIONS

3.0

i

""" "" b = 4 × 1 0 's cm

2.0

3~1019

lxlO 20

3x1020

lx1021

3x1021

1.102~poise

Fig. 9. Critical values of n for the existence of " a p p r o p r i a t e solutions", as a function of ~ for J = 2 . 0 × 1 0 - 8 / c m and 3 . 5 × 1 0 - 8 / c m and L = 4 0 km. For the values of n and ~ located in the hatched area no " a p p r o p r i a t e solutions" can be found.

n should be less than 2.6 for L = 40 km, and 2.7 for L = 60 km. Two curves in Fig. 9 show critical values of n as a function of ~s for J = 3.5 × 1 0 - S / c m and J = 2.0 × 10 S / c m when L = 40 km. Similarly, parameter points {n, ~ } from the two regions with oblique lines cannot give appropriate solutions. For J = 2.0 × 10 S / c m only if ~ > 3.55 × 102° poise the appropriate solution could be obwoo cm/s

"\

\ \

".

J:3.5×I() B n: 3

\

-6 10

\

"\ "\

\

"\

\ '

~-~

j:2.0x16gn=2

.....

\

\

'\

L

L : 6 x 10 c r ~

i

\

1G 7

10-8 1.x10 2°

L

I

,

3 ~ 1 0 2°

, , , , ,I 1 × 1 0 21

'

I

......

3 ~ 1 0 21

l x 1027poise

Fig. 10. Relations a m o n g Woo, L a n d ~s b y which " a p p r o p r i a t e solutions" c o u l d be o b t a i n e d for two l i m i t i n g conditions: (I) J = 3.5 × 10 S/cm, n = 3 (solid lines); (II) J = 2.0 x 10 a / c m , n = 2 (broken lines).

77

tained for any n, 2 < n < 3. However, for ~ls = 1 x 10 20 poise, n must be less than 2.7, and for ~s =- 3 × 10 ~9 poise, n must be < 2.3. Figure 10 shows the relations a m o n g Woo, L and ~ by which the appropriate solutions could be obtained for two limiting conditions: (I) Y = 3.5 x 1 0 - S / c m , n = 3 (in solid lines), and (II) J = 2.0 × 1 0 - S / c m , n -- 2 (in broken lines). All curves show that a lower ~s, the viscosity of the matrix requires a higher Woo, the velocity of the matrix at X = 0 and Z = 0, to obtain an appropriate solution. If the m a x i m u m of Woo is set to be at 3><10 6 c m / s , the m i n i m u m of ~ is 4 x 1 0 2 o poise for (II) and L = 60 km. To sum up, for accepted values of the main parameters concerned, i.e. 2 < n < 3, 2.0 × 1 0 - 8 , / c m _ 0 and u r < 0. It is obvious that melt-focussing occurs only when l ur I > U. If • is assumed to be constant, as done by Spiegelman and McKenzie [4], I Ur [ increases only as ~ increases. They found that only for 7l~ > 10 22 poise the melt-focussing can occur. Melting provides a continuous process and leads to an upward increase of the porosity of the

78

matrix, which can decrease the critical value of rIs for focussing the melt towards the center. I n the present paper this value of rls has been f o u n d to be as low as 10 2o poise. While the melt flows towards the axis of symmetry it accelerates upwards in response to the vertical gradient of the piezometric pressure which is the sum of c o n t r i b u t i o n s from the d e f o r m a t i o n s of the matrix a n d the b u o y a n c y of the melt. Both c o n t r i b u t i o n s , by focussing and accelerating the melt, can explain the m e c h a n i s m of melt extraction from partially m o l t e n regions a n d the formation of a m a g m a c h a m b e r at its top. It must be emphasized that the c o n t r i b u t i o n of the deformations of the matrix to the total piezometric pressure is d e p e n d e n t on depth a n d plays a less imp o r t a n t role w h e n a p p r o a c h i n g the u p p e r b o u n d a r y . But it can not be neglected for the whole region. Ribe a n d Smooke [11] c o n c l u d e d that the gradients of pressure associated with the d e f o r m a t i o n of the matrix are everywhere small relative to the b u o y a n c y Aog of the melt a n d the rate of the melt extraction from a p l u m e is not controlled by the viscosity of the matrix. Their conclusion was d r a w n by a first-order a p p r o x i m a tion in the lateral e x p a n s i o n of the matrix velocity. However, effects associated with the d e f o r m a t i o n of the matrix are expressed by the second-order partial derivatives of the velocity of the matrix a n d then the second-order terms in the e x p a n s i o n c a n n o t be neglected. The effects of the d e f o r m a t i o n of the matrix are most i m p o r t a n t for p r o d u c i n g the horizontal gradient of the piezometric pressure, b u t for the vertical gradient it competes with the b u o y a n c y of the melt. The d e f o r m a t i o n s of the matrix can be divided into two parts: (1) volume conserving shear deformation, a n d (2) volume changes of the matrix from c o m p a c t i o n or expansion. The latter part results from a n i n h o m o g e n e o u s d i s t r i b u t i o n of the porosity a n d a n o n - z e r o rate of mass transfer from the matrix to the melt, F. Models in which the porosity was assumed c o n s t a n t a n d melting was ignored, as suggested by Spiegelman a n d M c K e n z i e [4] have neglected the c o m p a c t i o n or expansion. However, in our results the effect on the convergence of the melt streamlines due to the volume changes of the matrix has been f o u n d to be greater than that due to the shear d e f o r m a t i o n .

LI Y I N T I N G E T AL.

7. Conclusions The results presented in this paper show that a matrix viscosity of 10 2o poise still allows meltfocussing. Thus, the n a r r o w zone of mid-ocean volcanism is explained by our model in a physical a n d m u c h more general way. A n o t h e r m a i n result that has come from this study is that the flow process b e n e a t h spreading centers is d o m i n a t e d by the rate of liquid supplied by the m e l t i n g p a r e n t rock, both processes occurring at the same time. The c o n d i t i o n s for the existence of the a p p r o p r i a t e solutions were easily fulfilled by i n t r o d u c i n g progressive melting into the model.

Acknowledgments We t h a n k Prof. L. Stegena a n d two u n k n o w n reviewers for their critical c o m m e n t s on this paper. The study was s u p p o r t e d by the Bundesminister i u m fiir F o r s c h u n g a n d Technologie ( B M F T ) and G E O M A R , Kiel.

References 1 D. McKenzie, The generation and compaction of partially molten rock, J. Petrol. 25, 713 765, 1984. 2 D.R. Scott and D.J. Stevenson, Magma solitons, Geophys. Res. Lett. 11, 1161 1164, 1984. 3 D.R. Scott and D.J. Stevenson, Magma ascent by porous flow, J. Geophys. Res. 91, 9283 9296, 1986. 4 M. Spiegelman and D. McKenzie, Simple 2-D models for melt extraction at mid-ocean ridges and island arcs, Earth Planet. Sci. Lett. 83, 137 152, 1987. 5 J.P. Morgan, Melt migration beneath mid-ocean spreading centers, Geophys. Res. Lett. 14, 1238-1241, 1987. 6 N.M. Rihe, Dynamical geochemistry of the Hawaiian plume, Earth Planet. Sci. Lett. 88, 37 46, 1988. 7 ki Yinting, R.O. Meissner and Xue En, The thermal and mechanical structure of a two-dimensional plume in the Earth's mantle, Phys. Earth Planet. Int. 33, 219-225, 1983. 8 D.E. Loper and F.D. Stacey, The dynamical and thermal structure of deep mantle plumes, Phys. Earth Planet. Int. 33, 304-317, 1983. 9 D. McKenzie and M.J. Bickle, The volume and composition of melt generated by extension of the lithosphere, J. Petrol. 29, 625 679, 1988. 10 N.M. Ribe, The generation and composition of partial melts in the earth's mantle, Earth Planet. Sci. Len. 73. 361 376, 1985. 11 N.M. Ribe and M. Smooke, A stagnation point flow model for melt extraction from a mantle plume, J. Geophys. Res. 92, 6437-6443, 1987.