- Email: [email protected]
The dynamics of plume-ridge interaction, 1: Ridge-centered plumes
EPSL ELSEVIER
Earth and Planetary Science Letters 134 (1995) 155-168
The dynamics of plume-ridge interaction, 1: Ridge-centered plumes N.M. Ribe 1, U.R. Christensen, J. Theiging Institut ffir Geophysik, Universiti~t Gfittingen, Herzberger Landstrasse 180, 37075 G6ttingen, Germany
Received 12 January 1995; revised 25 May 1995; accepted 6 June 1995
Abstract
We investigate the dynamics of mantle plumes rising beneath mid-ocean ridges. To clarify the physics of this process, we examine first a simple 'lubrication theory' model in which a point source (analogous to a plume 'stem') of volume flux Q located directly beneath the ridge releases buoyant fluid into a viscous comer flow driven by a velocity boundary condition u(x) = U tanh(x/d), where U is the half-spreading rate and the 'gap width' between the diverging plates is ~ 5d. Numerical solutions of the differential equation governing the plume head thickness S ( x , y ) show how the width W of the plume head along the ridge depends on Q,U,d, and o- =- gAp~48,1, where Ap is the density deficit of the plume and ~/ is the viscosity. In the geophysically relevant 'narrow gap' limit (0//0-) 1/4 >> d, W ~ (Q//U)I/2IIb 0"053, where Mb = Q0-//U 2 is the 'buoyancy number'. Numerical solutions of a more realistic 3D convection model with strongly temperature and pressure-dependent viscosity obey a nearly identical scaling law, and show no evidence that W is increased by 'upslope' flow of plume material toward the ridge along the sloping base of the rheological lithosphere. To apply our model to Iceland, we incorporate into it a melting parameterization that allows prediction of the excess crustal thickness produced by melting in the plume head. This extended model shows that the observed depth anomalies along the Mid-Atlantic Ridge near Iceland cannot be explained by a hot (temperature contrast AT ~ 250°C) and narrow (radius ~ 60 km) ridge-centered plume. Instead, the anomalies are consistent with a much cooler and broader upwelling. 1. Introduction
Many of the earth's hotspots lie on or very near mid-ocean ridges. Examples include the Iceland, Azores, Afar and Galapagos hotspots. The strong correlation of hotspot locations with ocean ridges was first noted by Morgan [1]. At about the same time, extensive surveying and sampling of the ocean floor revealed that the association of hotspots with ridges is not merely geometrical, but physical and
1 Permanent address: Department of Geology and Geophysics, Yale University, POB 208109, New Haven, CT 06520-8109, USA Elsevier Science B.V. SSD! 0012-821X(95)00116-6
chemical as well. Near Iceland, for example, the seafloor topography and the L a / S m ratio of dredged basalts are anomalously high for several hundred kilometers along the Mid-Atlantic Ridge [2,3], and similar anomalies exist near other hotspots located on ridges [4]. These anomalies are thought to be due to the lateral spreading of hot plume material beneath the ridge [3]. The present paper is a study of the dynamics of plume-ridge interaction for the special case of a ' ridge-centered' plume rising directly beneath a ridge. Our principal goal is to determine how far buoyant plume material can spread along the strike of a ridge before being swept 'downstream' by the diverging
156
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
plates. Our study consists of three parts. We first examine a simple 'lubrication theory' model that reveals the fundamental dynamics of plume-ridge interaction, and we compare the predictions of this model with the scaling laws of Feighner and Richards [5] and Feighner et al. [6]. We then study a more realistic 3D numerical convection model with temperature and pressure-dependent rheology. Finally, we apply the convection model to the Iceland hotspot, the classic example of a ridge-centered plume.
plex argument [9]. The quantity UzF(x,z) is simply the streamfunction for the corner flow in the halfspace z > 0. Eq. 1 is a nonlinear advection/diffusion equation describing a balance between injection of plume fluid, its advection by the corner flow, and its buoyancy-driven ' self-spreading'. We first use Eq. 1 to determine the length scales for the thickness S and width (normal to the spreading direction) L(x) of the plume head. By balancing the self-spreading and source terms we obtain o-S4/L 2 ~
2. Lubrication theory model
S~
In this section we make use of the theory of flow in thin layers to formulate a simple model for plume-ridge interaction. The idea behind the model is that buoyant plume material ascending from great depth through a narrow conduit will spread laterally when it reaches the base of the diverging plates, forming a layer that is thin compared to its lateral extent. Similar models have already been employed successfully to describe the geometry of swells produced by plumes far from ridges [7,8]. Our model envisions a passive 'corner flow' driven in a halfspace z > 0 of fluid with constant viscosity 7/ and density p by a lateral velocity u(x) = U tanh ( x / d ) imposed at the surface (z = 0). Here U is the half-spreading rate. Moreover, because u(__+2.5d) = _0.99U, the 'gap width' between the diverging plates is ~ 5d. Within the corner flow, a point source (analogous to the top of the plume 'stem') of volume flux Q located at ( x , y ) = (0,0) (directly beneath the ridge) releases buoyant fluid with viscosity 77 and density p - zip, which spreads along the upper surface of the halfspace as a layer (the plume 'head') of thickness S(x,y). The differential equation governing S can be derived using lubrication theory (Appendix), and is
U(FS) x =
0"~7254
-~ Q 6 ( X ) 6 ( y )
2 (1 _ S_ /x) F ( x , S ) = -~ Im ~O -~ + ~ d + ~--d
(2)
where 'Ira' denotes the imaginary part and O(p)--( d / d p ) l n F ( p ) is the digamma function of a c o m -
- So
(3)
The scale for L can now be determined by balancing the advection and source terms in Eq. 1. The result depends on one's distance from the ridge. Far from the ridge (x >> d), F = 1, which implies U(FS) x US/L. Balancing this with the source term (order Q/L2), we obtain
L
~ Q3/4U-lo'l/4
=L o
(4)
This lengthscale is identical to that derived for a plume beneath a surface moving at uniform velocity U [81. At the ridge (x = 0), however, the 'waist width' [5] L(0) - W is controlled by two additional lengthscales. Because S x = 0 at the ridge by symmetry, the advection term there is
USFx(O,S ) = - ~ b
~ 7rd +
(5)
This expression has two distinct limiting forms, depending on the ratio of the characteristic layer thickness S O to the gap width d. In the 'narrow gap' limit So/d >> 1, USFx(O,S) = 2U/Tr, which when balanced with the source term (order Q / W 2) yields
Q )1/2
(1)
where o-= gAp/48rl, Vn is the horizontal Laplacian, and
a l L 2, or
W~(~
-Ls
(6)
In the 'wide gap' limit So/d << 1, by contrast, USFx(O,S) = US~d, which when balanced with the source term yields W ~ Q3/8U-1/2dl/2or
1/8 _ Lw
(7)
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
We now examine the solutions of the full lubrication equation 1. To avoid a singularity, we replace the point source by an equivalent Gaussian source distribution with a small but finite radius a << W, viz., Q6(x)6( y) --, (Q/Tra2)exp[-(x z + y 2 ) / a 2 ] . Nondimensionalization of Eq. 1 shows that the solutions depend only on the dimensionless groups Sold and Lo/d. Fig. 1 shows the solution for So/d--- 1 and Lo/d = 10, obtained using a finite-difference technique with explicit time stepping [7] and a grid spacing Ax = Ay = 0.02L 0. The increasing width of the plume head away from the ridge, which is typical of the numerical solutions, reflects the balance between the widening effect of self-spreading and the narrowing effect of advection. Our next task is to determine a scaling law for the waist width W. Dimensional analysis shows that W/L s can depend only on So/d and Lo/d. However, this dependence takes a particularly simple form in the 'narrow gap' limit Sold >> 1, which is the relevant limit for the earth (next section). Our scaling analysis for this limit has already shown that W does not depend on the gap width d. This in turn implies that W/L s depends only on the quotient of So/d and Lo/d, or equivalently on the 'buoyancy number' [5]:
n b = v2 -
To
(8)
The dependence of W/L s on H b can be determined by solving Eq. 1 numerically for different values of //b" To avoid numerical inaccuracy associated with the difficulty of resolving the sharp front at the edge of the plume, we define the waist width as the value of y for which the scaled plume thickness is 0.1 (i.e., SolS(O,W)=0.1). The results are shown in Fig. 2. The two points with H b = 1 were obtained for So/d = Lo/d = 15 and 30; their near coincidence validates our assumption that W does not depend on d. The points in Fig. 2 lie within 1.5% of a regression line given by
( a }1/2
/-/0.053
W = 1.65~
(9)
¢
The theory outlined above is readily extended to the case of a plume head with viscosity r/p spreading
157
in an ambient fluid with a higher viscosity 7/0. To determine approximately how such a viscosity contrast modifies the waist width, we model the plume head as a thin layer with viscosity %, uniform thickness S, and infinite lateral extent, overlying a fluid halfspace with viscosity T/0. The corner flow driven in this two-layer system by a surface velocity u = Utanh(x/d) can be found analytically using a sine transform. By determining the new form of the advection term in Eq. 1, evaluating it at (x,z)= (O,S), and requiring that it balance the source term, we find that the generalization of Eq. 9 is
W~1.65( U I ~ ) l/2II °°53
(10)
where 3" = */o/% and cosh k + y sinh k l('y) = 2 f o (1 + y2)cosh2k + 23, sinhZk + (1 - y2)(1 + 2k 2) dk
(11) The function 1( 7 ) is simply the ratio of the vertical velocity at (x,z)= (0,S) for the two-layer ( y ~ 1) comer flow to the same velocity for an isoviscous (3' = 1) corner flow. 1(3") is shown in Fig. 3. The predictions of our lubrication theory model can be compared directly with the laboratory measurements of Feighner and Richards [5]. These authors estimated a scaling law for W by assuming that W/L s depends only on 3' and not on //bHowever, a re-analysis of their data that allows W/L s to depend on both y and H b yields the following improved version of their law:
( Q ]~/z W ~ 1.64[ ~ )
3'0.39//0.039
For 3' = H b = 1, this result agrees to within less than 1% with our Eq. 10. Moreover, both laws predict that the waist width W should increase as either the viscosity contrast 3" or the buoyancy number //b increases. However, our Eq. 10 predicts a slightly stronger dependence on //b and a weaker dependence on 3"(W ~ 3'0.17 for large 3") than does Eq. 12. Although the reason for this discrepancy is not clear, the difference is not large enough to be of much consequence for geophysical applications. Using an isoviscous model with chemically buoy-
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
158
I
~ 0.5
~
,
I
~
I
f
,
f
f
I
,
,
,
,
I
~
_
0"0 0.0
0.5
1.0
1.5
x/L o Fig. 1. Scaled plume thickness S o l S ( x / L o , Y / L o ) for the lubrication model, obtained by solving numerically Eq. 1 with So/d = 1 and L o / d = 10. Point source is at lower left corner.
ant tracer particles, Feighner et al. [6] determined yet a third scaling law for the waist width of ridgecentered plumes: Q )1/2 W 1.28 ~ (13)
3. 3Dconvectionmodel
This law differs in two ways from those given by Eq. 10 and 12: the multiplicative constant is about 20% smaller, and there is no dependence on the buoyancy number. This suggests that the law of Feighner et al. [6] is probably less reliable than either our law (Eq.
We now turn to the study of a more realistic 3D convection model, which is sketched in Fig. 4. A rectangular box of depth h ( = 400 km) and adjustable length and width is filled with a fluid whose viscosity ~ varies as a function of depth and temper-
(
10) or our improved version of Feighner and Richards' [5] law (Eq. 12).
1 ,
,
,
,
. . . .
I
i
i
i
,
.
.
.
.
.
.
.
.
II
,
,
,
,
,
,,
, , , ,
,J
1
.
0.1
.
.
.
.
.
.
.
I
1
.
.
.
.
.
.
.
.
10
IIb
Fig. 2. Normalized waist width W / L s as a function of buoyancy number H b. Points shown were obtained by numerical solution of Eq. 1. Solid line is defined by Eq. 9.
0.1 =
........
~ 10
........ 100
Fig. 3. F u n c t i o n 1(3,) defined b y Eq. 11.
159
N.M. Ribe et aL /Earth and Planetary Science Letters 134 (1995) 155-168
ature. The symmetry of the problem allows us to limit our calculation to the region x > 0,y > 0 (note that z now increases upward). The flow E within the box consists of two parts: a passive (corner flow) component Ep generated by the motion of the upper surface, and an active component u~ driven by the buoyancy of a hot thermal plume. The passive component is forced by the same boundary condition used for the lubrication theory model, viz., u(x,y,h) = U tanh (x/d), and also satisfies the sidewall outflow condition u(x+,y,z) = u+(z), where u+ is determined as described below. The active (plume) component of the flow is generated by imposing at the bottom boundary the condition
T - Bi( w) OzT = Tr + aTp( x, y)
e+pv RT
Up = - VX ~ 0 ~ + VX VX th,2 and Ea = VX VX thE2+ VX X2
RTr
]
(
~ ° = tanh d + D c o s ( T r z / 2 )
(15)
Vx=Wx=Tx--0 ,-._ / =
(17)
where ~o(X,Z) is a streamfunction, ~bl(x,y,z) and ~b2(x,y,z) are poloidal scalars, X(x,y,z) is a toroidal scalar, and 2, ~ and 2 are unit vectors. The streamfunction qt0 is chosen a priori to satisfy all the mechanical boundary conditions shown in Fig. 3, but does not in general satisfy the governing equations. The poloidal scalar ~bI is the 'correction' required to make Ep satisfy the governing equations. The choice of ~0 is arbitrary within the constraints imposed by the boundary conditions. We use the expression
e+ ooghV] u - u.(z)
(16)
(14)
where T~ is a constant reference temperature, Bi is an effective thermal resistance that depends on the vertical velocity w, ATp = AT exp(r2/a 2) is a temperature anomaly of radius a and amplitude AT, and r is the lateral distance from the center of the plume. The horizontal velocity components u and v are zero at the bottom of the box, while the condition O-z~= 0 there permits vertical flow through this boundary. More discussion of these boundary conditions can be found in [8]. Finally, we assume that the viscosity ~7 depends on temperature and pressure according to T] -~- 7~r exp
where E is the activation energy, V is the activation volume, and ~r is the reference viscosity at the bottom of the box. We use reduced values of the activation energy and volume that mimic the effect of nonlinear theology [10]. We obtain numerical solutions using a hybrid spectral/finite difference technique [8,11], which solves directly the equations of viscous flow for depth-dependent viscosity and iterates for the effects of lateral viscosity variations. The passive and active components of the flow are written as
)Jz
u+dz
(18)
u - U tanh (x/d) = v = w = T = 0
/
/
z
/
x
I
u - u+(z) :
x=x_
x=X+vx =Wx=Tx=0 T-BiT z-Tp(x,y)= u=V=Ozz =0
Fig. 4. Sketch of the 3D convectionmodel. A rectangular box of depth h is filled with a variableviscosityfluid. The passive(ridge-like) component of the flow is driven by an imposed horizontal velocity u = U tanh(x/d), where U is the half-spreadingrate and d is the characteristic 'gap width' of the ridge. The active(buoyancy-driven)componentof the flow is generatedby a temperatureanomaly To(x,y) on the bottom of the box. Additional boundary conditions required are shown and discussed further in the text.
160
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
where D (set to 0.5) is arbitrarily as long as the value of the tanh function at x = x+ is close to unity for all depths. The shear flow profile u + ( z ) is obtained at each time step by solving the equation d(~+(z)d~z)dz
=0
(19)
where ~ + ( z ) is the profile of laterally averaged viscosity at x = x + and u+(0) = u + ( h ) U = O. This procedure ensures that the vertical profile of velocity in the fluid leaving the box is dynamically consistent with the vertical viscosity profile. We found it necessary to limit the viscosity to a cutoff value */cut= 2.5~7r The implied lithospheric viscosity of 2.5 × 1021 Pa s is quite low, although it is nevertheless two orders of magnitude higher than the minimum viscosity beneath the plate, and sufficient to ensure the mechanical coherence of the thermal boundary layer. The low cutoff is necessary to limit lateral viscosity contrasts between the hot material under the ridge and the cold adjacent plate to less than 1000. For higher lateral contrasts, the numerical method effectively overestimates the viscosity beneath the ridge. The numerical solutions described below were obtained using a grid spacing /ix = ziy = A z = h / 6 4 ( = 6.25 km) in a box with 0 < x < 2h and 0 < y < 2h. Resolution tests with different grids, as well as comparison of results of 2D test cases with those obtained using a finite-element code, indicate that the maximum error in local velocity values is of the order of 10%. The error in integrated properties, such as the waist width, is of the order of a few percent. For each numerical run, time stepping was carried out until the flow reached a steady state, which exists for all parameter values we tested. The values assumed for the model input parameters are listed in Table 1. As in the previous section, our principal goal is to determine a scaling law for the waist width W. To this end, we performed a total of 23 numerical experiments (Table 2). The volume flux of the plume was adjusted by varying the plume radius a while holding the excess plume temperature A T = 300°C constant in most runs. The volume flux Q was calculated from the relation 1 B O = - - [ w 6 T dS (20) A T Js
Po a z i T
where w is the vertical velocity, /~T is the excess plume temperature (relative to a solution with the plume absent), B is the buoyancy flux, and the integral is over the bottom surface of the box. The range of values 32.5 < Q / K h < 397 corresponds to buoyancy fluxes 350 kg s - 1 < B < 4400 kg s - 1, and the range 135 < U h / r < 900 to half-spreading rates 0.85 cm y r - l < U < 5.7 cm yr -1. The gap width was d = 0.03h for all experiments except one. The parameter tr was calculated using a density contrast zip = Po a z i T and a nominal viscosity 7/= 1019 Pa s, which in the model is intermediate between the low viscosity in the plume head and the higher ambient viscosity at z = h and T = Tr. Fig. 5 shows the temperature field for solution 3 (Table 2) in the planes x = 0 and y = 0, and Fig. 6 shows for the same solution the topography anomaly P°a
6~- (po~-Pw)
f:
P°aAT 6r dz = - ( po _ Pw) S
(21)
where Pw is the density of seawater. The quantity 3ff is just the topography anomaly produced by the plume if isostatic compensation is strictly local. In this section we ignore the contribution to the topography caused by an enhanced rate of crustal production above the plume. Note that 6ff is similar in form to the function S ( x , y ) of Fig. 1, although the sharp drop at the edge is now missing because the plume head is more diffuse. We now use our numerical results to determine a scaling law for the waist width W. Again we define W by S ( O , W ) = 0.1S 0, where S O and S are defined by Eq. 3 and 21 respectively. First, we note that in the numerical runs (Table 2) S o / d is in the range 5.2-9.6, so that all runs are close to the 'narrow gap' limit S o l d >> 1. This suggests that the finite value of d has negligible influence and that the appropriate scale for W is L s - ( Q / U ) 1/2 (Eq. 6). The validity of our assumption is supported by run 7.2 of Table 2: doubling d to 0.06 increases W / L s by only 1.5%. Fig. 7 shows the scaled waist width W / L s as a function of the buoyancy number H b for the numerical runs in Table 2. The points lie within 5% of a regression line given by
w 207( )lJ2.o052
22,
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
The form of this law agrees closely with that of the analogous law for the lubrication theory model (Eq. 9). The difference between the multiplicative constants in the two laws is only apparent, and is due to the different ways Q is defined in the two models (see the discussion section). Is the scatter of the points about the regression
161
line in Fig. 7 due to some additional systematic effect? One candidate is 'upslope' flow of low-density plume material toward the ridge along the sloping base of the rheological lithosphere. Because the slope is controlled by the plate speed, we anticipate that the importance of upslope flow will be measured by a Peclet number ULJK =-Pes. However, a plot
Table 1 Notation
Parameter a d B g E h Lo
Meaning plume radius ridge width buoyancy flux gravitational acceleration activation energy fluid d e p t h -
~ QZ/SU-1/2dl/~o'I/s
= ULo/~
_
T, U V W a
lib Po Pc p~,
reference t e m p e r a t u r e half spreading rate activation volume plume head waist width t h e r m a l expansion coefficient excess plume t e m p e r a t u r e t o p o g r a p h y anomaly plume density anomaly t e m p e r a t u r e change across plume viscosity reference viscosity t h e r m a l diffusivity buoyancy n u m b e r ( - Qo'/U 2) reference mantle density crustal density seawater density
a
- gAp/48r 1
6( Ap AT ~/ r/~
m m
volume flux plume head thickness
So
6T
Units m m kg S - 1 m s -2 J m m
(Q/U) 112
L,
S
9.81 2 × 10 s 4 x 10 s
~- Q 3 / 4 U - l o 1 / 4
Lw
Pe~
Value
t 1620 K for Iceland model. * Variable for Iceland model.
m 3 s -1 m m
1597 t 4 x 10 -s m 3.5 x 10 -s
300 ~ 10 ~1 8 x 10 - r 3.3 x 103 2.8 x 103 1.0 x 10 a
K m s -1 m3 K -x K m kg m -3 K Pa s Pa s m s S-1
kg m -3 kg m -3 kg m -3 m - i s -1
N.M. Ribe et a t / E a r t h and Planetary Science Letters 134 (1995) 155-168
162
of WH~°°52/L s vs. Pe s shows no systematic trend. We conclude that the waist width is not significantly affected by upslope flow, and that the scatter in Fig. 7 is probably due to numerical error. To determine whether the viscosity contrast between the plume head and the ambient mantle influences the waist width, we ran three additional cases (6.1, 7.1 and 8.1 in Table 2) for which the maximum excess temperature of the plume AT was reduced to 150°C and the plume radius a was increased to maintain the flux Q constant. The maximum viscosity contrast between the plume head and the ambient mantle is approximately 11 for AT = 300°C and 3.6 for A T = 150°C. The scaled waist width WH~°°52/L S differs by less than 3% for the two
cases, suggesting that W does not depend strongly on the viscosity contrast. We close this section by examining the scaling for the along-ridge topography profile 8¢(y). Fig. 8a shows the (dimensional) topography profiles for the runs in Table 2. Only the portion y > 3a of each profile is shown, because the thinlayer scaling breaks down near the plume stem. The topography profiles exhibit a large range of amplitudes and characteristic widths due to the large variations in plume flux and spreading rate among our model runs. However, if we now scale y by the waist width W and 6~ by poSo a A T / ( Po -- Pw), all the topography profiles collapse approximately onto a single curve (Fig. 8b). Even though all profiles must pass through the point ( y / W , S / S o) = (0.5,0.1)
Table 2 Variable viscosity runs
Run 1 2 3 4 5
6 6.1
7 7.1 7.2 8 8.1 9 10 11 12 13 14 15 16 17 18 19
Uh/,~ a/h 135 135 135 135 200 200 200 200 200 200 200 200 300 300 300 435 435 435 635 635 635 635 900
Q/,~h Lo/h So/h II~ UL,/,¢
WIL ,
32.5 50.7 83.3 121 32.5 51.1 55.1 84.9 89.4 85.1 139 142 51.1 85.5 143 85.6 145 241 85.1 145 247 397 247
2.58 2.58 2.70 2.96 2.59 2.57 2.57 2.63 2.58 2.67 2.85 2.66 2.49 2.48 2.58 2.45 2.45 2.61 2.48 2.41 2.41 2.69 2.37
0.I00 0.112 0.127 0.140 0.100 0.112 0.144 0.127 0.163 0.127 0.144 0.184 0.112 0.127 0.144 0.127 0.144 0.163 0.127 0.144 0.163 0.184 0.163
Runs 6.1, 7.1 and 8.1: AT = 150°C. Run 7.2: d = 0.06.
.490 .613 .786 .945 .403 .505 .525 .652 .669 .652 .834 .842 .413 .533 .691 .444 .576 .744 .366 .478 .624 .790 .524
.155 .173 .196 .215 .155 .173 .210 .197 .237 .197 .223
.266 .173 .197 .224 .197 .225 .255 .197 .225 .257 .289 .257
102 160 260 370 46.3 71.8 39.1 121 63.0 121 198 100 32.3 54.1 90.7 25.8 43.3 71.8 12.0 20.5 35.0 55.7 17.4
66 83 106 127 81 101 105 130 134 130 167 168 124 160 207 193 251 323 232 304 396 502 472
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
,
P
,
I
500
I
i
I
163
I
J
i
I
I
,
, _ _
/ -1 O0 -
400 /200
/
E ¢-- -2oo
300
-
/
/
s°°j ~ f
E v.
"o
-300
:>,
(a)
-400
I
I
100
200
'
0
'
I
f
300
400
200-
~500"
100 -
~
500
x (km) o
'
0 -100 -
S
I
200
'
I
300
400
'
500
x (km)
~3s°~
Fig. 6. Topography anomaly (m) for solution 3 (Table 2), assuming local isostatic compensation of the thermal anomaly in the mantle (crustal contribution ignored).
E V' r-- -200 Q.. 03 "0
-3oo
(b)
-400 0
I 100
'
I 200
'
1 300
'
I 400
500
y (km) Fig. 5. Temperature (°C) for solution 3 (Table 2) in the planes (a) x = 0 and (b) y = 0.
by definition, the good overall agreement confirms that the plume thickness scales approximately as S O - - ( a / / o - ) 1/4 even when the viscosity is variable.
4. Models
100
for the Iceland
hotspot
We now apply our variable viscosity model to Iceland, the best known example of a hotspot associated with a ridge-centered plume [12]. The MidAtlantic Ridge near Iceland exhibits pronounced anomalies in bathymetry, petrology, geochemistry and isotopic ratios relative to 'normal' ridge segments far from plumes. Fig. 9, for example, shows the topography along the ridge (circles) as a function
of distance from Iceland. Also shown is a symmetrized representation of the data (solid line), obtained by projecting the raw data onto a great circle with pole at 36.7°E, -14.0°S, smoothing using a + 100 km moving average, and then averaging the northern and southern branches about a presumed center of symmetry at 19.9°W, 65.6°N. Our aim in this section is to estimate the values of the maximum temperature contrast AT and the buoyancy flux B for a ridge-centered plume beneath Iceland by trying to fit the observed topography
4
I
I
I
t
l i l , l
0
_j~3-
2
.
.
.
.
.
.
.
.
i
10
0
,
,
100
FIb Fig. 7. Scaled waist width W / L s for the runs of Table 2 (integer run numbers only) as a function of buoyancy number H b. The straight line is defined by Eq. 22.
164
N.M. Ribe et al. / Earth and Planetary ScienceLetters 134 (1995) 155-168
anomaly (Fig. 9). This anomaly is the sum of two parts: a dynamic component supported by flow in the plume head, and a static component due to the excess crustal thickness generated by pressure-release melting in the plume. Strictly speaking, the dynamic topography should be calculated from the normal stress ~zz at z = h. However, the imposed (arbitrary) variation of the surface velocity near the ridge would create large artifacts. We therefore calculate the dynamic topography approximately using Eq. 21, where 6T is the excess plume temperature obtained by disregarding the effect of melting. However, the integration is carried out only from the
400
I
I
I
I
E
"-"200' ¢<3
I
600
400
y (km) 0.5
I
,
I
,
I
,
I
,
0.4-
o0.3 CO U) 0 . 2 -
0.1
0.0
,
l
,
i
0.2
0.4
0.6
0.8
i
i
i
i
i
,
i
i
i
r
o o
4oo-
0 E ~ -400o
-~ -a0o-1200
-
-l~o
o
d i s t a n c e a l o n g ridge (km)
Fig. 9. Bathymetryalong the Mid-Atlantic Ridge as a function of distance from Iceland (©, data from [13] and R. Kingsley [pets. commun., 1994]). Negative distances are south of Iceland. The smoothed and symmetrized curve shown was obtained as described in the text.
(a)
200
8o0
1.0
y/W Fig. 8. (a) Along-ridgetopographyprofiles ~ ( y ) (portion y > 3 a only) for the runs with integer numbers in Table 2. (b) Rescaled versions of the same profiles.
surface to a depth of 100 or 200 km to avoid an unrealistically large contribution from the hot plume stem. We ignored the additional contributions to the dynamic topography from the lower density of the unmelted residuum and temperature reduction due to absorption of latent heat of melting, which have opposite signs and are much smaller in magnitude than the static (crustal) topography. Our method for calculating the crustal topography consists of six steps: (1) Determine the melt fraction X ( x , y , z ) at all points in the model box, using the melting parameterization of McKenzie and Bickle [14]. The surface z = h of our 3D model is taken to be the Moho, and the pressure there is calculated from the load of the overlying crust. The model temperature is interpreted as potential temperature and is corrected for the effect of latent heat absorption during melting using the method of Watson and McKenzie [15], but assuming that the entropy of fusion is AS = 400 J k g - 1 oC - 1 [16]. (2) Calculate the local melt production rate F ( x , y , z ) = po E. V X , where ~ ( x , y , z ) is the model velocity. (3) Calculate the total melt production rate q ( y ) per unit ridge length by integrating F over vertical planes: q ( y ) = f(1 - X ) - I F d x dz.
N.M. Ribe et al. // Earth and Planetary Science Letters 134 (1995) 155-168 I
i
I
i
I
10-
v
oQ.
o
0
~ ...........
0
200
- ...........
400
,
600
800
distance along ridge (km) Fig. 10. Predicted topography anomaly (light line) for a plume with AT = 263°Cand B = 1390 kg s- 1, comparedto the smoothed observed topography(heavy line) from Fig. 9. Dynamictopography due to excess plume temperaturebetween the surface and the 100 km depth is shown by the dashed line.
(4) Determine the crustal thickness c(y) by equating q(y) to the flux of new crust carried away by the diverging plates [16]. (5) Repeat steps (1-4) iteratively to correct for the effect of the extra pressure imposed on the Moho by the crust, which tends to suppress melting beneath the ridge. (6) Use an Airy isostatic model [17] to transform c(y) to an equivalent topography anomaly relative to an assumed reference depth of 2900 m [12,17] for 'normal' ridge far from the plume. To apply our model to Iceland, we choose a half-spreading rate U = 2.73 × 10 -1° m s -1 (0.85 cm yr 1; [12]) and a mantle (potential) temperature Tr = 1324°C. Our model without a plume then predicts a 'normal' crustal thickness c(~) = 5.5 km, an appropriate average value for a slow-spreading ridge [16]. For our initial modeling attempt, we assume the values AT = 263°C and B = 1390 kg s -1 (Q = 45.7 m 3 s -1) estimated by Schilling [12] for the Iceland plume. This value of B is similar to that (B = 1400 kg s -1) estimated by Sleep [18]. In our model, these values of AT and B correspond to a plume radius a = 62 km. The size of the model box and the grid is the same as in the previous section. Fig. 10 shows the dynamic (dashed line) and the total (light line) topography predicted by the model, together with the
165
smoothed topography (heavy line) from Fig. 9. Close to the plume, the topography anomaly is dominated by the crustal contribution (light line minus dashed line). At y = 0, the predicted topography is ~ 10 km, which corresponds to a crustal thickness of 68 km (the crustal thickness would be 320 km if the damping effect of the extra crustal pressure were ignored). This far exceeds the crustal thickness of ~ 20-34 km predicted for the same AT but purely passive upwelling [14,19]. Due to the active plume upwelling, the upward flux of source rock to shallow depth under the ridge at y = 0 is an order of magnitude higher than the flux for passive spreading. At a few hundred kilometers distance from the plume not much excess melting occurs, although the temperature below the ridge is significantly higher than ambient mantle temperature (compare Fig. 5a with parameter values similar to those for the Iceland model). This is because the plume stem rises almost vertically to shallow depth, where it loses a great deal of its potential melt, and then already depleted rock spreads almost horizontally along the strike of the ridge. In any case, Fig. 10 clearly shows that the topography anomaly predicted by this model is inconsistent with the observations. Fig. 11 shows the topography anomaly predicted for a much cooler and broader plume with AT = 93°C
i
I
~
I
i
I
3-
8' 1
o
i
i
i
400
800
1200
distance along ridge (km) Fig. 11. Same as Fig. 10, but for a plume with AT=93°C, a = 300 km, and B = 1370 kg s -1. Additional topography data ( ~ ) from [19]. Dynamic topography obtained by integration to the 200 km depth.
166
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
and radius a = 300 km. In order to achieve an acceptable fit to the topography at great distances from Iceland, we increased the mantle potential temperature outside the plume to 1347°C compared to 1324°C for normal upper mantle. This simulates in a crude way a very broad thermal halo of the plume or a generally anomalous upper mantle beneath the North Atlantic. Compared to the temperature of the halo, the excess temperature of the plume is 70°C. The buoyancy flux for this plume is B = 1370 kg s -1, the volume flux is Q = 193 m 3 s -1, and the mean upwelling velocity is 2 cm yr -1. The numerical solution was calculated in a box of 800 km in length and 1200 km in width on a grid of 129 × 161 X 65 points. The predicted topography anomaly is now much closer to the observed one, although it is still too narrow in the center. An even wider and cooler upwelling could probably lead to a better fit, but limitations of computer memory prevented us from obtaining a numerical solution in the larger model box that would be required for this.
5. Discussion
The principal goal of this study has been to determine the scaling law for the 'waist width' W for a plume rising beneath an ocean ridge. We found that both a simple lubrication theory model and a more realistic 3D convection model predict similar scaling laws (Eq. 10 and 22). This indicates that the essential dynamics of plume-ridge interaction are largely independent of whether the buoyancy of the plume is thermal or chemical in origin. Nevertheless, there are two relatively minor differences between the two laws. The first is that the lubrication theory model predicts a dependence of W on the mantleplume viscosity contrast "y, whereas no such dependence is evident in the 3D convection model. However, neither model is really capable of resolving this dependence properly: the lubrication theory model incorporates the viscosity contrast only in an approximate way, and numerical difficulties prevent us from using large lateral viscosity contrasts in the 3D model. The second difference between the two laws is that the waist width W is about 25% higher for the 3D model (Eq. 22) than for the lubrication theory model (Eq. 10). However, this discrepancy is only
apparent, because the volume flux Q (unlike the buoyancy flux B) is not uniquely defined for a thermal plume. We chose to define Q as B / p o t~AT, where AT is the maximum temperature contrast across the plume. However, it is probably more appropriate to replace AT in this definition by a smaller average temperature difference ATav, which would reduce the difference between the two results. Another effect that could explain some of the difference between the two results is conductive heat loss to the surface in the variable viscosity model. Our variable viscosity numerical results provide no evidence that the waist width is significantly affected by ' upslope' flow of plume material toward the ridge along the sloping base of the rheological lithosphere. This result may seem surprising until one considers that whereas upslope flow tends to increase W, the 'freezing' of plume material into the lithosphere near the ridge tends to decrease it. Moreover, the magnitude of both effects increases as the slope of the base of the lithosphere increases. We therefore suggest that the absence of evidence for upslope flow in our numerical experiments is due to a balance between the opposing effects of upslope flow and freezing. We do however anticipate that upslope flow may play a greater role in the dynamics of off-ridge plumes. The most surprising result of this study, however, is that the observed anomalous topography of the Mid-Atlantic Ridge suggests the presence of a very cool (AT< 100°C), wide (radius a > 300 km) and slow (vertical velocity ~ 2 cm yr 1) upwelling beneath Iceland. By contrast, a hot (AT = 263°C) and narrow (a = 62 km) plume of the type hypothesized by Schilling [12] predicts a topography anomaly that is far too narrow and a rate of crust production directly over the plume stem that is far too high. Whatever values of a and AT are assumed, the static component of the topography due to excess crustal thickness exceeds the dynamic component due to flow in the plume head. The topography anomaly our model predicts thus depends on the assumptions made in calculating the excess crustal thickness. The most important of these assumptions is that melt generated beneath the ridge is extracted sufficiently rapidly that the ridge-parallel component of its velocity is zero. If this were not so, melt generated near the plume might be extracted further
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
away from it, thereby widening the topography anomaly and making a 'narrow plume' model more viable. The importance of along-ridge melt migration can be estimated by means of a simple calculation. The velocity of the melt at any point (O,y,z) directly beneath the ridge is u'f = v~ + (w + wf)2, where v and w are the lateral and vertical components of velocity of the unmelted matrix and wf is the velocity of buoyancy-driven vertical melt migration relative to the matrix. If melt migration obeys Darcy's Law for flow in a porous medium, wf is given by
k4~gApf wf-
~b~?~--~-
(23)
where k6 is the permeability of the matrix, ¢ is the porosity, Apf is the difference between the densities of the matrix and the melt, and ~/f is the viscosity of the melt. Because the matrix has a component of velocity in the ~ direction, a parcel of melt ascending from depth h - z to the surface will travel a distance Ay ~ (h - z)v/(w + wf) along the ridge. As an example, consider the point (x,y,z)= (0 km,100 km,275 km) in our 'narrow plume' numerical solution where the along-ridge velocity v has its maximum value. At this point, v = 1.6 × 1 0 - 9 m s -1 and w = 1.0 × 1 0 - 9 m s 1. Moreover, by using the values [20] k6 = 1 0 - 1 4 m 2, A p f = 600 kg m 3, r/f = 1 Pa s and ~b = 0.01 we find wf = 5.9 × 10 9 m s -1. Thus Ay ~ 30 km. Lateral melt migration over such small distances is incapable of widening the topography anomaly predicted by our 'narrow plume' model (Fig. 10) sufficiently to match the observations. Our disregard of lateral melt migration is therefore valid to a first order. We close by noting that the low excess temperature of the Iceland plume in our model (AT = 93°C) is consistent with its low buoyancy flux. The buoyancy flux of the Iceland plume is smaller than that of the Hawaiian plume by a factor of 3 - 5 [8,18]. In the case of Hawaii, an excess temperature of 250-300°C is in agreement with the observed geometry of the topographic swell [8], and is required to explain the production rates and composition of melts [15]. O1son et al. [23] have shown that the rate at which plumes cool and broaden by diffusion and entrainment as they rise through the mantle depends inversely on the buoyancy flux. If we assume that both
167
the Hawaiian plume and the Iceland plume rise from the core-mantle boundary, where they start with the same excess temperature, the theory of Olson et al. predicts that the Hawaiian plume would be scarcely affected by diffusive cooling and broadening. For the Iceland plume, however, the effect would be quite significant, in agreement with our results. By implication, plumes that are much weaker than the Iceland plume cannot possibly rise from the core-mantle boundary.
Acknowledgements This research was supported by NSF grant EAR 92-04644 and a fellowship from the Alexander yon Humboldt-Stiflung (Germany) to N.M.R., and by grants Ch-77/4 and Ch-77/8 from the Deutsche Forschungsgemeinschafi to U.R.C. We thank J. Bown, L. Fleitout, R. Kingsley, J.G. Schilling and R. White for their generous assistance, M. Feighner and M. Richards for making their work available to us before publication, and C. Langmuir, D. McKenzie and an anonymous referee for helpful reviews. [CL]
Appendix A Lubrication theory equation The general differential equation governing the steady thickness S(x,y) of a thin buoyant fluid layer spreading in an ambient flow whose horizontal velocity is uve = u~} + v~ can be obtained using lubrication theory [7], and is
~7. foS~dz = 0-~7234 + Q6( x) 6(y)
(24)
where Vnz is the horizontal Laplacian and ~r and Q are defined in the text. For the case of a 2D corner flow described by a streamfunction ~ ( x , z )
V f o dZ=OxfSudz--Ox[ ,(x,S)] An analytical expression for ~ ( x , boundary condition ~z(X,O)= U determined by writing ~ = ~z(°)(x,0) = U and ~z(1)(x,0) =
(25)
z) that satisfies the tanh(x/d) can be ~(oh/t (1), where U[1 - tanh(x/d)].
168
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
The streamfunction qt(0)- 2 . r r - l U z tan-l(x/z) is just the standard corner flow solution [21], and 1/f(1) is the correction required by the modified velocity near x = 0. The biharmonic equation governing ~(1) can be solved by applying a sine transform and invoking formula 3.987.1 of [22], yielding q,(1) = _~Udf °~ sinkx [ 1kd 2
2 sinh(~kd/2) 7r ]z
× exp( - kz)dk
(26)
This integral can be evaluated using Eq. 3.941.1 and 4.131.3 of [22], whereupon we find that the first part cancels ~(0). The final result, after evaluation at z=S, is 2
(1
S
/x)
~ ( x,S) = -~USIm tp ~ + zrd + zr---d
(27)
where 'Im' denotes the imaginary part and ~O(z)-= (d/dz)ln F(z) is the digamma function of a complex argument [9].
References [1] W.J. Morgan, Plate motions and deep mantle convection, Mem. Geol. Soc. Am. 132, 7-22, 1972. [2] M. Talwani, C.C. Windisch and M.G. Langseth, Jr., Reykjanes Ridge crest: A detailed geophysical study, J. Geophys. Res. 76, 473-517, 1971. [3] J.G. Schilling, Iceland mantle plume: Geochemical evidence along Reykjanes Ridge, Nature 242, 565-571, 1,973. [4] J.G. Schilling, Upper mantle heterogeneities and dynamics, Nature 314, 62-67, 1985. [5] M.A. Feighner and M.A. Richards, The fluid dynamics of plume-ridge and plume-plate interactions: an experimental investigation, Earth Planet. Sci. Lett. 129, 171-182, 1995. [6] M.A. Feighner, L.H. Kellogg and B.J. Travis, Numerical modeling of chemically buoyant mantle plumes at spreading ridges, Geophys. Res. Lett. 22, 715-718, 1995. [7] P. Olson, Hot spots, swells and mantle plumes, in: Magma
Transport and Storage, M.P. Ryan, ed., pp. 33-51, Wiley, New York, 1990. [8] N.M. Ribe and U. Christensen, Three-dimensional modeling of plume-lithosphere interaction, J. Geophys. Res. 99, 669682, 1994. [9] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1046 pp., Dover, New York, 1965. [10] U. Christensen, Convection with pressure and temperature dependent non-Newtonian rheology, Geophys. J. R. Astron. Soc. 77, 242-284, 1984. [11] U. Christensen and H. Harder, 3D convection with variable viscosity, Geophys. J. Int. 104, 213-226, 1991. [12] J.G. Schilling, Fluxes and excess temperatures of mantle plumes inferred from their interaction with migrating mid-ocean ridges, Nature 352, 397-403, 1991. [13] J.G. Schilling, M. Zajac, R. Evans, T. Johnston, W. White, J.D. Devine and R. Kingsley, Petrologic and geochemical variations along the mid-Atlantic ridge from 29°N to 73°N, Am. J. Sci. 283, 510-586, 1983. [14] D. McKenzie and M. Bickle, The volume and composition of melt generated by extension of the lithosphere, J. Petrol. 29, 625-679, 1988. [15] S. Watson and D. McKenzie, Melt generation in plumes: a study of Hawaiian volcanism, J. Petrol. 32, 501-537, 1991. [16] J.W. Bown and R.S. White, Variation with spreading rate of oceanic crustal thickness and geochemistry, Earth Planet. Sci. Lett. 121, 435-449, 1994. [17] G. Schubert and D. Sandwell, Crustal volumes of the continents and of oceanic and continental submarine plateaus, Earth Planet. Sci. Lett. 92, 234-246, 1989. [18] N.H. Sleep, Hotspots and mantle plumes: some phenomenology, J. Geophys. Res. 95, 6715-6736, 1990. [19] E.M. Klein and C.H. Langmuir, Global correlation of ocean ridge basalt chemistry with axial depth and crustal thickness, J. Geophys. Res. 92, 8089-8115, 1987. [20] M.J. Cordery and J. Phipps Morgan, Convection and melting at mid-ocean ridges, J. Geophys. Res. 98, 19477-19503, 1993. [21] G.K. Batchelor, An Introduction to Fluid Dynamics, 615 pp., Cambridge University Press, Cambridge, 1967. [22] I.M. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, 1160 pp., Academic Press, New York, 1965. [23] P. Olson, G. Schubert and C. Anderson, Structure of axisymmetric mantle plumes, J. Geophys. Res. 98, 6829-6844, 1993.
Earth and Planetary Science Letters 134 (1995) 155-168
The dynamics of plume-ridge interaction, 1: Ridge-centered plumes N.M. Ribe 1, U.R. Christensen, J. Theiging Institut ffir Geophysik, Universiti~t Gfittingen, Herzberger Landstrasse 180, 37075 G6ttingen, Germany
Received 12 January 1995; revised 25 May 1995; accepted 6 June 1995
Abstract
We investigate the dynamics of mantle plumes rising beneath mid-ocean ridges. To clarify the physics of this process, we examine first a simple 'lubrication theory' model in which a point source (analogous to a plume 'stem') of volume flux Q located directly beneath the ridge releases buoyant fluid into a viscous comer flow driven by a velocity boundary condition u(x) = U tanh(x/d), where U is the half-spreading rate and the 'gap width' between the diverging plates is ~ 5d. Numerical solutions of the differential equation governing the plume head thickness S ( x , y ) show how the width W of the plume head along the ridge depends on Q,U,d, and o- =- gAp~48,1, where Ap is the density deficit of the plume and ~/ is the viscosity. In the geophysically relevant 'narrow gap' limit (0//0-) 1/4 >> d, W ~ (Q//U)I/2IIb 0"053, where Mb = Q0-//U 2 is the 'buoyancy number'. Numerical solutions of a more realistic 3D convection model with strongly temperature and pressure-dependent viscosity obey a nearly identical scaling law, and show no evidence that W is increased by 'upslope' flow of plume material toward the ridge along the sloping base of the rheological lithosphere. To apply our model to Iceland, we incorporate into it a melting parameterization that allows prediction of the excess crustal thickness produced by melting in the plume head. This extended model shows that the observed depth anomalies along the Mid-Atlantic Ridge near Iceland cannot be explained by a hot (temperature contrast AT ~ 250°C) and narrow (radius ~ 60 km) ridge-centered plume. Instead, the anomalies are consistent with a much cooler and broader upwelling. 1. Introduction
Many of the earth's hotspots lie on or very near mid-ocean ridges. Examples include the Iceland, Azores, Afar and Galapagos hotspots. The strong correlation of hotspot locations with ocean ridges was first noted by Morgan [1]. At about the same time, extensive surveying and sampling of the ocean floor revealed that the association of hotspots with ridges is not merely geometrical, but physical and
1 Permanent address: Department of Geology and Geophysics, Yale University, POB 208109, New Haven, CT 06520-8109, USA Elsevier Science B.V. SSD! 0012-821X(95)00116-6
chemical as well. Near Iceland, for example, the seafloor topography and the L a / S m ratio of dredged basalts are anomalously high for several hundred kilometers along the Mid-Atlantic Ridge [2,3], and similar anomalies exist near other hotspots located on ridges [4]. These anomalies are thought to be due to the lateral spreading of hot plume material beneath the ridge [3]. The present paper is a study of the dynamics of plume-ridge interaction for the special case of a ' ridge-centered' plume rising directly beneath a ridge. Our principal goal is to determine how far buoyant plume material can spread along the strike of a ridge before being swept 'downstream' by the diverging
156
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
plates. Our study consists of three parts. We first examine a simple 'lubrication theory' model that reveals the fundamental dynamics of plume-ridge interaction, and we compare the predictions of this model with the scaling laws of Feighner and Richards [5] and Feighner et al. [6]. We then study a more realistic 3D numerical convection model with temperature and pressure-dependent rheology. Finally, we apply the convection model to the Iceland hotspot, the classic example of a ridge-centered plume.
plex argument [9]. The quantity UzF(x,z) is simply the streamfunction for the corner flow in the halfspace z > 0. Eq. 1 is a nonlinear advection/diffusion equation describing a balance between injection of plume fluid, its advection by the corner flow, and its buoyancy-driven ' self-spreading'. We first use Eq. 1 to determine the length scales for the thickness S and width (normal to the spreading direction) L(x) of the plume head. By balancing the self-spreading and source terms we obtain o-S4/L 2 ~
2. Lubrication theory model
S~
In this section we make use of the theory of flow in thin layers to formulate a simple model for plume-ridge interaction. The idea behind the model is that buoyant plume material ascending from great depth through a narrow conduit will spread laterally when it reaches the base of the diverging plates, forming a layer that is thin compared to its lateral extent. Similar models have already been employed successfully to describe the geometry of swells produced by plumes far from ridges [7,8]. Our model envisions a passive 'corner flow' driven in a halfspace z > 0 of fluid with constant viscosity 7/ and density p by a lateral velocity u(x) = U tanh ( x / d ) imposed at the surface (z = 0). Here U is the half-spreading rate. Moreover, because u(__+2.5d) = _0.99U, the 'gap width' between the diverging plates is ~ 5d. Within the corner flow, a point source (analogous to the top of the plume 'stem') of volume flux Q located at ( x , y ) = (0,0) (directly beneath the ridge) releases buoyant fluid with viscosity 77 and density p - zip, which spreads along the upper surface of the halfspace as a layer (the plume 'head') of thickness S(x,y). The differential equation governing S can be derived using lubrication theory (Appendix), and is
U(FS) x =
0"~7254
-~ Q 6 ( X ) 6 ( y )
2 (1 _ S_ /x) F ( x , S ) = -~ Im ~O -~ + ~ d + ~--d
(2)
where 'Ira' denotes the imaginary part and O(p)--( d / d p ) l n F ( p ) is the digamma function of a c o m -
- So
(3)
The scale for L can now be determined by balancing the advection and source terms in Eq. 1. The result depends on one's distance from the ridge. Far from the ridge (x >> d), F = 1, which implies U(FS) x US/L. Balancing this with the source term (order Q/L2), we obtain
L
~ Q3/4U-lo'l/4
=L o
(4)
This lengthscale is identical to that derived for a plume beneath a surface moving at uniform velocity U [81. At the ridge (x = 0), however, the 'waist width' [5] L(0) - W is controlled by two additional lengthscales. Because S x = 0 at the ridge by symmetry, the advection term there is
USFx(O,S ) = - ~ b
~ 7rd +
(5)
This expression has two distinct limiting forms, depending on the ratio of the characteristic layer thickness S O to the gap width d. In the 'narrow gap' limit So/d >> 1, USFx(O,S) = 2U/Tr, which when balanced with the source term (order Q / W 2) yields
Q )1/2
(1)
where o-= gAp/48rl, Vn is the horizontal Laplacian, and
a l L 2, or
W~(~
-Ls
(6)
In the 'wide gap' limit So/d << 1, by contrast, USFx(O,S) = US~d, which when balanced with the source term yields W ~ Q3/8U-1/2dl/2or
1/8 _ Lw
(7)
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
We now examine the solutions of the full lubrication equation 1. To avoid a singularity, we replace the point source by an equivalent Gaussian source distribution with a small but finite radius a << W, viz., Q6(x)6( y) --, (Q/Tra2)exp[-(x z + y 2 ) / a 2 ] . Nondimensionalization of Eq. 1 shows that the solutions depend only on the dimensionless groups Sold and Lo/d. Fig. 1 shows the solution for So/d--- 1 and Lo/d = 10, obtained using a finite-difference technique with explicit time stepping [7] and a grid spacing Ax = Ay = 0.02L 0. The increasing width of the plume head away from the ridge, which is typical of the numerical solutions, reflects the balance between the widening effect of self-spreading and the narrowing effect of advection. Our next task is to determine a scaling law for the waist width W. Dimensional analysis shows that W/L s can depend only on So/d and Lo/d. However, this dependence takes a particularly simple form in the 'narrow gap' limit Sold >> 1, which is the relevant limit for the earth (next section). Our scaling analysis for this limit has already shown that W does not depend on the gap width d. This in turn implies that W/L s depends only on the quotient of So/d and Lo/d, or equivalently on the 'buoyancy number' [5]:
n b = v2 -
To
(8)
The dependence of W/L s on H b can be determined by solving Eq. 1 numerically for different values of //b" To avoid numerical inaccuracy associated with the difficulty of resolving the sharp front at the edge of the plume, we define the waist width as the value of y for which the scaled plume thickness is 0.1 (i.e., SolS(O,W)=0.1). The results are shown in Fig. 2. The two points with H b = 1 were obtained for So/d = Lo/d = 15 and 30; their near coincidence validates our assumption that W does not depend on d. The points in Fig. 2 lie within 1.5% of a regression line given by
( a }1/2
/-/0.053
W = 1.65~
(9)
¢
The theory outlined above is readily extended to the case of a plume head with viscosity r/p spreading
157
in an ambient fluid with a higher viscosity 7/0. To determine approximately how such a viscosity contrast modifies the waist width, we model the plume head as a thin layer with viscosity %, uniform thickness S, and infinite lateral extent, overlying a fluid halfspace with viscosity T/0. The corner flow driven in this two-layer system by a surface velocity u = Utanh(x/d) can be found analytically using a sine transform. By determining the new form of the advection term in Eq. 1, evaluating it at (x,z)= (O,S), and requiring that it balance the source term, we find that the generalization of Eq. 9 is
W~1.65( U I ~ ) l/2II °°53
(10)
where 3" = */o/% and cosh k + y sinh k l('y) = 2 f o (1 + y2)cosh2k + 23, sinhZk + (1 - y2)(1 + 2k 2) dk
(11) The function 1( 7 ) is simply the ratio of the vertical velocity at (x,z)= (0,S) for the two-layer ( y ~ 1) comer flow to the same velocity for an isoviscous (3' = 1) corner flow. 1(3") is shown in Fig. 3. The predictions of our lubrication theory model can be compared directly with the laboratory measurements of Feighner and Richards [5]. These authors estimated a scaling law for W by assuming that W/L s depends only on 3' and not on //bHowever, a re-analysis of their data that allows W/L s to depend on both y and H b yields the following improved version of their law:
( Q ]~/z W ~ 1.64[ ~ )
3'0.39//0.039
For 3' = H b = 1, this result agrees to within less than 1% with our Eq. 10. Moreover, both laws predict that the waist width W should increase as either the viscosity contrast 3" or the buoyancy number //b increases. However, our Eq. 10 predicts a slightly stronger dependence on //b and a weaker dependence on 3"(W ~ 3'0.17 for large 3") than does Eq. 12. Although the reason for this discrepancy is not clear, the difference is not large enough to be of much consequence for geophysical applications. Using an isoviscous model with chemically buoy-
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
158
I
~ 0.5
~
,
I
~
I
f
,
f
f
I
,
,
,
,
I
~
_
0"0 0.0
0.5
1.0
1.5
x/L o Fig. 1. Scaled plume thickness S o l S ( x / L o , Y / L o ) for the lubrication model, obtained by solving numerically Eq. 1 with So/d = 1 and L o / d = 10. Point source is at lower left corner.
ant tracer particles, Feighner et al. [6] determined yet a third scaling law for the waist width of ridgecentered plumes: Q )1/2 W 1.28 ~ (13)
3. 3Dconvectionmodel
This law differs in two ways from those given by Eq. 10 and 12: the multiplicative constant is about 20% smaller, and there is no dependence on the buoyancy number. This suggests that the law of Feighner et al. [6] is probably less reliable than either our law (Eq.
We now turn to the study of a more realistic 3D convection model, which is sketched in Fig. 4. A rectangular box of depth h ( = 400 km) and adjustable length and width is filled with a fluid whose viscosity ~ varies as a function of depth and temper-
(
10) or our improved version of Feighner and Richards' [5] law (Eq. 12).
1 ,
,
,
,
. . . .
I
i
i
i
,
.
.
.
.
.
.
.
.
II
,
,
,
,
,
,,
, , , ,
,J
1
.
0.1
.
.
.
.
.
.
.
I
1
.
.
.
.
.
.
.
.
10
IIb
Fig. 2. Normalized waist width W / L s as a function of buoyancy number H b. Points shown were obtained by numerical solution of Eq. 1. Solid line is defined by Eq. 9.
0.1 =
........
~ 10
........ 100
Fig. 3. F u n c t i o n 1(3,) defined b y Eq. 11.
159
N.M. Ribe et aL /Earth and Planetary Science Letters 134 (1995) 155-168
ature. The symmetry of the problem allows us to limit our calculation to the region x > 0,y > 0 (note that z now increases upward). The flow E within the box consists of two parts: a passive (corner flow) component Ep generated by the motion of the upper surface, and an active component u~ driven by the buoyancy of a hot thermal plume. The passive component is forced by the same boundary condition used for the lubrication theory model, viz., u(x,y,h) = U tanh (x/d), and also satisfies the sidewall outflow condition u(x+,y,z) = u+(z), where u+ is determined as described below. The active (plume) component of the flow is generated by imposing at the bottom boundary the condition
T - Bi( w) OzT = Tr + aTp( x, y)
e+pv RT
Up = - VX ~ 0 ~ + VX VX th,2 and Ea = VX VX thE2+ VX X2
RTr
]
(
~ ° = tanh d + D c o s ( T r z / 2 )
(15)
Vx=Wx=Tx--0 ,-._ / =
(17)
where ~o(X,Z) is a streamfunction, ~bl(x,y,z) and ~b2(x,y,z) are poloidal scalars, X(x,y,z) is a toroidal scalar, and 2, ~ and 2 are unit vectors. The streamfunction qt0 is chosen a priori to satisfy all the mechanical boundary conditions shown in Fig. 3, but does not in general satisfy the governing equations. The poloidal scalar ~bI is the 'correction' required to make Ep satisfy the governing equations. The choice of ~0 is arbitrary within the constraints imposed by the boundary conditions. We use the expression
e+ ooghV] u - u.(z)
(16)
(14)
where T~ is a constant reference temperature, Bi is an effective thermal resistance that depends on the vertical velocity w, ATp = AT exp(r2/a 2) is a temperature anomaly of radius a and amplitude AT, and r is the lateral distance from the center of the plume. The horizontal velocity components u and v are zero at the bottom of the box, while the condition O-z~= 0 there permits vertical flow through this boundary. More discussion of these boundary conditions can be found in [8]. Finally, we assume that the viscosity ~7 depends on temperature and pressure according to T] -~- 7~r exp
where E is the activation energy, V is the activation volume, and ~r is the reference viscosity at the bottom of the box. We use reduced values of the activation energy and volume that mimic the effect of nonlinear theology [10]. We obtain numerical solutions using a hybrid spectral/finite difference technique [8,11], which solves directly the equations of viscous flow for depth-dependent viscosity and iterates for the effects of lateral viscosity variations. The passive and active components of the flow are written as
)Jz
u+dz
(18)
u - U tanh (x/d) = v = w = T = 0
/
/
z
/
x
I
u - u+(z) :
x=x_
x=X+vx =Wx=Tx=0 T-BiT z-Tp(x,y)= u=V=Ozz =0
Fig. 4. Sketch of the 3D convectionmodel. A rectangular box of depth h is filled with a variableviscosityfluid. The passive(ridge-like) component of the flow is driven by an imposed horizontal velocity u = U tanh(x/d), where U is the half-spreadingrate and d is the characteristic 'gap width' of the ridge. The active(buoyancy-driven)componentof the flow is generatedby a temperatureanomaly To(x,y) on the bottom of the box. Additional boundary conditions required are shown and discussed further in the text.
160
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
where D (set to 0.5) is arbitrarily as long as the value of the tanh function at x = x+ is close to unity for all depths. The shear flow profile u + ( z ) is obtained at each time step by solving the equation d(~+(z)d~z)dz
=0
(19)
where ~ + ( z ) is the profile of laterally averaged viscosity at x = x + and u+(0) = u + ( h ) U = O. This procedure ensures that the vertical profile of velocity in the fluid leaving the box is dynamically consistent with the vertical viscosity profile. We found it necessary to limit the viscosity to a cutoff value */cut= 2.5~7r The implied lithospheric viscosity of 2.5 × 1021 Pa s is quite low, although it is nevertheless two orders of magnitude higher than the minimum viscosity beneath the plate, and sufficient to ensure the mechanical coherence of the thermal boundary layer. The low cutoff is necessary to limit lateral viscosity contrasts between the hot material under the ridge and the cold adjacent plate to less than 1000. For higher lateral contrasts, the numerical method effectively overestimates the viscosity beneath the ridge. The numerical solutions described below were obtained using a grid spacing /ix = ziy = A z = h / 6 4 ( = 6.25 km) in a box with 0 < x < 2h and 0 < y < 2h. Resolution tests with different grids, as well as comparison of results of 2D test cases with those obtained using a finite-element code, indicate that the maximum error in local velocity values is of the order of 10%. The error in integrated properties, such as the waist width, is of the order of a few percent. For each numerical run, time stepping was carried out until the flow reached a steady state, which exists for all parameter values we tested. The values assumed for the model input parameters are listed in Table 1. As in the previous section, our principal goal is to determine a scaling law for the waist width W. To this end, we performed a total of 23 numerical experiments (Table 2). The volume flux of the plume was adjusted by varying the plume radius a while holding the excess plume temperature A T = 300°C constant in most runs. The volume flux Q was calculated from the relation 1 B O = - - [ w 6 T dS (20) A T Js
Po a z i T
where w is the vertical velocity, /~T is the excess plume temperature (relative to a solution with the plume absent), B is the buoyancy flux, and the integral is over the bottom surface of the box. The range of values 32.5 < Q / K h < 397 corresponds to buoyancy fluxes 350 kg s - 1 < B < 4400 kg s - 1, and the range 135 < U h / r < 900 to half-spreading rates 0.85 cm y r - l < U < 5.7 cm yr -1. The gap width was d = 0.03h for all experiments except one. The parameter tr was calculated using a density contrast zip = Po a z i T and a nominal viscosity 7/= 1019 Pa s, which in the model is intermediate between the low viscosity in the plume head and the higher ambient viscosity at z = h and T = Tr. Fig. 5 shows the temperature field for solution 3 (Table 2) in the planes x = 0 and y = 0, and Fig. 6 shows for the same solution the topography anomaly P°a
6~- (po~-Pw)
f:
P°aAT 6r dz = - ( po _ Pw) S
(21)
where Pw is the density of seawater. The quantity 3ff is just the topography anomaly produced by the plume if isostatic compensation is strictly local. In this section we ignore the contribution to the topography caused by an enhanced rate of crustal production above the plume. Note that 6ff is similar in form to the function S ( x , y ) of Fig. 1, although the sharp drop at the edge is now missing because the plume head is more diffuse. We now use our numerical results to determine a scaling law for the waist width W. Again we define W by S ( O , W ) = 0.1S 0, where S O and S are defined by Eq. 3 and 21 respectively. First, we note that in the numerical runs (Table 2) S o / d is in the range 5.2-9.6, so that all runs are close to the 'narrow gap' limit S o l d >> 1. This suggests that the finite value of d has negligible influence and that the appropriate scale for W is L s - ( Q / U ) 1/2 (Eq. 6). The validity of our assumption is supported by run 7.2 of Table 2: doubling d to 0.06 increases W / L s by only 1.5%. Fig. 7 shows the scaled waist width W / L s as a function of the buoyancy number H b for the numerical runs in Table 2. The points lie within 5% of a regression line given by
w 207( )lJ2.o052
22,
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
The form of this law agrees closely with that of the analogous law for the lubrication theory model (Eq. 9). The difference between the multiplicative constants in the two laws is only apparent, and is due to the different ways Q is defined in the two models (see the discussion section). Is the scatter of the points about the regression
161
line in Fig. 7 due to some additional systematic effect? One candidate is 'upslope' flow of low-density plume material toward the ridge along the sloping base of the rheological lithosphere. Because the slope is controlled by the plate speed, we anticipate that the importance of upslope flow will be measured by a Peclet number ULJK =-Pes. However, a plot
Table 1 Notation
Parameter a d B g E h Lo
Meaning plume radius ridge width buoyancy flux gravitational acceleration activation energy fluid d e p t h -
~ QZ/SU-1/2dl/~o'I/s
= ULo/~
_
T, U V W a
lib Po Pc p~,
reference t e m p e r a t u r e half spreading rate activation volume plume head waist width t h e r m a l expansion coefficient excess plume t e m p e r a t u r e t o p o g r a p h y anomaly plume density anomaly t e m p e r a t u r e change across plume viscosity reference viscosity t h e r m a l diffusivity buoyancy n u m b e r ( - Qo'/U 2) reference mantle density crustal density seawater density
a
- gAp/48r 1
6( Ap AT ~/ r/~
m m
volume flux plume head thickness
So
6T
Units m m kg S - 1 m s -2 J m m
(Q/U) 112
L,
S
9.81 2 × 10 s 4 x 10 s
~- Q 3 / 4 U - l o 1 / 4
Lw
Pe~
Value
t 1620 K for Iceland model. * Variable for Iceland model.
m 3 s -1 m m
1597 t 4 x 10 -s m 3.5 x 10 -s
300 ~ 10 ~1 8 x 10 - r 3.3 x 103 2.8 x 103 1.0 x 10 a
K m s -1 m3 K -x K m kg m -3 K Pa s Pa s m s S-1
kg m -3 kg m -3 kg m -3 m - i s -1
N.M. Ribe et a t / E a r t h and Planetary Science Letters 134 (1995) 155-168
162
of WH~°°52/L s vs. Pe s shows no systematic trend. We conclude that the waist width is not significantly affected by upslope flow, and that the scatter in Fig. 7 is probably due to numerical error. To determine whether the viscosity contrast between the plume head and the ambient mantle influences the waist width, we ran three additional cases (6.1, 7.1 and 8.1 in Table 2) for which the maximum excess temperature of the plume AT was reduced to 150°C and the plume radius a was increased to maintain the flux Q constant. The maximum viscosity contrast between the plume head and the ambient mantle is approximately 11 for AT = 300°C and 3.6 for A T = 150°C. The scaled waist width WH~°°52/L S differs by less than 3% for the two
cases, suggesting that W does not depend strongly on the viscosity contrast. We close this section by examining the scaling for the along-ridge topography profile 8¢(y). Fig. 8a shows the (dimensional) topography profiles for the runs in Table 2. Only the portion y > 3a of each profile is shown, because the thinlayer scaling breaks down near the plume stem. The topography profiles exhibit a large range of amplitudes and characteristic widths due to the large variations in plume flux and spreading rate among our model runs. However, if we now scale y by the waist width W and 6~ by poSo a A T / ( Po -- Pw), all the topography profiles collapse approximately onto a single curve (Fig. 8b). Even though all profiles must pass through the point ( y / W , S / S o) = (0.5,0.1)
Table 2 Variable viscosity runs
Run 1 2 3 4 5
6 6.1
7 7.1 7.2 8 8.1 9 10 11 12 13 14 15 16 17 18 19
Uh/,~ a/h 135 135 135 135 200 200 200 200 200 200 200 200 300 300 300 435 435 435 635 635 635 635 900
Q/,~h Lo/h So/h II~ UL,/,¢
WIL ,
32.5 50.7 83.3 121 32.5 51.1 55.1 84.9 89.4 85.1 139 142 51.1 85.5 143 85.6 145 241 85.1 145 247 397 247
2.58 2.58 2.70 2.96 2.59 2.57 2.57 2.63 2.58 2.67 2.85 2.66 2.49 2.48 2.58 2.45 2.45 2.61 2.48 2.41 2.41 2.69 2.37
0.I00 0.112 0.127 0.140 0.100 0.112 0.144 0.127 0.163 0.127 0.144 0.184 0.112 0.127 0.144 0.127 0.144 0.163 0.127 0.144 0.163 0.184 0.163
Runs 6.1, 7.1 and 8.1: AT = 150°C. Run 7.2: d = 0.06.
.490 .613 .786 .945 .403 .505 .525 .652 .669 .652 .834 .842 .413 .533 .691 .444 .576 .744 .366 .478 .624 .790 .524
.155 .173 .196 .215 .155 .173 .210 .197 .237 .197 .223
.266 .173 .197 .224 .197 .225 .255 .197 .225 .257 .289 .257
102 160 260 370 46.3 71.8 39.1 121 63.0 121 198 100 32.3 54.1 90.7 25.8 43.3 71.8 12.0 20.5 35.0 55.7 17.4
66 83 106 127 81 101 105 130 134 130 167 168 124 160 207 193 251 323 232 304 396 502 472
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
,
P
,
I
500
I
i
I
163
I
J
i
I
I
,
, _ _
/ -1 O0 -
400 /200
/
E ¢-- -2oo
300
-
/
/
s°°j ~ f
E v.
"o
-300
:>,
(a)
-400
I
I
100
200
'
0
'
I
f
300
400
200-
~500"
100 -
~
500
x (km) o
'
0 -100 -
S
I
200
'
I
300
400
'
500
x (km)
~3s°~
Fig. 6. Topography anomaly (m) for solution 3 (Table 2), assuming local isostatic compensation of the thermal anomaly in the mantle (crustal contribution ignored).
E V' r-- -200 Q.. 03 "0
-3oo
(b)
-400 0
I 100
'
I 200
'
1 300
'
I 400
500
y (km) Fig. 5. Temperature (°C) for solution 3 (Table 2) in the planes (a) x = 0 and (b) y = 0.
by definition, the good overall agreement confirms that the plume thickness scales approximately as S O - - ( a / / o - ) 1/4 even when the viscosity is variable.
4. Models
100
for the Iceland
hotspot
We now apply our variable viscosity model to Iceland, the best known example of a hotspot associated with a ridge-centered plume [12]. The MidAtlantic Ridge near Iceland exhibits pronounced anomalies in bathymetry, petrology, geochemistry and isotopic ratios relative to 'normal' ridge segments far from plumes. Fig. 9, for example, shows the topography along the ridge (circles) as a function
of distance from Iceland. Also shown is a symmetrized representation of the data (solid line), obtained by projecting the raw data onto a great circle with pole at 36.7°E, -14.0°S, smoothing using a + 100 km moving average, and then averaging the northern and southern branches about a presumed center of symmetry at 19.9°W, 65.6°N. Our aim in this section is to estimate the values of the maximum temperature contrast AT and the buoyancy flux B for a ridge-centered plume beneath Iceland by trying to fit the observed topography
4
I
I
I
t
l i l , l
0
_j~3-
2
.
.
.
.
.
.
.
.
i
10
0
,
,
100
FIb Fig. 7. Scaled waist width W / L s for the runs of Table 2 (integer run numbers only) as a function of buoyancy number H b. The straight line is defined by Eq. 22.
164
N.M. Ribe et al. / Earth and Planetary ScienceLetters 134 (1995) 155-168
anomaly (Fig. 9). This anomaly is the sum of two parts: a dynamic component supported by flow in the plume head, and a static component due to the excess crustal thickness generated by pressure-release melting in the plume. Strictly speaking, the dynamic topography should be calculated from the normal stress ~zz at z = h. However, the imposed (arbitrary) variation of the surface velocity near the ridge would create large artifacts. We therefore calculate the dynamic topography approximately using Eq. 21, where 6T is the excess plume temperature obtained by disregarding the effect of melting. However, the integration is carried out only from the
400
I
I
I
I
E
"-"200' ¢<3
I
600
400
y (km) 0.5
I
,
I
,
I
,
I
,
0.4-
o0.3 CO U) 0 . 2 -
0.1
0.0
,
l
,
i
0.2
0.4
0.6
0.8
i
i
i
i
i
,
i
i
i
r
o o
4oo-
0 E ~ -400o
-~ -a0o-1200
-
-l~o
o
d i s t a n c e a l o n g ridge (km)
Fig. 9. Bathymetryalong the Mid-Atlantic Ridge as a function of distance from Iceland (©, data from [13] and R. Kingsley [pets. commun., 1994]). Negative distances are south of Iceland. The smoothed and symmetrized curve shown was obtained as described in the text.
(a)
200
8o0
1.0
y/W Fig. 8. (a) Along-ridgetopographyprofiles ~ ( y ) (portion y > 3 a only) for the runs with integer numbers in Table 2. (b) Rescaled versions of the same profiles.
surface to a depth of 100 or 200 km to avoid an unrealistically large contribution from the hot plume stem. We ignored the additional contributions to the dynamic topography from the lower density of the unmelted residuum and temperature reduction due to absorption of latent heat of melting, which have opposite signs and are much smaller in magnitude than the static (crustal) topography. Our method for calculating the crustal topography consists of six steps: (1) Determine the melt fraction X ( x , y , z ) at all points in the model box, using the melting parameterization of McKenzie and Bickle [14]. The surface z = h of our 3D model is taken to be the Moho, and the pressure there is calculated from the load of the overlying crust. The model temperature is interpreted as potential temperature and is corrected for the effect of latent heat absorption during melting using the method of Watson and McKenzie [15], but assuming that the entropy of fusion is AS = 400 J k g - 1 oC - 1 [16]. (2) Calculate the local melt production rate F ( x , y , z ) = po E. V X , where ~ ( x , y , z ) is the model velocity. (3) Calculate the total melt production rate q ( y ) per unit ridge length by integrating F over vertical planes: q ( y ) = f(1 - X ) - I F d x dz.
N.M. Ribe et al. // Earth and Planetary Science Letters 134 (1995) 155-168 I
i
I
i
I
10-
v
oQ.
o
0
~ ...........
0
200
- ...........
400
,
600
800
distance along ridge (km) Fig. 10. Predicted topography anomaly (light line) for a plume with AT = 263°Cand B = 1390 kg s- 1, comparedto the smoothed observed topography(heavy line) from Fig. 9. Dynamictopography due to excess plume temperaturebetween the surface and the 100 km depth is shown by the dashed line.
(4) Determine the crustal thickness c(y) by equating q(y) to the flux of new crust carried away by the diverging plates [16]. (5) Repeat steps (1-4) iteratively to correct for the effect of the extra pressure imposed on the Moho by the crust, which tends to suppress melting beneath the ridge. (6) Use an Airy isostatic model [17] to transform c(y) to an equivalent topography anomaly relative to an assumed reference depth of 2900 m [12,17] for 'normal' ridge far from the plume. To apply our model to Iceland, we choose a half-spreading rate U = 2.73 × 10 -1° m s -1 (0.85 cm yr 1; [12]) and a mantle (potential) temperature Tr = 1324°C. Our model without a plume then predicts a 'normal' crustal thickness c(~) = 5.5 km, an appropriate average value for a slow-spreading ridge [16]. For our initial modeling attempt, we assume the values AT = 263°C and B = 1390 kg s -1 (Q = 45.7 m 3 s -1) estimated by Schilling [12] for the Iceland plume. This value of B is similar to that (B = 1400 kg s -1) estimated by Sleep [18]. In our model, these values of AT and B correspond to a plume radius a = 62 km. The size of the model box and the grid is the same as in the previous section. Fig. 10 shows the dynamic (dashed line) and the total (light line) topography predicted by the model, together with the
165
smoothed topography (heavy line) from Fig. 9. Close to the plume, the topography anomaly is dominated by the crustal contribution (light line minus dashed line). At y = 0, the predicted topography is ~ 10 km, which corresponds to a crustal thickness of 68 km (the crustal thickness would be 320 km if the damping effect of the extra crustal pressure were ignored). This far exceeds the crustal thickness of ~ 20-34 km predicted for the same AT but purely passive upwelling [14,19]. Due to the active plume upwelling, the upward flux of source rock to shallow depth under the ridge at y = 0 is an order of magnitude higher than the flux for passive spreading. At a few hundred kilometers distance from the plume not much excess melting occurs, although the temperature below the ridge is significantly higher than ambient mantle temperature (compare Fig. 5a with parameter values similar to those for the Iceland model). This is because the plume stem rises almost vertically to shallow depth, where it loses a great deal of its potential melt, and then already depleted rock spreads almost horizontally along the strike of the ridge. In any case, Fig. 10 clearly shows that the topography anomaly predicted by this model is inconsistent with the observations. Fig. 11 shows the topography anomaly predicted for a much cooler and broader plume with AT = 93°C
i
I
~
I
i
I
3-
8' 1
o
i
i
i
400
800
1200
distance along ridge (km) Fig. 11. Same as Fig. 10, but for a plume with AT=93°C, a = 300 km, and B = 1370 kg s -1. Additional topography data ( ~ ) from [19]. Dynamic topography obtained by integration to the 200 km depth.
166
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
and radius a = 300 km. In order to achieve an acceptable fit to the topography at great distances from Iceland, we increased the mantle potential temperature outside the plume to 1347°C compared to 1324°C for normal upper mantle. This simulates in a crude way a very broad thermal halo of the plume or a generally anomalous upper mantle beneath the North Atlantic. Compared to the temperature of the halo, the excess temperature of the plume is 70°C. The buoyancy flux for this plume is B = 1370 kg s -1, the volume flux is Q = 193 m 3 s -1, and the mean upwelling velocity is 2 cm yr -1. The numerical solution was calculated in a box of 800 km in length and 1200 km in width on a grid of 129 × 161 X 65 points. The predicted topography anomaly is now much closer to the observed one, although it is still too narrow in the center. An even wider and cooler upwelling could probably lead to a better fit, but limitations of computer memory prevented us from obtaining a numerical solution in the larger model box that would be required for this.
5. Discussion
The principal goal of this study has been to determine the scaling law for the 'waist width' W for a plume rising beneath an ocean ridge. We found that both a simple lubrication theory model and a more realistic 3D convection model predict similar scaling laws (Eq. 10 and 22). This indicates that the essential dynamics of plume-ridge interaction are largely independent of whether the buoyancy of the plume is thermal or chemical in origin. Nevertheless, there are two relatively minor differences between the two laws. The first is that the lubrication theory model predicts a dependence of W on the mantleplume viscosity contrast "y, whereas no such dependence is evident in the 3D convection model. However, neither model is really capable of resolving this dependence properly: the lubrication theory model incorporates the viscosity contrast only in an approximate way, and numerical difficulties prevent us from using large lateral viscosity contrasts in the 3D model. The second difference between the two laws is that the waist width W is about 25% higher for the 3D model (Eq. 22) than for the lubrication theory model (Eq. 10). However, this discrepancy is only
apparent, because the volume flux Q (unlike the buoyancy flux B) is not uniquely defined for a thermal plume. We chose to define Q as B / p o t~AT, where AT is the maximum temperature contrast across the plume. However, it is probably more appropriate to replace AT in this definition by a smaller average temperature difference ATav, which would reduce the difference between the two results. Another effect that could explain some of the difference between the two results is conductive heat loss to the surface in the variable viscosity model. Our variable viscosity numerical results provide no evidence that the waist width is significantly affected by ' upslope' flow of plume material toward the ridge along the sloping base of the rheological lithosphere. This result may seem surprising until one considers that whereas upslope flow tends to increase W, the 'freezing' of plume material into the lithosphere near the ridge tends to decrease it. Moreover, the magnitude of both effects increases as the slope of the base of the lithosphere increases. We therefore suggest that the absence of evidence for upslope flow in our numerical experiments is due to a balance between the opposing effects of upslope flow and freezing. We do however anticipate that upslope flow may play a greater role in the dynamics of off-ridge plumes. The most surprising result of this study, however, is that the observed anomalous topography of the Mid-Atlantic Ridge suggests the presence of a very cool (AT< 100°C), wide (radius a > 300 km) and slow (vertical velocity ~ 2 cm yr 1) upwelling beneath Iceland. By contrast, a hot (AT = 263°C) and narrow (a = 62 km) plume of the type hypothesized by Schilling [12] predicts a topography anomaly that is far too narrow and a rate of crust production directly over the plume stem that is far too high. Whatever values of a and AT are assumed, the static component of the topography due to excess crustal thickness exceeds the dynamic component due to flow in the plume head. The topography anomaly our model predicts thus depends on the assumptions made in calculating the excess crustal thickness. The most important of these assumptions is that melt generated beneath the ridge is extracted sufficiently rapidly that the ridge-parallel component of its velocity is zero. If this were not so, melt generated near the plume might be extracted further
N.M. Ribe et al. /Earth and Planetary Science Letters 134 (1995) 155-168
away from it, thereby widening the topography anomaly and making a 'narrow plume' model more viable. The importance of along-ridge melt migration can be estimated by means of a simple calculation. The velocity of the melt at any point (O,y,z) directly beneath the ridge is u'f = v~ + (w + wf)2, where v and w are the lateral and vertical components of velocity of the unmelted matrix and wf is the velocity of buoyancy-driven vertical melt migration relative to the matrix. If melt migration obeys Darcy's Law for flow in a porous medium, wf is given by
k4~gApf wf-
~b~?~--~-
(23)
where k6 is the permeability of the matrix, ¢ is the porosity, Apf is the difference between the densities of the matrix and the melt, and ~/f is the viscosity of the melt. Because the matrix has a component of velocity in the ~ direction, a parcel of melt ascending from depth h - z to the surface will travel a distance Ay ~ (h - z)v/(w + wf) along the ridge. As an example, consider the point (x,y,z)= (0 km,100 km,275 km) in our 'narrow plume' numerical solution where the along-ridge velocity v has its maximum value. At this point, v = 1.6 × 1 0 - 9 m s -1 and w = 1.0 × 1 0 - 9 m s 1. Moreover, by using the values [20] k6 = 1 0 - 1 4 m 2, A p f = 600 kg m 3, r/f = 1 Pa s and ~b = 0.01 we find wf = 5.9 × 10 9 m s -1. Thus Ay ~ 30 km. Lateral melt migration over such small distances is incapable of widening the topography anomaly predicted by our 'narrow plume' model (Fig. 10) sufficiently to match the observations. Our disregard of lateral melt migration is therefore valid to a first order. We close by noting that the low excess temperature of the Iceland plume in our model (AT = 93°C) is consistent with its low buoyancy flux. The buoyancy flux of the Iceland plume is smaller than that of the Hawaiian plume by a factor of 3 - 5 [8,18]. In the case of Hawaii, an excess temperature of 250-300°C is in agreement with the observed geometry of the topographic swell [8], and is required to explain the production rates and composition of melts [15]. O1son et al. [23] have shown that the rate at which plumes cool and broaden by diffusion and entrainment as they rise through the mantle depends inversely on the buoyancy flux. If we assume that both
167
the Hawaiian plume and the Iceland plume rise from the core-mantle boundary, where they start with the same excess temperature, the theory of Olson et al. predicts that the Hawaiian plume would be scarcely affected by diffusive cooling and broadening. For the Iceland plume, however, the effect would be quite significant, in agreement with our results. By implication, plumes that are much weaker than the Iceland plume cannot possibly rise from the core-mantle boundary.
Acknowledgements This research was supported by NSF grant EAR 92-04644 and a fellowship from the Alexander yon Humboldt-Stiflung (Germany) to N.M.R., and by grants Ch-77/4 and Ch-77/8 from the Deutsche Forschungsgemeinschafi to U.R.C. We thank J. Bown, L. Fleitout, R. Kingsley, J.G. Schilling and R. White for their generous assistance, M. Feighner and M. Richards for making their work available to us before publication, and C. Langmuir, D. McKenzie and an anonymous referee for helpful reviews. [CL]
Appendix A Lubrication theory equation The general differential equation governing the steady thickness S(x,y) of a thin buoyant fluid layer spreading in an ambient flow whose horizontal velocity is uve = u~} + v~ can be obtained using lubrication theory [7], and is
~7. foS~dz = 0-~7234 + Q6( x) 6(y)
(24)
where Vnz is the horizontal Laplacian and ~r and Q are defined in the text. For the case of a 2D corner flow described by a streamfunction ~ ( x , z )
V f o dZ=OxfSudz--Ox[ ,(x,S)] An analytical expression for ~ ( x , boundary condition ~z(X,O)= U determined by writing ~ = ~z(°)(x,0) = U and ~z(1)(x,0) =
(25)
z) that satisfies the tanh(x/d) can be ~(oh/t (1), where U[1 - tanh(x/d)].
168
N.M. Ribe et al. / Earth and Planetary Science Letters 134 (1995) 155-168
The streamfunction qt(0)- 2 . r r - l U z tan-l(x/z) is just the standard corner flow solution [21], and 1/f(1) is the correction required by the modified velocity near x = 0. The biharmonic equation governing ~(1) can be solved by applying a sine transform and invoking formula 3.987.1 of [22], yielding q,(1) = _~Udf °~ sinkx [ 1kd 2
2 sinh(~kd/2) 7r ]z
× exp( - kz)dk
(26)
This integral can be evaluated using Eq. 3.941.1 and 4.131.3 of [22], whereupon we find that the first part cancels ~(0). The final result, after evaluation at z=S, is 2
(1
S
/x)
~ ( x,S) = -~USIm tp ~ + zrd + zr---d
(27)
where 'Im' denotes the imaginary part and ~O(z)-= (d/dz)ln F(z) is the digamma function of a complex argument [9].
References [1] W.J. Morgan, Plate motions and deep mantle convection, Mem. Geol. Soc. Am. 132, 7-22, 1972. [2] M. Talwani, C.C. Windisch and M.G. Langseth, Jr., Reykjanes Ridge crest: A detailed geophysical study, J. Geophys. Res. 76, 473-517, 1971. [3] J.G. Schilling, Iceland mantle plume: Geochemical evidence along Reykjanes Ridge, Nature 242, 565-571, 1,973. [4] J.G. Schilling, Upper mantle heterogeneities and dynamics, Nature 314, 62-67, 1985. [5] M.A. Feighner and M.A. Richards, The fluid dynamics of plume-ridge and plume-plate interactions: an experimental investigation, Earth Planet. Sci. Lett. 129, 171-182, 1995. [6] M.A. Feighner, L.H. Kellogg and B.J. Travis, Numerical modeling of chemically buoyant mantle plumes at spreading ridges, Geophys. Res. Lett. 22, 715-718, 1995. [7] P. Olson, Hot spots, swells and mantle plumes, in: Magma
Transport and Storage, M.P. Ryan, ed., pp. 33-51, Wiley, New York, 1990. [8] N.M. Ribe and U. Christensen, Three-dimensional modeling of plume-lithosphere interaction, J. Geophys. Res. 99, 669682, 1994. [9] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1046 pp., Dover, New York, 1965. [10] U. Christensen, Convection with pressure and temperature dependent non-Newtonian rheology, Geophys. J. R. Astron. Soc. 77, 242-284, 1984. [11] U. Christensen and H. Harder, 3D convection with variable viscosity, Geophys. J. Int. 104, 213-226, 1991. [12] J.G. Schilling, Fluxes and excess temperatures of mantle plumes inferred from their interaction with migrating mid-ocean ridges, Nature 352, 397-403, 1991. [13] J.G. Schilling, M. Zajac, R. Evans, T. Johnston, W. White, J.D. Devine and R. Kingsley, Petrologic and geochemical variations along the mid-Atlantic ridge from 29°N to 73°N, Am. J. Sci. 283, 510-586, 1983. [14] D. McKenzie and M. Bickle, The volume and composition of melt generated by extension of the lithosphere, J. Petrol. 29, 625-679, 1988. [15] S. Watson and D. McKenzie, Melt generation in plumes: a study of Hawaiian volcanism, J. Petrol. 32, 501-537, 1991. [16] J.W. Bown and R.S. White, Variation with spreading rate of oceanic crustal thickness and geochemistry, Earth Planet. Sci. Lett. 121, 435-449, 1994. [17] G. Schubert and D. Sandwell, Crustal volumes of the continents and of oceanic and continental submarine plateaus, Earth Planet. Sci. Lett. 92, 234-246, 1989. [18] N.H. Sleep, Hotspots and mantle plumes: some phenomenology, J. Geophys. Res. 95, 6715-6736, 1990. [19] E.M. Klein and C.H. Langmuir, Global correlation of ocean ridge basalt chemistry with axial depth and crustal thickness, J. Geophys. Res. 92, 8089-8115, 1987. [20] M.J. Cordery and J. Phipps Morgan, Convection and melting at mid-ocean ridges, J. Geophys. Res. 98, 19477-19503, 1993. [21] G.K. Batchelor, An Introduction to Fluid Dynamics, 615 pp., Cambridge University Press, Cambridge, 1967. [22] I.M. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, 1160 pp., Academic Press, New York, 1965. [23] P. Olson, G. Schubert and C. Anderson, Structure of axisymmetric mantle plumes, J. Geophys. Res. 98, 6829-6844, 1993.